← Back to team overview

yade-users team mailing list archive

[Question #703024]: Particles in the middle of the cylinder disappeared after adding the thermal engine

 

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

Hi,everyone!

I have completed a relatively good fluid-solid coupling simulation before. Recently, I wanted to add temperature effect, so I added thermalengine to the original script. However, there was a problem as described in the title: I opened the 3D view and found that the particles in the middle of the cylinder had disappeared, so I could not get the desired results.

I do not know what happened,thank you for any possible help or suggestion..

------------------code-------------------------
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys

#readParamsFromTable(Sy=4e6)
#from yade.params import table
#global Sy
Sy=4e6
damp=0.4

key='_triax_base_'
young=66.2e9
name='JCFPM_triax'
poisson=0.522
name='cylinder'
preStress=-60e6
strainRate = -10
OUT=str(Sy)+'_JCFPM_triax'
r1=0.005
r2=0.008

nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05
A=0.25*3.14*width*width

allx,ally,allz,r=np.genfromtxt('packing-cylinder.spheres', unpack=True)
mnx=min(allx)*1.01 #same as aabbExtrema, get the xyz of the lowest corner, multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*1.01
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, get the xyz of the top corner, multiply by factor to increase it to let the wall surround all the spheres.
mxy=max(ally)*1.01
mxz=max(allz)*1.01

O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=8e6,cohesion=40e6,frictionAngle=radians(25),label='sphere'))
O.materials.append(JCFpmMat(type=1,label='wall1',young=young,poisson=poisson,frictionAngle=radians(25)))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall2'))

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0,material='wall2')
wallIds=O.bodies.append(walls)

bottom,top=Vector3(0,0,0),Vector3(0,0,0.05)
sp=pack.SpherePack()
pred=pack.inCylinder(bottom,top,radius)
sp=pack.randomDensePack(pred,radius=0.0005,rRelFuzz=0,returnSpherePack=True,memoizeDb='/tmp/triax.sqlite')
spheres=sp.toSimulation(color=(0.9,0.8,0.6),material='sphere')

bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
 s.shape.color = (1,0,0)
 s.state.blockedDOFs = 'xyzXYZ'
 #s.state.vel = (0,0,-vel)
for s in top:
 s.shape.color = Vector3(0,1,0)
 s.state.blockedDOFs = 'xyzXYZ'
 s.state.vel = (0,0,vel)

facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
	for h in range(nh):
		v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
		v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
		v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
		f1 = facet((v1,v2,v3),color=(0,0,1),material='wall1')
		f2 = facet((v1,v3,v4),color=(0,0,1),material='wall1')
		facets.extend((f1,f2))

O.bodies.append(facets)
mass = O.bodies[7].state.mass
for f in facets:
	f.state.mass = mass
	f.state.blockedDOFs = 'XYZz'



def lateral():
 elatTot=0.0
 nTot=0
 for b in O.bodies:
 	x=b.state.refPos[0]
 	y=b.state.refPos[1]
 	d=sqrt(pow(x,2)+pow(y,2))
 	if d > r1 and d < r2:
 		b.shape.color=(0.6,0.5,0.15)
 		xnew=b.state.pos[0]
 		ynew=b.state.pos[1]
 		dnew=sqrt(pow(xnew,2)+pow(ynew,2))
 		elat=(dnew-d)/d
 		elatTot+=elat
 		nTot+=1
 elat_avg=elatTot/nTot
 return elat_avg

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.0125,0.0125,0.025)

def addForces():
	for f in facets:
		n = f.shape.normal
		a = f.shape.area
		O.forces.addF(f.id,preStress*a*n)

def stopIfDamaged(maxEps=0.001):
	extremum = max(abs(s) for s in plot.data['s'])
	s = abs(plot.data['s'][-1])
	e = abs(plot.data['e'][-1])
	if O.iter < 1000 or e < maxEps:
		return
	if abs(s)/abs(extremum) < 0.5 :
		print('Simulation finished')
		presentcohesive_count = 0
		for i in O.interactions:
        		if hasattr(i.phys, 'isCohesive'):
            			if i.phys.isCohesive == True:
                			presentcohesive_count+=1
		print('the number of cohesive bond now is:',presentcohesive_count)
		print('Max stress and strain:',extremum,e)
		O.wait()

O.dt=.01*utils.PWaveTimeStep()
newton=NewtonIntegrator(damping=damp)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
	),

	PyRunner(iterPeriod=1,command="addForces()"),
	FlowEngine(dead=1,label="flow"),
	ThermalEngine(dead=1,label='thermal'),
	newton,
	PyRunner(command='plotAddData()',iterPeriod=1000),
	PyRunner(iterPeriod=1000,command='stopIfDamaged()'),

]
def plotAddData():
	elat_avg=lateral()
	Qin = flow.getBoundaryFlux(4)
	perm = abs(Qin) * flow.viscosity * height / (A * Sy)
	f1 = sum(O.forces.f(b.id)[2] for b in top)
	f2 = sum(O.forces.f(b.id)[2] for b in bot)
	f = .5*(f2-f1)
	s = f/(pi*.25*width*width) 
	e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
	plot.addData(
		elateral = elat_avg,
		i = O.iter,
		s = -s,
		e = -e,
		tc = interactionLaw.nbTensCracks,
		sc = interactionLaw.nbShearCracks,
		permeability = perm
	)
	plot.saveDataTxt(OUT)
O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
cohesiveCount = 0
for i in O.interactions:
	if hasattr(i.phys, 'isCohesive'):
		if i.phys.isCohesive == True:
			cohesiveCount+=1
print('the origin total number of cohesive bond is:',cohesiveCount)

flow.debug=False
flow.permeabilityMap = False
flow.defTolerance=-1
flow.meshUpdateInterval=-1
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=1
flow.viscosity= 0.001
flow.decoupleForces =  False
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,Sy]
flow.dead=0
flow.emulateAction()

flow.bndCondIsTemperature=[1,1,1,1,1,1] 
flow.thermalEngine=True
flow.thermalBndCondValue=[100,100,100,100,100,100] 
flow.tZero=25
flow.tZero=25

thermal.dead=1
thermal.conduction=True
thermal.fluidConduction=True
thermal.debug=0

thermal.thermoMech=True
thermal.solidThermoMech = True
thermal.fluidThermoMech = True

thermal.advection=True
thermal.useKernMethod=1
thermal.bndCondIsTemperature=[1,1,1,1,1,1]
thermal.thermalBndCondValue=[25,25,25,25,25,25]
thermal.fluidK = 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 = 2640
thermal.tsSafetyFactor = 0.8 #0.01
thermal.uniformReynolds =10
thermal.minimumThermalCondDist=0

thermal.dead=0

plot.plots = { 'e':('s',), 'elateral':('s'),}
plot.plot()
O.run()
########################

Best wishes!

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