Stability of advection modeling by ThermalEngine

Asked by Zoheir Khademian

Hello,

I am working on the advection through a spherical packing. The problem is when thermal advection is on in the model, after about 30000 steps or about 2 s the body and cell temperatures exceed the target temp (by several orders of magnitudes) and thus the model becomes unstable. When I reduce the timestep to 1e-5 (from 1e-4), the same problem comes up but after 8 s into the simulation.

So, my question is the origin of this instability if it is not merely me missing something in the model. If not, is it correct that there may be some polyhedral mesh with bad geometries somewhere in the model away from the heat source so it takes a couple thousand steps before the flow reaches these bad geometries and model crashes? If this is the case, how can I identify these "bad" geometries and remove them from calculations? Any other reasons for this delayed instability?

Here is the MWE:

from yade import pack, ymport, plot, utils, export, timing
import numpy as np

young=5e6

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)

identifier = '-thm_coupling'

if not os.path.exists('VTK'+identifier):
 os.mkdir('VTK'+identifier)

if not os.path.exists('txt'+identifier):
 os.mkdir('txt'+identifier)

O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=.33,num=200,seed=11)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(

  maxMultiplier=0,
  finalMaxMultiplier=0,
  thickness = 0,
  goal1=-50000,
  goal2=-50000,
  goal3=-50000,
  max_vel=0.1,
  stressMask = 7,
  internalCompaction=False,
  dead=True,
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 VTKRecorder(iterPeriod=500,fileName='VTK'+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=1,label='VTKrec'),
 NewtonIntegrator(damping=0.5)
]
O.dt=PWaveTimeStep()*.3
O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
triax.dead=False

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.1 and abs(triax.goal3-triax.meanStress)/-triax.goal3<0.001:
    break
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.dynamic=False
triax.internalCompaction=False
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
print(minX,minY,minZ)
print(maxX,maxY,maxZ)
dz=maxZ-minZ
dy=maxY-minY
dx=maxX-minX

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 10
flow.meshUpdateInterval=5
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,10]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]
flow.tZero=25

flow.dead=0

thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,0,1]
thermal.thermalBndCondValue=[0,0,0,0,0,45]
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2700
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0

timing.reset()
O.dt=1e-4
O.dynDt=False
thermal.dead=0
flow.emulateAction()

def bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))

from yade import plot

def MaxTempInModel():
 MaxFluidTemp=0
 for i in range(0,flow.nCells()):
  MaxFluidTemp=max(MaxFluidTemp,flow.getCellTemperature(i))
 MaxTemp=0
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   MaxTemp=max(MaxTemp,b.state.temp)
 print("######### MAX SOLID TEMPERATURE IS:",MaxTemp,"##########")
 print("######### MAX FLUID TEMPERATURE IS:",MaxFluidTemp,"##########")
O.engines=O.engines+[PyRunner(iterPeriod=100,command='MaxTempInModel()',label='MaxTempInModel')]

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature(((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))),

  t=O.time,
  i = O.iter,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)
 #plot.saveDataTxt('Record.txt',vars=('t','i','p','ftemp1','bodyOfIntTemp'))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

def pressureField():
 flow.saveVtk('VTK'+identifier+'/',withBoundaries=True)

O.engines=O.engines+[PyRunner(iterPeriod=2000,command='pressureField()')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))}

plot.plot(subPlots=False)
def ColorScaler():
 for s in O.bodies:
  s.shape.color=scalarOnColorScale(s.state.temp,25,45)

O.engines=O.engines+[PyRunner(command='ColorScaler()',iterPeriod=20)]
ColorScaler()

O.run(100000,0)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Hey Zoheir,

>So, my question is the origin of this instability if it is not merely me missing something in the model.

You are applying constant 45 degrees to 1 of the 6 particle boundaries, yet you have no heat-flux for all 6 boundaries of the flow advection calculation <- that seems suspect. I would apply more traditional boundary conditions a bit to confirm that you can reach stable states for your setup:

flow.bndCondIsTemperature=[0,0,0,0,1,1]
flow.thermalBndCondValue=[0,0,0,0,25,45]
thermal.bndCondIsTemperature=[0,0,0,0,1,1]
thermal.thermalBndCondValue=[0,0,0,0,25,45]

In fact, it brings me to the main question: which example script are you modelling this after? What changes have you made in comparison to the example script? Which change results in the instability?

>>If not, is it correct that there may be some polyhedral mesh with bad geometries somewhere in the model away from the heat source so it takes a couple thousand steps before the flow reaches these bad geometries and model crashes? If this is the case, how can I identify these "bad" geometries and remove them from calculations?

You can guarantee a clean mesh by setting relFuzz=0. This will set a monodisperse cubic packing that you can be sure will have a clean triangulation.

That said, the triangulation is not actually an issue in the pore-finite volume (PFV) scheme. Inherent to the method we don't actually care about tetrahedral shape - it is simply used for geometric reductions. The one issue we may run into is when two tetrahedra centers occur on top of one another (centers are bary centers in PFV). This usually only occurs in regular packings, and yade is designed to handle these situations elegantly. You can take some manual control with thermal.minimumFluidCondDist[1]. However, you should also be able to figure out if that is the issue by simply turning on and off thermal.fluidConduction.

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.minimumFluidCondDist

