Incorrect porosity after predicate filtering

Hello,

> sp.load('packtest_rMean0.0012_sig-100.txt')

please provide a MWE [1].

read the documentation [2].
Probably you should use volume argument. Otherwise it assumes volume of enclosing axis aligned box.
So basically you have same volume, but less particles, so the resulting higher porosity is reasonable.

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.utils.html#yade._utils.porosity

Hello Jan,
thank you for the answer.
Regarding to point [1], here is the "get to the right porosity" Script:

readParamsFromTable(rMean=.0012, sig3=-100)
from yade import pack, plot
from yade.params.table import *

mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1)
compFricDegree=30
damp=0.2
young=5e6
stabilityThreshold=0.01
targetPorosity=.4

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2500,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,rMean=rMean,seed=1)
sp.makeCloud(mn,mx,rMean=.5*rMean,seed=1)
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

newton=NewtonIntegrator(damping=damp)

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,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

triax.goal1=triax.goal2=triax.goal3=sig3

while 1:
 O.run(1000, True)
 unb=unbalancedForce()
 print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
 if unb<stabilityThreshold and abs(sig3-triax.meanStress)/sig3<0.001:
  print('Compression equilibrium reached')
  break

import sys
while triax.porosity>targetPorosity:
 compFricDegree*=.95
 setContactFriction(radians(compFricDegree))
 print('\r Friction: ',compFricDegree,' porosity:',triax.porosity,)
 sys.stdout.flush()
 O.run(250,1)

if triax.porosity<targetPorosity:
  sp.fromSimulation()
  sp.save(f'packtest_rMean{rMean}_sig{sig3}_n_{triax.porosity}.txt')
  print('Pack ist exportiert')
  O.pause()

waitIfBatch()

---------------------------------------------------------------------
Regarding point [2]:

I have tried the volume argument and unfortunately get different results:
-------------------------------------------------------------------------------
For script 3 (Predicate is a box).
I get a similar porosity when I insert the volume of the predicate into the porosity function.
.5672 without argument
.56898 with Predicate Volume as argument

from yade import pack

sp=SpherePack()
sp.load('packtest_rMean0.0012_sig-100.txt')

print('target porosity is 0.38542')

pred=pack.inAlignedBox(minAABB=(-.01,-.01,-.01), maxAABB=(.01, .01, .01))
O.bodies.append(geom.facetBox((0, 0, 0), (.01, .01, .01), wallMask=31))

spf=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True)
spf.toSimulation()
print('porosity without argument',porosity())
print('porosity with argument volume',porosity(volume=(.02*.02*.02)))

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

For script 2 (Predicate is a cylinder).
I get a lower porosity if I insert the volume of the predicate into the porosity function. Nevertheless, this is much too high.

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

from yade import pack

sp=SpherePack()
sp.load('packtest_rMean0.0012_sig-100.txt')
print('target porosity is 0.38542')

O.bodies.append(geom.facetCylinder(center=(0, 0, .01), radius=.01, height=.02,segmentsNumber=50, wallMask=6))
pred=pack.inCylinder(centerBottom=(0,0,0), centerTop=(0,0,.02), radius=.01)

spf2=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True)
spf2.toSimulation()
print('porosity without argument',porosity())
print('porosity with argument', porosity(volume=((pi*.01**2)*.02)))

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

Unfortunately, at this point I do not know how to proceed, since now only the proof has been given that the porosity function uses the volume of the predicate, which it is supposed to do.
Why there is a larger deviation with the cylinder predicate, I cannot explain myself.

I would still like to find a way to work with the porosity function, since I would like to map the pore number over the stress later in the oedometer experiment.

I conclude that the filtered section of the imported sphere packing has a smaller total particle volume in relation to the new (predicate) volume than the total particle volume in relation to the final compressed particle cloud (old volume) from [1].

Since in [1] the particle cloud was compressed isotropically, I assumed that the porosity is equally distributed over the entire volume and the cropped part accordingly also has the same porosity as the entire sphere packing volume.

How do I manage to get the same porosity in my predicate as it is in the old volume (imported pack)?

Thanks in advance,
Lukas

