Question #701642 on Yade changed: https://answers.launchpad.net/yade/+question/701642
Description changed to: 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() [1]https://yade- dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.thermalBndFlux -- 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@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp