Today we will learn how to build a script that simulates heat conduction through a spherical packing and compares the numerical values to Fourier’s analytical solution.
Day 3 - Thermal Hands-on part 1¶
We know where to start, let’s import the necessary libraries and set our variables:
from yade import pack
from yade import timing
import numpy as np
num_spheres=1000
young=1e6
rad=0.003
mn,mx=Vector3(0,0,0),Vector3(1.0,0.008,0.008)
These are all recognizable variables from previous hands-on sessions. Next, we append our materials and walls as we’ve done in the past:
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'))
Here we see that we are appending a new type of sphere packing called regularOrtho. As the name suggests, this creates a regular orthogonal packing which will be useful for ensuring that randomness doesn’t affect our comparison to the analytic conduction solution later.
Next, we need to construct our engines list as usual:
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(label="flow"),
ThermalEngine(label='thermal'),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
VTKRecorder(iterPeriod=500,fileName='VTK'+timeStr+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=1,label='VTKrec'),
NewtonIntegrator(damping=0.2)
]
Most of this should look familiar based on our previous hands-on sessions. But we see two important components including FlowEngine and ThermalEngine. These two engines rely intimately on one another for simulating THM processes, and thus ThermalEngine cannot be used without FlowEngine. We instantiate both of these engines without setting any parameters so that we can do so in detail in following steps.
We are only interested in validating the thermal conduction scheme in Yade, so we need to turn many default functionalities off starting with body dynamics:
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.dynamic = False
b.dynamic is a body parameter which tells Yade to consider it for force calculations or not. Setting this value to false ensures that these spheres will not move during the entirety of our simulation.
Next, we set our thermal parameters:
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.particleDensity = 2600 # kg/m^3
thermal.particleT0 = 400 # K
thermal.particleCp = 710 #J(kg K)
thermal.particleK = 2. #W/(mK)
thermal.particleAlpha = 11.6e-3
thermal.useKernMethod = False
The full set of available ThermalEngine parameters and all their specific details can be found here inside our Class Reference. We see that we need to ensure many of the functionalities are set to False for the basic conduction example here. Next, we set our boundary conditions in the same way we learned how to set boundary conditions during the previous FlowEngine hands-on session. Meanwhile, the initial temperature of the particles is set with particleT0. Finally, we set the basic thermal conduction parameters such as the particle density (particleDensity), thermal conductivity (particleK), heat capacity (particleCp), and diffusivity (particleAlpha).
Now we need to employ the FlowEngine for one step so that it can identify the boundaries for our ThermalEngine. We do not require the FlowEngine beyond this step because we are not simulating any fluid fluxes in the present conduction example:
O.dt=1.
O.dynDt=False
flow.updateTriangulation=True
flow.dead=0
flow.emulateAction()
flow.dead=1
Here we see that we are forcing FlowEngine to update the triangulation in a fake timestep with flow.emulateAction. Once this is done, we reset the FlowEngine to dead=1 So that we do not waste computational effort calculating pressure fields.
Gathering field data¶
Since we are comparing our numerical conduction to an analytical scheme, we need a way to obtain field data from arbitrary coordinates. Here is an example of one way to do so:
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
Where we simply feed it arbitrary coordinates and it will return the closest body with which we can extract physical quantities such as temperature, velocity, etc.
Let’s use this function to grab 10 bodies along the x-axis for us to track during the simulation:
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]
Additionally, we need a way to compute the analytical solution. Here is the solution to the heat equation for a uniform initial temperature condition and boundary conditions at 0 K:
k = 2
Cp = 710
rho = 2600
alpha = 6.*k/(np.pi*Cp*rho)
def analyticalHeatSolution(x,t,u0,L,alpha):
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((-alpha*(n*np.pi/L)**2)*t)
return solution
Where x is the x coordainte along the x-axis, t is the time of measurement, u0 is the initial temperature of the rod, L is the length of the rod, and k is the thermal diffusivity of the rod. alpha is an effective thermal diffusivity which scales the discrete elements to cubical continuum elements.
Finally, we need to collect and plot the data during the simulation. The temperature can be obtained via the body state. And you have the bodies of interest set in axisBodies. Using the information from previous hands-on sessions, fill out the following template to collect data for
def history():
plot.addData(
t = O.time,
i = O.iter,
temp1 = ________,
temp2 = ________,
temp3 = ________,
tempAnalytic1 = analyticalHeatSolution(________),
tempAnalytic2 = analyticalHeatSolution(________),
tempAnalytic3 = analyticalHeatSolution(________)
)
plot.saveDataTxt('conductionAnalyticalComparison.txt',vars=('t','i','temp1','temp2','temp3','tempAnalytic1','tempAnalytic2','tempAnalytic3'))
O.engines=O.engines+[PyRunner(iterPeriod=500,command='history()',label='recorder')]
Use the lessons we learned from previous hands-on sessions to:
plot the comparison between the numerical temperature and the analytical temperature.
ensure that our VTKRecorder is also collecting and printing files for paraview.
start the simulation.
Example script¶
Please find a full script located in the examples folder