p.s. turning various physics (thermal.advection, conduction, fluidConduction) on and off is a good way to start locating the issue. You can also turn debug on to get some helpful messages. You can *not* break it :)

Let me know what you come up with,

Robert

Revision history for this message
Zoheir Khademian (zkhademian) said :
#2

Hi Robert and thanks for the response.

>>What changes have you made in comparison to the example script?
1- Turning internal compaction off in the triaxial engine and compacted the sample by inwards-moving walls.
2- Reducing the particle density from 2600e10 to 2600 to run the triaixal engine.
3- Setting bodies dynamic false after compaction so there is no movement.
4- Turning off the thermo-mechanical analyses in thermal engine.
5- flow.bndCondIsPressure=[0,0,0,0,1,1] ; flow.bndCondValue=[0,0,0,0,0,10]
6- changing rRelFuzz from 0.333 to 0.33 #Small changes in the packing can make the model unstable like changing (mn,mx) a bit.

After making these changes, when I run, the model maximum body and fluid temp exceed the boundary temp after 2 s.

>> I would apply more traditional boundary conditions a bit to confirm that you can reach stable states for your setup:
I changed the boundary conditions to the traditional ones as you suggested but it did not make any difference in stability.

>>You can take some manual control with thermal.minimumFluidCondDist[1].
I tried assigning different values (0.0015 to 0.015) to minimumFluidCondDist[1] but it did not make any difference.

>>However, you should also be able to figure out if that is the issue by simply turning on and off thermal.fluidConduction.
I turned off thermal.fluidConduction and it solved the instability. I checked the centers of neighboring cells and some of them had pretty close X or Y or Z coordinates (less than 1e-5). Is this the issue? If so, could you give me some hint on how to avoid this issue.

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.minimumFluidCondDist

Revision history for this message
Robert Caulk (rcaulk) said :
#3

Hey Zoheir,

Thanks for the good feedback. Can you precise to me the example script that you are basing this on? We have 4 example scripts in the repository. I guess you are using flowScenario.py since you had b.dynamic=False before as well as a pressure boundary condition?

>> If so, could you give me some hint on how to avoid this issue.

We've discussed in the past how to estimate the stable timestep dynamically [1]. Did you explore that already for your case? I started adding that feature to ThermalEngine [2][3], but it is still experimental.

Let's try to debug this through some kind of trial and error, I think it will identify the problem most rapidly. We can consider the factors discussed in [1] and hold various components constant. Since you know it is in the fluid conduction, then you just need to consider the parameters associated with fluid conduction such as thermal.fluidConductionAreaFactor[4] and thermal.fluidK. As we discuss in [1], increasing the diffusion coefficient [5] is going to decrease the maximum allowable timestep. You have pythonic control over all 3 parameters in that diffusion equation which means you should be able to determine the exact cause of the instability through trial and error. I would first start by holding them all constant at extraordinary values i.e. force that diffusion constant to be extremely low to test if the system is perfectly stable. Then open up each parameter back to your desired values indivdiually. BTW to be strict about the test you need to set flow.meshUpdateInterval=-1 and b.dynamic=False - since those two components adjust geometric quantities used for the diffusion coefficient and which evolve with time during the simulation.

Another easy test would be to not run compaction (basically just holding bodies static) and see if the simulation remains stable.

Thanks for testing this stuff, let me know what you come up with.

[1]https://answers.launchpad.net/yade/+question/695542
[2]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L177
[3]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.tsSafetyFactor
[4]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.fluidConductionAreaFactor
[5]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L544

Thanks,

Robert

Revision history for this message
Zoheir Khademian (zkhademian) said :
#4

Thanks so much Robert for the great suggestions. I forgot to mention that the example I am basing the MWE on is [1] but I added the wall compaction and made the particles dynamic false after compaction.

I am going through the trial and error approach you suggested, but I had a question about what you said earlier in this thread: "The one issue we may run into is when two tetrahedra centers occur on top of one another". How can this be an issue? I mean how the alignment of two meshes on the top of each other can contribute to the instability of TH? Could you direct me to the section in the source code that may be sensitive to this specific configuration of meshes?

I will get back to you about the results of the trial and error exercise.

Thank again

Zoheir

[1] https://gitlab.com/yade-dev/trunk/-/blob/07085490a7a47f24aa962f914b6ba062ea11a026/examples/ThermalEngine/thermoHydroMechanical_coupling.py

Revision history for this message
Robert Caulk (rcaulk) said :
#5

Hey Zoheir,

The calculation of the diffusion coefficient for the fluid conduction scheme is inversely proportional to the distance between the centers i.e. dividing by zero is a no-no.

>>Could you direct me to the section in the source code that may be sensitive to this specific configuration of meshes?

[1]

Cheers,

Robert

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L544

Revision history for this message
Robert Caulk (rcaulk) said :
#6

>> I forgot to mention that the example I am basing the MWE on is [1] but I added the wall compaction and made the particles dynamic false after compaction.

If it is true, you forgot to mention you changed the meshUpdateInterval.

