Porosity overlapping particles

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

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

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

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

[1] https://www.yade-dem.org/wiki/Howtoask

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

> 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

> 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

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

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

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

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

[1] https://www.sciencedirect.com/science/article/abs/pii/S0032591010004481

Hi Bruno,

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

Best,
Mithu