Interparticle spacing

Hello,

I would like to compute the spacing between adjacent individual particles (spheres) in a given packing (i.e. interparticle spacing). To do that, I thought about getting the length of the sides from the Delaunay triangulation process in the FlowEngine. So I created a packing of spheres, then run the flow engine to do the triangulation, then I want to get the length of the side. My code is shown below and I am using Yadedaily in running this code. My questions are:

1- Is there a more efficient approach to get the interparticle spacing?
2- I found two functions in Yade documentation [1]: getVertices and printVertices, are those vertices belongs to the Delaunay triangulation (i.e. the vertices of the tetrahedra that forms the mesh for flow modeling)?
3- What are the values that I get from getVertices function? are they the Ids of the cells?
4- printVertices gives a txt file with six columns: id, x, y, z, alpha, fictitious. What does alpha and id mean?

Thank you,
Othman

-----------

P=1 #Pa
visc=1e-3 #Pa.sec - taken from Catalano, Chareyre, Bathelemy (2014)
density=1000 #kg/m3
g=9.81 #m/s2

########## create walls ##########

mnx=0
mny=0
mnz=0
mxx=50.8/1000
mxy=50.8/1000
mxz=50.8/1000

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0)
wallIds=O.bodies.append(walls)
########## spheres ##########
sp=pack.SpherePack()
sp.toSimulation()
Height=max(utils.aabbDim())
dP=P/Height #Pa/m
##check how to fix the particles
#for i in sp:
# body= O.bodies[i]
# body.state.blockedDOFs='xyzXYZ'

print ('porosity = ', utils.porosity())

########## flow engine ##########

flow=FlowEngine()
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=visc
flow.bndCondIsPressure=[0,0,0,0,0,1]
flow.bndCondValue=[0,0,0,0,0,P]
flow.boundaryUseMaxMin=[1,1,1,1,1,1]
O.dt=1e-6
O.dynDt=False
flow.emulateAction()

#get verticies positions
flow.printVertices()

n=flow.nCells()
x=[]
for i in range (n):
l=flow.getVertices(i)
x.append(l)

print (x)

Question information

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

>1- Is there a more efficient approach to get the interparticle spacing?

Yes. Iterate on interactions and compute distance between particle centers.

>2- I found two functions in Yade documentation [1]: getVertices and printVertices, are those vertices belongs to the Delaunay triangulation (i.e. the vertices of the tetrahedra that forms the mesh for flow modeling)?

Probably yes, but you are linking to the entire doc. Better to link to individual functions like this [1] so we know what you are referencing.

>3- What are the values that I get from getVertices function? are they the Ids of the cells?

getVertices gives you the 4 vertex ids comprising the argument cell.

>4- printVertices gives a txt file with six columns: id, x, y, z, alpha, fictitious. What does alpha and id

id - vertex id
alpha - if the vertex is alpha boundary or not. Special boundary type. If you dont know what this means then you don't need to worry about it.

 Revision history for this message Othman Sh (othman-sh) said on 2021-01-25: #2

Hi Robert,

Thanks for your answers. Would you please elaborate on #1 (Iterate on interactions and compute distance between particle centers)? how do I get the ids of adjacent/interacting bodies?

Thanks,
Othman

 Revision history for this message Robert Caulk (rcaulk) said on 2021-01-25: #3

You would iterate on the interaction container and query the constituent body ids:

separations =[]

for i in O.interactions:
id1,id2 = i.id1,i.id2
x = np.linalg.norm(O.bodies[id1].state.pos - O.bodies[id2].state.pos)
separations.append(x)

print(separations)

Cheers,

Robert

 Revision history for this message Othman Sh (othman-sh) said on 2021-01-25: #4

Thanks Robert. I did the code below but strangely I am not getting any outcome. The separation array is still empty!.

I am copying what I get in the terminal

