Extract a cylinder from a cloud of spheres

Asked by Othman Alshareedah on 2019-01-07

Hello all,

I'm new in Yade and in using python in general. I'm trying to generate a cylinder full of spheres that are packed to a specific porosity. Then I will use that cylinder in another script to apply uniaxial compression. So to achieve that, first I create a cloud of spheres in a box and then I apply triaxial compaction until I get the desired porosity. The issue is that, after I get the desired porosity, how can I extract (or trim/crop) a cylinder from the box that is full of compacted spheres??

I would greatly appreciate if anyone can help me to solve this problem. Below is my code.

Thank you
Othman

############################ Material properties #############################
sigmaIso=-1e3
frictangel=.52
targetp=0.7 #this is the targeted porosity
pervconc=O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = radians(frictangel), density=1920))

############################ generate loose packing #############################

from yade import pack, qt, plot, export
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.1,.1,.23),rMean=.0045,rRelFuzz=.015,periodic=True)

sp.toSimulation()
yade.qt.View()

############################ Triaxial compaction #############################
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   PeriTriaxController(label='triax',
      goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
      dynCell=True,maxStrainRate=(100,100,100),
      maxUnbalanced=.1,relStressTol=1e-3,
      doneHook='compactionFinished()'
   ),
   NewtonIntegrator(damping=.2),
   PyRunner(command='P()',iterPeriod=5000)
]
O.dt=.5*PWaveTimeStep()
## stop conditions

def P():

 global frictangel
 if utils.porosity()>targetp:
  frictangel = frictangel
  frictangel = .95* frictangel #decreasing friction angle
  setContactFriction(radians(frictangel))
  print ('porosity = ', utils.porosity(), 'friction angle', frictangel)
 if utils.porosity()<targetp:
  O.pause()
  print 'Finished'

def compactionFinished():
   # set the current cell configuration to be the reference one
   O.cell.trsf=Matrix3.Identity
   # change control type: keep constant confinement in x,y, 20% compression in z
   triax.goal=(sigmaIso,sigmaIso,-.2)
   triax.stressMask=3
   # allow faster deformation along x,y to better maintain stresses
   triax.maxStrainRate=(1.,1.,.1)
   # next time, call triaxFinished instead of compactionFinished
   triax.doneHook='triaxFinished()'
   # do not wait for stabilization before calling triaxFinished
   triax.maxUnbalanced=10

def triaxFinished():
   print 'Finished'
   O.pause()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2019-01-08
Last query:
2019-01-08
Last reply:
2019-01-08
Robert Caulk (rcaulk) said : #1

Welcome to Yade!

The Yade doc will be your best friend [1] (modules and class reference are particularly useful).

It seems you are looking for filterSpherePack [2]. If you want an example of how to compact a cube and then filter it, that is exactly what randomDensePack does [3] (but it doesn't reach a desired porosity in randomDensePack).

Cheers,

Robert

[1]https://yade-dem.org/doc/index-toctree.html#
[2]https://yade-dem.org/doc/yade.pack.html#yade.pack.filterSpherePack
[3]https://yade-dem.org/doc/_modules/yade/pack.html#randomDensePack

Othman Alshareedah (othman-sh) said : #2

Hi Robert,

Thanks for your answers. I tried randomDensePack for a cylinder but it only give a porosity of about .65 while I need porosity~.2-.3.

Also, I tried using filterSpherePack as below, but nothing appeared in my simulation viewer. Also I don't know if I'm filtering the original cloud or the compacted one. My goal is to crop a cylinder, delete the rest of the speheres and export the data to a txt file. Can you please help with that. Below is my code.

Thank you
Othman

############################ Material properties #############################
sigmaIso=-1e3
frictangel=.52
targetp=0.74 #this is the targeted porosity
pervconc=O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = radians(frictangel), density=1920))

############################ generate loose packing #############################

from yade import pack, qt, plot, export
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.4,.4,.6),rMean=.015,rRelFuzz=.015,periodic=True)

sp.toSimulation()
yade.qt.View()

############################ Triaxial compaction #############################
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   PeriTriaxController(label='triax',
      goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
      dynCell=True,maxStrainRate=(100,100,100),
      maxUnbalanced=.1,relStressTol=1e-3,
      doneHook='compactionFinished()'
   ),
   NewtonIntegrator(damping=.2),
   PyRunner(command='P()',iterPeriod=5000)
]
O.dt=.5*PWaveTimeStep()
## stop conditions

def P():

 global frictangel
 if utils.porosity()>targetp:
  frictangel = frictangel
  frictangel = .95* frictangel #decreasing friction angle
  setContactFriction(radians(frictangel))
  print ('porosity = ', utils.porosity(), 'friction angle', frictangel)
 if utils.porosity()<targetp:
  O.pause()
  print 'Finished'

def compactionFinished():
   # set the current cell configuration to be the reference one
   O.cell.trsf=Matrix3.Identity
   # change control type: keep constant confinement in x,y, 20% compression in z
   triax.goal=(sigmaIso,sigmaIso,-.2)
   triax.stressMask=3
   # allow faster deformation along x,y to better maintain stresses
   triax.maxStrainRate=(1.,1.,.1)
   # next time, call triaxFinished instead of compactionFinished
   triax.doneHook='triaxFinished()'
   # do not wait for stabilization before calling triaxFinished
   triax.maxUnbalanced=10

def triaxFinished():
   print 'Finished'
   O.pause()

############################ Extracting a cylinder #############################

ss=pack.filterSpherePack(pack.inCylinder((.05,.05,.05),(.15,.15,.2),.02), sp, returnSpherePack=True)
ss.toSimulation()

Best Robert Caulk (rcaulk) said : #3

>>I tried randomDensePack for a cylinder but it only give a porosity of about .65 while I need porosity~.2-.3

randomDensePack is not currently capable of returning a packing with a desired porosity.

>>Can you please help with that.

Following the method shown in [1], your code might want to look something like this:

## porosity reached above##
pred =pack.inCylinder((.05,.05,.05),(.15,.15,.2),.02)
sp=SpherePack()
sp.fromSimulation()
O.switchScene()
spFiltered = filterSpherePack(pred,sp,material=material,returnSpherePack=returnSpherePack)
spFiltered.toSimulation()
export.textExt('name.spheres','x_y_z_r')

[1]https://yade-dem.org/doc/_modules/yade/pack.html#randomDensePack

Othman Alshareedah (othman-sh) said : #4

Robert,

Thank you so much. Problem solved.