Question #698948 on Yade changed: https://answers.launchpad.net/yade/+question/698948
Zoheir Khademian posted a new comment: Thanks Robert! I went through a trial and error process using a more traditional BC as you suggested: low.bndCondIsTemperature=[0,0,0,0,1,1] flow.thermalBndCondValue=[0,0,0,0,25,45] thermal.bndCondIsTemperature=[0,0,0,0,1,1] thermal.thermalBndCondValue=[0,0,0,0,25,45] and avoid mesh updating by meshUpdateInterval=-1 - First, I turned off the compaction and ran the model. Temperature ranges seemed to be within expected range 25-45 K. Turning on the compaction made the model unstable in a way that solid and fluid temp goes above the assigned boundary condition (45 K). - I tried to minimize the particle-pore and pore-pore resistivities by assigning very small numbers to thermal.fluidConductionAreaFactor[1] and thermal.fluidK[2]. I also made sure that the distance between cells wont be zero by assigning thermal.minimumFluidCondDist[3] of 0.001. I tried different numbers but the particle temp still goes up to 1e3 very quickly while cell temp stays below 45 K. - I also tried avoiding fictious cells by thermal.ignoreFictiousConduction[4]=True but still no chance. I also tried using HertzMethod [5] for area calculation. The cell temp was in the range but particle temp went to negative. - I also limited the Reynolds (thermal.uniformReynolds[6]) to 10 to make sure Nusselt number calc is not the reason. [1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.fluidConductionAreaFactor [2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.fluidK [3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.minimumFluidCondDist [4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.ignoreFictiousConduction [5] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.useHertzMethod [6] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ThermalEngine.uniformReynolds Here is the MWE for reference: ########################################### from yade import pack, ymport, plot, utils, export, timing import numpy as np young=5e6 mn,mx=Vector3(0,0,0),Vector3(0.05,0.05,0.05) identifier = '-thm_coupling' if not os.path.exists('VTK'+identifier): os.mkdir('VTK'+identifier) if not os.path.exists('txt'+identifier): os.mkdir('txt'+identifier) 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=.333,num=200,seed=11) sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres') triax=TriaxialStressController( maxMultiplier=0, finalMaxMultiplier=0, thickness = 0, goal1=-50000, goal2=-50000, goal3=-50000, max_vel=0.1, stressMask = 7, internalCompaction=False, dead=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'), VTKRecorder(iterPeriod=500,fileName='VTK'+identifier+'/spheres-',recorders=['spheres','thermal','intr'],dead=0,label='VTKrec'), NewtonIntegrator(damping=0.5) ] O.dt=PWaveTimeStep()*.3 O.step() ss2sc.interactionDetectionFactor=-1 is2aabb.aabbEnlargeFactor=-1 triax.dead=False while 1: O.run(1000, True) unb=unbalancedForce() print('unbalanced force:',unb,' mean stress: ',triax.meanStress) if unb<0.1 and abs(triax.goal3-triax.meanStress)/-triax.goal3<0.001: break for b in O.bodies: if isinstance(b.shape,Sphere): b.dynamic=False triax.internalCompaction=False maxY=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) maxX=max([b.state.pos[0]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) minY=min([b.state.pos[1]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) minX=min([b.state.pos[0]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) minZ=min([b.state.pos[2]-b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) maxZ=max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]) print(minX,minY,minZ) print(maxX,maxY,maxZ) dz=maxZ-minZ dy=maxY-minY dx=maxX-minX 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 flow.viscosity= 0.001 flow.decoupleForces = False flow.bndCondIsPressure=[0,0,0,0,1,1] flow.bndCondValue=[0,0,0,0,0,10] ## Thermal Stuff flow.bndCondIsTemperature=[0,0,0,0,1,1] flow.thermalEngine=True flow.thermalBndCondValue=[0,0,0,0,25,45] flow.tZero=25 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,1,1] thermal.thermalBndCondValue=[0,0,0,0,25,45] thermal.fluidK = 1e-6#0.650 thermal.fluidBeta = 2e-5 # 0.0002 thermal.particleT0 = 25 thermal.particleK = 2.0 thermal.particleCp = 710 thermal.particleAlpha = 3.0e-5 thermal.particleDensity = 2700 thermal.tsSafetyFactor = 0 #0.01 thermal.uniformReynolds =10 thermal.minimumThermalCondDist=0 thermal.minimumFluidCondDist=1e-3 thermal.fluidConductionAreaFactor=1e-6 #thermal.ignoreFictiousConduction=True timing.reset() O.dt=1e-3 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((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.)) from yade import plot def MaxTempInModel(): MaxFluidTemp=0 for i in range(0,flow.nCells()): MaxFluidTemp=max(MaxFluidTemp,flow.getCellTemperature(i)) MaxTemp=0 for b in O.bodies: if isinstance(b.shape,Sphere): MaxTemp=max(MaxTemp,b.state.temp) print("######### MAX SOLID TEMPERATURE IS:",MaxTemp,"##########") print("######### MAX FLUID TEMPERATURE IS:",MaxFluidTemp,"##########") O.engines=O.engines+[PyRunner(iterPeriod=100,command='MaxTempInModel()',label='MaxTempInModel')] def history(): plot.addData( ftemp1=flow.getPoreTemperature(((maxZ-dz/2.),(maxZ-dz/2.),(maxZ-dz/2.))), t=O.time, i = O.iter, bodyOfIntTemp = O.bodies[bodyOfInterest.id].state.temp, ) #plot.saveDataTxt('Record.txt',vars=('t','i','p','ftemp1','bodyOfIntTemp')) O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')] def pressureField(): flow.saveVtk('VTKPressure'+identifier+'/',withBoundaries=True) O.engines=O.engines+[PyRunner(iterPeriod=500,command='pressureField()')] 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=20)] ColorScaler() O.run(2000,0) -- 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