Packing with two PSD

Asked by Othman Sh on 2021-03-11

Hello,

I am modeling compaction of two materials "modeled as spheres" with two particle size distribution. Material 1 sphere sizes is much smaller than material 2 as seen in the code below. When I run this code for ~20 seconds, spheres pass out the facets. If I comment this line "spFilter2.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)", the code works perfectly!.

I have asked similar question [1] and was told to reduce O.dt, and use large stiffness for the facets. I have tried using O.dt=.01*utils.PWaveTimeStep() and a young = 1e20 for Mat 3 (facet materials) and I ran the code for few minuets. The smaller size spheres still escape the facet.

 Can anyone please explain why this problem occurs and how to solve it?

Thanks,
Othman

[1] https://answers.launchpad.net/yade/+question/695929
Note: I am using yadedaily to run the code below.

-----------------------
from yade import pack, export
import numpy as np

Mat1=O.materials.append(FrictMat(young = 1e6, poisson = 0.45,frictionAngle = 0.00349066, density=2340,label='portlandite'))
Mat2=O.materials.append(FrictMat(young = 7.2e10, poisson = 0.17,frictionAngle = 0.558, density=2650,label='sand'))
Mat3=O.materials.append(FrictMat(young = 1e14, poisson = 0.1))
SG1=2.34 ##portlandite
SG2=2.65 ##ASTM quartz sand

##cylinder dimensions
radiuscyl=(500e-6/2)
heightcyl=610e-6
##center of cylinder
cx=0
cy=0
cz=0
##Initial cube dimensions ###
mnx=cx-(radiuscyl*1.1)
mny=cy-(radiuscyl*1.1)
mnz=0
mxx=cx+(radiuscyl*1.1)
mxy=cy+(radiuscyl*1.1)
mxz=heightcyl*1.1
volume=(pi*radiuscyl**2)*heightcyl

############################ spheres #############################
sp1=pack.SpherePack()

sizes1=1e-6*np.array([15,15.76,22.94,25]) #modified to avoid spheres escaping facet due to large size difference, dt issues
passing1=[0,0.75,0.9,1]

sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes1,psdCumm=passing1,num=5811)

#### cylinder extraction
pred=pack.inCylinder((cx,cy,cz),(cx,cy,heightcyl),radiuscyl)
spFilter1=filterSpherePack(pred,sp1,Material=Mat1, returnSpherePack=True)
print (len (spFilter1))

spFilter1.toSimulation(color=(0.533, 0.803, 0.780),material=Mat1)

mass1=utils.getSpheresMass()
#### sizes and distribution are from gradation curve of sand #######
sp2=pack.SpherePack()
sizes2=1e-6*np.array([149,150,300,425,600,850]) #Diameters of portlandite
passing2=[0,0.02,0.23,0.675,0.98,1]
sp2.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=51)

#### cylinder extraction
pred=pack.inCylinder((cx,cy,cz),(cx,cy,heightcyl),radiuscyl)
spFilter2=filterSpherePack(pred,sp2,Material=Mat2, returnSpherePack=True)
print (len (spFilter2))

spFilter2.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)
mass2=utils.getSpheresMass()-mass1

print ('mass 1 mass2 in g',mass1*1000, mass2*1000)
print ('mass2/mass1 ', mass2/mass1)

total_mass=utils.getSpheresMass()
############################ facets #############################

facets=geom.facetCylinder((cx,cy,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=150,wallMask=4,color=(0.156, 0.164, 0.152),material=Mat3)

cylinder=O.bodies.append(facets)
yade.qt.View()

##creating disks

d1=geom.facetCylinder((cx,cy,heightcyl),radiuscyl*0.99,0,segmentsNumber=150,wallMask=1,color=(0.156, 0.164, 0.152),material=Mat3)
d2=geom.facetCylinder((cy,cx,cz),radiuscyl*0.99,0,segmentsNumber=150,wallMask=1,color=(0.156, 0.164, 0.152),material=Mat3)

disk1IDs= O.bodies.append(d1)
disk2IDs= O.bodies.append(d2)

for i in disk1IDs:
 body= O.bodies[i]
 body.state.vel = (0,0,-1)

for n in disk2IDs:
 body= O.bodies[n]
 body.state.vel = (0,0,0)

############################# compaction #############################
O.dt=.1*utils.PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),
 PyRunner(iterPeriod=2500,command='force()',initRun=True),

]
O.run()

