some questions about conductive heat flux between DEM particles

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,

What is "bigger?" Please quantify your observations i.e. "when I set
*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 [1]. 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 [1] to see if that is
increasing or decreasing with your boundary flux estimate.

Cheers,

Robert

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/Thermal.cpp#L490

On Tue, May 3, 2022 at 1:56 PM Lin 1999 <
<email address hidden>> wrote:

> New question #701642 on Yade:
> https://answers.launchpad.net/yade/+question/701642
>
> 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[1].
> 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 yade import pack,plot,timing,os
> 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
> O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(frictionangle),density=density,label='spheres'))
>
>
> O.materials.append(FrictMat(young=wallyoung,poisson=wallpoisson,frictionAngle=radians(wallfrictionangle),density=walldensity,label='walls'))
>
> 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[2],
> wall_top_id = wallIds[3],
> wall_left_id = wallIds[0],
> wall_right_id = wallIds[1],
> wall_back_id = wallIds[4],
> wall_front_id = wallIds[5],
> thickness=0.00001,
> internalCompaction=False,
> stressMask=7,
> goal1=compress,
> goal2=compress,
> goal3=compress,
> #dead=True,
> )
>
> newton=NewtonIntegrator(damping=damp)
>
> ThermalEngine = ThermalEngine(dead=1,label='thermal');
>
> 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,
> FlowEngine(dead=1,label="flow",multithread=False), #
> ThermalEngine,
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> #quick
> newton
> ]
>
> triax.dead=False
> while 1:
> O.run(1000, True)
> unb=unbalancedForce()
> s1=triax.stress(triax.wall_right_id)[0]
> s2=triax.stress(triax.wall_top_id)[1]
> s3=triax.stress(triax.wall_front_id)[2]
> volume =
> (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
> 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)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/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)[0],'s22',-triax.stress(triax.wall_top_id)[1],'s33',-triax.stress(triax.wall_front_id)[2],'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
>
> flow.dead=0
>
> thermal.dead=0
> thermal.debug=False
> thermal.conduction=True
> thermal.thermoMech=False
> thermal.advection=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()
> flow.dead=1
>
> 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[4]
> disx=triax.width
> disy=triax.height
> disz=triax.depth
>
> keff=(thermal.thermalBndFlux[4])*float(disz)/(float(disx)*float(disy)*deltaT)
> if thermal.thermalBndFlux[4]+thermal.thermalBndFlux[5] < -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')]
>
> def addPlotData():
> volume =
> (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*(O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2])
> Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if
> isinstance(b.shape,Sphere))/volume
> plot.addData(
> t=O.time,
> i = O.iter,
> porosity=Porosity,
> flux4=thermal.thermalBndFlux[4],
> flux5=thermal.thermalBndFlux[5],
> disx=triax.width,
> disy=triax.height,
> disz=triax.depth,
> )
>
> O.engines=O.engines+[PyRunner(iterPeriod=5000,command='addPlotData()',label='recorder')]
> plot.plots={'t':(('flux4','k-'),('flux5','r-'))}
>
> plot.plot()
>
> O.run()
>
> #pylab.show()
>
>
>
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

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

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

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

Changing status to avoid janitor cleanup.