It is important to keep track of all changes in comparison to the example script, this way we can determine the culprit.

Cheers,

Robert

Revision history for this message
Zoheir Khademian (zkhademian) said :
#7

Thanks Robert!
I went through a trial and error process using a more traditional BC as you suggested:

low.bndCondIsTemperature=[0,0,0,0,1,1]
flow.thermalBndCondValue=[0,0,0,0,25,45]
thermal.bndCondIsTemperature=[0,0,0,0,1,1]
thermal.thermalBndCondValue=[0,0,0,0,25,45]

and avoid mesh updating by meshUpdateInterval=-1

- First, I turned off the compaction and ran the model. Temperature ranges seemed to be within expected range 25-45 K. Turning on the compaction made the model unstable in a way that solid and fluid temp goes above the assigned boundary condition (45 K).

- I tried to minimize the particle-pore and pore-pore resistivities by assigning very small numbers to thermal.fluidConductionAreaFactor[1] and thermal.fluidK[2]. I also made sure that the distance between cells wont be zero by assigning thermal.minimumFluidCondDist[3] of 0.001. I tried different numbers but the particle temp still goes up to 1e3 very quickly while cell temp stays below 45 K.
- I also tried avoiding fictious cells by thermal.ignoreFictiousConduction[4]=True but still no chance. I also tried using HertzMethod [5] for area calculation. The cell temp was in the range but particle temp went to negative.
- I also limited the Reynolds (thermal.uniformReynolds[6]) to 10 to make sure Nusselt number calc is not the reason.

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.fluidConductionAreaFactor
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.fluidK
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.minimumFluidCondDist
[4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.ignoreFictiousConduction
[5] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.useHertzMethod
[6] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.uniformReynolds
Here is the MWE for reference:
###########################################
from yade import pack, ymport, plot, utils, export, timing
import numpy as np

young=5e6

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05)

identifier = '-thm_coupling'

if not os.path.exists('VTK'+identifier):
 os.mkdir('VTK'+identifier)

if not os.path.exists('txt'+identifier):
 os.mkdir('txt'+identifier)

O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=.333,num=200,seed=11)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(

  maxMultiplier=0,
  finalMaxMultiplier=0,
  thickness = 0,
  goal1=-50000,
  goal2=-50000,
  goal3=-50000,
  max_vel=0.1,
  stressMask = 7,
  internalCompaction=False,
  dead=True,
)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
 triax,
 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 VTKRecorder(iterPeriod=500,fileName='VTK'+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=0,label='VTKrec'),
 NewtonIntegrator(damping=0.5)
]
O.dt=PWaveTimeStep()*.3
O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
triax.dead=False

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.1 and abs(triax.goal3-triax.meanStress)/-triax.goal3<0.001:
    break
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.dynamic=False
triax.internalCompaction=False
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
print(minX,minY,minZ)
print(maxX,maxY,maxZ)
dz=maxZ-minZ
dy=maxY-minY
dx=maxX-minX

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 10
flow.meshUpdateInterval=-1
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,10]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,1,1]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,25,45]
flow.tZero=25

flow.dead=0

thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,1,1]
thermal.thermalBndCondValue=[0,0,0,0,25,45]
thermal.fluidK = 1e-6#0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2700
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0
thermal.minimumFluidCondDist=1e-3
thermal.fluidConductionAreaFactor=1e-6
#thermal.ignoreFictiousConduction=True
timing.reset()
O.dt=1e-3
O.dynDt=False
thermal.dead=0
flow.emulateAction()

def bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))

from yade import plot

def MaxTempInModel():
 MaxFluidTemp=0
 for i in range(0,flow.nCells()):
  MaxFluidTemp=max(MaxFluidTemp,flow.getCellTemperature(i))
 MaxTemp=0
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   MaxTemp=max(MaxTemp,b.state.temp)
 print("######### MAX SOLID TEMPERATURE IS:",MaxTemp,"##########")
 print("######### MAX FLUID TEMPERATURE IS:",MaxFluidTemp,"##########")
O.engines=O.engines+[PyRunner(iterPeriod=100,command='MaxTempInModel()',label='MaxTempInModel')]

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature(((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))),

  t=O.time,
  i = O.iter,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)
 #plot.saveDataTxt('Record.txt',vars=('t','i','p','ftemp1','bodyOfIntTemp'))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

def pressureField():
 flow.saveVtk('VTKPressure'+identifier+'/',withBoundaries=True)

O.engines=O.engines+[PyRunner(iterPeriod=500,command='pressureField()')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))}

plot.plot(subPlots=False)
def ColorScaler():
 for s in O.bodies:
  s.shape.color=scalarOnColorScale(s.state.temp,25,45)

O.engines=O.engines+[PyRunner(command='ColorScaler()',iterPeriod=20)]
ColorScaler()

O.run(2000,0)

Revision history for this message
Robert Caulk (rcaulk) said :
#8

Hey zoheir,

I just ran your MWE out 160 seconds, 170k iterations. It is stable, temperatures do not exceed the boundary conditions. Am I missing something?

Which version of yade are you running?

Cheers,

Robert

Revision history for this message
Zoheir Khademian (zkhademian) said :
#9

