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.
[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/ThermalEngine/flowScenario.py
The MWE follows:
Thanks - Zoheir
from yade import pack, ymport
from yade import timing
import numpy as np
import shutil
young=1e9
rad=0.003
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
r = rad
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=radians(3),density=2600,label='spheres'))
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,material='spheres'),
sphere((0.000778022,0.000753066,0.00145943),0.0000875846,material='spheres'),
sphere((0.000888722,0.000810072,0.00133622),0.0000875846,material='spheres'),
sphere((0.000903823,0.000647224,0.00139896),0.0000875846,material='spheres'),
sphere((0.00105629,0.000841192,0.00137666),0.0000875846,material='spheres'),
sphere((0.00101452,0.00070423,0.00127576),0.0000875846,material='spheres'),
sphere((0.00084695,0.00067311,0.00123532),0.0000875846,material='spheres'),
sphere((0.000999422,0.000867078,0.00121301),0.0000875846,material='spheres'),
sphere((0.00103358,0.00108745,0.0011483),0.0000815788,material='spheres'),
sphere((0.000902032,0.00103675,0.00153707),0.0000735711,material='spheres'),
sphere((0.000986413,0.00091621,0.00153692),0.0000735711,material='spheres'),
sphere((0.00111325,0.000915283,0.00146233),0.0000735711,material='spheres'),
sphere((0.00108202,0.000666186,0.00154242),0.0000665643,material='spheres'),
sphere((0.00111005,0.000740653,0.00132599),0.0000665643,material='spheres'),
sphere((0.00080908,0.00109404,0.0014756),0.0000499888,material='spheres'),
sphere((0.000751015,0.00103522,0.00153186),0.0000499888,material='spheres'),
sphere((0.000670623,0.00108745,0.00150348),0.0000499888,material='spheres'),
sphere((0.000687866,0.00108202,0.00126569),0.0000459731,material='spheres'),
sphere((0.000703081,0.00107221,0.00135584),0.0000459731,material='spheres'),
sphere((0.00110443,0.00083307,0.00156066),0.0000420579,material='spheres'),
sphere((0.00112293,0.000765198,0.00151454),0.0000420579,material='spheres'),
sphere((0.00106475,0.000819695,0.0014877),0.0000420579,material='spheres'),
sphere((0.0011389,0.000814544,0.00144832),0.0000420579,material='spheres'),
sphere((0.00112143,0.000746581,0.00109922),0.0000439814,material='spheres'),
sphere((0.00105308,0.000800634,0.00111123),0.0000439814,material='spheres'),
sphere((0.00110668,0.000870009,0.00110401),0.0000439814,material='spheres'),
sphere((0.00112389,0.000807084,0.00116302),0.0000439814,material='spheres'),
sphere((0.00103833,0.000924061,0.00111603),0.0000439814,material='spheres'),
sphere((0.000792443,0.000977976,0.00151641),0.0000220555,material='spheres'),
sphere((0.000760407,0.000979304,0.00148612),0.0000220555,material='spheres'),
sphere((0.000836178,0.000937735,0.0014684),0.0000220555,material='spheres'),
sphere((0.000775855,0.00092188,0.00146726),0.0000220555,material='spheres'),
sphere((0.000798293,0.000958519,0.00147726),0.0000220555,material='spheres'),
sphere((0.00073797,0.000942664,0.00147612),0.0000220555,material='spheres'),
sphere((0.000766257,0.000959847,0.00144696),0.0000220555,material='spheres'),
sphere((0.000728372,0.000980632,0.00145582),0.0000220555,material='spheres'),
sphere((0.00081374,0.000901096,0.0014584),0.0000220555,material='spheres'),
sphere((0.000842027,0.000918279,0.00142924),0.0000220555,material='spheres'),
sphere((0.000781704,0.000902424,0.00142811),0.0000220555,material='spheres'),
sphere((0.000804142,0.000939063,0.0014381),0.0000220555,material='spheres'),
sphere((0.000772106,0.000940391,0.00140781),0.0000220555,material='spheres'),
sphere((0.000734221,0.000961175,0.00141667),0.0000220555,material='spheres'),
sphere((0.000809991,0.000919607,0.00139895),0.0000220555,material='spheres'),
sphere((0.00113905,0.00100891,0.00119864),0.0000235576,material='spheres'),
sphere((0.000698534,0.000815074,0.00126433),0.0000765702,material='spheres'),
sphere((0.000695882,0.00071517,0.00138037),0.0000765702,material='spheres'),
sphere((0.000670978,0.000701228,0.00116568),0.0000765702,material='spheres'),
sphere((0.000683566,0.000664341,0.00153398),0.0000481936,material='spheres'),
sphere((0.00113412,0.000910585,0.00138428),0.000040366,material='spheres'),
sphere((0.00112752,0.00102435,0.0013913),0.000040366,material='spheres'),
sphere((0.0011135,0.000914201,0.0013063),0.000040366,material='spheres'),
sphere((0.00100741,0.000944425,0.00133576),0.000040366,material='spheres'),
sphere((0.00108155,0.000966389,0.00135897),0.000040366,material='spheres'),
sphere((0.0011069,0.00102797,0.00131333),0.000040366,material='spheres'),
sphere((0.00102898,0.00102219,0.00133366),0.000040366,material='spheres'),
sphere((0.00107495,0.00108015,0.00136599),0.000040366,material='spheres'),
sphere((0.00112847,0.00110573,0.00131123),0.000040366,material='spheres'),
sphere((0.00113885,0.000975778,0.00126066),0.000040366,material='spheres'),
sphere((0.00106093,0.000970005,0.001281),0.000040366,material='spheres'),
sphere((0.00105433,0.00108377,0.00128802),0.000040366,material='spheres'),
sphere((0.000905936,0.000944606,0.00110321),0.000105624,material='spheres'),
sphere((0.000756408,0.00108773,0.00114542),0.000105624,material='spheres'),
sphere((0.00074381,0.000861183,0.00150084),0.0000322503,material='spheres'),
sphere((0.000794579,0.000887874,0.00153035),0.0000322503,material='spheres'),
sphere((0.000828512,0.000858476,0.00142805),0.0000322503,material='spheres'),
sphere((0.000884233,0.00108217,0.00130317),0.000033125,material='spheres'),
sphere((0.000901453,0.00111232,0.00135959),0.000033125,material='spheres'),
sphere((0.000690365,0.00103048,0.00144934),0.0000328425,material='spheres'),
sphere((0.000744769,0.00106632,0.00144098),0.0000328425,material='spheres'),
sphere((0.000777308,0.000901361,0.00119494),0.000045,material='spheres'),
sphere((0.000928126,0.00101227,0.00131929),0.000045,material='spheres'),
sphere((0.000844602,0.000967877,0.00126048),0.000045,material='spheres'),
sphere((0.00102741,0.000992656,0.00113012),0.000045,material='spheres'),
sphere((0.000758144,0.000949585,0.00128654),0.000045,material='spheres'),
sphere((0.00100301,0.00111497,0.00140571),0.000045,material='spheres'),
sphere((0.000824715,0.00100449,0.0014298),0.000045,material='spheres'),
sphere((0.00101191,0.000708242,0.00142522),0.000045,material='spheres'),
sphere((0.000932121,0.000920137,0.00140155),0.000045,material='spheres'),
sphere((0.000815047,0.00105458,0.00135521),0.000045,material='spheres'),
])
ThermalEngine = ThermalEngine(dead=1,label='thermal');
newton=NewtonIntegrator(damping=0.2)
intRadius = 1
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
),
FlowEngine(dead=1,label="flow",multithread=False),
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.dead=0
flow.defTolerance=-1 #0.3
flow.meshUpdateInterval=-1
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.maxKdivKmean=1
flow.minKdivmean=0.0001;
thermal.dead=0
thermal.debug=False
thermal.fluidConduction=True
thermal.ignoreFictiousConduction=True
thermal.conduction=True
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
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
O.dt=0.1e-5
O.dynDt=False
O.run(1,1)
flow.dead=0
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)
from yade import plot
def history():
print(bodyOfInterest.state.temp)
plot.addData(
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)
