# Pore size distribution

Dear all,

I want to obtain the pore size distribution of a sphere packing.

Does anyone have any suggestions for this purpose?

Thank you so much in advance.

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2022-03-17: #1

Hello,

- what is "a sphere packing"? a MWE [1] would be excellent.
- what is "pore" w.r.t this sphere packing? For normal sphere packings, there is just connected "void net" and no "real pores".
- anything else, which would help somebody like us (having absolutely no experience with your problems, your ideas, etc.) to understand more the problem

Cheers
Jan

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-17: #2

Dear Jan;

Thank you very much for your kind reply. I am sorry for my vague question.

I am a new user of YADE. A part of my PhD project is investigation of the effects of compaction on pore structure and soil hydraulic properties. To this aim, we want to follow the DEM approach proposed by ”Mahmoodlu, M. G., Raoof, A., Sweijen, T., & Van Genuchten, M. T. (2016). Effects of sand compaction and mixing on pore structure and the unsaturated soil hydraulic properties. Vadose Zone Journal, 15(8).” which is done using YADE.

I am trying to reproduce the compaction process of a granular material under oedometric conditions based on my desired particle size distribution.

But at first, I need to have a function to obtain the pore size distribution of the media.

For sphere packing I am using the makeCloud function. For example, sphere packing in this code [1].

I am writing the pore size distribution definition for your consideration. “The pore size distribution is the relative abundance of each pore size in a representative volume of soil. It can be represented with a function f(r), which has a value proportional to the combined volume of all pores whose effective radius is within an infinitesimal range centered on r.”

For example in the mentioned paper, a regular triangulation method was used to extract the geometry of the pore structure which referred to the [1]. They obtained the pore body and pore throat.

 Revision history for this message Jan Stránský (honzik) said on 2022-03-18: #3

Hello,

The article describes the pore evaluation of pores in "Calculation of Pore Structure" section well.

As I have no experience with this topic, I don't know if there are some predefined tools in Yade for this.
But it should not be difficult with ordinary Python.

Cheers
Jan

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-19: #4

Dear Jan;

They used this approach for obtaining permeability. I do not know if I can obtain the pore size distribution from that or not.

I hope there will be a predefined tool in yade.

Cheers
Soheil

 Revision history for this message Jan Stránský (honzik) said on 2022-03-19: #5

> They used this approach for obtaining permeability. I do not know if I can obtain the pore size distribution from that or not.

They also defined pores.
You (Delaunay?) tetrehdralize your packing, centers of spheres define vertices of tetrahedrons.
Each tetrahedron represents one pore.
Pore (size) is defined as inscribed sphere (size) between four "vertex" spheres (each sphere in on vertex).
Once you have the pores and sizes, size distribution is easy.

> I hope there will be a predefined tool in yade.

If there is no such tool, you do your job and still hope in this, you are free to publish your code :-)

Cheers
Jan

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-21: #6

Dear Jan;

Thank you very much for your valuable information.

I will work on it. But I am not good enough in Python, also I have a short time for this task.

There is a useful image analysis tool (PoreSpy) which is based on python. It has a predefined function for obtaining the pore size distribution.

Can I export Yade images and use this tool? In which formats?

I would really appreciate it if you could provide me some information regarding this matter.

Thank you so much in advance.

cheers,
Soheil

 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-21: #7

Hello,

I will assume Jan has found the quantity that you search (his definition seems reasonable to me): the inscribed sphere radius between 4 vertices of a Delaunay triangulation. TwoPhaseFlowEngine allows you to check that value [5]. However, this is not implemented in FlowEngine (i.e. oedometer.py). If you need this function in FlowEngine, you can make a request for it on gitlab [6].

FlowEngine does, however, give you the ability to print all the vertices [4], which you can post-process (using python) in the same way as [5] to get the inscribed sphere radius distribution.

>>There is a useful image analysis tool (PoreSpy) which is based on python

I am unfamiliar with PoreSpy, but someone else in this forum was trying to do something similar by converting a Yade simulation to pixelated image so he could do pore geometry interpretation [1].

Cheers,

Robert

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-21: #8

Dear Robert;

Thank you very much for your valuable information.

I am sorry for that, It is just because I am fresh in yade.

Yes, Jan read the article which is very valuable for me. I am really grateful for your kind help.

Unfortunately, there is no specific definition for “pore”. I mentioned that article because they used the yade for their simulation.

I wanted to know if anyone knows about any predefined function for this purpose or not. Because I could not find any through the document and website.

I will ask for this function in FlowEngine on gitlab.

Many thanks for your useful recommendations.

Cheers,
Soheil

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-21: #9