stress=[]
def force():
 f1= [O.forces.f(i)[2] for i in disk1IDs]
 f=np.mean(f1)
 s=f/(pi*radiuscyl**2)
 stress.append(s)
 thickness=(O.bodies[disk1IDs[1]].state.pos[2])-(O.bodies[disk2IDs[1]].state.pos[2])
 packing_density=utils.getSpheresVolume()/(thickness*pi*radiuscyl**2)
 print("stress, thickness, packing density ",s,thickness,packing_density)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Othman Sh
Solved:
2021-03-18
Last query:
2021-03-18
Last reply:
2021-03-16
Jan Stránský (honzik) said : #1

Hello,

> young = 1e20 for Mat 3
> Mat3=O.materials.append(FrictMat(young = 1e14, poisson = 0.1))

please be consistent.

> sp1=pack.SpherePack()
> sp2=pack.SpherePack()

use ONE SpherePack, calling makeCloud 2x.
The current way, the two sphere packs are independent, knowing nothing about each other, so the resulting bodies may be (most likely are) overlapping, which may cause some small particles "firing" from large particles, giving them high velocity, higher than needed to detect sphere-facet interaction.
Using one sphere pack, particles are guaranteed to be non-overlapping.

Also, using >0 thick boundary, e.g. box instead of facet, should help.

cheers
Jan

Othman Sh (othman-sh) said : #2

Hi Jan,

I followed your advise and used one SpherePack and 2 makeCloud to avoid overlapping. I also used box instead of cylinder. My new code is copied below. I reduced the O.dt to 0.005*utils.PWaveTimeStep() and the young modulus of the facets to 1e20. I have the following questions:

1- I still have the same issue of spheres passing through the facet. What could be the issue?
2- When you say "box instead of facet", are you referring to geom.facetBox or aabbWalls?
3- How to check if there is sphere overlap at the initial packing (i.e. iteration= 0)?

Thanks a lot,
Othman

---------------------
from yade import pack
import numpy as np

Mat1=O.materials.append(FrictMat(young = 1e6, poisson = 0.45,frictionAngle = 0.00349066, density=2340,label='portlandite'))
Mat2=O.materials.append(FrictMat(young = 7.2e10, poisson = 0.17,frictionAngle = 0.558, density=2650,label='sand')) #update properties for sand
Mat3=O.materials.append(FrictMat(young = 1e20, poisson = 0.1))

##cylinder dimensions
cube_dimensions=100e-6 #micron

##Initial cube dimensions ###
mnx=0
mny=0
mnz=0
mxx=cube_dimensions
mxy=cube_dimensions
mxz=cube_dimensions

box_center=(mxx/2,mxy/2,mxz/2)
box_extents=(mxx/2,mxy/2,mxz/2)
plate_center=(mxx/2,mxy/2,1.05*mxz/2)
plate_extents=(mxx/2,mxy/2,0.98*mxz/2)
############################ spheres #############################
sp1=pack.SpherePack()

sizes1=1e-6*np.array([9,9.16,15.76,22.94,25])
passing1=[0,0.5,0.75,0.9,1]
sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes1,psdCumm=passing1,num=5811)
sp1.toSimulation(color=(0.533, 0.803, 0.780),material=Mat1)

mass1=utils.getSpheresMass()

sizes2=1e-6*np.array([149,150,300,425,600,850]) #Diameters of portlandite
passing2=[0,0.02,0.23,0.675,0.98,1]
sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=51)
sp1.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)

mass2=utils.getSpheresMass()-mass1

print ('mass 1 mass2 in g',mass1*1000, mass2*1000)
print ('mass2/mass1 ', mass2/mass1)

total_mass=utils.getSpheresMass()

############################# box and plate #############################
yade.qt.View()
box=O.bodies.append(geom.facetBox(box_center,box_extents,wallMask=31,material=Mat3)) #open from one side
plate1IDs=O.bodies.append(geom.facetBox(plate_center,plate_extents,material=Mat3,wallMask=32))

for i in plate1IDs:
 body= O.bodies[i]
 body.state.vel = (0,0,-10)

############################## compaction #############################
O.dt=.005*utils.PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),

]
O.run()

Jan Stránský (honzik) said : #3

Hello,

> 1- I still have the same issue of spheres passing through the facet. What could be the issue?
> 2- When you say "box instead of facet", are you referring to geom.facetBox or aabbWalls?

I did not meant box as a shape of the boundary, but utils.box as a bodies forming the boundary.
Facets have zero thickness, so small particles can "jump" through the facet without contact detection.
This should not be the case for utils.box, which would have non-zero thickness and such "jump" would be more difficult.