Hi Lukas,

When I worked with filterSpherePack function, I made some observations, which might be helpful.

> Since in [1] the particle cloud was compressed isotropically, I assumed that the porosity is equally distributed over the entire volume

I would expect that the packing’s porosity is lowest in the center of the packing. At the boundaries there will be empty spaces, which increase the porosity close to the packing’s boundary.

> and the cropped part accordingly also has the same porosity as the entire sphere packing volume.

When you crop the packing, spheres which are partially inside and outside of the predicate will be deleted. This will cause an increase in empty space and thus increase porosity. How much porosity will increase due to this effect, will depend on your packing (sphere sizes) and predicate.

> How do I manage to get the same porosity in my predicate as it is in the old volume (imported pack)?

Good question. One option would be to add a facet cylinder and to grow particles until you meet the target porosity.

Hope it helps. Regards,
Bettina

> here is the "get to the right porosity" Script:

sorry, I run it for a while without any convergence.
Can you just share the resulting .txt file?

> I get a lower porosity if I insert the volume of the predicate into the porosity function. Nevertheless, this is much too high.

please provide specific values

Otherwise follow Bettina's answer, the boundary effects (both for the cube and the filtered cylinder case) were my first guesses.

Cheers
Jan

Hey Bettina and Jan,

>Can you just share the resulting .txt file?

here is the saved .txt Sphere Pack: https://drive.google.com/file/d/13l0QlEVK_3OOR0oGTQh6O5CiE9DXSiy8/view?usp=sharing

> please provide specific values

target porosity is .385
box predicate porosity without volume argument in porosity function is .5672
box predicate porosity with volume argument in porosity function is .5689

cylinder predicate porosity without volume argument in porosity function is .6647
cylinder predicate porosity with volume argument in porosity function is .5821

> One option would be to add a facet cylinder and to grow particles until you meet the target porosity.

for understanding:
the filtering process is skipped as the porosity of my particle cloud is performed by enlarging the spheres in the same cylinder vessel in which the oedometer simulation is subsequently performed?

So far, I wanted to achieve the setting of the correct porosity by "external" compaction, because I can not estimate how much a different particle size has influence on the material calibration by oedometer and direct shear test. In addition, settlement should be calculated in the future.

> Otherwise follow Bettina's answer, the boundary effects (both for the cube and the filtered cylinder case) were my first guesses.

It is probably not very helpful to compress the particle cloud with PeriTriaxController to the desired pore number, because mainly the deletion of overlapping spheres in the edge region of the predicate are the cause for the too high porosity?

Can oedometer and direct shear tests be described and modelled/simulated with periodic boundary conditions and if so, would this be a solution for the boundary effect problem when it comes to having the right porosity?

Best regards,
Lukas

>> One option would be to add a facet cylinder and to grow particles until you meet the target porosity.
> for understanding:
> the filtering process is skipped as the porosity of my particle cloud is performed by enlarging the spheres in the same cylinder vessel in which the oedometer simulation is subsequently performed?

yes

> So far, I wanted to achieve the setting of the correct porosity by "external" compaction, because I can not estimate how much a different particle size has influence on the material calibration by oedometer and direct shear test. In addition, settlement should be calculated in the future.

particle growth is just one option.
You can e.g. perform simple gravity deposition or directly oedometric compression.

> It is probably not very helpful to compress the particle cloud with PeriTriaxController to the desired pore number, because mainly the deletion of overlapping spheres in the edge region of the predicate are the cause for the too high porosity?

seems reasonable to me (according to so far results)

> Can oedometer and direct shear tests be described and modelled/simulated with periodic boundary conditions and if so, would this be a solution for the boundary effect problem when it comes to having the right porosity?

for oedometric case, yes and yes

for direct shear test, in **my opinion**, no and no

Cheers
Jan

Hello,
I tried a new approach that allowed me to achieve the target porosity, at least for a cuboid predicate.

For a cylindrical predicate, with the same approach, there is a higher porosity at the beginning, but it slowly decreases with increasing iteration (I did not run the simulation completely due to time constraints).

