# Energy conservation in thermal and flow Engines

I am evaluating the conservation of energy within a cubic box of spheres with voids filled by a fluid. There is no flow in the model and thermal boundary condition forces transfer of heat at only one boundary. The initial temperature of solid and fluid is 45 K and boundary temperature of both solid and fluid is fixed at 25 K. I calculate the change in the internal energy of solid and fluid phases and compare it with the total thermal flux * time.

The difference between them implies "thermalBndFlux" excludes calculation of flux of energy by the fluid cells. Is this true?
If so, how can I calculate the boundary thermal flux through fluid? I have difficulty finding the throat index relating two neighboring cells so that I can use "getDiffusionCoeff" and cell temperatures to estimate thermal flux between the cells.

Thanks-Zoheir

Here is a MWE:

from yade import pack, ymport, plot, utils, export, timing
import numpy as np
t0=45
young=5e6

mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1)

O.materials.append(FrictMat(young=young*100,poisson=0.5,frictionAngle=0,density=2600,label='walls'))

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

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=0.333,num=400,seed=11)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = 0,
internalCompaction=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,
NewtonIntegrator(damping=0.5)
]

O.step()

ss2sc.interactionDetectionFactor=-1# double check
is2aabb.aabbEnlargeFactor=-1#double check

tri_pressure = 1000
triax.goal1=triax.goal2=triax.goal3=-tri_pressure
while 1:
O.run(1000, True)
unb=unbalancedForce()
print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
if unb<0.1 and abs(-tri_pressure-triax.meanStress)/tri_pressure<0.001:
break

triax.internalCompaction=False
#-------------------------------------------------
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.dynamic=False # mechanically static

#-------------------------------------------------
flow.debug=False
flow.permeabilityMap = False
#flow.pZero = 10
flow.meshUpdateInterval=1
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=-1e-5#for calculating or assigning conductivity; negetive allows you to assign conductivity
flow.viscosity= 0.001
flow.decoupleForces = False# double check
flow.bndCondIsPressure=[0,0,0,0,0,0]
flow.bndCondValue=[0,0,0,0,0,0]

## Thermal Stuff
flow.bndCondIsTemperature=[0,0,0,0,0,1]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,25]
flow.tZero=t0

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.useKernMethod=False
thermal.bndCondIsTemperature=[0,0,0,0,0,1]
thermal.thermalBndCondValue=[0,0,0,0,0,25]
thermal.fluidK = 0.650
thermal.fluidBeta = 2e-5 # 0.0002
thermal.particleT0 = t0
thermal.particleK = 2.0
thermal.particleCp = 710
thermal.particleAlpha = 3.0e-5
thermal.particleDensity = 2600.
thermal.tsSafetyFactor = 0 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0
#-------------------------------------------------

timing.reset()
O.dt=1.
O.dynDt=False
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(0.025,0.025,0.025)
#-------------------------------------------------

def history():
ftemp1=flow.getPoreTemperature((0.025,0.025,0.025)),
p=flow.getPorePressure((0.025,0.025,0.025)),
t=O.time,
i = O.iter,
bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp,
)

O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]

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=100)]

f= open("RecordedDataTwoPhase.txt", 'w+')
f.close()
def Recorder():
Q=(abs(thermal.thermalBndFlux[0])*O.dt+abs(thermal.thermalBndFlux[1])*O.dt+abs(thermal.thermalBndFlux[2])*O.dt+abs(thermal.thermalBndFlux[3])*O.dt+abs(thermal.thermalBndFlux[4])*O.dt+abs(thermal.thermalBndFlux[5])*O.dt)
f= open("RecordedDataTwoPhase.txt", 'a')
f.write('%f %f' % (Q,O.iter))
f.write('\n')
f.close()

O.engines=O.engines+[PyRunner(iterPeriod=1,command='Recorder()',label='fluxSum')]
#-------------------------------------------------
#Finding Initial IE
cellsHit = []
numPoints = 40
mnx=0.
mxx=0.1
mny=0.
mxy=0.1
mnz=0.
mxz=0.1
xs = np.linspace(0.95*mnx,1.05*mxx,numPoints)
ys = np.linspace(0.95*mny,1.05*mxy,numPoints)
zs = np.linspace(0.95*mnz,1.05*mxz,numPoints)

