specific porosity and radius of spheres in a cloud

Hi Wojtek!

You have two solutions:
1/ from target porosity and target radius you can deduce the number of
particles analyticaly, then you give this number to makeCloud()
2/ instead of growing the particles, you give them the target radius and
you compress them by moving the walls [1], downside: the final size of
the box is changing, and the compaction is not that homogeneous

Bruno

Hi Bruno!

1. I use a slightly different way. First I determine the target porosity and the number of particles, and next I calculate the theoretical volume of the bed. It turned out that the final volume is bigger as the initial volume of the cloud (given by min and max corners). After several attempts I found the optimal size of the initial volume. The defect is that the average radius is in every case is a little different. For this reason I want to know how to calculate the average radius in YADE.

2. I tried likewise to resolve my task in this way. The problem is that the calculation time is very long and in fact I don't understand fully the python script. I tried to connect the trick with precisely defined target porosity and the example with triax test. I afraid that some parts may be without sense. For example I still don't know what is the triax.goal.

If you have few minutes - check my script below and give me some remarks how to improve it, please. If not, I will try do it alone.

Regards,

Wojtek.

-------
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 = 1000

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

#liczba iteracji:
maxIter = 150000

#plik:
file = 'p06_25'

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

#czytanie danych o rozkladzie zloza:
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()

#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 ' 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 = 5

#--------------------------------------------------------------------------------------------
#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)

#kasowanie danych do wizualizacji:
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(file+'_'+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()

triax.goal1 = triax.goal2 = triax.goal3 = 10000

import sys #this is only for the flush() below
while triax.porosity > targetPorosity:
  compFricDegree = 0.95*compFricDegree
  setContactFriction(radians(compFricDegree))
  print "\r Friction:",compFricDegree
  print "\r Porosity:",triax.porosity
  sys.stdout.flush()
  # while we run steps, triax will tend to grow particles as the packing
  # keeps shrinking as a consequence of decreasing friction. Consequently
  # porosity will decrease
  O.run(500,1)

from yade import plot, utils

print "\r --------"
print "\r Final friction:",compFricDegree
print "\r Final porosity:",triax.porosity

from yade import export
export.text(file+'.yade')

O.pause()

Hi Bruno!

I sent you wrong script - this is my last try:

#--------------------------------------------------------------------------------------------
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

#porowatosc docelow:
targetPorosity = 0.41

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

#liczba iteracji:
maxIter = 150000

#plik:
file = 'p06_25'

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

#czytanie danych o rozkladzie zloza:
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()

#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 ' 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 = 5 # 0.449

#--------------------------------------------------------------------------------------------
#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()

#---------------------------------------------------------------------------------
triax = TriaxialStressController(
  thickness = 0,
  stressMask = 7,
  internalCompaction = True,
  )

#--------------------------------------------------------------------------------------------
#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()]
 ),
 #silnik intgracji rownan ruchu Newtona:
 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(file+'_'+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 = 0.4*PWaveTimeStep()

#???
triax.goal1 = triax.goal2 = triax.goal3 = 10000

import sys
while triax.porosity > targetPorosity:
  spFrAngle = 0.95*spFrAngle
  setContactFriction(radians(spFrAngle))
  print "\r Friction:",spFrAngle
  print "\r Porosity:",triax.porosity
  sys.stdout.flush() #???
  O.run(500,1)

from yade import plot, utils

print "\r --------"
print "\r Final friction:",spFrAngle
print "\r Final porosity:",triax.porosity

from yade import export
export.text(file+'.yade')

O.pause()

#--------------------------------------------------------------------------------------------

Regards,

Wojtek

For 1/, I don't understand why you don't get exactly what you calculate
theoreticaly. Usually it works. There must be something wrong somewhere.
How do you calculate the theoretical volume? Then how do you calculate
the average radius?
I don't see it very helpful at this point to think about how yade
calculates average radius. The short answer is: it does not.
Average radius is not well defined for arbitrary particle size
distributions.

Missing in my post above, the [1] link:
[1]
https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py#L94

Bruno

Hi,

About 2, unfortunately Murphy's law says that you usually do not get what you expect until you fully understand your script.. So, try to achieve this goal, we're here to try to answer - precise, ideally ;-) - questions !

For example, triax.goal1/2/3 is generally a python variable reflecting the corresponding attribute of a TriaxialStressController (or similar) engine used in a simulation, see [1] . This engine being labelled 'triax', see [2] and examples if you do not know what a label is in Yade.

Here, in the 2d version of your script, you define indeed (directly, not by label) such an Engine as 'triax', but you do not use it in your simulation. (Because it is not in the O.engines list definition)

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.goal1 You will need to understand also stressMask variable to fully understand goal variable. If the corresponding doc is not enough, give a look to Bruno's link in answer #4, few lines above l.94)
[2] https://yade-dem.org/doc/user.html#labeling-things

Hi Bruno,

And if you define the size

mn = Vector3(0,0,0)
mx = Vector3(length,height,width)

Is it true that this is the volume inside which the cloud is generated on the beginning? But if you increase the size of particles, some of them will be protrude from this volume after some time. In consequence the volume containing all spheres will be bigger on the end compare to the beginning. Is my thinking correct?

Thank you for help.

Regards,

Wojtek.

Hi Jérôm!

thank you for the answer. You're right of course that I should understand all parts of the script. This is my aim for the future. Thank you for links.

I will try to resolve my problem with the first method - looks easier.

Regards,

Wojtek.

On 20/10/14 21:36, Wojciech Sobieski wrote:
> Is it true that this is the volume inside which the cloud is generated
> on the beginning? But if you increase the size of particles, some of
> them will be protrude from this volume after some time.

If the volume is bounded by rigid boundaries, then it is always the same
(exception for the very small overlap that can occure at contacts but it
is irrelevant to the present discussion).
As soon a particle would protrude it is pushed back to an inner part of
the box.

B

Hi Bruno,

I found a bug in my calculations. Now all looks good.

Thank you for help.

Regards,

Wojtek.