yade-users team mailing list archive

[Lin 1999 <question701642@xxxxxxxxxxxxxxxxxxxxx>]

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.