Hey Robert,

I ran it in both Yadedaily 20201215-4510~714a723~bionic1 Using python version: 3.6.9
and "eraseBodyThermalEngine" [1] trunk version.

Of 10 times of running the model, I get 8 stable results and 2 unstable results (after 10000 iterations)
If I change the sphere size to rMean=0.001 from 0.0015, all tries lead to instability. Please reduce the mean sphere size to 0.001 and let me know if you can still get stable results.

[1]https://gitlab.com/yade-dev/trunk/-/tree/eraseBodyThermalEngine

Revision history for this message
Robert Caulk (rcaulk) said :
#10

Hey Zoheir,

I have no way of testing that yadedaily version, and it is not recommended to be compiling old branches of yade. That was a branch that was merged into the master, so it is best to stick to the master as you get updates. But, I can confirm that running the latest master version is stable in all the condtions you asked me to test. Please update your sources and recompile:

cd ../trunk
git remote update
git pull upstream master
cd ../build
make -jX install # X being the numberof cores on your computer

Let me know if this fixes it, sorry for the delay,

Robert

Revision history for this message
Zoheir Khademian (zkhademian) said :
#11

Hi Robert,

I updated my Yade version to yade-2021.01a and tried the script above with rMean=0.001. The particle temperature still goes up to 1000 K after 20000 iterations. I even tried the script on other machines but got the same results. Please note that with rMean=0.0015, the temp range is OK (under 45 K) but once I change the mean size to rMean=0.001, the particle temperature exceeds the boundary temp and goes to 1000 K.

Would you be able to try this one more time with rMean=0.001 on your machine? I would really appreciate your help

Zoheir

Revision history for this message
Robert Caulk (rcaulk) said :
#12

Hey Zoheir,

The version that you have installed is actually an *earlier* version than the branch you were using. Best to start using a *later* version so as to keep up with bug fixes etc. For example, I am using the latest, it is 52f0c5a7, but installing it is as simple as following these instructions [1].

>>Would you be able to try this one more time with rMean=0.001 on your machine? I would really appreciate your help

Yes, it works when I decrease the timestep, as expected. The magnitude of timestep decrease depends on the seed and packing - we are sampling a uniform particle distribution [0.001+/-0.00033] 100 times. So sometimes we may get a tiny sphere other times we may not. That has significant impact on stability as I will teach to you below. (btw, we do not have the same seed values since we are on separate computers.) So you can change seed to demonstrate that process.

I think it will help you for me to review to you directly the fluid conduction scheme:

All particles are triangulated to create a connected set of tetrahedra. Each tetrahedron has 4 neighbors, and with those 4 neighbors, it shares 1 facet each. The fluid conduction is computed using the fluid area (A) of the incident facet divided by the distance between the two neighboring tetrahedra centers (L). Fluid area is simply the area of the facet not consumed by a sphere. Please refer to our paper for visualization of these geometries. So we have:

q_12 = k_12 * A_12/L_12 * (T1 - T2)

As we have discussed at length, the time step depends on this "diffusion coefficient" kA/L. As that increases, the allowable timestep decreases. So what does this mean? If the area is big, the timestep is small. If L is small, the timestep is small, if k is big, the timestep is small. We can refer to A/L as a characteristic length scale of the problem.

When you start adding a wide distribution of spheres (rRelFuzz=0.33), you start to play with these geometries in the triangulation indirectly. Which means you are adjusting your maximum allowable timestep. So if you are constrained by the timestep of 1e-3, but you insist on decreasing the characteristic length scale of the fluid conduction scheme (decreasing sphere size), then you need to make some sacrifices elsewhere. At least, you can demonstrate this is the root of your instability by manually controlling the diffusion coefficient via python:

thermal.minimumFluidCondDist=rMean
thermal.fluidConductionAreaFactor=0.1

This means you have a fluid area that is anomalously large compared to the rest of your triangulation, or a length that is anomalously small. I leave it to you to play with these factors to determine which one it is.

By the way, depending on the sensitivity of the triangulation (it can be quite sensitive when you have a wide range of vertex weights (sphere radii) ), you need to isolate the problem from any sort of external stochasticity such as parallelized compaction such as you are showing here. You have a slightly different triangulation each time (even using seed in make cloud). Thus, best practice is to export the packing and then import the identical packing each time.

Cheers,

Robert

[1]https://yade-dem.org/doc/installation.html#source-code

Revision history for this message
Robert Caulk (rcaulk) said :
#13
Revision history for this message
Zoheir Khademian (zkhademian) said :
#14

Hey Robert,
Sorry for my delays. I try to give your comments enough thoughts and do investigation before replying.

I tried both released version [1] and development version [2] and still observed the instabilities meaning the particle temp goes above the max boundary temp while cell temp is still within the range.

Regarding the diffusivity coefficient kA/L, in the MWE I uploaded last, I have reduced kA/L by reducing fluidk and fluidConductionAreaFactor to arbitrarily low values. I also tried your values (thermal.minimumFluidCondDist=rMean and thermal.fluidConductionAreaFactor=0.1) but still the particle temp is unstable.