> 3- How to check if there is sphere overlap at the initial packing (i.e. iteration= 0)?

O.step()
print(len(i for i in O.interactions))

cheers
Jan

Jan Stránský (honzik) said : #4

> sp1.toSimulation(color=(0.533, 0.803, 0.780),material=Mat1)
> ...
> sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=51)
> sp1.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)

do sp1.toSimulation just once, this way you have much more spheres in the packing then desired.
So you still have overlapping spheres and now the situation is even worse :-D

### a MWE showing this
sp = SpherePack([
    ((0,0,0),1),
    ((2,0,0),1),
])
print(len(sp)) # 2
sp.toSimulation()
print(len(O.bodies)) # 2
sp.toSimulation()
print(len(O.bodies)) # 4 !!!
###

cheers
Jan

Othman Sh (othman-sh) said : #5

Hi Jan,

I have updated my code to use utils.box instead of facet. My question is how to apply force on the spheres without using a "plate" facet? I am currently using geom.facetBox to make a flat surface "i.e. plate" and move it downward to apply force, but still smaller spheres can pass through it.

Also, how to make a packing with 2 psdSizes and psdCumm while using one packSpherepack() and one toSimulation to avoid initial overlaping of spheres?

Thanks,
Othman

---------
from yade import pack
import numpy as np

Mat1=O.materials.append(FrictMat(young = 1e6, poisson = 0.45,frictionAngle = 0.00349066, density=2340,label='portlandite'))
Mat2=O.materials.append(FrictMat(young = 7.2e10, poisson = 0.17,frictionAngle = 0.558, density=2650,label='sand'))
Mat3=O.materials.append(FrictMat(young = 1e20, poisson = 0.1))

cube_dimensions=100e-6 #micron

##Initial cube dimensions ###
mnx=0
mny=0
mnz=0
mxx=cube_dimensions
mxy=cube_dimensions
mxz=cube_dimensions

box_center=(mxx/2,mxy/2,mxz/2)
box_extents=(mxx/2,mxy/2,mxz/2)
plate_center=(mxx/2,mxy/2,1.05*mxz/2)
plate_extents=(mxx/2,mxy/2,0.98*mxz/2)
############################ spheres #############################
sp1=pack.SpherePack()

sizes1=1e-6*np.array([9,9.16,15.76,22.94,25])
passing1=[0,0.5,0.75,0.9,1]
sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes1,psdCumm=passing1,num=5811)
sp1.toSimulation(color=(0.533, 0.803, 0.780),material=Mat1)

sizes2=1e-6*np.array([149,150,300,425,600,850]) #Diameters of portlandite
passing2=[0,0.02,0.23,0.675,0.98,1]
sp1.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=51)
sp1.toSimulation(color=(0.152, 0.368, 0.988),material=Mat2)

############################# box and plate #############################
yade.qt.View()

box=O.bodies.append(utils.box(box_center,box_extents,mask=31,material=Mat3,wire=True))
plate1IDs=O.bodies.append(geom.facetBox(plate_center,plate_extents,material=Mat3,wallMask=32))

for i in plate1IDs:
 body= O.bodies[i]
 body.state.vel = (0,0,-10)

############################## compaction #############################
O.dt=.005*utils.PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),

]
#O.run()

Hi,

> how to make a packing with 2 psdSizes and psdCumm while using one packSpherepack()

sp.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes1,psdCumm=passing1,num=N1)
sp.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes2,psdCumm=passing2,num=N2)

HTH
Bruno

Jan Stránský (honzik) said : #7

Hi Jan,

> My question is how to apply force on the spheres without using a "plate" facet?

just use "plate" box instead of "plate" facet

> I am currently using geom.facetBox ... but still smaller spheres can pass through it.

use utils.box(es) instead.
As discussed, there are two un-investigated "problems" in your simulation.
1) "walls" are zero-thickness facets. IMO, using boxes with finite thickness should help
2) you have overlapping packing, giving (all, but especially the smaller ones) particles unrealistically high velocities, which helps them to "escape".
I cannot say which (if any) modification will help more, you have to try..

> Also, how to make a packing with 2 psdSizes and psdCumm while using one packSpherepack() and one toSimulation to avoid initial overlaping of spheres?

the same way as now, just delete the first sp.toSimulation()
If you worry about material, you can assign it afterwards according to particle sizes (if the two PSDs are "distinct") or according to len(sp) before the second makeCloud

cheers
Jan

Othman Sh (othman-sh) said : #8

Thank you Bruno and Jan.