# yade-users team mailing list archive

## [Question #695989]: Packing with two PSD

```New question #695989 on Yade:

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

Note: I am using yadedaily to run the code below.

-----------------------
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
heightcyl=610e-6
##center of cylinder
cx=0
cy=0
cz=0
##Initial cube dimensions ###
mnz=0
mxz=heightcyl*1.1

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

cylinder=O.bodies.append(facets)

##creating disks

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)
stress.append(s)
thickness=(O.bodies[disk1IDs[1]].state.pos[2])-(O.bodies[disk2IDs[1]].state.pos[2])