Interparticle spacing

Asked by Othman Sh on 2021-01-24

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

[1] https://yade-dem.org/doc/Yade.pdf
-----------

from yade import pack

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

########## 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)
yade.qt.View()
########## spheres ##########
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=sp_radius)
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:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2021-01-27
Last query:
2021-01-27
Last reply:
2021-01-27
Robert Caulk (rcaulk) said : #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.

[1]http://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.getVertices

Othman Sh (othman-sh) said : #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

Robert Caulk (rcaulk) said : #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

Othman Sh (othman-sh) said : #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

Welcome to Yade 20201215-4510~714a723~focal1
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
--------------------------------

from yade import pack
import numpy as np

sp_radius=2.36/1000

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)

yade.qt.View()
########## spheres ##########
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=sp_radius)
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)

Jérôme Duriez (jduriez) said : #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()

Jan Stránský (honzik) said : #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?
Or Yade's TesselationWrapper [1]?

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TesselationWrapper

Othman Sh (othman-sh) said : #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

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.printVertices
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngineT.getCellVolume
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TesselationWrapper.computeVolumes

Best Jan Stránský (honzik) said : #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
from yade import pack
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

[4] https://en.wikipedia.org/wiki/Delaunay_triangulation#Relationship_with_the_Voronoi_diagram
[5] https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html

Robert Caulk (rcaulk) said : #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.

Othman Sh (othman-sh) said : #10

Thank you Jan and Robert.

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

Othman Sh (othman-sh) said : #11

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