Finally I looked at the thermal.cpp and noticed that applyTempDeltaToSolids [3] adds pore temp differences to the solid temp for no apparent reason (particle temp is updated here [4]). I deactivated applyTempDeltaToSolids and recompiled the source. I ran the MWE and the results are stable now.

Let me know if applyTempDeltaToSolids must not be deactivated or I am missing something here

Thanks
Zoheir

[1] https://launchpad.net/yade/+download

[2] https://gitlab.com/yade-dev/trunk

[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L94

[4] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L605

Revision history for this message
Robert Caulk (rcaulk) said :
#15

Hey Zoheir,

This is quite odd: applyTempDeltaToSolids is deactivated by default: the
user needs to manually set thermal.delT>0 for that function to activate
[1]. As you see in [1], thermal.delT is set to 0 if user does nothing. So
compiling without the line you indicated will change nothing in the output
unless your user script manually activated thermal.delT>0. I don’t see you
setting that value in any of the MWEs that you provided here. So this
cannot be the issue (unless you are using a different script to the one
provided here where you use delT).

Can you please find the unstable packing, then export those sphere radii
and location to text file, the copy paste that file to this thread so I can
replicate exactly your problem. Also please copy paste the identical script
that goes with running that packing.

Then please type here the exact output of yade -v where yade is the
executable you are using (it may be named differently depending on how you
compile or which package you download).

Cheers,

Robert

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.hpp#L133

Le mar. 19 oct. 2021 à 23:01, Zoheir Khademian <
<email address hidden>> a écrit :

> Question #698948 on Yade changed:
> https://answers.launchpad.net/yade/+question/698948
>
> Zoheir Khademian posted a new comment:
> Hey Robert,
> Sorry for my delays. I try to give your comments enough thoughts and do
> investigation before replying.
>
> I tried both released version [1] and development version [2] and still
> observed the instabilities meaning the particle temp goes above the max
> boundary temp while cell temp is still within the range.
>
>
> Regarding the diffusivity coefficient kA/L, in the MWE I uploaded last, I
> have reduced kA/L by reducing fluidk and fluidConductionAreaFactor to
> arbitrarily low values. I also tried your values
> (thermal.minimumFluidCondDist=rMean and
> thermal.fluidConductionAreaFactor=0.1) but still the particle temp is
> unstable.
>
> Finally I looked at the thermal.cpp and noticed that
> applyTempDeltaToSolids [3] adds pore temp differences to the solid temp
> for no apparent reason (particle temp is updated here [4]). I
> deactivated applyTempDeltaToSolids and recompiled the source. I ran the
> MWE and the results are stable now.
>
> Let me know if applyTempDeltaToSolids must not be deactivated or I am
> missing something here
>
> Thanks
> Zoheir
>
>
> [1] https://launchpad.net/yade/+download
>
> [2] https://gitlab.com/yade-dev/trunk
>
> [3] https://gitlab.com/yade-
> dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L94
>
> [4] https://gitlab.com/yade-
> dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L605
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Robert Caulk (rcaulk) said :
#16

Also, please indicate your version of Ubuntu.

Revision history for this message
Zoheir Khademian (zkhademian) said (last edit ):
#17

Hi Robert,
>>This is quite odd: applyTempDeltaToSolids is deactivated by default: the
user needs to manually set thermal.delT>0 for that function to activate

I think the problem is delT is being used as a variable here as well [1], so during fluid-fluid conduction calculations, delT becomes positive and it activates applyTempDeltaToSolids automatically. I deleted detT variable here [1] and directly used "cell->info().temp() - neighborCell->info().temp()" for conductionEnergy calculation here[2]. Then recompiled and the results are stable.

Let me know if this does not makes sense to you.

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L530
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L545

In any case, here is the list of particles (all the same size) and the MWE to run the model.

Thanks - Zoheir

#format x_y_z_r
0.0177741 0.0293044 0.0178294 0.0015
0.032225 0.0197924 0.0178472 0.0015
0.0322266 0.0291739 0.0178539 0.0015
0.0279236 0.032249 0.0261029 0.0015
0.0322118 0.0302264 0.0320853 0.0015
0.0272061 0.0316643 0.0320534 0.0015
0.029361 0.0178245 0.0290514 0.0015
0.0223492 0.0237293 0.0212975 0.0015
0.0322662 0.0177355 0.029063 0.0015
0.0249392 0.0297701 0.0178535 0.0015
0.0278041 0.0253036 0.0320824 0.0015
0.0209031 0.0239399 0.0259082 0.0015
0.029996 0.0247971 0.0227321 0.0015
0.0231086 0.0224061 0.0178694 0.0015
0.0237889 0.026771 0.0280437 0.0015
0.0271796 0.0296372 0.0298686 0.0015
0.0293275 0.0207069 0.0294883 0.0015
0.029728 0.0307252 0.0238577 0.0015
0.029455 0.0177438 0.03204 0.0015
0.0220291 0.0251347 0.0320876 0.0015
0.0228457 0.0177657 0.0205171 0.0015
0.0247976 0.020158 0.0187476 0.0015
0.0221528 0.017793 0.0297048 0.0015
0.0177809 0.0178246 0.0320527 0.0015
0.0178286 0.0204492 0.0192357 0.0015
0.0191019 0.0272907 0.0227524 0.0015
0.023182 0.0266745 0.0241058 0.0015
0.0228497 0.0296931 0.0320761 0.0015
0.0203349 0.0184782 0.0277075 0.0015
0.0178115 0.0232242 0.0181458 0.0015
0.0236951 0.0187262 0.0273578 0.0015
0.017815 0.0178142 0.0290756 0.0015
0.0268862 0.0177484 0.0275311 0.0015
0.0293104 0.0293306 0.0320274 0.0015
0.0178224 0.0228895 0.0229322 0.0015
0.0229362 0.0182461 0.0233847 0.0015
0.0234661 0.0204362 0.0296891 0.0015
0.0308867 0.0228868 0.0282905 0.0015
0.0322099 0.0250677 0.0296246 0.0015
0.032281 0.0178211 0.0320448 0.0015
0.0322146 0.0204259 0.0245334 0.0015
0.0322253 0.0322509 0.0239898 0.0015
0.0279165 0.0228391 0.0279433 0.0015
0.0234598 0.0322574 0.0219894 0.0015
0.0202221 0.0215691 0.017819 0.0015
0.0292976 0.025252 0.0289125 0.0015
0.0178179 0.0322929 0.0207288 0.0015
0.0219679 0.0204761 0.0198587 0.0015
0.0219003 0.0322297 0.0244164 0.0015
0.0188566 0.0202206 0.0221634 0.0015
0.0276402 0.0177882 0.0232129 0.0015
0.0178161 0.0322694 0.0269258 0.0015
0.0295711 0.0280749 0.0290676 0.0015
0.0265758 0.0178165 0.0320975 0.0015
0.0321235 0.0295753 0.0253156 0.0015
0.0203762 0.0177426 0.0220801 0.0015
0.0292814 0.0231458 0.0251311 0.0015
0.0299286 0.0268797 0.0186235 0.0015
0.0288255 0.0202145 0.0243078 0.0015
0.0178049 0.0276979 0.0202425 0.0015
0.0322722 0.0256006 0.0178161 0.0015
0.0223008 0.0202683 0.0253642 0.0015
0.0255165 0.0266773 0.0206327 0.0015
0.0199502 0.0216179 0.0289218 0.0015
0.030126 0.0177624 0.0178698 0.0015
0.0254739 0.0197561 0.0231361 0.0015
0.0322638 0.0228734 0.0228553 0.0015
0.0265408 0.022283 0.0178242 0.0015
0.0254871 0.0259749 0.0258332 0.0015
0.0322915 0.0227056 0.0178061 0.0015
0.0319731 0.0202977 0.0214805 0.0015
0.0235946 0.0323112 0.0178428 0.0015
0.0250886 0.0203689 0.0321208 0.0015
0.0311899 0.0205841 0.0320759 0.0015
0.0293867 0.0323455 0.0178564 0.0015
0.017819 0.0319405 0.0320753 0.0015
0.0230816 0.0291808 0.0260803 0.0015
0.0322467 0.0233608 0.0320651 0.0015
0.027053 0.0279131 0.0274663 0.0015
0.0191108 0.0322681 0.0295038 0.0015
0.0298198 0.0309172 0.0294285 0.0015
0.0265891 0.027943 0.0320729 0.0015
0.0177748 0.0178067 0.0235421 0.0015
0.0226982 0.0292137 0.0291665 0.0015
0.0292801 0.0177808 0.0258092 0.0015
0.025091 0.0177539 0.0296773 0.0015
0.0322576 0.0321652 0.0178683 0.0015
0.0303121 0.0204239 0.0267807 0.0015
0.0263177 0.0252864 0.0286105 0.0015
0.0220065 0.032246 0.0292096 0.0015
0.0204973 0.0258364 0.0205884 0.0015
0.0228527 0.0303463 0.0198524 0.0015
0.0266579 0.029408 0.0204424 0.0015
0.0199349 0.0270645 0.0318538 0.0015
0.0207689 0.0322449 0.020817 0.0015
0.0177812 0.0248037 0.0206708 0.0015
0.0273755 0.0289692 0.0232743 0.0015
0.0232237 0.0229252 0.0245358 0.0015
0.0265182 0.0322593 0.0177995 0.0015
0.0195693 0.0213444 0.0248116 0.0015
0.0243073 0.0322565 0.0320575 0.0015
0.0236199 0.0178027 0.0321506 0.0015
0.0211121 0.0268755 0.0261824 0.0015
0.0322276 0.0276056 0.0203509 0.0015
0.0322252 0.0177805 0.0199996 0.0015
0.0229195 0.0277541 0.0178142 0.0015
0.0206611 0.0177697 0.0321174 0.0015
0.0272737 0.0322313 0.0232492 0.0015
0.0267728 0.0226499 0.0314472 0.0015
0.0177887 0.0256031 0.0246928 0.0015
0.0273974 0.0323373 0.0289861 0.0015
0.0295738 0.0219256 0.0222416 0.0015
0.0213171 0.0322629 0.0320744 0.0015
0.0206219 0.02938 0.0180342 0.0015
0.0184764 0.0280415 0.0267233 0.0015
0.0201655 0.0183195 0.0249349 0.0015
0.0239698 0.0272352 0.030993 0.0015
0.0187189 0.0227822 0.0321331 0.0015
0.0207101 0.0182969 0.0182606 0.0015
0.0259197 0.0300618 0.025721 0.0015
0.0235804 0.0177613 0.0177803 0.0015
0.0248667 0.0246848 0.0320612 0.0015
0.0191357 0.0322515 0.0233922 0.0015
0.0231734 0.0294893 0.0230771 0.0015
0.0322124 0.0178675 0.026089 0.0015
0.0321633 0.027037 0.0273948 0.0015
0.0302813 0.0217504 0.0195146 0.0015
0.0177592 0.0289692 0.0320756 0.0015
0.0263988 0.0204391 0.0294554 0.0015
0.029371 0.0296149 0.0182665 0.0015
0.020224 0.0298721 0.0308167 0.0015
0.017831 0.0194442 0.0263944 0.0015
0.0265363 0.0177864 0.0178349 0.0015
0.0203174 0.0275839 0.0289467 0.0015
0.0319785 0.0246858 0.0205816 0.0015
0.0177702 0.0298453 0.0221658 0.0015
0.0178166 0.0223711 0.0269294 0.0015
0.0300936 0.0322618 0.0320543 0.0015
0.0322475 0.0205852 0.0293348 0.0015
0.0322126 0.0297129 0.0223999 0.0015
0.0177868 0.0265145 0.02991 0.0015
0.0296558 0.0240756 0.0178264 0.0015
0.0231045 0.0251686 0.0188838 0.0015
0.0253767 0.0238163 0.0199777 0.0015
0.0276649 0.0261025 0.0239278 0.0015
0.0283892 0.0266703 0.0210978 0.0015
0.0322364 0.0298972 0.0282113 0.0015
0.0269041 0.0204408 0.0265127 0.0015
0.0295258 0.0293948 0.0212187 0.0015
0.0245398 0.0231523 0.0295862 0.0015
0.0178015 0.0261738 0.0177968 0.0015
0.026021 0.0322515 0.0206297 0.0015
0.0301032 0.0277903 0.0236731 0.0015
0.0322264 0.0278804 0.0303835 0.0015
0.0229634 0.0223551 0.0320091 0.0015
0.0307124 0.0259274 0.032048 0.0015
0.0177696 0.0301867 0.0249686 0.0015
0.0293013 0.0296542 0.0265791 0.0015
0.0283612 0.0200495 0.0181101 0.0015
0.0322158 0.0258758 0.0244254 0.0015
0.0256938 0.0230807 0.0260243 0.0015
0.0322614 0.0322903 0.0299867 0.0015
0.0271081 0.0279631 0.0178916 0.0015
0.0282175 0.0203153 0.0321007 0.0015
0.0178086 0.0204051 0.0305596 0.0015
0.0248368 0.0246643 0.0227298 0.0015
0.0207381 0.0246782 0.0294686 0.0015
0.0178079 0.0178202 0.0206456 0.0015
0.0296785 0.0259988 0.0260802 0.0015
0.0178222 0.0178226 0.0177958 0.0015
0.0322144 0.0232181 0.0257478 0.0015
0.0247606 0.0312801 0.0293061 0.0015
0.0269235 0.0223401 0.0234643 0.0015
0.0269121 0.0251292 0.0178316 0.0015
0.0242797 0.0215566 0.0213311 0.0015
0.0253235 0.0177877 0.0251117 0.0015
0.0235819 0.0321655 0.0267623 0.0015
0.0265927 0.0187846 0.0206085 0.0015
0.019697 0.0225662 0.0205688 0.0015
0.0321963 0.0178108 0.0230955 0.0015
0.0216489 0.0212353 0.0226878 0.0015
0.0225049 0.022332 0.0276505 0.0015
0.017783 0.0322525 0.0178251 0.0015
0.0295761 0.0232292 0.0309606 0.0015
0.020183 0.0246431 0.0232522 0.0015
0.0186716 0.025176 0.0274007 0.0015
0.0205473 0.0292455 0.0244887 0.0015
0.0178014 0.0296022 0.0291546 0.0015
0.0282425 0.0237411 0.0203954 0.0015
0.0206179 0.020632 0.031634 0.0015
0.0320111 0.0322371 0.0208441 0.0015
0.0206879 0.0322691 0.0178517 0.0015
0.0203723 0.0247275 0.0178577 0.0015
0.0227773 0.0275868 0.0207949 0.0015
0.0290156 0.0322597 0.02082 0.0015
0.0205505 0.030286 0.0273229 0.0015
0.0203694 0.0293112 0.0209934 0.0015
0.0177824 0.0235551 0.0296056 0.0015
0.030871 0.0322143 0.0267549 0.0015
0.0293854 0.017816 0.0207962 0.0015

# Here is the MWE for running the model. Please change the address for the location of the txt files (compacted.txt and BinSizeHP.txt).

from yade import pack, ymport, plot, utils, export, timing

import numpy as np

young=5e6

identifier = '-thm_coupling'

if not os.path.exists('VTK'+identifier):
 os.mkdir('VTK'+identifier)

if not os.path.exists('txt'+identifier):
 os.mkdir('txt'+identifier)

O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600,label='walls'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))

