Issues in using psd to create a cloud of spheres

Asked by Othman Sh on 2019-04-19

Hi all,

I have a code that generate a cloud of spheres and compress the cloud to a certain porosity. It works well when I make the cloud with specific radius and rRelFuzz. However, when I try to make the cloud with specific size distribution using psdSizes and psdCumm, the code run but the spheres will freeze and not compact. Can you please help and let me know what went wrong? My code is below.

Thank you,
Othman

----------

from yade import pack

targetp = .3 ##specify the targeted porosity

##specimen geometry
radiuscyl=.05
heightcyl=0.3

# material parameters

O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = atan(.2), density=1920))

############################ spheres #############################
sp=pack.SpherePack()
psdSizes=[0,.00238,.00475,.00952,.0127]
psdCumm=[0,0.0173,.0289,.9469,1]
sp.makeCloud((0,0,0),(.3,.3,1.5),psdSizes=psdSizes,psdCumm=psdCumm)
#sp.makeCloud((0,0,0),(.3,.3,1.5),rMean=.00475) #if I uncomment this line and comment the above 3 lines, the spheres will not move
#### cylinder extraction
pred=pack.inCylinder((.2,.2,0),(.2,.2,heightcyl),radiuscyl)
spFilter=filterSpherePack(pred,sp,Material=Material, returnSpherePack=True)
spFilter.toSimulation()
############################ facets #############################

facets=geom.facetCylinder((.2,.2,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=50,wallMask=4)

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

##creating disks

d1=geom.facetCylinder((.2,.2,heightcyl),radiuscyl-.002,0,segmentsNumber=50,wallMask=1)
d2=geom.facetCylinder((.2,.2,0),radiuscyl-.002,0,segmentsNumber=50,wallMask=1)

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

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

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

############################ compaction #############################
O.dt=.5*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),

 PyRunner(iterPeriod=500,command='stop()'),

]
O.run()
# reset interaction detection enlargement
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0

def stop():
   if utils.porosity()<targetp:
 O.pause()
 print 'Finished'
 for i in disk1IDs: O.bodies.erase(i)
 for i in disk2IDs: O.bodies.erase(i)
 for i in cylinder: O.bodies.erase(i)
 print ('unbalanced forces = ', utils.unbalancedForce())
 print ('cylinder dimensions = ', utils.aabbDim())
 print ('porosity = ', utils.porosity())
 print ('No. of spheres = ', len (O.bodies))

Question information

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

Hello,

Thanks for the MWE. The particles are moving, but the simulation runs slower for two reasons:

1) time step is a few orders of magnitude smaller because you introduce much smaller particles with the PSD.
2) slower simulation speed because you introduce more particles with the PSD.

Cheers,

Robert

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

Hello,

> but the spheres will freeze and not compact

please always try to get also numerical evidence of such statements. If you print position of particles, you will see that they DO move and compact. Then the right question would be, why do they compact so slowly.

I you print minimal particle size and O.dt before O.run:
##
print min(b.shape.radius for b in O.bodies)
print O.dt
##
in the "working" case I got
0.00475
4.65403051129e-07
while for the "frozen" case it was
2.93328060983e-05
2.87401630659e-09
so the "frozen" case runs ~200x slower than the "working" case

The solution is to modify the PSD

cheers
Jan

Robert Caulk (rcaulk) said : #3

Jan's advice about gathering numerical evidence will help you far beyond this forum.

>The solution is to modify the PSD.

Is there ever a valid reason for using '0' as a minimum value in the Yade PSD generator? It seems like that's just asking the RNG to give you the smallest timestep possible. Maybe we should add a warning in the documentation?

Othman Sh (othman-sh) said : #4

Thank you Jan and Robert for the explanations and recommendations. I modified the PSD to solve this issue.

Regards,

Othman