Voronoi volumes and local density - Introduction

In a typical static random packing of the spheres, the particles can occupy approximately 60% to 65% of the free volume. During flow, this packing fraction can be decreased by several percent, since the particles must have room to rearrange. Tracking these changes can provide important insights into the dynamics of the random packing.

Ideally we would like to look at these density fluctuations on a scale of several particle diameters. However, due to the discreteness of the problem, getting an accurate local measurement of density is difficult, and naive methods of counting up the number of particles in a small region can have errors larger than the scale of interest.

Voronoi cells

Voronoi cells for a section of packing

We have therefore made use of Voronoi cells. For each particle in a granular packing, we can define a cell around it which consists of all the space which is closer to that particle than any other. Voronoi cells take the form of irregular polyhedra, the sides of which are the planes which are perpendicular bisectors between a particle and its neighbors. Dividing a cell volume by the sphere volume gives an accurate measurement of local packing fraction at a scale of a single particle. The following diagram shows Voronoi cells computed for a section of random packing between horizontal planar walls eight particle diameters apart.

Voronoi cells
	at an oblique corner
Voronoi cells at an oblique corner

A voronoi
	cell approximation to a sphere
An approximation to a sphere made by intersecting many random planes

Our algorithm

Several methods for computing Voronoi tesselations are available, and one is built into the popular mathematical package Matlab, which makes use of the dual Delaunay triangulation. For my group's purposes, I created a plane-cutting algorithm in C++ which would rapidly calculate Voronoi cells directly. The algorithm starts with a large simple polyhedron encompassing the particle. For each neighboring particle we need to make a planar cut of this polyhedron. We make use of the following algorithm:

  1. Pick an arbitrary vertex.
  2. Trace from vertex to vertex until an edge intersecting the cutting plane is found. If none exists, then the plane does not intersect the polyhedron.
  3. Trace from this vertex to find all intersected edges.
  4. Recompute a face.
  5. Delete all remaining points that were cut.

Such an algorithm lends itself well to cases where we wish to handle complicated geometries, such as this oblique corner boundary shown above. The algorithm has also been use to look at a conical funnel. The algorithm runs very fast, handling tens of thousands of Voronoi cells per second. It can also compute very complicated cells, such as this approximation to a sphere, created by cutting a polyhedron by many random planes all of unit distance away from the origin.

Density comparison for Spot Model and DEM

The following two sequences of images track small changes in local density for the spot model and DEM simulations using Voronoi cells. Before flow starts, the packing fraction is approximately 63% (colored cyan), and during flow, this drops to around 60% (colored dark blue) or 57% (colored yellow). Density decreases to 50% or lower are colored red and can be seen in the vicinity of the orifice.

DEM Simulation
Montage of density plots for the DEM Simulation
Spot Model
Montage of density plots for the spot model

While the spot model provides a very good match to DEM in many aspects, this sequence of plots highlights a discrepancy between the two simulations. We see that the two simulations agree on how fast the density decrease propagates through the packing. However, the current method of spot motion, based on simple random walks, predicts that the largest density decrease should be in the center of the flow, where the velocity is highest, but it appears that the largest density decrease is on either side of the center, in the regions of highest shear. This suggests the possibility of finding better, more physically realistic spot rules that take into account forces and shear at a mesoscopic scale. To investigate this possibility we have recently been carrying out a systematic study of these quantities in a variety of granular flows, with very interesting results.

See also a comparison movie between the spot model (left) and DEM (right).