O.saveTmp('Init')

_sp = ymport.textExt('compacted.txt',format='x_y_z_r', material='spheres')
sp = O.bodies.append(_sp)
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
print ("minX:",minX,"maxX:",maxX,"minY:",minY,"maxY:",maxY,"minZ:",minZ,"maxZ:",maxZ)

Hsize=(minX,minY,minZ,maxX,maxY,maxZ)
output='\t%g\t%g\t%g\t%g\t%g\t%g'%tuple(Hsize[i] for i in range(6))
with open('BinSizeHP.txt', 'w') as f:
    for item in Hsize:
        f.write("%s\n" % item)

O.loadTmp('Init')

with open('BinSizeHP.txt') as f:
  BSP2 = f.readlines()
f.close()

mn,mx=Vector3(float(BSP2[0]),float(BSP2[1]),float(BSP2[2])),Vector3(float(BSP2[3]),float(BSP2[4]),float(BSP2[5])) # corners of the initial packing

walls=aabbWalls([mn,mx],material='walls',thickness=0)#xx*.1)
wallIds=O.bodies.append(walls)
_sp = ymport.textExt('compacted.txt',format='x_y_z_r', material='spheres')
sp = O.bodies.append(_sp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1,label='is2aabb'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1,label='ss2sc'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),
 #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),

 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 VTKRecorder(iterPeriod=500,fileName='VTK'+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=0,label='VTKrec'),
 NewtonIntegrator(damping=0.5)
]
O.dt=PWaveTimeStep()*.1

