Unstable advection simulaton for gas flow

Hey Zoheir,

> The problem is the model becomes unstable.

What does that mean? Is there an error? If so, what is the error? If not, what is happening?

> I used this example [1] and only modified the fluid properties to represent vapor.

Are you sure? When I use those parameters and decrease the timestep proportionally to viscosity, the model is stable. Please provide an MWE [1] ;).

> I used fluid density of 0.013 kg/m3 and heat capacity of 1996 J/kgK and viscosity of 0.00013 pa.s.

The stability (read maximum allowable timestep) of the coupled thermo-hydro system depends on the viscosity indeed. In fact, assuming the pressure field remains identical, a back of the napkin scale analysis demonstrates that the time step is directly proportional to viscosity. So in this case you should only need to decrease the timestep by 1 order of magnitude. As I mention above, it is in fact stable by decreasing timestep by 1 order of magnitude. If you change the particle size, however, you are playing with a quadratic relationship between particle size and maximum allowable timestep.

> I reduced the timestep by two orders of magnitude and still no luck.

Hence why I suspect you changed particle sizes or pressure field.

> One more question: How can I specify different thermal properties (like conductivity) to different spheres in the model?

I am happy to answer this question but please open a separate thread. The goal here is to maintain an organized knowledge base that is easy for future yade users to search and parse information. Thanks for helping us build this knowledge base :-). If you are interested in the guidelines that we adhere to on launchpad, you are welcome to review them here [2]

[1]https://www.yade-dem.org/wiki/Howtoask

Cheers,

Robert

Thanks, Robert

>What does that mean? Is there an error?
There is no error but temperature of spheres goes up to 1e200 in a few timestep. I

>Are you sure? When I use those parameters and decrease the timestep proportionally to viscosity, the model is stable. Please provide an MWE [1] ;).

I put the MWE here, which is basically this example [1] but viscosity, fluid density, and heat capacity are reduced to match those of gas. I also reduced the time step from 1e-4 to 1e-6. If you run this MWE, you see that both fluid and solid temperature go up to 1e200 after afew timesteps.

> I am happy to answer this question but please open a separate thread.
Sorry, I will open another question

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/ThermalEngine/flowScenario.py

from yade import pack, ymport
from yade import timing
import numpy as np
import shutil
num_spheres=1000# number of spheres
young=1e9
rad=0.003

mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05) # 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

identifier = '-flowScenario'

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('5cmEdge_1mm.spheres', 'x_y_z_r',color=(0.1,0.1,0.9), material='spheres'))

print('num bodies ', len(O.bodies))

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,
)

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),
 triax,

 newton
]

#goal = -1e5
#triax.goal1=triax.goal2=triax.goal3=goal

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.00013#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 = 0.013#997
flow.fluidCp = 1996#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
timing.reset()

O.dt=0.1e-5#0.1e-3
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

#bodyOfInterest = bodyByPos(15.998e-3,0.0230911,19.5934e-3)
bodyOfInterest = bodyByPos(0.025,0.025,0.025)

# find 10 bodies along x axis
axis = np.linspace(mn[0], mx[0], num=5)
axisBodies = [None] * len(axis)
axisTrue = np.zeros(len(axis))
for i,x in enumerate(axis):
 axisBodies[i] = bodyByPos(x, mx[1]/2, mx[2]/2)
 axisTrue[i] = axisBodies[i].state.pos[0]

print("found body of interest at", bodyOfInterest.state.pos)

from yade import plot

## a function saving variables
def history():
 plot.addData(
  ftemp1=flow.getPoreTemperature((0.024,0.023,0.02545)),
  p=flow.getPorePressure((0.025,0.025,0.025)),
  t=O.time,
  i = O.iter,
  temp1 = axisBodies[0].state.temp,
  temp2 = axisBodies[1].state.temp,
  temp3 = axisBodies[2].state.temp,
  temp4 = axisBodies[3].state.temp,
  temp5 = axisBodies[4].state.temp,
  bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp
)

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

def endFlux():
 if O.time >= 30:
  flux = 0
  n=utils.porosity()
  for i in flow.getBoundaryVel(1):
   flux +=i[0]*i[3]/n # area * velocity / porosity (dividing by porosity because flow engine is computing the darcy velocity)
  massFlux = flux * 997

  K = abs(flow.getBoundaryFlux(1))*(flow.viscosity*0.5)/(0.5**2.*(10.-0))
  d=8e-3 # sphere diameter
  Kc = d**2/180. * (n**3.)/(1.-n)**2

  print('Permeability', K, 'kozeny', Kc)
  print('outlet flux(with vels only):',massFlux,'compared to CFD = 0.004724 kg/s')
  print('sim paused')
  O.pause()

O.engines=O.engines+[PyRunner(iterPeriod=10,command='endFlux()')]

from yade import plot

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

print("starting thermal sim")
O.run()

Ok, thank you for the MWE.

Yes I see now, your density value is being used automatically for the calculation of the Reynolds number, which is used to estimate a Nusselt number, which is used to estimate the heat flux between the fluid and the particles. So you are unstable in the conduction domain in this case.

The current Nusselt estimation is empirical for water in dense granular packing. For simulating a gas, you will need to either assign the Reynolds number manually (thermal.uniformReynolds) or modify the source code to generate a proper local Nusselt number for your gas.

Cheers,

Robert

Thanks, Robert.

Could you give some hints or direct me to a reference on the relationship between timestep and heat transfer properties? What factors effects maximum allowable timestep?

Zoheir

If you only seek the maximum time step for the conduction, you would consider the finite difference scheme currently employed and identify the thermal resistivity that destabilizes the scheme. You need to consider both resistivities associated with particle-particle and particle-pore. In your current case, the Reynolds number (Nusselt) contributes to the total thermal resistivity between the particle-pore. Obviously, these thermal resistivities are highlighted in our paper [1].

In other words, you need to set up a finite difference approximation for a single reservoir (particle) connected to other particles and other pores, which is highlighted in the paper [1]. The consideration of stability is simple, you consider this reservoir at 0 temperature initial state, and perturb it by delta T. You can then use the resulting finite-difference formulation to identify the timestep that prevents the next temperature step from oscillating below 0.

Of course, this is only considering the stability of the finite-difference conduction scheme. In your current case, you are also considering advective heat transfer, which is two-way coupled with the conduction scheme. Further, you are simulating fluid fluxes so there is a maximum timestep associated with that implicit scheme - much more computationally intensive to compute directly, but we highlight a method in the appendix of our other paper [2]. So things can become complicated if you do not know which regime you are in, or if your problem crosses regimes.

[1]https://www.sciencedirect.com/science/article/abs/pii/S0045782520304771
[2]https://www.sciencedirect.com/science/article/abs/pii/S0010465519303340?via%3Dihub

Thanks, Robert. This solved my problem.