Particles in the middle of the cylinder disappeared after adding the thermal engine

Hi Ziyu,

I tried to run your example but it requires loading sphere packing from an external file:

>OSError: packing-cylinder.spheres not found.

Cheers,
Karol

Hi Karol,

Sorry for ignore this error,in fact,the packing-cylinder.spheres is generated by following script(I don't know how to upload external files, so I can only provide this method...)

######
from yade import pack,export,ymport

name='cylinder'

bottom,top=Vector3(0,0,0),Vector3(0,0,0.05)
radius=0.0125
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))

export.text('packing-'+name+'.spheres')
######

Thank you for your kindly help again!

Hi,

I don't know whether such a coupling can be used for your problem. I noticed that the example from tutorial [1] uses extremely high density. Maybe someone else could provide a better answer...

Cheers,
Karol

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/ThermalEngine/thermoHydroMechanical_coupling.py

Hi,Karol,

I noticed the density problem in the example you mentioned. Although I thought the density value in the example was unreasonable, I changed the density in my script so that it was consistent with the example. However, this did not solve the problem. The particles in the middle still disappeared..

Anyway, thank you for your attention and kindly help!

Best thing for us to help you is if you could tell us "when I change
*this* in the example thermal script, the particles disappear in the
middle."If you cannot do this, I cannot help you.

On Sun, Sep 4, 2022 at 2:47 PM Ziyu Wang <
<email address hidden>> wrote:

> Question #703024 on Yade changed:
> https://answers.launchpad.net/yade/+question/703024
>
> Ziyu Wang gave more information on the question:
> Hi,Karol,
>
> I noticed the density problem in the example you mentioned. Although I
> thought the density value in the example was unreasonable, I changed the
> density in my script so that it was consistent with the example.
> However, this did not solve the problem. The particles in the middle
> still disappeared..
>
> Anyway, thank you for your attention and kindly help!
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Robert,

In fact, I think there may be a misunderstanding between us, or I misunderstood your meaning(Sorry for that..)

My description in 1# is to add the Thermalengine and related parameters on the basis of the fluid-solid coupling code I wrote myself. The parameter values refer to the examples, but do not involve the modification of the example script.

I'll put the fluid solid coupling code that was successfully run before below. I don't know if it's what you need..

Best wishes.

######
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys
Sy=4e6
damp=0.4

key='_triax_base_'
young=66.2e9
name='JCFPM_triax'
poisson=0.522
name='cylinder'
preStress=-90e6
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)

spheres=O.bodies.append(ymport.text('packing-'+name+'.spheres',scale=1,shift=Vector3(0,0,0),material='sphere',color=(0.25,0.25,0.25)))

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=.003*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=1000
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()

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

By the way,the packing file is generated by the method in #2.Thanks again!

This question was expired because it remained in the 'Open' state without activity for the last 15 days.