yade-users team mailing list archive

[Zoheir Khademian <question698948@xxxxxxxxxxxxxxxxxxxxx>]

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

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(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=.33,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=1,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=5
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,0,0] 
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0] 
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,0,1]
thermal.thermalBndCondValue=[0,0,0,0,0,45]
thermal.fluidK = 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


timing.reset()
O.dt=1e-4
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('VTK'+identifier+'/',withBoundaries=True)

O.engines=O.engines+[PyRunner(iterPeriod=2000,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(100000,0)

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