How to calculate the surface roughness of sphere particle packing?

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
      else:
         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)
O.bodies.append(sphs)

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))
      O.bodies.append((f1,f2))

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

cheers
Jan

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?

Thanks,
Xuesong

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 :-)

cheers
Jan

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.
https://drive.google.com/file/d/1Kl_uAiYYobjhhyx9_wKClCLnPxbDTtot/view?usp=sharing

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.
https://drive.google.com/file/d/10bkLM-_SYdskcw5FL64D1pbb0PVDtTVu/view?usp=sharing
Much thanks.