O.step()
#ss2sc.interactionDetectionFactor=-1
#is2aabb.aabbEnlargeFactor=-1

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.dynamic=False

maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
print(minX,minY,minZ)
print(maxX,maxY,maxZ)
dz=maxZ-minZ
dy=maxY-minY
dx=maxX-minX

flow.debug=False
# add flow
flow.permeabilityMap = False
flow.pZero = 10
flow.meshUpdateInterval=-1
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,0,10]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,1,1]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,25,45]
flow.tZero=25

flow.dead=0

thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,1,1]
thermal.thermalBndCondValue=[0,0,0,0,25,45]
thermal.fluidK = 0.65
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2700
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =-1
thermal.minimumThermalCondDist=0
thermal.minimumFluidCondDist=0.0015
thermal.fluidConductionAreaFactor=.1
timing.reset()
O.dt=1e-5
O.dynDt=False
thermal.dead=0
flow.emulateAction()

def bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))

from yade import plot

def MaxTempInModel():
 MaxFluidTemp=0
 for i in range(0,flow.nCells()):
  MaxFluidTemp=max(MaxFluidTemp,flow.getCellTemperature(i))
 MaxTemp=0
 for b in O.bodies:
  if isinstance(b.shape,Sphere):
   MaxTemp=max(MaxTemp,b.state.temp)
 print("######### MAX SOLID TEMPERATURE IS:",MaxTemp,"##########")
 print("######### MAX FLUID TEMPERATURE IS:",MaxFluidTemp,"##########")
