# Porosity overlapping particles

Asked by Mithushan Soundaranathan

Hi,
I have compacted spherical particles in a cylinder, and the particle are overlapping. I am struggling to measure the porosity of my compacys. I have tried utils.porosity, voxelporosity and also tried using the usual porosity Eq. (1- (rho_compacts/rho_t)). They are giving me either negative values upto (-60%) or to high values (without accounting for the overlap). I have looked through the forum for similar problems, I have not any solution which applicable for my simulation. I was wondering if anyone have experienced similar issues and how they solved it.

Best,
Mithu

Here is my code

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export

o = Omega()

# Physical parameters
fr = 0.54
rho = 1050
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 126000
kp = 12.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999

Tab_rad=0.005
Cyl_height=0.02

Comp_press1=1.4e8
Comp_force1=Comp_press1*cross_area

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression
sp=pack.SpherePack()
sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-8*Diameter),(3.0*Diameter,3.0*Diameter,28.0*Diameter), rMean=Diameter/2.0,num=527)
#cyl = pack.inCylinder((0,0,0),(0,0,40*Diameter),40*Diameter-0.006)
#sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)

######################################################################
#O.bodies.append(geom.facetBox((0,0,0), (4.0*Diameter,4.0*Diameter,4.0*Diameter), wallMask=63-32, material=mat1))
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

# Add engines
o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
#VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=10000),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]

def checkForce():
if O.iter < 1000000:
return
timing.reset()
if unbalancedForce() > 0.2:
return
# add plate at upper box side

highSphere = 0.0
for b in O.bodies:
if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
highSphere = b.state.pos[2]
else:
pass

O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
global plate
plate = O.bodies[-1]
plate.state.vel = (0, 0, -.3)
fCheck.command = 'unloadPlate()'

def unloadPlate():
if abs(O.forces.f(plate.id)[2]) > Comp_force1:
plate.state.vel *= -1
fCheck.command = 'stopUnloading()'

def stopUnloading():
if abs(O.forces.f(plate.id)[2]) == 0:
O.pause()

## Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
 Revision history for this message Jérôme Duriez (jduriez) said on 2021-05-18: #1

Hi,

It is somewhat unusual in YADE to have such overlaps that would give negative porosities, when computing those the usual way. Are you sure you want to have these overlaps in your simulations (with respect to your application) ?

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2021-05-18: #2

Hi,
I agree with Jérôme: this is a weird case.
First question is: how would you define the "correct" porosity with words/equations in such case?
How to compute it, is only a secondary question.
Regards
Bruno

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-05-24: #3

Hi,

I would define correct porosity as the total volume of void space in the compacts.

In the real world pharma compaction process these overlap exist.

Best,
Mithu

 Revision history for this message Jan Stránský (honzik) said on 2021-05-25: #4

Hello,

please provide a MWE [1] (W = working), here I got: NameError: name 'cross_area' is not defined

For this case, the running is not important at all, same for materials. A MWE is the compacted state (a few spheres with hard-coded coordinates, possibly walls/facets if they are important to define the volume).
Pre-defined spheres has the advantage, that you can set the spheres such that you know the real porosity.

> I have tried utils.porosity, voxelporosity and also tried using the usual porosity Eq. (1- (rho_compacts/rho_t)).

please provide how you use it. It is not present in the code, an option is that you do not use it correctly.

Cheers
Jan

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-05-25: #5

Hi Jan,

I defined cross_Area as:
cross_area=math.pi*(Tab_rad**2)

The code I used to calculate porosity with equation:
m_tot=particleMass*527,
V=math.pi*((0.5*utils.aabbDim[0])**2)*utils.aabbDim[2]
porosity=1-(m_tot/(V*rho))

 Revision history for this message Jan Stránský (honzik) said on 2021-05-25: #6

> I defined cross_Area as:
> cross_area=math.pi*(Tab_rad**2)

Ok, but I run the simulation for almost an hour and no signal that it wants to end.. it is really far from a MWE.
Please provide the final packing with boundary, that is the true MWE :-)

and the other porosities? this task sounds like ideal for voxelPorosity.. maybe it works only for box shape?

Cheers
Jan

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2021-05-25: #7

