# yade-users team mailing list archive

## [Question #676248]: Radius Expansion Method

```New question #676248 on Yade:

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

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

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'))
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,
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
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()

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

--