The approach is based on Bettina's hints that spheres from the sphere packing are inserted into the predicate; however, spheres that intersect the predicate are completely removed. This leads to an increased porosity, which is basically a boundary value problem.

The idea is this:
the predicate is chosen larger than the facet volume.
The predicate boundary is increased by the maximum existing sphere radius in the sphere packing in each spatial dimension by a factor of 2.

If the predicate is larger than the actual volume, 2 things happen:
a) spheres that lie in the predicate but intersect the facet volume are not deleted. This increases the number of spheres in the volume and thus decreases the porosity.
b) spheres with a sphere radius smaller than the maximum sphere radius can lie outside the facet volume (since the predicate is larger than it).

Now the following is done:
Since the spatial facet boundaries are known, all spheres are deleted from the filtered sphere packing that have their center (state.pos) outside these boundaries (condition b))
Spheres whose centers are within the boundaries but intersect the facets will slip into the volume when the simulation is started (condition a)).

The start of the simulation is somewhat chaotic, but after an equilibrium has been established, the actual simulation (e.g. oedometer test) can be started.

--------------------------------------------------------------------------------------------------------
generating sphere pack with target porosity:
-----------------------------------------------------------------------------------------------------

readParamsFromTable(factor=5, sig3=-10000)
from yade import pack, plot
from yade.params.table import *

targetPorosity = 0.4 #the porosity we want for the packing
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1) # corners of the initial packing
compFricDegree=20
psdSizes,psdCumm=[factor*.0003,factor*.0004,factor*.0006,factor*.0008],[0,.5,.5,1]

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1)
sp.toSimulation(material='spheres')
print('Anzahl Kugeln:', len(sp))

triax=TriaxialStressController(
 thickness = 0,
 stressMask = 7,
 internalCompaction=False, # If true the confining pressure is generated by growing particles
 goal1=sig3,
 goal2=sig3,
 goal3=sig3
)

newton=NewtonIntegrator(damping=damp)

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,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton,
 PyRunner(command='checkUnbalanced()', iterPeriod=1000, label='checker'),
 PyRunner(command='addPlotData()', iterPeriod=200),
]

def checkUnbalanced():
 unb=unbalancedForce()
 print('unbalanced Force: ', unb, ' mean stress: ', triax.meanStress)
 if O.iter> 1000 and unb<stabilityThreshold and abs((triax.goal1-triax.meanStress)/triax.goal1)<.001:
  checker.iterPeriod=500
  checker.command='reduceFric()'

def reduceFric():
 if triax.porosity>targetPorosity:
  compFricDegree=.95*compFricDegree
  setContactFriction(radians(comFricDegree))
  print('Friction: ', compFricDegree, 'Porosity: ',triax.porosity)
  #sys.stdout.flush()
 else:
  #sp.fromSimulation()
  #sp.save(f'STRESSED_Pack_factor_{factor}_porosity_{triax.porosity}.txt')
  print('Destress')
  checker.iterPeriod=150
  checker.command='destress()'

def destress():
 triax.goal1=sig3*.01
 triax.goal2=sig3*.01
 triax.goal3=sig3*.01
 if abs((triax.goal1-triax.meanStress)/triax.goal1)<.001 and triax.porosity<targetPorosity:
  print('Porosity Reached: ',triax.porosity)
  sp.fromSimulation()
  #sp.save(f'DESTRESSED_Pack_factor_{factor}_porosity_{triax.porosity}.txt')
  sp.save('testPack.txt')
  O.pause()

def addPlotData():
 plot.addData(unbalanced=unbalancedForce(), Porosity=triax.porosity, Sigma=triax.goal1, MeanStress=triax.meanStress, i=O.iter)

plot.plots = {'i': ('unbalanced'), 'i ': ('Porosity'), 'i ': ('Sigma', 'MeanStress')}
plot.plot()

O.run()

waitIfBatch()

-----------------------------------------------------------------
filter sphere pack to prediacte
-----------------------------------------------------------------

from yade import pack, plot
import numpy as np

#from yade.params.table import *
height=.02
#radiusoed=.0356825
radiusoed=.03
overlap=True
Cyl=False

sp=SpherePack()

