Pore size distribution

Hello,

please read [1] and provide more information. Namely:
- 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

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

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.

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/FluidCouplingPFV/oedometer.py

Hello,

thanks for more information.
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

Dear Jan;

Thank you very much for your reply.

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

> 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

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

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

Thanks Robert Caulk, that solved my question.

> 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

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(dead=1,label="flow")
FlowEngine.printVertices()

I am getting following error message:

"ArgumentError: Python argument types in
    FlowEngineT.printVertices()
did not match C++ signature:
    printVertices(yade::TemplateFlowEngine_FlowEngineT<yade::FlowCellInfo_FlowEngineT, yade::FlowVertexInfo_FlowEngineT, yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> >, yade::CGT::FlowBoundingSphereLinSolv<yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> >, yade::CGT::FlowBoundingSphere<yade::CGT::_Tesselation<yade::CGT::TriangulationTypes<yade::FlowVertexInfo_FlowEngineT, yade::FlowCellInfo_FlowEngineT> > > > > {lvalue})
"

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

[4]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.printVertices

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

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

Generating pores "home-made"
Cheers
Jan

###
import scipy.spatial

# https://www2.mps.mpg.de/homes/daly/CSDS/t4h/tetra.htm
# https://math.stackexchange.com/questions/2820212/circumradius-of-a-tetrahedron
# 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]
    radii = [b.shape.radius 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
    radius = (center - vertices[0]).norm() - radii[0]
    return center,radius

# testing packing, colored red
O.bodies.append((
    sphere((1.0,2.0,3.0),1.0,color=(1,0,0),mask=0b001),
    sphere((4.0,2.0,3.0),1.6,color=(1,0,0),mask=0b001),
    sphere((2.5,4.0,3.0),1.2,color=(1,0,0),mask=0b001),
    sphere((2.6,3.0,5.0),1.1,color=(1,0,0),mask=0b001),
    sphere((5.0,4.0,5.0),1.4,color=(1,0,0),mask=0b001),
    sphere((5.0,5.0,2.5),1.3,color=(1,0,0),mask=0b001),
    sphere((6.5,3.0,3.0),1.2,color=(1,0,0),mask=0b001),
    sphere((7.4,4.8,4.0),1.2,color=(1,0,0),mask=0b001),
))

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:
    O.bodies.append(sphere(c,r,color=(0,1,1),mask=0b010))

# 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((
        facet((v1,v2,v3),mask=0b100),
        facet((v1,v2,v4),mask=0b100),
        facet((v1,v3,v4),mask=0b100),
        facet((v2,v3,v4),mask=0b100),
    ))

# improve sphere display quality
yade.qt.Gl1_Sphere.quality = 2
###

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