# 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
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)*O.dt+abs(thermal.thermalBndFlux)*O.dt+abs(thermal.thermalBndFlux)*O.dt+abs(thermal.thermalBndFlux)*O.dt+abs(thermal.thermalBndFlux)*O.dt+abs(thermal.thermalBndFlux)*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)
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  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 . So you should be able to get the heat flux from:

Q = flow.getBoundaryFlux * fluidCp * fluidRho * O.dt * flow.thermalBndCndValue

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  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 .

>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