sp.load('testPack.txt')
maxR=max(s[1] for s in sp)
minR=min(s[1] for s in sp)
meanR=np.average([s[1] for s in sp])
if Cyl==False:
 O.bodies.append(geom.facetBox(center=(radiusoed,radiusoed,height/2),extents=(radiusoed,radiusoed,height/2,),wallMask=63))
 if overlap==False:
  pred=pack.inAlignedBox(minAABB=(0,0,0),maxAABB=(radiusoed*2,radiusoed*2,height))
 else:
  pred=pack.inAlignedBox(minAABB=(0-maxR,0-maxR,0-maxR),maxAABB=(radiusoed*2+maxR,radiusoed*2+maxR,height+maxR))

else:
 O.bodies.append(geom.facetCylinder(center=(0, 0, height/2), radius=radiusoed, height=height,segmentsNumber=2000, wallMask=7)) #oder6
 if overlap==False:
  pred=pack.inCylinder(centerBottom=(0,0,0), centerTop=(0,0,height), radius=radiusoed)
 else:
  pred=pack.inCylinder(centerBottom=(0,0,0-maxR/4), centerTop=(0,0,height+maxR/4), radius=radiusoed+maxR/3)

sp=pack.filterSpherePack(predicate=pred,spherePack=sp,returnSpherePack=True)
sp.toSimulation()

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.shape.color=scalarOnColorScale(x=b.shape.radius, xmin=minR, xmax=maxR)

print('Anzahl Körper vor löschen: ',len(O.bodies))
ListKugeln=[]
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  ListKugeln.append(b.id)

ListDel=[]

if Cyl==False:
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   if b.state.pos[0]<0:
    ListDel.append(b.id)
   if b.state.pos[0]>radiusoed*2:
    ListDel.append(b.id)
   if b.state.pos[1]<0:
    ListDel.append(b.id)
   if b.state.pos[1]>radiusoed*2:
    ListDel.append(b.id)
   if b.state.pos[2]<0:
    ListDel.append(b.id)
   if b.state.pos[2]>height:
    ListDel.append(b.id)
else:
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   if b.state.pos[2]<0:
    ListDel.append(b.id)
   if b.state.pos[2]>height:
    ListDel.append(b.id)
   if sqrt(b.state.pos[0]**2+b.state.pos[1]**2)>radiusoed:
    ListDel.append(b.id)

ListKugeln=list(set(ListKugeln))
ListDel=list(set(ListDel))
for x in ListDel:
 O.bodies.erase(x)

print('Anzahl aller Kugeln: ', len(ListKugeln))
print('Anzahl gelöschter Kugeln: ',len(ListDel))
print('Anzahl alle Kugeln nach Bereinigung: ', len(ListKugeln)-len(ListDel))

print('Porosity after filter: ', porosity())

O.engines = [
        ForceResetter(),

        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(

                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.9),

        # the label creates an automatic variable referring to this engine
        # we use it below to change its attributes from the functions called
        PyRunner(command='checkUnbalanced()', iterPeriod=20, label='checker'),
        PyRunner(command='addPlotData()', iterPeriod=20),
]

O.dt =PWaveTimeStep()

def checkUnbalanced():
 if unbalancedForce()<.005:
  print('Porenanteil n : ',porosity())
  print('Porenzahl e: ', porosity()/(1-porosity()))
  ListLast=[]
  for b in O.bodies:
   if isinstance(b.shape, Sphere):
    ListLast.append(b.id)
  print('Anzahl Kugeln nach Bewegung: ',len(ListLast))
  print('Kugeln die verloren gegangen sind: ',len(ListKugeln)-len(ListDel)-len(ListLast))
  O.pause()

def addPlotData():
 plot.addData(unbalanced=unbalancedForce(),Porosity=porosity(),i=O.iter)

plot.plots = {'i': ('unbalanced'), 'i ': ('Porosity')}
plot.plot()

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

Note that at the beginning of the filter predicate script there are boolean variables that can be used to set the predicate form (Cyl=False for cuboid, Cyl=True for cylinder).
In addition, the particle-facet overlap can also be switched on and off here (overlap=True).

best regards,
Lukas