How to calculate the surface roughness of sphere particle packing?

Asked by gaoxuesong on 2019-05-22

    I have got the sphere particle packing using Yade. I want to get the top surface roughness which is the height deviation to an mean value of the particle on the surface. The definition is shown in the picture as linked.

    A possible way i can think of is convert the particle packing into a stl format file and then obtain the surface (particle packing surface) geometrical information by some method. But i don't know how to conduct it.
   I hope someone can give me some advice.


Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Last query:
Last reply:
Jan Stránský (honzik) said : #1

Hi Xuesong,

nice task.
Since the packing is porous and in reality there is no "top surface", triangulation of particles would nor help..

A few options (depending on your needs and requirements):
- enlarge the particles "a bit". Then you could get "real surface" either analytically or by triangulation
- using some kind of a potential function of the spheres and assume the surface at certain level of this potential (see below)
- ...

### a MWE for the "potential" approach
# method params (you can play with the values)
level = .98 # closer to 1 -> closer to actual surface
factor = 2 # higher -> closer to actual surface
# dimensions of the packing
dim = 10.
# using only [dd,dim-dd] interval
dd = 2.5
# createa grid of the surface
n = 50

# "potential" of a body. =1 inside or on the surface, decreases to 0 for increasing distance from surface
def potential1(x,y,z,body,factor=1.):
   xyz = Vector3(x,y,z)
   p = body.state.pos
   r = body.shape.radius
   d = (xyz-p).norm() - r # distance from surface
   f = factor/r
   if d < 0: # inside
      ret = 1
   else: # outside
      ret = exp(-f*d) # could be other functions, like 1/x or so...
   return ret

# sum of potentials of all bodies. Very inoptimal (!) using all bodies, a few "boudary bodies" would be enough
def potential(x,y,z,factor=1.):
   return sum(potential1(x,y,z,b,factor) for b in O.bodies)

# for given x,y find z such that potential=level
# The closer is level to 1, the closer is the final z to actual surface
# using bisection method with initial boundaries z1,z2
def findZ(level,x,y,z1,z2,factor=1.):
   for i in range(20): # TODO a more sohisticated end condition
      z = .5*(z1+z2)
      l = potential(x,y,z,factor)
      if l < level:
         z2 = z
         z1 = z
   return z

# create some packing
pred = pack.inAlignedBox((0,0,0),(dim,dim,dim))
sphs = pack.randomDensePack(pred,1,rRelFuzz=.5,spheresInCell=100)

surf = [[None for y in range(n+1)] for x in range(n+1)]
for ix in range(n+1):
   x = dd + ix*(dim-2*dd)/n
   for iy in range(n+1):
      y = dd + iy*(dim-2*dd)/n
      z = findZ(level,x,y,.5*dim,2*dim,factor)
      surf[ix][iy] = (x,y,z)
# create facets to see the result
for ix in range(n):
   for iy in range(n):
      v1 = surf[ix+0][iy+0]
      v2 = surf[ix+1][iy+0]
      v3 = surf[ix+0][iy+1]
      v4 = surf[ix+1][iy+1]
      f1 = facet((v1,v2,v3))
      f2 = facet((v2,v3,v4))

# export to vtk for Paraview
from yade import export
vtk = export.VTKExporter("surf")


gaoxuesong (260582472-9) said : #2

Hi, Jan,

    Much thanks. Your method as well as code reallys works. I am very impressed by this method. Could you tell me more about it? Like the math behind it.

    I am a bit confused by the potential of one body. For the potential function outside a body, what does the factor/r really mean? Why the radius is set as denominator?


Jan Stránský (honzik) said : #3

I think there are many ways how to define the potential function, this is just a quick form, you can try to experiment.

The basic idea is that inside the body, the potential is 1 and with increasing distance, it tends to zero. So, even summed from all bodies, it is slightly less than 1 around the "boundary" (assuming approaching "from outside")
Decreasing exponential function is used for this purpose.
Factor controls the initial slope at the particle boundary (faster tends to zero for higher factor).
radius maybe should be in the numerator to make smaller effect of smaller particles.. Or maybe radius should not be present as all.. As said, it was just a quick try to make it work somehow.. Have a try and let us and future readers know what worked best :-)


gaoxuesong (260582472-9) said : #4

Hi. Sorry for late test.
Based on my experiments about this method, i find the parameter 'factor' is very import. Larger potential factor, more accurate the result. Here i put an image about how the 'factor' parameter affects the results. The 'factor' set as 8 is already enough for my case. Below is the link.

Also for the potential function, the parameter of 'r' can be removed. But the 'factor' should be increased to a larger value. In my case, the radius is in the order of microns and the accuracy becomes poor when i leave out the 'r'. Because the 'factor' is magnified by 1e6 equivalently. So if one wants to neglect the 'r', a larger 'factor' should be used. Taking the 'r' as denominator makes it easy to pick a 'factor' without consideration of the order of the particle size. The following result corresponds to the case where the 'r' is left out and the 'factor' is set as 1e7 instead.
Much thanks.