# Unstable advection simulaton for gas flow

I am trying to simulate advection with a hot gas flowing though a cold packing. I used this example [1] and only modified the fluid properties to represent vapor. I used fluid density of 0.013 kg/m3 and heat capacity of 1996 J/kgK and viscosity of 0.00013 pa.s. The problem is the model becomes unstable.

I reduced the timestep by two orders of magnitude and still no luck. So, I was wondering if there is any other parameters I should modify to make the model stable.

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

Thanks - Zoheir

[1] https:/

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Zoheir Khademian

- Solved:
- 2021-02-18

- Last query:
- 2021-02-18

- Last reply:
- 2021-02-16

Robert Caulk (rcaulk) said : | #1 |

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]

Robert Caulk (rcaulk) said : | #2 |

Cheers,

Robert

Zoheir Khademian (zkhademian) said : | #3 |

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:/

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(

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.

# macro diffusivity

identifier = '-flowScenario'

O.materials.

O.materials.

walls=aabbWalls

wallIds=

sp = O.bodies.

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

triax=TriaxialS

maxMultiplier=

finalMaxMultip

thickness = 0,

stressMask = 7,

internalCompac

)

ThermalEngine = ThermalEngine(

newton=

intRadius = 1

O.engines=[

ForceResetter(),

InsertionSortC

InteractionLoop(

[Ig2_

[Ip2_

[Law2_

),

FlowEngine(

ThermalEngine, GlobalStiffness

triax,

newton

]

#goal = -1e5

#triax.

for b in O.bodies:

if isinstance(b.shape, Sphere):

b.dynamic=False # mechanically static

flow.dead=0

flow.defToleran

flow.meshUpdate

flow.useSolver=4

flow.permeabili

flow.viscosity= 0.00013#0.001

flow.bndCondIsP

flow.bndCondVal

flow.thermalEng

flow.debug=False

flow.fluidRho = 0.013#997

flow.fluidCp = 1996#4181.7

flow.getCHOLMOD

flow.bndCondIsT

flow.thermalEng

flow.thermalBnd

flow.tZero=t0

flow.pZero=0

flow.maxKdivKmean=1

flow.minKdivmea

thermal.dead=0

thermal.debug=False

thermal.

thermal.

thermal.

thermal.

thermal.

thermal.

thermal.

thermal.

thermal.

thermal.fluidK = 0.6069 #0.650

thermal.

thermal.particleT0 = t0

thermal.

thermal.particleK = thermalCond

thermal.particleCp = heatCap

thermal.

#thermal.

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(

for b in O.bodies:

if isinstance(b.shape, Sphere):

dist = b.state.pos - Vector3(x,y,z)

if np.linalg.

cDist = dist

cBody = b

print('found closest body ', cBody.id, ' at ', cBody.state.pos)

return cBody

#bodyOfInterest = bodyByPos(

bodyOfInterest = bodyByPos(

# 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[

print("found body of interest at", bodyOfInterest.

from yade import plot

## a function saving variables

def history():

plot.addData(

ftemp1=

p=flow.

t=O.time,

i = O.iter,

temp1 = axisBodies[

temp2 = axisBodies[

temp3 = axisBodies[

temp4 = axisBodies[

temp5 = axisBodies[

bodyOfIntTemp = O.bodies[

)

O.engines=

def endFlux():

if O.time >= 30:

flux = 0

n=utils.

for i in flow.getBoundar

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.

d=8e-3 # sphere diameter

Kc = d**2/180. * (n**3.)/(1.-n)**2

print(

print('outlet flux(with vels only):'

print('sim paused')

O.pause()

O.engines=

from yade import plot

plot.plots=

plot.plot()

O.saveTmp()

print("starting thermal sim")

O.run()

Robert Caulk (rcaulk) said : | #4 |

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.

Cheers,

Robert

Zoheir Khademian (zkhademian) said : | #5 |

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

Robert Caulk (rcaulk) said : | #6 |

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:/

[2]https:/

Zoheir Khademian (zkhademian) said : | #7 |

Thanks, Robert. This solved my problem.