O.engines=O.engines+[PyRunner(iterPeriod=100,command='MaxTempInModel()',label='MaxTempInModel')]

def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature(((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))),

  t=O.time,
  i = O.iter,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)
 #plot.saveDataTxt('Record.txt',vars=('t','i','p','ftemp1','bodyOfIntTemp'))

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

def pressureField():
 flow.saveVtk('VTKPressure'+identifier+'/',withBoundaries=True)

O.engines=O.engines+[PyRunner(iterPeriod=500,command='pressureField()')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))}

plot.plot(subPlots=False)
def ColorScaler():
 for s in O.bodies:
  s.shape.color=scalarOnColorScale(s.state.temp,25,45)

O.engines=O.engines+[PyRunner(command='ColorScaler()',iterPeriod=20)]
ColorScaler()

O.run(25000,0)

Revision history for this message
Zoheir Khademian (zkhademian) said :
#18

BTW, my Ubuntu version is: 20.04.3

Revision history for this message
Best Robert Caulk (rcaulk) said :
#19

>>I think the problem is delT is being used as a variable here as well [1], so during fluid-fluid conduction calculations, delT becomes positive and it activates applyTempDeltaToSolids automatically

Yes, that is a big bug. Great spot, Zoheir. I am pushing the fix upstream immediately, here [1] you can follow the progress.

This is quite interesting, as it would seem that in the code's current state, it would only work when delT is always a negative value, which depends entirely on how the facet list is created with respect to the boundary conditions. Anyways, this is a necessary fix - thank you.

Cheers,

Robert

[1]https://gitlab.com/yade-dev/trunk/-/merge_requests/786

Revision history for this message
Zoheir Khademian (zkhademian) said :
#20

I am glad we were able to fix this after all.

Thanks

Zoheir

Revision history for this message
Zoheir Khademian (zkhademian) said :
#21

Thanks Robert Caulk, that solved my question.