Energy conservation in thermal and flow Engines

Asked by Zoheir Khademian on 2020-10-30

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'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))

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,
 stressMask = 7,
 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,
 FlowEngine(dead=1,label="flow",multithread=False),
 ThermalEngine(dead=1,label='thermal'),
 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
triax.stressMask=7
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
# add flow
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

flow.dead=0
thermal.dead=1

thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0
thermal.thermoMech=False
thermal.solidThermoMech = False
thermal.fluidThermoMech = False
thermal.advection=True
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
thermal.dead=0
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)
#-------------------------------------------------
from yade import plot

def history():
 plot.addData(
  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")
lines = infile.readlines()
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:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-30
Last reply:
2020-10-31

This question was reopened

Robert Caulk (rcaulk) said : #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

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.thermalBndFlux
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getBoundaryFlux

Robert Caulk (rcaulk) said : #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

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.thermalBndFlux

Zoheir Khademian (zkhademian) said : #3

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

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.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

Can you help with this problem?

Provide an answer of your own, or ask Zoheir Khademian for more information if necessary.

To post a message you must log in.