Thanks Robert Caulk, that solved my question.

 Revision history for this message Jan Stránský (honzik) said on 2022-03-21: #10

> But I am not good enough in Python, also I have a short time for this task.

time to improve :-)
short time is not the best conditions, thought..

> There is a useful image analysis tool (PoreSpy) which is based on python. It has a predefined function for obtaining the pore size distribution.
> Can I export Yade images and use this tool? In which formats?

Yes, you it is possible to do it like this, but DO NOT do it like this, it is too screwed.
You have data and you can directly get what you want, data -> result.
Do not go direction data -> image (already some loss) -> image analysis (uncertain results, further loss, ...) -> result.

I will try to give some MWE soon
Cheers
Jan

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2022-03-21: #11

Hi Robert and Jan,

I am working on similar problem, I thought it would be better to continue in this question instead of opening a new.
I am trying implement the solution suggested by Robert in #7 by using print all the vertices [4].
I have compacted spherical particles into a cylinder shape, and want estimates/measure pore size in the compact. I was running the suggestion by print all vertices.

By including following MWE after my existing code:

FlowEngine.printVertices()

I am getting following error message:

"ArgumentError: Python argument types in
FlowEngineT.printVertices()
did not match C++ signature:
"

My question how should implement this code to extract the vertices and type of the void/pore space?

 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-21: #12

Hello,

>>I am working on similar problem, I thought it would be better to continue in this question instead of opening a new.

It is better to open a new thread. Especially considering this one has been "solved" already. So please open a new one, and in doing so please attach a full MWE (See [1]). Thank you for helping keep this knowledge base organized and sustainable, and thanks for giving us the necessary information we need to help you.

Cheers,

Robert

 Revision history for this message Jan Stránský (honzik) said on 2022-03-21: #13

Cheers
Jan

###
import scipy.spatial

# https://www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm
# filter out "degenerate" tetrahedrons
def delaunay2poreTetras(delaunay,limit=0.1):
vertices = delaunay.points
tetras = delaunay.simplices
def quality(tetra):
vs = v0,v1,v2,v3 = [Vector3(vertices[i]) for i in tetra]
d1 = v1 - v0
d2 = v2 - v0
d3 = v3 - v0
v = d1.dot(d2.cross(d3)) / 6
v = abs(v) # volume
a = d1.norm()
b = d2.norm()
c = d3.norm()
A = (v2-v3).norm()
B = (v1-v3).norm()
C = (v1-v2).norm()
aa,bb,cc = a*A, b*B, c*C
den = (aa+bb+cc)*(aa+bb-cc)*(aa-bb+cc)*(-aa+bb+cc)
den = abs(den)
r = sqrt(den)/(24*v) # radius of circum sphere
q = (9*sqrt(3)/8*v)**(1/3.)/r # quality in range [0,1]
return q
return [t for t in tetras if quality(t) > limit]

# some home-made inefficient (but working :-) iterative approach
def tetra2pore(tetra):
bodies = [O.bodies[int(i)] for i in tetra]
vertices = [b.state.pos for b in bodies]
center = 0.25*sum(vertices,Vector3.Zero) # initial guess, center of mass of tetrahedron
for i in range(100): # TODO stop condition
ds = d1,d2,d3,d4 = [(center-v).norm()-r for v,r in zip(vertices,radii)] # distances from "touching"
dMin,dMax = min(ds),max(ds)
iMin,iMax = ds.index(dMin), ds.index(dMax)
vMin,vMax = vertices[iMin],vertices[iMax]
dirMin = (center - vMin).normalized()
dirMax = (center - vMax).normalized()
# based on min and max distance "push" the center a bit
dd = (dMax - dMin) * 0.1
center += dirMin * dd - dirMax * dd

# testing packing, colored red
O.bodies.append((
))

centers = [b.state.pos for b in O.bodies]
delaunay = scipy.spatial.Delaunay(centers) # tetrahedralization
tetras = delaunay2poreTetras(delaunay) # only "good" tetrahedrons
pores = [tetra2pore(t) for t in tetras] # pores out of tetrahedrons

# pores, colored cyan
for c,r in pores:

# tetrahedrons as facets, for visualization
for ii in tetras:
vs = [delaunay.points[i] for i in ii]
v1,v2,v3,v4 = [Vector3(list(v)) for v in vs]
O.bodies.append((
))

# improve sphere display quality
###

 Revision history for this message Soheil Safari (soheilsafari) said on 2022-03-22: #14

Dear Jan;

That is great. Thank you very much for your valuable time spent on this useful code.

I really appreciate your kind effort in solving my issue.

Best regards,
Soheil