yade-users team mailing list archive

[Robert Caulk <question701642@xxxxxxxxxxxxxxxxxxxxx>]

Question #701642 on Yade changed:
https://answers.launchpad.net/yade/+question/701642

    Status: Open => Answered

Robert Caulk proposed the following answer:
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 <
question701642@xxxxxxxxxxxxxxxxxxxxx> 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp
>

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.