Yade

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

Link existing bug

Revision history for this message

Jan Stránský (honzik) said on 2021-03-11:

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

Revision history for this message

Othman Sh (othman-sh) said on 2021-03-15:

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

Hi Jan,

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

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

mass2=utils.getSpheresMass()-mass1

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

total_mass=utils.getSpheresMass()

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

Revision history for this message

Jan Stránský (honzik) said on 2021-03-15:

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

Revision history for this message

Jan Stránský (honzik) said on 2021-03-15:

> 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

Revision history for this message

Othman Sh (othman-sh) said on 2021-03-16:

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

cube_dimensions=100e-6 #micron

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

############################# 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 Jan,

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))