Welcome to ask.woodem.org. You may post when you login through your GitHub account.

porosity calculation for periodic boundaries

0 votes
asked Oct 18, 2018 by MaxWiebicke (500 points)
edited Oct 18, 2018 by MaxWiebicke


is there a way to precisely calculate the porosity in a box with periodic boundaries? I tried the commands 


but it wouldn't take my input cell. In your examples you use an AlignedBox3 for the calculation.

Another question is how you compute the porosity in the box -- do you consider the overlap at the contacts?

Right now, I am calculating the porosity with, which does not consider the overlap:

vs=sum([p.shape.volume for p in S.dem.par])
porosity = (v-vs)/v


1 Answer

0 votes
answered Oct 23, 2018 by eudoxos (48,570 points)
Hi Max, this is correct, there is no function to consider overlaps. A few reasons for that: (a) no-one has written it yet, (b) it will work only for shape couples that are implemented and their number is combinatorial (e.g. sphere-sphere but not sphere-capsule, unless also implemented etc), (c) you will have to neglect 3-overlaps as those are really messy in terms of the geometry involved, IIRC. What you could do is just write a simple function to compute sphere overlap volume (takes 2 radii and distance, returns volume ) and then sum that over all contacts, subtracting from *vs*. I might also add it to the codebase, that is for consideration. What do you think? v
commented Oct 24, 2018 by MaxWiebicke (500 points)

Yes, that's what I thought.
I'll try to do it for spheres, after I am back in Dresden in 2 weeks. I am in Sydney working with Itai right now ;)

commented Oct 26, 2018 by MaxWiebicke (500 points)

Hey Vaclav, I did write a small function and tested it on the stress range in which I perform my simulations. the influence of the overlap on the solid volume is minimal -- if overlaps get larger it might be of importance however.

this is the small python script I used:

def porosityOverlap(S):
# considering overlap at the contacts
vs=sum([p.shape.volume for p in S.dem.par])    
# for every contact we need to substract the overlap once
vOall = 0
for C in S.dem.con:
    rA = S.dem.par[C.pA.id].shape.radius
    rB = S.dem.par[C.pB.id].shape.radius
    # calculate the distance from overlap uN 
    #   otherwise problems at the periodic boundaries
    dAB = rA + rB + C.geom.uN
    #calculate overlap
    # https://en.wikipedia.org/wiki/Spherical_cap
    vO = math.pi/(12*dAB) * (rA+rB-dAB)**2 * (dAB**2 + 2*dAB*(rA+rB) - 3*(rA-rB)**2 )
    vOall += vO

vs = vs - vOall
return (v-vs)/v