Running thermal engine with different initial temperatures

Asked by Zoheir Khademian on 2020-07-17

I am trying to re-assign different temperature values after one cycle of thermal engine by "body.state.temp". However, the conduction scheme stops working and newly assigned body temperature does not change by the boundary conditions. Here is the script:

#*************************************************************************
# Copyright (C) 2019 by Robert Caulk *
# <email address hidden> *
# *
# This program is free software; it is licensed under the terms of the *
# GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************/
#
# Script demonstrating the use of ThermalEngine by comparing conduction
# scheme to analytical solution to Fourier (rod cooling with constant
# boundary conditions). See details in:
#
# Caulk, R., Scholtes, L., Kraczek, M., Chareyre, B. (In Print) A
# pore-scale Thermo-Hydro-Mechanical coupled model for particulate systems.
# Computer Methods in Applied Mechanics and Engineering. Accepted July 2020.
#

from yade import pack
from yade import timing
import numpy as np
import shutil
timeStr = time.strftime('%m-%d-%Y')
num_spheres=1000# number of spheres
young=1e6
rad=0.003

mn,mx=Vector3(0,0,0),Vector3(1.0,0.008,0.008) # corners of the initial packing

thermalCond = 2. #W/(mK)
heatCap = 710. #J(kg K)
t0 = 400. #K

r = rad
k = 2*2.0*r # 2*k*r
Cp = 710.
rho = 2600.
D = 2.*r
m = 4./3.*np.pi*r**2/rho
# macro diffusivity
thermalDiff = 6.*k/(D*np.pi*Cp*rho)

identifier = '-conductionVerification'

if not os.path.exists('VTK'+timeStr+identifier):
 os.mkdir('VTK'+timeStr+identifier)
else:
 shutil.rmtree('VTK'+timeStr+identifier)
 os.mkdir('VTK'+timeStr+identifier)

if not os.path.exists('txt'+timeStr+identifier):
 os.mkdir('txt'+timeStr+identifier)
else:
 shutil.rmtree('txt'+timeStr+identifier)
 os.mkdir('txt'+timeStr+identifier)

shutil.copyfile(sys.argv[0],'txt'+timeStr+identifier+'/'+sys.argv[0])

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)

O.bodies.append(pack.regularOrtho(pack.inAlignedBox(mn,mx),radius=rad,gap=-1e-8,material='spheres'))

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

ThermalEngine = ThermalEngine(dead=1,label='thermal');

newton=NewtonIntegrator(gravity=(0,0,-10), 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,
 VTKRecorder(iterPeriod=500,fileName='VTK'+timeStr+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=1,label='VTKrec'),
 newton
]

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

#for b in O.bodies:
 #if isinstance(b.shape, Sphere):
  #b.dynamic=False

# we only need flow engine to detect boundaries, there is no flow computed for this
flow.dead=0
flow.defTolerance=-1 #0.3
flow.meshUpdateInterval=-1 #200
flow.useSolver=4
flow.permeabilityFactor=-0.7e-7 #1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,0,0,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.debug=False
flow.fluidRho = 1000
flow.fluidCp = 4184
flow.fluidK = 0.650
flow.bndCondIsTemperature=[1,1,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]

flow.tZero=t0
flow.pZero=0
thermal.dead=0
thermal.conduction=True
thermal.thermoMech=False
thermal.advection=False
thermal.fluidThermoMech = False
thermal.solidThermoMech = False
thermal.fluidConduction= False

thermal.bndCondIsTemperature=[1,1,0,0,0,0]
thermal.thermalBndCondValue=[0,0,0,0,0,0]
thermal.tsSafetyFactor=0
thermal.particleDensity=2600
thermal.particleT0=t0
thermal.particleCp=heatCap
thermal.particleK=thermalCond
thermal.particleAlpha =11.6e-3
thermal.useKernMethod=False

timing.reset()
#ThermalEngine.dead=0

flow.updateTriangulation=True
O.dt=1.
O.dynDt=False

O.run(1,1)
flow.dead=1

for b in O.bodies:
 if isinstance(b.shape,Sphere) and b.state.pos[1]>.5:
  b.state.temp=700.
 else:
  b.state.temp=100.

#triax.goal2=-11000

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

# solution to the heat equation for constant initial condition , BCs=0, and using series for approx
def analyticalHeatSolution(x,t,u0,L,k):
 ns = np.linspace(1,1000,1000)
 solution = 0
 for i,n in enumerate(ns):
  integral = (-2./L)*u0*L*(np.cos(n*np.pi)-1.) / (n*np.pi)
  solution += integral * np.sin(n*np.pi*x/L)*np.exp((-k*(n*np.pi/L)**2)*t)
 return solution

# find 10 bodies along x axis
axis = np.linspace(mn[0], mx[0], num=11)
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]
np.savetxt('txt'+timeStr+identifier+'/xdata.txt',axisTrue)
print("Axis length used for analy ", max(axisTrue)-min(axisTrue))
from yade import plot

## a function saving variables
def history():
 plot.addData(
  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,
  temp6 = axisBodies[5].state.temp,
  temp7 = axisBodies[6].state.temp,
  temp8 = axisBodies[7].state.temp,
  temp9 = axisBodies[8].state.temp,
  temp10 = axisBodies[9].state.temp,
  temp11 = axisBodies[10].state.temp,
)

O.engines=O.engines+[PyRunner(iterPeriod=500,command='history()',label='recorder')]
##make nice animations:
#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk(withBoundaries=False)')]
VTKrec.dead=0
from yade import plot

plot.plots={'t':(('temp4','k-'),('temp3','r-'))} #
#plot.plots={'t':(('bndFlux1','k-'),('bndFlux2','r-'))} #
plot.plot()
O.saveTmp()
O.timingEnabled=1
from yade import timing
print("starting oedometer simulation")
O.run(20000,1)
timing.stats()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Zoheir Khademian
Solved:
2020-07-22
Last query:
2020-07-22
Last reply:
2020-07-21
Robert Caulk (rcaulk) said : #1

> the conduction scheme stops working

I'm sorry but we cannot help you without more information. In what way does it "stop working"? Please read closely [1] and elaborate the details necessary.

>newly assigned body temperature does not change by the boundary conditions

Again, please provide details for why you say this.

Cheers,

Robert

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

Zoheir Khademian (zkhademian) said : #2

Thanks Robert for the explanation. I assigned temperature values of 700 for the right half and 100 for the left half of the pack as shown above.

When I run the script after initializing , the temperature of spheres does not change as if there is no conduction between spheres. It seems I am fixing the temperature but I only wanted to initialize them and let the model predict final temperature values. Am I missing some thing?Is there a way to initialize temperature values instead of fixing them?

Robert Caulk (rcaulk) said : #3

Why did you remove the walls?

Why did you add gravity?

Have you checked to see if your particles are still touching despite these two changes?

Robert Caulk (rcaulk) said : #4

And why did you add dynamics?

If you are focused on a single problem, such as how to add a gradient of temperature across the specimen. It is best to make only changes relative to that first, then start adding your other desired modifications.

Zoheir Khademian (zkhademian) said : #5

Thank you Robert. I cleaned up the script and everything worked as expected.