Using python version: 3.8.5 (default, Jul 28 2020, 12:59:40)
[GCC 9.3.0]
TCP python prompt on localhost:9000, auth cookie `sycsuk'
XMLRPC info provider on http://localhost:21000
Running script inter.py
porosity = 0.7337488803402193
number of interactions = 3776
[]
[[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]]

My code is below.

Thanks
Othman
--------------------------------

import numpy as np

mnx=0
mny=0
mnz=0
mxx=50.8/1000
mxy=50.8/1000
mxz=50.8/1000

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)

########## spheres ##########
sp=pack.SpherePack()
sp.toSimulation()

print ('porosity = ', utils.porosity())

O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[
Ig2_Sphere_Sphere_ScGeom(),
Ig2_Facet_Sphere_ScGeom(),
],
[

Ip2_FrictMat_FrictMat_FrictPhys(),
],
[

Law2_ScGeom_FrictPhys_CundallStrack(),
],
),

NewtonIntegrator(gravity=(0,0,0),damping=.4),

]
O.step() #creating interactions
print('number of interactions = ',len(O.interactions))

separations =[]

for i in O.interactions:
id1,id2 = i.id1,i.id2
x = np.linalg.norm(O.bodies[id1].state.pos - O.bodies[id2].state.pos)
separations.append(x)

print(separations)

 Revision history for this message Jérôme Duriez (jduriez) said on 2021-01-26: #5

Maybe your interactions in O.interactions are not real. Non-real interactions are accounted for if you ask for len(O.interactions), but they are skipped in the for i in O.interactions: loop

Also note that for any Vector3 "vec" in yade terminal, like O.bodies[id1].state.pos - O.bodies[id2].state.pos, you have directly access to a .norm() function: vec.norm()

 Revision history for this message Jan Stránský (honzik) said on 2021-01-26: #6

Hello,

> I would like to compute the spacing between adjacent individual particles ... getting the length of the sides from the Delaunay triangulation process in the FlowEngine.

Isn't FlowEngine an overkill here?
Would "plain" Delaunay triangulation / Voronoi tessellation (using e.g. numpy) be an option?

cheers
Jan

 Revision history for this message Othman Sh (othman-sh) said on 2021-01-26: #7

Hello all,

As an update about using Delaunay triangulation of the FlowEngine, printVerticies [1] gave me a list of the vertices but they are not grouped for each tetrahedron (i.e Delaunay tetrahedon). Instead, they are listed based on the bodies ids. So I cannot use them to calculate the distances between the adjacent particle centers because I don't know if body #1 is actually the neighbor of body #2.

To solve this issue, I used getCellVolume [2] to compute the volume of the cells which I think it is the volume of the tetrahedrons. Then from each tetrahedron volume, I can get the length of the sides of that tetrahedron (volume= side length^3 / 6*sqrt(2))) which should be equal to the distance between the centers of the particles forming the tetrahedron. My question here is:

- Does "cell" in getCellVolume means the tetrahedron element? I just want to confirm if my approach is correct
- Is there any error in this approach that I am not aware of?

Jerome,
That seems right because I checked O.interactions.countReal() and this returned zero. Is this because there is no contact between the spheres? how to make the interactions real in this case? I want to count the interparticle spacing for the initial packing state (i.e. before running any simulation).

Jan,

That would be a good idea too. However, I couldn't find a function that can give the positions of the Delaunay tetrahedrons corners. There is computeVolumes [3] to compute the volume of Voronoi cells but this will not be useful because I think Voronoi cells are not the same as the Delaunay tetrahedron (right?), so I can't use that in calculating the interparticle spacing.

Thanks,
Othman

 Revision history for this message Jan Stránský (honzik) said on 2021-01-26: #8

> I couldn't find a function that can give the positions of the Delaunay tetrahedrons corners

position of Delaunay tetrahedron corners are centers of particles

> Voronoi cells are not the same as the Delaunay tetrahedron (right?)

right, they are dual to each other [4]

A MWE using scipy.spatial.Delaunay [5]:
###
from scipy.spatial import Delaunay
sp = pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.1,rRelFuzz=.5,num=10,seed=1)
sp.toSimulation()
pts = [b.state.pos for b in O.bodies]
#
t = Delaunay(pts)
for i1,i2,i3,i4 in t.simplices:
for j1,j2 in ((i1,i2),(i1,i3),(i1,i4),(i2,i3),(i2,i4),(i3,i4)):
p1,p2 = pts[j1],pts[j2]
print("spacing",j1,j2,(p1-p2).norm())
# TODO: remove duplicites?
# TODO: consider radius of particles?
###

cheers
Jan

 Revision history for this message Robert Caulk (rcaulk) said on 2021-01-27: #9

>- Is there any error in this approach that I am not aware of?

Yes. You are assuming all sides of the tetrahedron are equal.

 Revision history for this message Othman Sh (othman-sh) said on 2021-01-27: #10

Thank you Jan and Robert.

Jan, your MWE is great. This solved my problem.

 Revision history for this message Othman Sh (othman-sh) said on 2021-01-27: #11

Thanks Jan Stránský, that solved my question.