yade-users team mailing list archive

[Zoheir Khademian <question698948@xxxxxxxxxxxxxxxxxxxxx>]

Question #698948 on Yade changed:
https://answers.launchpad.net/yade/+question/698948

Zoheir Khademian posted a new comment:
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)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.