# some questions about conductive heat flux between DEM particles

Hello!
I am trying to use the thermal conduction scheme in packings of spherical particles. The top temperature of spheres are 5K and the bottom are 1K. When running the model, I can get the heat flux of the solid boundary by using thermalBndFlux.
However, when I increase the internal friction angle of the spheres, I find that the heat flux in heat balance at the solid boundary is bigger. I can't understand it. Actually I want to get the effective thermal conductivity of the 3D packed particles system using Fourier's law.

The principle is: keff=q*L/(A*deltaT)
keff(W/mk): Effective thermal conductivity;
q(W): The total heat flow out of the boundary plane in steady-state;
L(m): Depth of specimen;
A(m2): Cross sectional area of specimen;
deltaT(K): Temperature difference between upper and lower boundaries.

Under normal circumstances, the larger the porosity of the sample, the smaller the effective thermal conductivity, but I got the opposite result. Did I make any mistakes in the simulation?

Thanks
Lin

Here is a MWE:
from __future__ import print_function
from numpy import arange
import matplotlib;matplotlib.rc('axes',grid=True)
import itertools
import random
import numpy as np
import shutil
#import pylab

# Particles
density=2650
poisson=0.31
young=29e9
damp=0.25
number=2000
frictionangle=0
thermalCond = 3. #W/(mK)
heatCap = 2130. #J(kg K)
t0 = 1. #K
deltaT=4.0 #k
# Walls
walldensity=0
wallfrictionangle=0
wallpoisson=0.5
wallyoung=30000e9

# Compress
compress=-50000 #pa

# Corners of the initial packing
mn,mx=Vector3(0,0,0),Vector3(0.006,0.006,0.006)

# Materials

wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=0.00001,material='walls')) #

#PSD
psdSizes,psdCumm=[0.00025,0.0004,0.0006,0.00075],[0,0.15,0.85,1.0]
#pylab.plot(psdSizes,psdCumm,label='precribed mass PSD')

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.33,number,False,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
#pylab.plot(*sp.psd(bins=30,mass=True),label='PSD of (free) %d random spheres'%len(sp))
#pylab.legend()
sp.toSimulation(material='spheres')

O.trackEnergy=True

triax=TriaxialStressController(
wall_bottom_id = wallIds,
wall_top_id = wallIds,
wall_left_id = wallIds,
wall_right_id = wallIds,
wall_back_id = wallIds,
wall_front_id = wallIds,
thickness=0.00001,
internalCompaction=False,
goal1=compress,
goal2=compress,
goal3=compress,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys()],
[Law2_ScGeom_MindlinPhys_Mindlin()]

),
triax,
ThermalEngine, GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), #quick
newton
]

while 1:
O.run(1000, True)
unb=unbalancedForce()
s1=triax.stress(triax.wall_right_id)
s2=triax.stress(triax.wall_top_id)
s3=triax.stress(triax.wall_front_id)
volume = (O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)
Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume
mStress=-(triax.stress(triax.wall_right_id)+triax.stress(triax.wall_top_id)+triax.stress(triax.wall_front_id))/3.0,
#disx=triax.width,
#disy=triax.height,
#disz=triax.depth,
#print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
#if unb<0.1 and abs(triax.goal3-triax.meanStress)/-triax.goal3<0.001:
print('unbalanced force',unb,'meanstress',mStress,'s11',-triax.stress(triax.wall_right_id),'s22',-triax.stress(triax.wall_top_id),'s33',-triax.stress(triax.wall_front_id),'porosity',Porosity,'V',volume)
if unb<0.05 and abs(compress-s1)/(-compress)<0.01 and abs(compress-s2)/(-compress)<0.01 and abs(compress-s3)/(-compress)<0.01:
#if Porosity < 0.45:
break
for o in O.bodies:
o.dynamic=False

thermal.debug=False
thermal.conduction=True
thermal.thermoMech=False
thermal.fluidThermoMech = False
thermal.solidThermoMech = False
thermal.fluidConduction= False
thermal.bndCondIsTemperature=[0,0,0,0,1,1]
thermal.thermalBndCondValue=[0,0,0,0,1,5]
thermal.tsSafetyFactor=0
thermal.particleDensity=2650
thermal.particleT0=t0
thermal.particleCp=heatCap
thermal.particleK=thermalCond
thermal.particleAlpha =0
#thermal.useKernMethod=False
thermal.useHertzMethod=True

timing.reset()
flow.updateTriangulation=True
O.dt=0.1e-3
O.dynDt=False
flow.emulateAction()

def ColorScaler():
for s in O.bodies:
s.shape.color=scalarOnColorScale(s.state.temp,1,5)

O.engines=O.engines+[PyRunner(command='ColorScaler()',iterPeriod=100)]

