# 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.

Thanks,
Xuesong

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
gaoxuesong
Solved:
2019-06-04
Last query:
2019-06-04
2019-05-23
 Jan Stránský (honzik) said on 2019-05-23: #1

Hi Xuesong,

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

 gaoxuesong (260582472-9) said on 2019-05-23: #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?

Thanks,
Xuesong

 Jan Stránský (honzik) said on 2019-05-23: #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 :-)

cheers
Jan

 gaoxuesong (260582472-9) said on 2019-06-04: #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.