Radius Expansion Method

Hi,

I use the Radius Expansion Method (REM) to generate a granular bed with a specified particle distribution. I use a discrete form of the cumulative curve of particle distributions. However, after simulation the distribution of particles in the bed differs from the theoretical distribution (see: http://pracownicy.uwm.edu.pl/wojsob/pliki/tmp/yade-distribution.jpg). The effect is independent on the number of bins used to define the discrete distribution function and on the friction angle. I checked also that the initial cloud has correct distribution. Am I doing something wrong (see the code below)? Or maybe such effect it is a feature of the REM?

Best Regards from Poland,

Wojciech Sobieski

--------------------------------------------------------------------------------------------------------------------

from yade.pack import *
from yade import utils
from yade import export
from yade import plot
import sys
import os
import subprocess

file = 'p06_w1_10'
n_s = 10000
e_target = 0.413
mi = 0.5
seed = 5

n_band = sum(1 for line in open(file+'.txt'))
input = open(file+'.txt','r')
psdSizes = []
psdCumm = []
try:
  for line in input:
    psdSizes.append(float(line[0:16]))
    psdCumm.append(float(line[18:]))
finally:
  input.close()
psdCumm[0] = 0.0
psdCumm[n_band-1] = 1.0

subprocess.call('./calc_l.out '+str(file)+' '+str(n_s)+' '+str(e_target),shell=True)
input = open(str(file)+'.l_c','r')
try:
  for line in input:
    l = float(line[0:])
finally:
  input.close()

mn = Vector3(0,0,0)
mx = Vector3(l,2*l,l)

O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))
O.bodies.append(aabbWalls([mn,mx],thickness=0,material='walls'))
O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(mi),density=2600,label='spheres'))
sp = yade._packSpheres.SpherePack()
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=n_s,distributeMass=1,seed=seed)
O.bodies.append([sphere(s[0],s[1],material='spheres') for s in sp])

#-------------------------------------------------------------------------------------
# here the particle distribution is correct:
vtkExporter = export.VTKExporter(str(file)+'-init.vtk')
vtkExporter.exportSpheres(what=[('dist','b.state.pos.norm()')])

os.mkdir(str(file))

triax = TriaxialStressController(
   finalMaxMultiplier = 1.0001,
   maxMultiplier = 1.0001,
   stressMask = 7,
   internalCompaction = True,
   goal1 = 10000,goal2 = 10000,goal3 = 10000)

O.engines=[
 ForceResetter(),
 InsertionSortCollider(
   [Bo1_Sphere_Aabb(),
       Bo1_Box_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack()]),
 GlobalStiffnessTimeStepper(
   active=1,
   timeStepUpdateInterval=100,
   timestepSafetyCoefficient=0.8),
 triax,
 NewtonIntegrator(damping = 0.2,gravity=[0,0,0]),
 VTKRecorder(fileName=str(file)+'/cloud-',recorders=['all'],iterPeriod=10)
 ]

while triax.porosity > e_target:
  mi = 0.999*mi
  setContactFriction(radians(mi))
  print "\r Friction:",mi
  print "\r Porosity:",triax.porosity
  sys.stdout.flush()
  O.run(500,1)

#-------------------------------------------------------------------------------------
# here the particle distribution in not correct:
vtkExporter = export.VTKExporter(str(file)+'-final.vtk')
vtkExporter.exportSpheres(what=[('dist','b.state.pos.norm()')])

input = open('set_A1.porosity','a')
input.write('\r'+str(file)+' '+str(triax.porosity))
input.close()

export.text(str(file)+'.yade')
subprocess.call('./convert.out '+str(file)+' '+str(n_s), shell=True)
subprocess.call('clear', shell=True)
print "\r --------"
print "\r Final friction:",mi
print "\r Final porosity:",triax.porosity