yade-users team mailing list archive

[Othman Sh <question695989@xxxxxxxxxxxxxxxxxxxxx>]

New question #695989 on Yade:
https://answers.launchpad.net/yade/+question/695989

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)



-- 
You received this question notification because your team yade-users is
an answer contact for Yade.