a cloud with a target porosity

Hi,

I am a beginer in YADE. I want to prepare a cloud with defined distribution (table of psdSizes and psdCumm are from a file) and with defined porosity. Is it possible to obtain with a simple script like this below (no 1)? Why doesn't work the option "porosity="? The cloud has porosity 0.75 what is too big. I know from an experiment, that the bed has porosity equal to 0.41.

Alternatively, I try to resolve the problem with a trixial test, but the porosity is still too big (no 2) - it is about 0.46-0.47. Is in general the script no 2 prepared correctly for my problem?

Which way is better and what should be changed to resolve my problem?

I am using YADE 1.11.0-41-6f921d5~trusty.

I will be grateful for your help,

Wojtek

---------------- script no 1:

from yade import pack

l_x = 0.3 #[m]
l_y = 0.6 #[m]
l_z = 0.3 #[m]

file = 'p06_25'
number = 5000
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

corner1 = Vector3(0.0,0.0,0.0)
corner2 = Vector3(l_x,l_y,l_z)
sp = pack.SpherePack()
sp.makeCloud(corner1,corner2,psdSizes=psdSizes,psdCumm=psdCumm,num=number,porosity=0.42,distributeMass=True,periodic=False)
O.bodies.append([sphere(s[0],s[1],wire=False) for s in sp])

print ' particles = ',len(O.bodies)
print ' porosity = ',utils.porosity(l_x*l_y*l_z)

from yade import qt
qt.View()

---------------- script no 2 (sorry for Polish comments):

from yade import plot
from yade.pack import *

from yade import utils

O=Omega()

#--------------------------------------------------------------------------------------------
#rozmiar obszaru obliczeniowego:
length = 0.25
height = 0.5
width = 0.25

#grubosc scianek (musza to byc bryly a nie plaszczyzny):
thickness = 0.001

#liczba czastek:
number = 5000

#predkosc kompresji (scianki, ktora sie rusza):
v = 1.0

#liczba iteracji:
maxIter = 150000

#plik:
plik = 'p06_23'

#sprawdzenie liczby przedzialow:
n_band = sum(1 for line in open(plik+'.txt'))

#czytanie danych o rozkladzie zloza:
input = open(plik+'.txt','r')
psdSizes = []
psdCumm = []
try:
  for linijka in input:
    psdSizes.append(float(linijka[0:16]))
    psdCumm.append(float(linijka[18:]))
finally:
  input.close()

#korekta wartosci skrajjnych (musi byc 0 i 1):
psdCumm[0] = 0.0
psdCumm[n_band-1] = 1.0

#informacje o zadaniu:
os.system('clear')
print '|-------------------------------------|'
print '| Particle distribution Olsztyn 2014 |'
print '|-------------------------------------|'
print ''
print ' length = ', length
print ' height = ', height
print ' width = ', width
print ' thickness = ', thickness
print ' v = ', v
print ' maxIter = ', maxIter
print ' n_bands = ',n_band
print ' psdSizes, psdCumm '
print ''

#dane stali:
stDensity = 7860 #sprawdzone
stYoung = 2e11 #sprawdzone
stPoisson = 0.3 #sprawdzone
stFrAngle = 30 #

#dane szkla:
spDensity = 2600 #sprawdzone
spYoung = 50e9 #sprawdzone
spPoisson = 0.18 #sprawdzone 0.18-0.3
#spFrAngle = 15 # 0
#spFrAngle = 10 #
spFrAngle = 15 #

#--------------------------------------------------------------------------------------------
#definicja materialu scian:
O.materials.append(FrictMat(density=stDensity,young=stYoung,poisson=stPoisson,frictionAngle=radians(stFrAngle),label='wallMat'))

#definicja poszczegolnych bokow kostki tworzacej obszar obliczeniowy:
leftBox = box( center=(-thickness/2.0,(height)/2.0,0), extents=(thickness/2.0,5*(height/2.0+thickness),width/2.0) ,fixed=True,wire=True)
lowBox = box( center=(length/2.0,-thickness/2.0,0), extents=(length/2.0,thickness/2.0,width/2.0) ,fixed=True,wire=True)
rightBox = box( center=(length+thickness/2.0,height/2.0,0), extents=(thickness/2.0,5*(height/2.0+thickness),width/2.0) ,fixed=True,wire=True)
upBox = box( center=(length/2.0,height+thickness/2.0,0), extents=(length/2.0,thickness/2.0,width/2.0) ,fixed=True,wire=True)
behindBox = box( center=(length/2.0,height/2.0,-width/2.0-thickness/2.0), extents=(2.5*length/2.0,height/2.0+thickness,thickness/2.0), fixed=True,wire=True)
inFrontBox = box( center=(length/2.0,height/2.0,width/2.0+thickness/2.0), extents=(2.5*length/2.0,height/2.0+thickness,thickness/2.0), fixed=True,wire=True)

#dodawanie cial do sceny - mozna to uzyc w state.pos[]
O.bodies.append([leftBox,lowBox,rightBox,upBox,behindBox,inFrontBox])

#--------------------------------------------------------------------------------------------
#definicja materialu sfer:
O.materials.append(FrictMat(density=spDensity,young=spYoung,poisson=spPoisson,frictionAngle=radians(spFrAngle),label='spMat'))

#definiowanie pozycji chmury:
mincorner = Vector3(0,0.0,-width/2.0)
maxcorner = Vector3(length,height,width/2.0)

#tworzenie chmury:
sp = yade._packSpheres.SpherePack()
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)

#dodanie chmury do sceny:
O.bodies.append([sphere(s[0],s[1],material='spMat') for s in sp])

#--------------------------------------------------------------------------------------------
#informacja o liczbie sfer i o porowatosci:
print ' particles = ',len(O.bodies)
print ' porosity = ',utils.porosity(length*height*width)

#plot.resetData()

#--------------------------------------------------------------------------------------------
#definiowanie ustawien solwera:
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()]
 ),
 NewtonIntegrator(damping=.4),
 PyRunner(iterPeriod=1,command='letSave()'),
 PyRunner(command='letPlot()',iterPeriod=100),
        PyRunner(iterPeriod=100,command='print voxelPorosity(100,mincorner,O.bodies[3].state.pos)'),
 ]

def letSave():
   if O.iter == maxIter-1:
     from yade import export
     export.text(plik+'_'+str(number)+".yade")
     print "---------------------------"
     print "maxIter = ",maxIter
     print "mumber of particles = ",len(O.bodies)
     print "hight current = ", O.bodies[3].state.pos[1]-thickness/2.0
     print "porosity = ",utils.porosity(length*(O.bodies[3].state.pos[1]-thickness/2.0)*width)

def letPlot():
    plot.addData(step=O.iter,e=utils.porosity(length*(O.bodies[3].state.pos[1]-thickness/2.0)*width) )

plot.plots={'step':('e')}
plot.plot()

from yade import qt; qt.Controller();qt.View()

#definiowanie warunkow kompresji:
O.engines = O.engines+[KinemCTDEngine(compSpeed=v,sigma_save=(),temoin_save=(),targetSigma=40000.0,LOG=False)]
O.dt=.4*PWaveTimeStep()

#O.run(maxIter)
O.stopAtIter=maxIter

from yade import plot, utils