Yade

Dense packing with given distribution of radii

Asked by Olav Ragnvaldsen

Hi all,

I want to create a spherical dense packing of spheres, where the radii of the spheres are given by a Gaussian distribution. Since I want to assign the radii myself, I have seen on some questions (https://answers.launchpad.net/yade/+question/406741) that filterSpherePack should be the way to go. But I am unsure how I can make my spherePack dense after filling the predicate? In the question linked above it seems as I could do it using an isotropic compression, but when I try that all my spheres "explode". And since my body of spheres is also spherical, it seems to me that a compression from 6 different planes will only distort my sphere.

 Is there another way to achieve a densely packed spherical body with assigned radii? I have seen on https://answers.launchpad.net/yade/+question/211583
that it is suggested to change the source code for spherePack, how can this be done?

I am new with Yade so I greatly appreciate any help!

Cheers, Olav

Below is my script for creating the spherical loose-packed body using filterSpherePack:

from yade import pack
import numpy as np
## Creating a body of spheres of varying radii

# Using numpy to generate a Gaussian distribution of sphere diameters
mu, sigma=20,2 #mean and standard deviation
numParticles=1000
d=np.random.normal(mu,sigma,numParticles) #diameter array

#Sorting the diameters
psdSizes=sorted(d)
#psdCumm=np.arange(1/numParticles,1+1/numParticles,1/numParticles)
#Creating a vector of accumulated factors of the particles
#Here we have one particle assigned to each radius
#for some reason Yade says 1/numParticles=0, so 0.001 is used instead
psdCumm=np.arange(0.001,1+0.001,0.001)

## making cloud of spheres to insert into predicate
sp=pack.SpherePack()
#giving the spheres a radius assigned by Normal distribution
sp.makeCloud((0,0,0),(200,200,200),psdSizes=psdSizes,psdCumm=psdCumm
              ,num=numParticles,porosity=0.2)

## constructing predicate
pred=pack.inSphere(center=(0,0,0),radius=100)

## packing predicate with loosely packed spheres
spheres=pack.filterSpherePack(pred,sp)

## add spheres to simulation
O.bodies.append(spheres)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Olav Ragnvaldsen
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hi,

> But I am unsure how I can make my spherePack dense after filling the predicate?
> And since my body of spheres is also spherical, it seems to me that a compression from 6 different planes will only distort my sphere.

one the approach of randomDensePack is to create a densePacking and than "crop" it according to predicate

> it is suggested to change the source code for spherePack, how can this be done?

e.g.:

####################################
def myMakeCloud(mn,mx,mu,sigma,rMin,num,maxTry=1000,seed=0): # [1]
   import random
   import numpy as np
   random.seed(seed)
   mn = Vector3(mn)
   mx = Vector3(mx)
   size = mx-mn
   rnd = random.random
   ret = []
   radii = np.random.normal(mu,sigma,num)
   radii = map(lambda r: max(r,rMin), radii)
   for r in radii:
      for t in xrange(maxTry):
         c = Vector3([mn[axis] + r + (size[axis]-2*r)*rnd() for axis in (0,1,2)])
         packSize = len(ret)
         overlap = False
         for pc,pr in ret:
            if pow(pr+r,2) >= (pc-c).squaredNorm():
               overlap = True
               break
         if not overlap:
            ret.append((c,r))
            break
      if t == maxTry:
         raise RuntimeError, "WRNING: Exceeded {} tries to insert non-overlapping sphere to packing. Only {} spheres were added, although you requested {}.".format(maxTry,i,num)
   return pack.SpherePack(ret)

def myRandomDensePack(predicate,mu,sigma,rMin,material=-1,spheresInCell=300,color=None,returnSpherePack=None,seed=0): # [2]
   dim=predicate.dim()
   if max(dim)==float('inf'): raise RuntimeError("Infinite predicate and no dimension of packing requested.")
   center=predicate.center()
   fullDim=dim
   cloudPorosity=0.25
   beta,gamma=fullDim[1]/fullDim[0],fullDim[2]/fullDim[0]
   N100=spheresInCell/cloudPorosity
   x1=mu*(1/(beta*gamma)*N100*(4/3.)*pi)**(1/3.)
   y1,z1=beta*x1,gamma*x1; vol0=x1*y1*z1
   O.switchScene(); O.resetThisScene() ### !!
   O.periodic=True
   O.cell.setBox((x1,y1,z1))
   O.materials.append(FrictMat(young=3e10,density=2400))
   sp = myMakeCloud(Vector3().Zero,O.cell.refSize,mu,sigma,rMin,spheresInCell,seed=seed)
   O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=.05*mu),InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*mu,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=1,keepProportions=True),NewtonIntegrator(damping=.6)]
   O.materials.append(FrictMat(young=30e9,frictionAngle=.5,poisson=.3,density=1e3))
   for s in sp: O.bodies.append(utils.sphere(s[0],s[1]))
   O.dt=utils.PWaveTimeStep()
   O.run(); O.wait()
   sp=SpherePack(); sp.fromSimulation()
   O.switchScene() ### !!
   sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
   return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)

mu, sigma = 20,2
rMin = .25*mu
pred = pack.inSphere(center=(0,0,0),radius=100)
sphs = myRandomDensePack(pred,mu,sigma,rMin,spheresInCell=300)
O.bodies.append(sphs)
####################################

cheers
Jan

[1] https://answers.launchpad.net/yade/+question/667571
[2] https://github.com/yade/trunk/blob/master/py/pack/pack.py#L454

Revision history for this message
Olav Ragnvaldsen (olavrag) said : 		#2

Thank you, Jan! That solved my problem

-Olav