Why the sum of TW.volume is not equal to the volume of pack

Asked by Chu on 2020-06-07

Dear all,

I use TesselationWrapper to compute the volume of spheres, but I find that the sum of TW.volume is not equal to the volume of pack. Does the Voronoi cell of each sphere overlap each other?

I run the MWE and I got

volume=0.6067921587125377
TW_Volume=0.8510512980654049

Thanks for any suggestion.

####MWE####

import numpy as np
O.periodic=True
O.cell.setBox(1,1,1)
num_spheres=1000 # number of spheres
den_ball=2600 # density of particles from the experimental test
young=1e8
iso_pressure=-1e4

compFricDegree = 30 # initial contact friction during the confining phase

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

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),
             -1,0.333,num=num_spheres,
             periodic=True,seed=1)# final porosity should be 0.24
sp.toSimulation(material='spheres')

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),
 NewtonIntegrator(damping=0.2),
 PeriTriaxController(dynCell=True,maxUnbalanced=0.01,relStressTol=0.02,goal=(iso_pressure,iso_pressure,iso_pressure),stressMask=7,globUpdate=5,maxStrainRate=(.1,.1,.1),doneHook='print("Hydrostatic load reached."); O.pause();',label='triax'),
]
O.dt=PWaveTimeStep()
O.run();O.wait()

TW=TesselationWrapper()
TW.triangulate()
TW.setState()
TW.computeVolumes()

TW_Volume=0
for i in range(1000):
    TW_Volume += TW.volume(i)

print('volume='+str(O.cell.volume)+'\nTW_Volume='+str(TW_Volume))

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Chareyre
Solved:
2020-06-11
Last query:
2020-06-11
Last reply:
2020-06-08
Robert Caulk (rcaulk) said : #1

I suppose the periodic cell volume may not be estimated in the same way as the periodic tesselationwrapper. Do you notice the same inconsistency when you use non-periodic boundary conditions?

Chu (arcoubuntu) said : #2

Dear Robert,

Thank you for your suggestion. I used non-periodic boundary and I got TW_Volume=1.0451466106445113, when the volume of pack is 1.0. It seems better but error still exists.

In addition, I got different result when I set different young, such as

TW_Volume=1.045146610644511 #young=1e6
TW_Volume=1.000684823430585 #young=1e8
TW_Volume=1.000056824921447 #young=1e10

It's interesting, the larger the young, the smaller the error.

####MWE####

from yade import pack

num_spheres=1000# number of spheres
young=1e6
compFricDegree = 3 # initial contact friction during the confining phase
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

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,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=0.2)

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()],label="iloop"
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

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

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:
    break

TW=TesselationWrapper()
TW.triangulate()
TW.setState()
TW.computeVolumes()

TW_Volume=0
for b in O.bodies:
    if isinstance(b.shape,Sphere):
        TW_Volume += TW.volume(b.id)

print('TW_Volume='+str(TW_Volume))

Robert Caulk (rcaulk) said : #3

Thankyou, and then increasing the number of particles decreases or increases the discrepancy?

Robert Caulk (rcaulk) said : #4

Looking at the source code, it looks like the original author (Bruno?) added a FIXME [1] saying we still need to properly define volumes for spheres on the borders. Looking at the function where the volume is computed for each sphere [2], it looks like fictitious vertices are neglected in the calculation of the volumes. It makes sense since these vertices are usually placed extremely far from the packing in order to generate a boundary. Further, if a cell is fictitious (if it is a boundary condition cell) then it is also neglected from the calculation. Finally, I see nothing here to consider periodicity. This all means that any volume that you might expect to be accounted for along the boundary of a packing appears to not be calculated. I guess these edge geometries are computed in the saveVTK(withBoundaries=True) function, but it would need to be copied over to this algorithm.

At the end, you should probably just take a volume as a relative value and stick with one function or the other.

[1]https://gitlab.com/yade-dev/trunk/-/blob/master/lib/triangulation/Tesselation.ipp#L1251
[2]https://gitlab.com/yade-dev/trunk/-/blob/master/lib/triangulation/Tesselation.ipp#L1196

Best Chareyre (bruno-chareyre-9) said : #5

Hi,
They are expected differences.
For the periodic case the problem is just that you apply non-periodic
triangulation to a periodic problem.
Triangulation defines a bounding box, and you get the volume of this box.
It is obviously larger than a period since some bodies are overlapping the
split between periods.

For effect of Young (much smaller) the reason is the same, if the spheres
overlap the walls then the bounding box is also a bit larger.

FlowEngine implements a triangulation with more options, so that we can
impose where the walls are and get exact volume.

Regards
Bruno

Chu (arcoubuntu) said : #6

Thanks Chareyre, that solved my question.

Chu (arcoubuntu) said : #7

Hi Robert,
Thank you for your suggestions. You point out that TesselationWrapper shoudn't be used in non-periodic boundary conditions.

>increasing the number of particles decreases or increases the discrepancy?

When the volume of pack is 1.0
TW_Volume=1.05088660135025 #1000 spheres
TW_Volume=1.02397027562159 #10000 spheres
TW_Volume=1.01293294882138 #100000 spheres

Chu (arcoubuntu) said : #8

Hi Bruno,
Thank you for your suggestions. I got good result with FlowEnigne in periodic boundary conditions.

Robert Caulk (rcaulk) said : #9

 > You point out that TesselationWrapper shoudn't be used in non-periodic boundary conditions.

Nope. I said periodic conditions are not considered in TesselationWrapper.

Chu (arcoubuntu) said : #10

My fault, non-periodic -> periodic :-)