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.

Sand in different mechanical states

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 short-lived 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.


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 inter-particle contact forces from adjacent particles. The latter requires mapping of the particle contacts, and application of a contact model.

DEM's typically use a soft-body or rigid-body formulation. In the soft-body 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 rigid-body (or non-smooth) formulations particles are by definition not allowed to have overlapping volumes. Velocities are discontinuous through time (hence the non-smooth name), and the solution procedure is an iterative process.

The soft-body approach requires that the time step is sufficiently small to resolve the elastic (seismic) wave propagation through the particles. The rigid-body formulation allows for larger time steps, but accelerations are not found. I will describe the simpler soft-body approach in the following.

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 inter-particle contacts have been identified and characterized by their overlap, the resulting forces from the interaction can be resolved. One of the simplest models is linear-elastic (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, non-linear 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 soft-body type choices include ESyS-Particle, LIGGGHTS/LAMMPS and YADE. Rigid-body types include LMGC90, the Chrono Engine and the Bullet Physics Engine. I’m going to give an example using my own GPU-accelerated soft-body 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.

( download
#!/usr/bin/env python
Example of two particles colliding.
Place script in sphere/python folder,
and invoke with `python`

# Import the sphere module for setting up, running, and analyzing the
# experiment. We also need the numpy module when setting arrays in the sphere
# object.
import sphere
import numpy


# Create a sphere object with two preallocated particles and a simulation ID
SB =
sphere.sim(np = 2, sid = 'collision')

SB.radius[:] = 0.3 # set radii to 0.3 m

# Define the positions of the two particles
SB.x[0, :] = numpy.array([0.0, 0.0, 0.0]) # particle 1 (idx 0)
SB.x[1, :] = numpy.array([1.0, 0.0, 0.0]) # particle 2 (idx 1)

# The default velocity is [0,0,0]. Slam particle 1 into particle 2 by defining
# a positive x velocity for particle 1.
SB.vel[0, 0] = 1.0

# Set the world limits and the particle sorting grid. The particles need to stay
# within the world limits for the entire simulation, otherwise it will stop!
SB.initGridAndWorldsize(margin = 10.0)

# Define the temporal parameters, e.g. the total time (total) and the file
# output interval (file_dt), both in seconds
SB.initTemporal(total = 2.0, file_dt = 0.1)

# Using a 'dry' run, the sphere main program will display important parameters.
# sphere will end after displaying these values. = True)


# Start the simulation on the GPU from the sphere program


# Plot the system energy through time, image saved as collision-energy.png
sphere.visualize(SB.sid, method = 'energy')

# Render the particles using the built-in raytracer

# Write the output to VTK files for ParaView

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 built-in functions that generate packings of many particles in a single call. A two-body 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.