v = np.array([0,0,0])

O.step()
O.step()

PM1=0.;VV=0;PM2=0.;tt=0
for b in O.bodies:
if isinstance(b.shape, Sphere):
PM1+=(t0-b.state.temp)*b.state.mass*thermal.particleCp

cellId=0
for x,y,z in itertools.product(xs, ys, zs):
cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
if cellId in cellsHit: continue #this prevents counting a cell multiple times
cellsHit.append(cellId)
VV+=flow.getCellVolume((x,y,z))
PM2+=(1./flow.getCellInvVoidVolume(cellId) ) *flow.fluidRho*flow.fluidCp* (t0-flow.getPoreTemperature((x,y,z)))
#print("Initial Solid IE is",PM1, "Initial Fluid IE is",PM2,"Total Volume is:",VV)

#-------------------------------------------------
O.run(1000,wait=True)
#-------------------------------------------------
dQsolid=0.
for b in O.bodies:
if isinstance(b.shape, Sphere):
dQsolid+=(t0-b.state.temp)*b.state.mass*thermal.particleCp

VV=0;VVv=0.; CellNo=0; dQfluid=0.;cellsHit=[]
for x,y,z in itertools.product(xs, ys, zs):
cellId = flow.getCell(x,y,z) #getCell return cell id for a point position xyz
if cellId in cellsHit: continue #this prevents counting a cell multiple times
cellsHit.append(cellId)
if (t0-flow.getPoreTemperature((x,y,z)))>0.:
dQfluid+=(1./flow.getCellInvVoidVolume(cellId)) *flow.fluidRho*flow.fluidCp* (t0-flow.getPoreTemperature((x,y,z)))

dQtotal=dQsolid+dQfluid-PM1-PM2

FluxTotal=0.
infile = open("RecordedDataTwoPhase.txt","r")
infile.close()
ret=[]
for line in lines:
data = line.split()
FluxTotal+=float(data[0])
print("Fluid Internal Energy Change is", dQfluid-PM2,", Solid Internal Energy Change is:",dQsolid-PM1,", Total Internal Energy Change is",dQtotal)
print("Total flux is",FluxTotal, ", and error is",(FluxTotal-dQtotal)/dQtotal)

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-10-30
2020-10-31

## This question was reopened

 Robert Caulk (rcaulk) said on 2020-10-31: #1

Hello,

Nice question. In the future, please try to link to the Yade documentation when you reference functions such as thermalBndFlux and getDiffusionCoeff etc. Similar to how I link to functions below. It makes the discussion easier to answer and easier to follow.

>The difference between them implies "thermalBndFlux" excludes calculation of flux of energy by the fluid cells. Is this true?

Yes, that is true. thermalBndFlux [1] only computes the heat flux for the solid boundary.

>If so, how can I calculate the boundary thermal flux through fluid?

You can get the fluid boundary flux using flow.getBoundaryFlux [2]. So you should be able to get the heat flux from:

Q = flow.getBoundaryFlux[5] * fluidCp * fluidRho * O.dt * flow.thermalBndCndValue[5]

Let me know if it works.

Cheers,

Robert

 Robert Caulk (rcaulk) said on 2020-10-31: #2

I am realizing you have no advection in your model, and my proposed solution will only account for advective heat flux in the fluid phase.

Instead, thermal.thermalBndFlux [1] may not be including the heat flux due to fluid conduction. You should be able to verify this by running your model with thermal.fluidConduction=False. If my suspicion is correct, you should end up with consistent numbers.

As for adding the fluidConduction flux to the total thermalBndFlux, it will need to be done through the source, I will have a look today and try to make a merge request for you.

Cheers,

Robert

Hi Robert,

Thanks for the response.

>You should be able to verify this by running your model with thermal.fluidConduction=False. If my suspicion is correct, you should end up with consistent numbers.

Yes, I tried it and the numbers are the same, meaning the fluid conduction is not implemented in the thermal flux calculation in thermalBndFlux [1].

>As for adding the fluidConduction flux to the total thermalBndFlux, it will need to be done through the source, I will have a look today and try to make a merge request for you.

Thanks and looking forward to it

Zoheir