def checktemp():
q=thermal.thermalBndFlux
disx=triax.width
disy=triax.height
disz=triax.depth
keff=(thermal.thermalBndFlux)*float(disz)/(float(disx)*float(disy)*deltaT)
if thermal.thermalBndFlux+thermal.thermalBndFlux < -0.0005:
plot.saveDataTxt('noflow_Bndflux_fri0.txt.bz2')
else:
print(keff)
print("The test is finished!")
O.pause()
O.engines=O.engines+[PyRunner(iterPeriod=5000,command='checktemp()',label='check')]

volume = (O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)
Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if isinstance(b.shape,Sphere))/volume
t=O.time,
i = O.iter,
porosity=Porosity,
flux4=thermal.thermalBndFlux,
flux5=thermal.thermalBndFlux,
disx=triax.width,
disy=triax.height,
disz=triax.depth,
)
plot.plots={'t':(('flux4','k-'),('flux5','r-'))}

plot.plot()

O.run()

#pylab.show()

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Robert Caulk (rcaulk) said on 2022-05-03: #1

Hello,

>>I find that the heat flux in heat balance at the solid boundary is
bigger. Under normal circumstances, the larger the porosity of the sample,
the smaller the effective thermal conductivity,

*blank* friction angle I get *blank* porosity and *blank heatflux*. When I
increase the frictionangle by *blank* percent, the porosity changes by
*blank* percent and the heatflux increases by *blank* percent." This helps
you and us to understand your problem. Right now, I don't see any
indication that you are controlling porosity.

The heat flux through your specimen depends on the area of contact between
your particles. You have selected to use thermal.useHertzMethod() for
computing the area of contact, which depends on the normal force between
the spheres . The greater the contact force, the greater the area, the
greater the flux. In the end, you can check your total area of contact by
iterating on the interactions and repeating  to see if that is
increasing or decreasing with your boundary flux estimate.

Cheers,

Robert

On Tue, May 3, 2022 at 1:56 PM Lin 1999 <

