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'))
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)