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

porosity calculation for periodic boundaries

0 votes
14 views
asked Oct 18 by MaxWiebicke (400 points)
edited Oct 18 by MaxWiebicke

Hi,

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

woo._utils2.porosity(S.dem,None,box=S.cell)

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])
v=S.cell.volume
porosity = (v-vs)/v

 

1 Answer

0 votes
answered Oct 23 by eudoxos (47,770 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 by MaxWiebicke (400 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 by MaxWiebicke (400 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
v=S.cell.volume
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
...