Disclaimer: Your browser will need to render Javascript to properly display the equations below.
Granular material is a very common form of matter, both in nature and industry. It can be defined as material consisting of interacting, discrete particles. Common granular materials include gravels, sands and soils, ice bergs, asteroids, powders, seeds, and other foods. Over 75% of the raw materials that pass through industry are granular. This wide occurrence has driven the desire to understand the fundamental mechanics of the material.
Contrary to other common materials such as gases, liquids and solids, a general mathematical formulation of it’s behavior hasn’t yet been found. Granular material can, however, display states that somewhat resemble gases, fluids and solids.
The left image above shows sand at rest, behaving like a solid. The sand piles are motionless and static, even though they have surface slopes and the gravitational pull attempts to turn the pile into a flat layer. When the surface slopes exceed the angle of repose, the grains form an avalanche and flow downwards, much like a liquid. This is seen in the center image, where the accumulated grains at the top of a dune tumble down the lee side. When the individual grains are moving at high velocities, like inside the dust storm in the above right image, the grains are interacting by shortlived inelastic collisions, much different than the interactions in conventional gases. For a well written review of the curious states and the transformation of granular material, I can recommend the 1996 article by Heinrich M. Jaeger and others.
The lack of a general mathematical description of granular mechanics leaves laboratory analyses and numerical experiments as the main tools to predict the behavior. In the numerical approach, the key to understanding the overall material mechanics lies in understanding the contact mechanics between a pair of particles. In this post, I’ll explain the theory and show an example of a method commonly used to simulate granular material numerically.
Theory
The Discrete Element Method (DEM) is a numerical method that can be used to simulate the interaction of grains. Originally derived from Molecular Dynamics, it simulates particles as separate entities, and calculates their positions, velocities, and accelerations through time.
In the DEM, the time is discretized into small steps ($\Delta t$). For each time step, the entire network of contacts is resolved, and the resulting forces and torques for each particle are found. With these values at hand, the new linear and rotational accelerations can be found using Newton’s second law of the motion of solid bodies. If a particle with mass $m$ at a point in time experiences a sum of forces denoted $\vec{F}$, the resultant acceleration ($\vec{a}$) can be found by rearranging Newton’s second law: $$ \begin{align} \vec{F} = m \vec{a} \Rightarrow \vec{a} = \frac{\vec{F}}{m} \end{align} $$ The new velocity and position is found by integrating the above equation with regards to time. The simplest integration scheme in this regard is the Euler method:
where $\vec{v}$ and $\vec{p}$ are the velocity and position vectors of the particle, respectively. And that’s it; the kinematic properties of the particle have been predicted at an interval $\Delta t$ after the current time, assuming that the force, acceleration, and velocity was constant during the interval.
Now the question arises on how to determine the sum of forces $\vec{F}$ acting on a given particle. To perform this force summation, it is necessary to determine the external forces, such as gravity, and the interparticle contact forces from adjacent particles. The latter requires mapping of the particle contacts, and application of a contact model.
DEM's typically use a softbody or rigidbody formulation. In the softbody approach, particles are allowed to slightly overlap. The repulsive forces and torques are computed from this overlap. During the temporal integration, new accelerations, velocities and positions are found as described above.
In rigidbody (or nonsmooth) formulations particles are by definition not allowed to have overlapping volumes. Velocities are discontinuous through time (hence the nonsmooth name), and the solution procedure is an iterative process.
The softbody approach requires that the time step is sufficiently small to resolve the elastic (seismic) wave propagation through the particles. The rigidbody formulation allows for larger time steps, but accelerations are not found. I will describe the simpler softbody approach in the following.
Contact search
If the particles are treated as homogeneous discs (in 2D) or spheres (in 3D), the particle geometry is fully described by the center position ($\vec{p}$) and radius ($r$). The state of two particles, denoted with superscripts $i$ and $j$, can be determined by calculating the overlap:
The term $\vec{p}^i  \vec{p}^j$ is the distance between the particle centers. If the value of $\delta_n^{ij}$ is positive, the two particles are not in contact. If the value is $0$, the particles are exactly touching. If the value is negative, the particle areas (2D) or volumes (3D) are overlapping.
Contact interaction
Now that the interparticle contacts have been identified and characterized by their overlap, the resulting forces from the interaction can be resolved. One of the simplest models is linearelastic (Hookean) repulsion between particles $i$ and $j$:
The fraction is there to specify the direction of particle contact. The parameter $k_n$ is the spring coefficient. The force above is called the normal force, since it characterizes the repulsion acting on the particles, normal (perpendicular) to the contact plane. The contact model can, and usually is, extended by a more complex, nonlinear interaction law, as well as a tangential contact force acting perpendicular to the normal force.
Numerical codes for simulating granular interaction
When using DEM simulation programs, most of the mathematical and numerical theory described above is hidden from the user, whose main task involves simulation setup and analysis. Quite a few free and open source scientific DEM codes exist. Popular softbody type choices include ESySParticle, LIGGGHTS/LAMMPS and YADE. Rigidbody types include LMGC90, the Chrono Engine and the Bullet Physics Engine. I’m going to give an example using my own GPUaccelerated softbody particle dynamics code, sphere.
Below is a sample script for a simple DEM simulation with sphere. It creates a simulation object, sets the simulation parameters, runs the simulation for a prescribed length of time, and finally renders the particles as images, and shows how energy evolved through time.

The rendered images of the particles are placed in the ../img_out
folder. The
particles are colored corresponding to the sum of pressure they experience at
that point in time. Blue colors denote low or no pressure, and red colors denote
high pressure. The color scale boundary values can be set in the SB.render()
call.
Using the above template, as many particles as desired can be added. In the sphere Python API it is often desirable to use the builtin functions that generate packings of many particles in a single call. A twobody problem doesn’t even get close to saturating the many cores of modern GPU’s. The sphere program is meant for simulating thousands of particles simultaneously. See the python/ folder for more examples.