Re: [Question #695065]: ThermoHydro analyses

```Question #695065 on Yade changed:

Description changed to:
Hello,

I am trying to run a TH analysis with warm fluid going through a cold
pack of polydisperse particles using the provided example in [1].

When I use a previously compacted, poly-disperse-particle sample for
this analysis, the temperature of all particles remains constant at
first and then turns to "nan". I was wondering what needs to be adjusted
in the thermal and flow engines following the features of the granular
sample  under study.

The MWE follows:

Thanks - Zoheir

import numpy as np
import shutil

young=1e9

mn,mx=Vector3(0.0005944,0.0005596,0.00099759),Vector3(0.00119,0.00119,0.0016) # corners of the initial packing
thermalCond = 2. #W/(mK)
heatCap = 710. #J(kg K)
t0 = 333.15 #K

# micro properties
k = 2.0   # 2*k*r
Cp = 710.
rho = 2600.
D = 2.*r
m = 4./3.*np.pi*r**2/rho
# macro diffusivity

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

#sp = O.bodies.append(ymport.textExt('ShCube', 'x_y_z_r',color=(0.1,0.1,0.9), material='spheres'))
O.bodies.append([
sphere((0.000945595,0.000784186,0.00149986),0.0000875846),
sphere((0.000778022,0.000753066,0.00145943),0.0000875846),
sphere((0.000888722,0.000810072,0.00133622),0.0000875846),
sphere((0.000903823,0.000647224,0.00139896),0.0000875846),
sphere((0.00105629,0.000841192,0.00137666),0.0000875846),
sphere((0.00101452,0.00070423,0.00127576),0.0000875846),
sphere((0.00084695,0.00067311,0.00123532),0.0000875846),
sphere((0.000999422,0.000867078,0.00121301),0.0000875846),
sphere((0.00103358,0.00108745,0.0011483),0.0000815788),
sphere((0.000902032,0.00103675,0.00153707),0.0000735711),
sphere((0.000986413,0.00091621,0.00153692),0.0000735711),
sphere((0.00111325,0.000915283,0.00146233),0.0000735711),
sphere((0.00108202,0.000666186,0.00154242),0.0000665643),
sphere((0.00111005,0.000740653,0.00132599),0.0000665643),
sphere((0.00080908,0.00109404,0.0014756),0.0000499888),
sphere((0.000751015,0.00103522,0.00153186),0.0000499888),
sphere((0.000670623,0.00108745,0.00150348),0.0000499888),
sphere((0.000687866,0.00108202,0.00126569),0.0000459731),
sphere((0.000703081,0.00107221,0.00135584),0.0000459731),
sphere((0.00110443,0.00083307,0.00156066),0.0000420579),
sphere((0.00112293,0.000765198,0.00151454),0.0000420579),
sphere((0.00106475,0.000819695,0.0014877),0.0000420579),
sphere((0.0011389,0.000814544,0.00144832),0.0000420579),
sphere((0.00112143,0.000746581,0.00109922),0.0000439814),
sphere((0.00105308,0.000800634,0.00111123),0.0000439814),
sphere((0.00110668,0.000870009,0.00110401),0.0000439814),
sphere((0.00112389,0.000807084,0.00116302),0.0000439814),
sphere((0.00103833,0.000924061,0.00111603),0.0000439814),
sphere((0.000792443,0.000977976,0.00151641),0.0000220555),
sphere((0.000760407,0.000979304,0.00148612),0.0000220555),
sphere((0.000836178,0.000937735,0.0014684),0.0000220555),
sphere((0.000775855,0.00092188,0.00146726),0.0000220555),
sphere((0.000798293,0.000958519,0.00147726),0.0000220555),
sphere((0.00073797,0.000942664,0.00147612),0.0000220555),
sphere((0.000766257,0.000959847,0.00144696),0.0000220555),
sphere((0.000728372,0.000980632,0.00145582),0.0000220555),
sphere((0.00081374,0.000901096,0.0014584),0.0000220555),
sphere((0.000842027,0.000918279,0.00142924),0.0000220555),
sphere((0.000781704,0.000902424,0.00142811),0.0000220555),
sphere((0.000804142,0.000939063,0.0014381),0.0000220555),
sphere((0.000772106,0.000940391,0.00140781),0.0000220555),
sphere((0.000734221,0.000961175,0.00141667),0.0000220555),
sphere((0.000809991,0.000919607,0.00139895),0.0000220555),
sphere((0.00113905,0.00100891,0.00119864),0.0000235576),
sphere((0.000698534,0.000815074,0.00126433),0.0000765702),
sphere((0.000695882,0.00071517,0.00138037),0.0000765702),
sphere((0.000670978,0.000701228,0.00116568),0.0000765702),
sphere((0.000683566,0.000664341,0.00153398),0.0000481936),
sphere((0.00113412,0.000910585,0.00138428),0.000040366),
sphere((0.00112752,0.00102435,0.0013913),0.000040366),
sphere((0.0011135,0.000914201,0.0013063),0.000040366),
sphere((0.00100741,0.000944425,0.00133576),0.000040366),
sphere((0.00108155,0.000966389,0.00135897),0.000040366),
sphere((0.0011069,0.00102797,0.00131333),0.000040366),
sphere((0.00102898,0.00102219,0.00133366),0.000040366),
sphere((0.00107495,0.00108015,0.00136599),0.000040366),
sphere((0.00112847,0.00110573,0.00131123),0.000040366),
sphere((0.00113885,0.000975778,0.00126066),0.000040366),
sphere((0.00106093,0.000970005,0.001281),0.000040366),
sphere((0.00105433,0.00108377,0.00128802),0.000040366),
sphere((0.000905936,0.000944606,0.00110321),0.000105624),
sphere((0.000756408,0.00108773,0.00114542),0.000105624),
sphere((0.00074381,0.000861183,0.00150084),0.0000322503),
sphere((0.000794579,0.000887874,0.00153035),0.0000322503),
sphere((0.000828512,0.000858476,0.00142805),0.0000322503),
sphere((0.000884233,0.00108217,0.00130317),0.000033125),
sphere((0.000901453,0.00111232,0.00135959),0.000033125),
sphere((0.000690365,0.00103048,0.00144934),0.0000328425),
sphere((0.000744769,0.00106632,0.00144098),0.0000328425),
sphere((0.000777308,0.000901361,0.00119494),0.000045),
sphere((0.000928126,0.00101227,0.00131929),0.000045),
sphere((0.000844602,0.000967877,0.00126048),0.000045),
sphere((0.00102741,0.000992656,0.00113012),0.000045),
sphere((0.000758144,0.000949585,0.00128654),0.000045),
sphere((0.00100301,0.00111497,0.00140571),0.000045),
sphere((0.000824715,0.00100449,0.0014298),0.000045),
sphere((0.00101191,0.000708242,0.00142522),0.000045),
sphere((0.000932121,0.000920137,0.00140155),0.000045),
sphere((0.000815047,0.00105458,0.00135521),0.000045),
])

newton=NewtonIntegrator(damping=0.2)
O.engines=[
ForceResetter(),
InteractionLoop(
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
),
ThermalEngine,	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
newton
]

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

flow.defTolerance=-1#0.3
flow.meshUpdateInterval=200
flow.useSolver=4
flow.permeabilityFactor= 1
flow.viscosity= 0.001
flow.bndCondIsPressure=[1,1,0,0,0,0]
flow.bndCondValue=[10,0,0,0,0,0]
flow.thermalEngine=True
flow.debug=False
flow.fluidRho = 997
flow.fluidCp = 4181.7
#flow.getCHOLMODPerfTimings=True
flow.bndCondIsTemperature=[1,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[343.15,0,0,0,0,0]
flow.tZero=t0
flow.pZero=0
#flow.clampKValues=False
flow.maxKdivKmean=1
flow.minKdivmean=0.0001;

thermal.debug=False
thermal.fluidConduction=True
thermal.ignoreFictiousConduction=False#True
#thermal.minimumFluidCondDist=1.89e-5
thermal.conduction=True
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.bndCondIsTemperature=[0,0,0,0,0,0]
thermal.thermalBndCondValue=[0,0,0,0,0,0]
thermal.fluidK = 0.6069#0.650
thermal.fluidConductionAreaFactor=1.
thermal.particleT0 = t0
thermal.particleDensity=2600.
thermal.particleK = thermalCond
thermal.particleCp = heatCap
thermal.useKernMethod=True
#thermal.useHertzMethod=False
timing.reset()

O.dt=0.1e-5
O.dynDt=False

O.run(1,1)

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
print('found closest body ', cBody.id, ' at ', cBody.state.pos)
return cBody
# Finding a body at the X flow boundary
bodyOfInterest = bodyByPos(0.00035,0.001,0.001)

def history():
print(bodyOfInterest.state.temp)
ftemp1=flow.getPoreTemperature((0.00035,0.001,0.001)),
t=O.time,
i = O.iter,
bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp
)

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

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-'))} #
plot.plot()
O.saveTmp()
def ColorScaler():
for s in O.bodies:

s.shape.color=scalarOnColorScale(s.state.temp,333,343)

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

ColorScaler()

O.run(1000)

--