Stability of advection modeling by ThermalEngine
Hello,
I am working on the advection through a spherical packing. The problem is when thermal advection is on in the model, after about 30000 steps or about 2 s the body and cell temperatures exceed the target temp (by several orders of magnitudes) and thus the model becomes unstable. When I reduce the timestep to 1e-5 (from 1e-4), the same problem comes up but after 8 s into the simulation.
So, my question is the origin of this instability if it is not merely me missing something in the model. If not, is it correct that there may be some polyhedral mesh with bad geometries somewhere in the model away from the heat source so it takes a couple thousand steps before the flow reaches these bad geometries and model crashes? If this is the case, how can I identify these "bad" geometries and remove them from calculations? Any other reasons for this delayed instability?
Here is the MWE:
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
young=5e6
mn,mx=Vector3(
identifier = '-thm_coupling'
if not os.path.
os.mkdir(
if not os.path.
os.mkdir(
O.materials.
O.materials.
walls=aabbWalls
wallIds=
sp=pack.
sp.makeCloud(
sp.toSimulation
triax=TriaxialS
maxMultiplier=0,
finalMaxMulti
thickness = 0,
goal1=-50000,
goal2=-50000,
goal3=-50000,
max_vel=0.1,
stressMask = 7,
internalCompa
dead=True,
)
O.engines=[
ForceResetter(),
InsertionSortC
InteractionLoop(
[Ig2_
[Ip2_
[Law2_
),
#GlobalStiffne
triax,
FlowEngine(
ThermalEngine(
VTKRecorder(
NewtonIntegrat
]
O.dt=PWaveTimeS
O.step()
ss2sc.interacti
is2aabb.
triax.dead=False
while 1:
O.run(1000, True)
unb=unbalance
print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
if unb<0.1 and abs(triax.
break
for b in O.bodies:
if isinstance(
b.dynamic=False
triax.internalC
maxY=max(
maxX=max(
minY=min(
minX=min(
minZ=min(
maxZ=max(
print(minX,
print(maxX,
dz=maxZ-minZ
dy=maxY-minY
dx=maxX-minX
flow.debug=False
# add flow
flow.permeabili
flow.pZero = 10
flow.meshUpdate
flow.fluidBulkM
flow.useSolver=4
flow.permeabili
flow.viscosity= 0.001
flow.decoupleForces = False
flow.bndCondIsP
flow.bndCondVal
## Thermal Stuff
flow.bndCondIsT
flow.thermalEng
flow.thermalBnd
flow.tZero=25
flow.dead=0
thermal.dead=1
thermal.
thermal.
thermal.debug=0
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = 25
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.
thermal.
thermal.
thermal.
thermal.
timing.reset()
O.dt=1e-4
O.dynDt=False
thermal.dead=0
flow.emulateAct
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
return cBody
bodyOfInterest = bodyByPos(
from yade import plot
def MaxTempInModel():
MaxFluidTemp=0
for i in range(0,
MaxFluidTemp=
MaxTemp=0
for b in O.bodies:
if isinstance(
MaxTemp=
print("######### MAX SOLID TEMPERATURE IS:",MaxTemp,
print("######### MAX FLUID TEMPERATURE IS:",MaxFluidTe
O.engines=
def history():
plot.addData(
ftemp1=
t=O.time,
i = O.iter,
bodyOfIntTemp = O.bodies[
)
#plot.
O.engines=
def pressureField():
flow.saveVtk(
O.engines=
plot.plots=
plot.plot(
def ColorScaler():
for s in O.bodies:
s.shape.
O.engines=
ColorScaler()
O.run(100000,0)
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- Yade Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Robert Caulk
- Solved:
- Last query:
- Last reply: