the definition of porosity

Hi Jan,

Thank you for spending so much time with me. I don't know how to repay you.Thank you again.

You said Yade (DEM in general) considers particles as perfectly rigid. But when I run the following Oedometric test script(from https://yade-dem.org/doc/tutorial-examples.html).I found that the wall would be compressed into the upper ball, and in addition, the ball would overlap with the ball.So i think it is soft ball model.I may have a problem with my understanding. Please point out my mistake. Thank you.

Feng

####################################
# gravity deposition, continuing with oedometric test after stabilization
# shows also how to run parametric studies with yade-batch

# The components of the batch are:
# 1. table with parameters, one set of parameters per line (ccc.table)
# 2. readParamsFromTable which reads respective line from the parameter file
# 3. the simulation muse be run using yade-batch, not yade
#
# $ yade-batch --job-threads=1 03-oedometric-test.table 03-oedometric-test.py
#

# load parameters from file if run in batch
# default values are used if not run from batch
readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
# make rMean, rRelFuzz, maxLoad accessible directly as variables later
from yade.params.table import *

# create box with free top, and ceate loose packing inside the box
from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean,rRelFuzz=rRelFuzz)
sp.toSimulation()

O.engines=[
   ForceResetter(),
   # sphere, facet, wall
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      # the loading plate is a wall, we need to handle sphere+sphere, sphere+facet, sphere+wall
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
   # 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()',realPeriod=2,label='checker'),
]
O.dt=.5*PWaveTimeStep()

# the following checkUnbalanced, unloadPlate and stopUnloading functions are all called by the 'checker'
# (the last engine) one after another; this sequence defines progression of different stages of the
# simulation, as each of the functions, when the condition is satisfied, updates 'checker' to call
# the next function when it is run from within the simulation next time

# check whether the gravity deposition has already finished
# if so, add wall on the top of the packing and start the oedometric test
def checkUnbalanced():
   # at the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable
   if O.iter<5000: return
   # the rest will be run only if unbalanced is < .1 (stabilized packing)
   if unbalancedForce()>.1: return
   # add plate at the position on the top of the packing
   # the maximum finds the z-coordinate of the top of the topmost particle
   O.bodies.append(wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
   global plate # without this line, the plate variable would only exist inside this function
   plate=O.bodies[-1] # the last particles is the plate
   # Wall objects are "fixed" by default, i.e. not subject to forces
   # prescribing a velocity will therefore make it move at constant velocity (downwards)
   plate.state.vel=(0,0,-.1)
   # start plotting the data now, it was not interesting before
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   # next time, do not call this function anymore, but the next one (unloadPlate) instead
   checker.command='unloadPlate()'

def unloadPlate():
   # if the force on plate exceeds maximum load, start unloading
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      # next time, do not call this function anymore, but the next one (stopUnloading) instead
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      # O.tags can be used to retrieve unique identifiers of the simulation
      # if running in batch, subsequent simulation would overwrite each other's output files otherwise
      # d (or description) is simulation description (composed of parameter values)
      # while the id is composed of time and process number
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()

O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()

##############################

> I found that the wall would be compressed into the upper ball, and in addition, the ball would overlap with the ball

yes. In reality (if the particles were soft), there would be no overlap but the particles would deform (being in contact, but not overlapping).
In Yade, the particles are rigid and do overlap. Repulsive forces are somehow computed from this overlap.
Porosity is computed using the original particles volume. This way, the computation is very easy and fast.

From existing overlaps, you can relatively easily compute actual volume of overlaps and adjust the evaluation if needed.

Also please, if the question is about porosity, next time provide the code you use to compute porosity. There are two of them, [1] using the computation you described, the other [2] using voxel approximation (but for your case, computing actual overlaps is not that difficult and much more precise).

Jan

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

Hi, Jan

Ask another question.When I run the code above. I increased maxLoad=1e6 to maxLoad=1e8. Why do some little balls go through walls? Is it because the wall goes through the center of the ball? I don't really understand how this works.

Thanks，
Feng

> Is it because the wall goes through the center of the ball?

Yes. The displacement is proportional to force, increasing force more and more thus leads bodies to cross each other.
Bruno

Thanks Jan Stránský, that solved my question.

Thanks Bruno too.
Feng

Dear Bruno,

is there a contact model exists in YADE which can handle non-overlapping spheres? I mean, is any solution to increase the compressional force while the spheres are not cross each other.

Thanks in advance, Feri

Dear Feri,

Yes there is solution for that. Just have a look at the script below which describes a uniaxial test with the JCFPM. More specifically, look at the intR parameter and how it is used.

Luc

ps: next time, please open a dedicated question.

----

from yade import ymport, plot

################# SIMULATIONS DEFINED HERE

#### pre-existing packing
PACK='box_112_5659.spheres'

#### Simulation Control
rate=-0.01 #deformation rate
iterMax=60000 # maximum number of iterations
saveData=int(iterMax/1000) # data record interval
saveVTK=int(iterMax/5.) # saving output files for paraview
OUT='uniaxialCompressionTest_0.01'

#### Microproperties (interparticle parameters)
DENS=3000 # this one can be adjusted for different reasons (porosity of packing vs porosity of material / increase time step (no gravity -> no real effect on the result)

intR=1.2 # allows near neighbour interaction (can be adjusted for every packing / the bigger -> the more brittle / careful when intR is too large -> bonds can be created "over" particles) -> intR can be calibrated to reach a certain coordination number K (see calculation on line 115)
YOUNG=10e9 # this one controls the Young's modulus of the material
ALPHA=0.15 # this one controls the material Poisson's ratio of the material
TENS=10e6 # this one controls the tensile strength UTS of the material
COH=10e6 # this one controls the compressive strength UCS of the material, more precisely, the ratio UCS/UTS (from my experience: COH should be >= to TENS, >= 10*TENS for competent materials like concrete)
FRICT=20 # this one controls the slope of the failure envelope (effect mainly visible on triaxial compression tests)

#### example values
### granite -> needs K=13
#YOUNG=68e9
#ALPHA=0.333
#TENS=8e6
#COH=16e7
#FRICT=10

### Fontainebleau sandstone -> needs K=10
#YOUNG=50e9
#ALPHA=0.25
#TENS=45e5
#COH=45e6
#FRICT=18

#### claystone -> needs K=8
#YOUNG=10e9
#ALPHA=0.1
#TENS=12e6
#COH=12e6
#FRICT=7

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,tensileStrength=TENS,cohesion=COH,frictionAngle=radians(FRICT))

#### import pre-existing specimen
O.bodies.append(ymport.text(PACK,scale=1.,shift=Vector3(0,0,0),material=sphereMat))

R=0
Rmax=0
Rmin=1e6
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.7,0.5,0.3)
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
   if o.shape.radius<Rmin:
     Rmin=o.shape.radius
Rmean=R/nbSpheres

print('nbSpheres=',nbSpheres,' | Rmean=',Rmean, ' | Rmax/Rmin=', Rmax/Rmin)

#### help define boundary conditions (see utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# DEM loop + ENGINES DEFINED HERE

O.engines=[ForceResetter(),
           InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
           InteractionLoop(
               [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
               [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
               [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=False,Key=OUT,label='interactionLaw')]
           ),
           UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
           GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep()),
           NewtonIntegrator(damping=0.5,label='newton'),
           PyRunner(iterPeriod=saveData,initRun=True,command='recorder()',label='data'),
           #VTKRecorder(iterPeriod=1,initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks','bstresses'],Key=OUT,dead=1,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
        'eps':strainer.strain,
        'sigma':strainer.avgStress,
        'tc':interactionLaw.nbTensCracks,
        'sc':interactionLaw.nbShearCracks,
        'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
#### initializes the interaction detection factor
#SSgeom.interactionDetectionFactor=-1.
#Saabb.aabbEnlargeFactor=-1.

#### special treatment for tension tests: reinforcement of bonds close to boundaries
if rate>0:
    dim=aabbExtrema()
    layerSize=0.15 # size of the reinforced zone
    for o in O.bodies:
        if isinstance(o.shape,Sphere):
            if ( o.state.pos[longerAxis]<(dim[0][longerAxis]+layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) or ( o.state.pos[longerAxis]>(dim[1][longerAxis]-layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) :
                o.shape.color=(1,1,1)

#### coordination number calculation
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if not i.isReal : continue
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
      # FOR TENSION TEST: these lines reinforce the bonds near the boundaries to avoid rupture along the boundaries
      if O.bodies[i.id1].shape.color==(1,1,1) or O.bodies[i.id2].shape.color==(1,1,1) :
 i.phys.FnMax*=100
 i.phys.FsMax*=100
    if i.phys.isCohesive :
      numCohesivelinks+=1
print ("K=", 2.0*numCohesivelinks/nbSpheres)

#vtk.dead=0
#O.step()
#vtk.iterPeriod=saveVTK

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
O.run(iterMax)

Dear Luc,

thank you for your quick response. Sorry, it seemed to me, that in this way it is more easier to explain what I would like to do.
Can you offer a detailed description about JCFP model? Until this time I used only linear elastic or Hertz-Mindlin contact law to analyse motion of cohesionless particulates. What do you think, is JCFP usable for modeling sandstone particles?

Thank you, Feri

Hi,

JCFPM is suitable for modeling rocks [1] and fractured rocks [2].

What do you want to model exacly? -> Please open a dedicated question since we are very far from "
the definition of porosity"

Luc

[1] https://www.sciencedirect.com/science/article/abs/pii/S0022509612002268
[2] https://www.sciencedirect.com/science/article/pii/S1365160912000391

Hi Ferenc,
There's also this feature in the laws using interaction physics CohFrictPhys, and also in other more exotic functors (e.g. VirtualLubricationPhys).
But Luc is right, it looks like a new question.
Bruno

Dear Bruno and Luc,

thank you for your help, it is really appreciated! The problem was solved based on Luc's comment. Moreover after some research I found out that aabbEnlargeFactor is working with CohFrictMat as well, however I had to increse the Young's modulus of the particles. Unfortunately the simulation is slower now, but the spheres are not intersecting.

My question is related a bit to porosity, as in case high pressure, because of many intesections the built-in porosity function gave me wrong result (it was checked by exporting the confined geometry into CAD).

Thanks a lot!
Feri