yade-users team mailing list archive

[Jiannan Wang <question692194@xxxxxxxxxxxxxxxxxxxxx>]

New question #692194 on Yade:
https://answers.launchpad.net/yade/+question/692194

Hello,

I'm trying to explore the thermal expansion with Yade, similar to the paper [1]. But while I heated up one side of the boundary: thermal.thermalBndCondValue=[343.15,0,0,0,0,0], the temperature of the body  went up beyond the boundary temperature. I'm wondering where did I go wrong? Thanks. Here is the MWE: 
###################

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)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(30),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=0.0015,rRelFuzz=0.333,num=1000)
sp.toSimulation(color=(0.752, 0.752, 0.752),material='spheres')

triax=TriaxialStressController(
	maxMultiplier=1.+2e4/young,
	finalMaxMultiplier=1.+2e3/young,
	thickness = 0,
	stressMask = 7,
	internalCompaction=True,
)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,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'),
	NewtonIntegrator(damping=0.5)
]

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1

O.dt=0.1e-5
O.dynDt=False

tri_pressure = 10000
triax.goal1=triax.goal2=triax.goal3=-tri_pressure
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<0.01 and abs(-tri_pressure-triax.meanStress)/tri_pressure<0.001:
    break

triax.internalCompaction=False
for b in O.bodies:
	if isinstance(b.shape, Sphere):
		b.dynamic=False

flow.dead=0
flow.defTolerance=0.003
flow.meshUpdateInterval=100
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
flow.useSolver=4
flow.permeabilityFactor= 1
flow.viscosity= 0.001
flow.bndCondIsPressure=[0,0,0,0,0,0]
flow.bndCondValue=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.debug=False
flow.fluidBulkModulus=2.2e11
flow.fluidRho = 997
flow.fluidCp = 4181.7
flow.bndCondIsTemperature=[0,0,0,0,0,0]
flow.thermalEngine=True
flow.thermalBndCondValue=[0,0,0,0,0,0]
flow.tZero=333.15
flow.pZero=0

thermal.dead=0
thermal.debug=False
thermal.fluidConduction=True
thermal.ignoreFictiousConduction=True
thermal.conduction=True
thermal.thermoMech=True
thermal.solidThermoMech = True
thermal.fluidThermoMech = True
thermal.advection=True
thermal.bndCondIsTemperature=[1,0,0,0,0,0]
thermal.thermalBndCondValue=[343.15,0,0,0,0,0]
thermal.fluidK = 0.6069
thermal.fluidConductionAreaFactor=1.
thermal.uniformReynolds=10
thermal.particleT0 = 333.15
thermal.particleDensity=2600.
thermal.particleK = 2.
thermal.particleCp = 710.
thermal.tsSafetyFactor=0
thermal.useKernMethod=False
thermal.useHertzMethod=True

timing.reset()
O.dt=1e-5
O.run(1,1)

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(0.025,0.025,0.025)

from yade import plot

def history():
	plot.addData(
		ftemp1=flow.getPoreTemperature((0.025,0.025,0.025)),
		p=flow.getPorePressure((0.025,0.025,0.025)),
		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')]

plot.plots={'t':(('ftemp1','k-'),('bodyOfIntTemp','r-')), 't ':('p')} #

plot.plot(subPlots=False)
O.run()
################

Best regards
Jiannan



[1]: https://mercurylab.co.uk/dem8/wp-content/uploads/sites/4/2019/07/217.pdf

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