> New question #701642 on Yade:
>
> Hello!
> I am trying to use the thermal conduction scheme in packings of spherical
> particles. The top temperature of spheres are 5K and the bottom are 1K.
> When running the model, I can get the heat flux of the solid boundary by
> using thermalBndFlux.
> However, when I increase the internal friction angle of the spheres, I
> find that the heat flux in heat balance at the solid boundary is bigger. I
> can't understand it. Actually I want to get the effective thermal
> conductivity of the 3D packed particles system using Fourier's law.
>
> The principle is: keff=q*L/(A*deltaT)
> keff(W/mk): Effective thermal conductivity;
> q(W): The total heat flow out of the boundary plane in steady-state;
> L(m): Depth of specimen;
> A(m2): Cross sectional area of specimen;
> deltaT(K): Temperature difference between upper and lower boundaries.
>
> Under normal circumstances, the larger the porosity of the sample, the
> smaller the effective thermal conductivity, but I got the opposite result.
> Did I make any mistakes in the simulation?
>
> Thanks
> Lin
>
> Here is a MWE:
> from __future__ import print_function
> from numpy import arange
> import matplotlib;matplotlib.rc('axes',grid=True)
> import itertools
> import random
> import numpy as np
> import shutil
> #import pylab
>
> # Particles
> density=2650
> poisson=0.31
> young=29e9
> damp=0.25
> number=2000
> frictionangle=0
> thermalCond = 3. #W/(mK)
> heatCap = 2130. #J(kg K)
> t0 = 1. #K
> deltaT=4.0 #k
> # Walls
> walldensity=0
> wallfrictionangle=0
> wallpoisson=0.5
> wallyoung=30000e9
>
> # Compress
> compress=-50000 #pa
>
> # Corners of the initial packing
> mn,mx=Vector3(0,0,0),Vector3(0.006,0.006,0.006)
>
> # Materials
>
>
>
> wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=0.00001,material='walls'))
> #
>
> #PSD
> psdSizes,psdCumm=[0.00025,0.0004,0.0006,0.00075],[0,0.15,0.85,1.0]
> #pylab.plot(psdSizes,psdCumm,label='precribed mass PSD')
>
> sp=pack.SpherePack()
>
> sp.makeCloud(mn,mx,-1,0.33,number,False,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
> #pylab.plot(*sp.psd(bins=30,mass=True),label='PSD of (free) %d random
> spheres'%len(sp))
> #pylab.legend()
> sp.toSimulation(material='spheres')
>
> O.trackEnergy=True
>
> triax=TriaxialStressController(
> wall_bottom_id = wallIds,
> wall_top_id = wallIds,
> wall_left_id = wallIds,
> wall_right_id = wallIds,
> wall_back_id = wallIds,
> wall_front_id = wallIds,
> thickness=0.00001,
> internalCompaction=False,
> goal1=compress,
> goal2=compress,
> goal3=compress,
> )
>
> newton=NewtonIntegrator(damping=damp)
>
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_MindlinPhys()],
> [Law2_ScGeom_MindlinPhys_Mindlin()]
>
> ),
> triax,
> ThermalEngine,
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> #quick
> newton
> ]
>
> while 1:
> O.run(1000, True)
> unb=unbalancedForce()
> s1=triax.stress(triax.wall_right_id)
> s2=triax.stress(triax.wall_top_id)
> s3=triax.stress(triax.wall_front_id)
> volume =
> (O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)
> Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if
> isinstance(b.shape,Sphere))/volume
>
> mStress=-(triax.stress(triax.wall_right_id)+triax.stress(triax.wall_top_id)+triax.stress(triax.wall_front_id))/3.0,
> #disx=triax.width,
> #disy=triax.height,
> #disz=triax.depth,
> #print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
> #if unb<0.1 and abs(triax.goal3-triax.meanStress)/-triax.goal3<0.001:
> print('unbalanced
> force',unb,'meanstress',mStress,'s11',-triax.stress(triax.wall_right_id),'s22',-triax.stress(triax.wall_top_id),'s33',-triax.stress(triax.wall_front_id),'porosity',Porosity,'V',volume)
> if unb<0.05 and abs(compress-s1)/(-compress)<0.01 and
> abs(compress-s2)/(-compress)<0.01 and abs(compress-s3)/(-compress)<0.01:
> #if Porosity < 0.45:
> break
> for o in O.bodies:
> o.dynamic=False
>
>
> thermal.debug=False
> thermal.conduction=True
> thermal.thermoMech=False
> thermal.fluidThermoMech = False
> thermal.solidThermoMech = False
> thermal.fluidConduction= False
> thermal.bndCondIsTemperature=[0,0,0,0,1,1]
> thermal.thermalBndCondValue=[0,0,0,0,1,5]
> thermal.tsSafetyFactor=0
> thermal.particleDensity=2650
> thermal.particleT0=t0
> thermal.particleCp=heatCap
> thermal.particleK=thermalCond
> thermal.particleAlpha =0
> #thermal.useKernMethod=False
> thermal.useHertzMethod=True
>
> timing.reset()
> flow.updateTriangulation=True
> O.dt=0.1e-3
> O.dynDt=False
> flow.emulateAction()
>
> def ColorScaler():
> for s in O.bodies:
> s.shape.color=scalarOnColorScale(s.state.temp,1,5)
>
> O.engines=O.engines+[PyRunner(command='ColorScaler()',iterPeriod=100)]
>
> def checktemp():
> q=thermal.thermalBndFlux
> disx=triax.width
> disy=triax.height
> disz=triax.depth
>
> keff=(thermal.thermalBndFlux)*float(disz)/(float(disx)*float(disy)*deltaT)
> if thermal.thermalBndFlux+thermal.thermalBndFlux < -0.0005:
> plot.saveDataTxt('noflow_Bndflux_fri0.txt.bz2')
> else:
> print(keff)
> print("The test is finished!")
> O.pause()
>
> O.engines=O.engines+[PyRunner(iterPeriod=5000,command='checktemp()',label='check')]
>
> volume =
> (O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)*(O.bodies.state.pos-O.bodies.state.pos)
> Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if
> isinstance(b.shape,Sphere))/volume
> t=O.time,
> i = O.iter,
> porosity=Porosity,
> flux4=thermal.thermalBndFlux,
> flux5=thermal.thermalBndFlux,
> disx=triax.width,
> disy=triax.height,
> disz=triax.depth,
> )
>
> plot.plots={'t':(('flux4','k-'),('flux5','r-'))}
>
> plot.plot()
>
> O.run()
>
> #pylab.show()
>
>
>
>
> --
>
> _______________________________________________
> Post to : <email address hidden>
>

 Revision history for this message Lin 1999 (sarcasticli) said on 2022-05-05: #2

I am Sorry, I didn't make my question clear. For example, when I set the frictionangle of particles to 0, the porosity of the specimen is 0.395. I use thermal.thermalBndFlux to get the heat flux through thermal boundary. When the specimen enters the thermal stable state, the heat flux through thermal boundary is about 0.080W. When I adjust the frictionangle to 45, the porosity is 0.447. At this time, the heat flux through thermal boundary is 0.095W. The heat flux is 18.8% higher than before. Usually, the larger the porosity of the sample, the worse the heat conduction, but this result does not appear in this simulation. This problem has been bothering me.

 Revision history for this message Launchpad Janitor (janitor) said on 2022-05-20: #3

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

 Revision history for this message Robert Caulk (rcaulk) said on 2022-05-20: #4

Changing status to avoid janitor cleanup.