> I would define correct porosity as the total volume of void space in the compacts.

Right. But what what would be the correct one _in a model_ when particles overlap?

> In the real world pharma compaction process these overlap exist.

This is confusing me. I use to teach that overlap do not exist in the real world, and that the numerical overlaps are just because we don't display the deformed shapes.
1 million\$ question is: how do we know porosity when we don't know the actual deformed shapes? It is up to you to decide.
With hydrogels, for instance, we know that the volume of each particle is constant during the deformation, it makes the problem easy: solid volume is constant (and then voxelPorosity is wrong). Does it apply in your case, I don't know, but in general I would say voxelPorosity is not a good approach for highly deformable particles.

Bruno

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-05-25: #8

Hi Jan,

try this cod instead, which is much more faster,

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]

from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export

o = Omega()

# Physical parameters
fr = 0.54
rho = 1050
Diameter = 0.0012
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 126000
kp = 12.0*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

o.dt = 1.0e-7
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho

Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Tab_rad=0.005
Cyl_height=0.045
#Comp_press=0.22e8

cross_area=math.pi*(Tab_rad**2)
#Comp_force=Comp_press*cross_area

Comp_press1=1.4e8
Comp_force1=Comp_press1*cross_area

#*************************************

# Add material
mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression

sp=pack.SpherePack()
sp.makeCloud((-3.0*Diameter,-3.0*Diameter,-18*Diameter),(3.0*Diameter,3.0*Diameter,18.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=527)
#cyl = pack.inCylinder((0,0,0),(0,0,40*Diameter),40*Diameter-0.006)
#sp = pack.filterSpherePack(cyl,sp,True,material=mat1)
sp.toSimulation(material=mat1)

######################################################################
#O.bodies.append(geom.facetBox((0,0,0), (4.0*Diameter,4.0*Diameter,4.0*Diameter), wallMask=63-32, material=mat1))
walls=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))

# Add engines
o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
VTKRecorder(fileName='vtk-',recorders=['all'],iterPeriod=95000),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]

def checkForce():
# 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 < 10000:
return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
highSphere = 0.0
for b in O.bodies:
if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
highSphere = b.state.pos[2]
else:
pass

O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1, material=mat1))
# without this line, the plate variable would only exist inside this
# function
global plate
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, -5)
# start plotting the data now, it was not interesting before
# next time, do not call this function anymore, but the next one
# (unloadPlate) instead
fCheck.command = 'unloadPlate()'

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

def stopUnloading():
if abs(O.forces.f(plate.id)[2]) == 0:
# 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()

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-05-25: #9

Hi Bruno,

Sorry, my fault. You are absolutely right. I was meaning to say deformation, not much overlap in the compacts.

The right model, I not sure. It is tricky, not sure which method gives the "correct" porosity. An idea would be to divided volumes of particles with total volume of the compacts.

Best,
Mithu

 Revision history for this message Jan Stránský (honzik) said on 2021-05-25: #10

Consider your physics, that should be the first criterion.

From purely mathmetical point of view, you can:
1)
assume that not more than two spheres are overlapping, then from the knowledge of sizes of all spheres and penetration depths of all interactions it should be an easy task to compute total solid volume
2)
compute the solid volume accurately as a volume of union of spheres. Have a look e.g. at paper "Computing the Volume of a Union of Balls: A Certified Algorithm"

Cheers
Jan

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2021-05-25: #11

If you have no clear idea, I think best by default is to assume 4/3*πR^3 is still a good approximation for the volumes of the deformed spheres (isochoric deformation). No need to go compute intersections since they are not physical intersections.

If it gives negative values then the modelling approach in itself is questionable and there is no magic trick to fix it.
In some advanced models [1] we tried to factor local porosity in the pairwise interactions, so that porosity would never go negative.
It hasn't been maintained unfortunately, but the concepts could be re-implemented.

Bruno

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2021-05-26: #12

Hi Bruno,

Thank you very much, do you mean to use Eq.4 in the paper?

Best,
Mithu

## Can you help with this problem?

Provide an answer of your own, or ask Mithushan Soundaranathan for more information if necessary.

To post a message you must log in.