Discrete Element Method (DEM): Simulating Granular Materials

1. What Makes Granular Materials Hard

A pile of sand has ~10¹⁰ grains. Even a small laboratory sample of 10 000 grains cannot be described by classical continuum mechanics (which treats materials as differentiable fields) because:

Contacts are unilateral: grains push but never pull. Forces are compressive-only, creating asymmetric, history-dependent behaviour.
Friction is directional and rate-independent: the Coulomb friction model relates tangential to normal force, not to velocity.
No thermodynamic equilibrium: granular systems are athermal — random thermal energy at room temperature is negligible compared to grain weight. The system is always in a non-equilibrium, metastable state.
Jammed states: shaken grains flow; stopped grains jam rigidly. The jamming transition is not well described by standard phase-transition theory.

The Discrete Element Method (DEM), introduced by Cundall and Strack in 1979, treats every grain as a rigid body and computes all contact forces explicitly. The collective grain behaviour emerges from millions of individual two-body interactions — no constitutive law for the bulk material is assumed.

DEM vs SPH vs FEM for granular flows: DEM is accurate but expensive (particle count limited to ~10⁶–10⁷). Smooth Particle Hydrodynamics (SPH) can handle 10⁸ particles but uses fluid assumptions. Continuum FEM with a Drucker–Prager or μ(I) rheology is fastest but misses grain-scale effects like force chains and size segregation. Choice depends on problem scale and required detail.

2. Contact Detection

Before computing forces, we must determine which pairs of particles are in contact. For spheres (the predominant DEM grain shape), contact is simple:

// Sphere-sphere contact detection Particles i and j have positions xᵢ, xⱼ and radii rᵢ, rⱼ Distance: d = |xᵢ - xⱼ| Overlap: δ = rᵢ + rⱼ - d Contact exists iff δ > 0 (grains overlap) Contact normal unit vector: n = (xⱼ - xᵢ) / d (points from i toward j)

Non-spherical grains (ellipsoids, polyhedral rocks, rounded cubes) require more complex contact detection — GJK or SAT algorithms — and are 10–100× more expensive per contact pair, which is why DEM studies overwhelmingly use spheres.

Broad-Phase: Finding Candidate Pairs

Checking every pair of N particles costs O(N²). With N = 100 000 this is 10¹⁰ checks per timestep — impossible. A spatial hash grid or bounding volume hierarchy (BVH) reduces this to O(N) with a small constant:

// Uniform grid spatial hash (for uniform grain sizes) cell_size = 2 × max_radius For each particle i: cell = floor(xᵢ / cell_size) // 3D integer cell index hash(cell) → bucket // hash to a 1D array Neighbour search: check particle against all particles in the 3×3×3 = 27 surrounding cells. Average particles per cell ≈ packing fraction × 27 ≈ 15 Contact candidates reduced from O(N²) → O(15N) per frame

The grid must be rebuilt each timestep (particles move). For highly dynamic simulations, a linked-list-per-cell structure that supports O(1) insert/delete is preferred over rebuilding from scratch.

3. Contact Force Models

Linear Spring-Dashpot (LSD)

The simplest contact model represents each contact as a spring (elastic repulsion proportional to overlap) and a dashpot (viscous damping proportional to relative velocity):

// Normal contact force (compression only, no tension) Overlap: δₙ = rᵢ + rⱼ - |xᵢ - xⱼ| (positive when in contact) Relative velocity along normal: vₙ = (vᵢ - vⱼ) \cdot n Normal force: Fₙ = kₙ \cdot δₙ - γₙ \cdot vₙ (spring - dashpot) kₙ = normal spring stiffness [N/m] γₙ = normal damping coefficient [N\cdots/m] n = contact normal unit vector Apply to particle i: force += Fₙ \cdot n Apply to particle j: force -= Fₙ \cdot n (Newton's 3rd law)

The spring stiffness kₙ is a critical parameter. It must be large enough that overlaps remain physically small (δ ≪ r), but larger kₙ requires smaller timesteps (stability condition: Δt ≤ π × √(m/kₙ)). In practice kₙ is chosen to yield overlaps of ~0.5–2% of particle radius.

Hertz–Mindlin Contact Model

The more physically accurate model derives from Hertz contact theory, which predicts that the contact area grows with load and the force scales as δ^(3/2) rather than δ linearly:

// Hertz (non-linear) normal force model Effective radius: R* = (rᵢ \times rⱼ) / (rᵢ + rⱼ) Effective modulus: E* = ((1-νᵢ²)/Eᵢ + (1-νⱼ²)/Eⱼ)⁻¹ Stiffness (load-dependent): kₙ_hertz = (4/3) \times E* \times \sqrtR* \times \sqrtδₙ (varies with δₙ) Hertz force: Fₙ = kₙ_hertz \times δₙ = (4/3) \times E* \times \sqrtR* \times δₙ^(3/2) The non-linear force-overlap relation captures the real physics of elastic deformation at the grain contact surface.

Hertz–Mindlin requires material parameters (Young's modulus E, Poisson's ratio ν) that can be measured experimentally, making calibration to real materials straightforward. LSD requires choosing kₙ somewhat arbitrarily.

4. Time Integration

DEM uses Newton's second law for both translational and rotational motion of each particle:

// Equations of motion for particle i Translational: mᵢ \cdot ẍᵢ = Σ Fᵢ_contact + mᵢ \cdot g Rotational: Iᵢ \cdot ω̇ᵢ = Σ Tᵢ_contact mᵢ = mass (sphere: (4/3)π rᵢ³ ρ) Iᵢ = moment of inertia (sphere: (2/5) mᵢ rᵢ²) g = gravitational acceleration Tᵢ = torque from tangential contact forces Integration (explicit Euler or Velocity Verlet): vᵢ(t+Δt) = vᵢ(t) + (Fᵢ/mᵢ) \cdot Δt xᵢ(t+Δt) = xᵢ(t) + vᵢ(t+Δt) \cdot Δt (symplectic Euler) ωᵢ(t+Δt) = ωᵢ(t) + (Tᵢ/Iᵢ) \cdot Δt

DEM timesteps are tiny: typically Δt ≈ 10⁻⁶–10⁻⁵ seconds (far smaller than the simulation duration of seconds to minutes). The critical timestep is set by the natural oscillation period of the spring-dashpot system. Violating it causes numerical instability — overlaps grow without bound.

Critical timestep: Δt_crit = π × √(m_min / kₙ_max). For sand grains of diameter 1 mm (m ≈ 10⁻⁶ kg) with kₙ = 10⁵ N/m, Δt_crit ≈ 10⁻⁵ s. Simulating 1 second of physical time needs ~100 000 timesteps. At 10 000 particles with 5 contacts each, that is 5×10⁹ force calculations per second of simulation time.

5. Friction, Rolling, and Cohesion

Coulomb Tangential Friction

A grain resists sliding tangentially. The classic Coulomb friction model limits the tangential force to a fraction of the normal force:

// Tangential (shear) force with Coulomb friction bound Track tangential displacement δₜ: integrate relative tangential velocity δₜ(t+Δt) = δₜ(t) + vₜ_rel \cdot Δt Tangential spring force: Fₜ_spring = kₜ \cdot δₜ (kₜ \approx (2/7) kₙ for Hertz-Mindlin) Coulomb limit: if |Fₜ_spring| > μ \times Fₙ: Fₜ = μ \times Fₙ \times sign(Fₜ_spring) // sliding friction reset δₜ = Fₜ / kₜ // trim spring to avoid accumulation else: Fₜ = Fₜ_spring // static friction μ = coefficient of sliding friction (typically 0.3-0.7)

Rolling Resistance

Perfect spheres roll without rolling resistance — they would form unrealistically flat piles. A rolling friction torque resists rotation at the contact, parametrised by a rolling friction coefficient μᵣ ≪ μ:

// Rolling resistance torque (contact-limited models) Torque magnitude: |Tᵣ| \leq μᵣ \times rᵢ \times Fₙ Direction: opposes ω_relative at contact (relative angular velocity projected onto contact plane)

Cohesion (van der Waals / Capillary)

Fine powders (d < 100 μm) experience significant inter-particle attraction from van der Waals forces. Wet sand gains cohesion from liquid bridges (capillary adhesion). These are modelled by adding an attractive force component, using a Johnson–Kendall–Roberts (JKR) or Derjaguin–Müller–Toporov (DMT) contact mechanics extension, or a simplified "cohesion energy density" approach where grains slightly outside contact still attract.

6. Emergent Phenomena: Force Chains & Angle of Repose

Force Chains

Strike a sand pile and the impact force does not spread uniformly — it travels along force chains: narrow, branching paths of heavily loaded grains surrounded by nearly unloaded "spectator" grains. The heterogeneity is striking: the most-loaded grains carry 10× the average load.

DEM simulations predicted force chains before they were imaged experimentally. Photoelastic disk experiments (using birefringent discs that show stress optically under polarised light) confirmed DEM predictions in the 1990s. Force chains are critical to understanding stress transmission in granular media and the sudden onset of flow (avalanche).

Angle of Repose

Pour sand into a pile: it builds to a maximum slope angle — the angle of repose θᵣ — and stops. Add more, and small avalanches maintain that slope. DEM naturally produces this:

Coulomb friction coefficient μ → tan(θᵣ) ≈ μ (for uniform spheres on a flat surface)
Rolling friction increases θᵣ (angular particles have higher effective rolling friction)
Grain size distributions affect θᵣ (wider distributions → higher θᵣ due to interlocking)

Empirical relationship: θᵣ \approx arctan(μ_eff) μ_eff = effective friction accounting for grain shape, size distribution, and rolling resistance For smooth glass spheres: θᵣ \approx 22°-25° (μ \approx 0.4-0.45) For dry sand: θᵣ \approx 30°-35° (includes shape interlocking) For wet sand: θᵣ \approx 45°+ (capillary cohesion) For angular gravel: θᵣ \approx 40°-45°

Size Segregation (Brazil-Nut Effect)

Shake a container of mixed grain sizes: larger grains rise to the top — the Brazil-nut effect. DEM simulations reproduce this through a convection/void-filling mechanism: small grains fall into gaps below large grains during vibration, ratcheting large grains upward. The effect reverses (inverse Brazil-nut) for dense large grains, and is sensitive to grain size ratio, density ratio, and vibration parameters.

7. Performance: Spatial Hashing & GPU DEM

CPU DEM Performance

A well-optimised CPU DEM code handles ~10⁵–10⁶ particles at interactive speeds (a few timesteps per second for slow dynamics). Key optimisations:

Verlet neighbour lists: Rebuild contact candidates only every few steps. Check every timestep if particles have moved more than a "skin depth" since last rebuild.
Linked-list-per-cell: O(1) insert/remove from spatial grid as particles move.
SIMD: Vectorise force calculations over batches of contact pairs using AVX2/AVX-512.
OpenMP: Parallelize force computation over particles (beware race conditions on shared force arrays — use thread-local accumulators or atomic operations).

GPU DEM

GPU DEM pushes to 10⁷–10⁸ particles by executing contact force kernels with one GPU thread per contact pair. Challenges unique to GPU DEM:

Atomic force accumulation: Multiple threads simultaneously add forces to the same particle. Use CUDA atomicAdd or split into phases where each thread "owns" a particle.
GPU spatial hashing: Sort particles by cell key using GPU radix sort (CUB library). Then scan to find cell boundaries. Rebuild every timestep but each step is only ~1 ms on a modern GPU.
Memory layout: AoS (Array of Structs) is cache-inefficient for GPU. SoA (Struct of Arrays) packs x-coordinates together, y-coordinates together — far better coalescing.

// CUDA kernel structure for contact force (pseudocode) __global__ void computeForces(float* pos, float* vel, float* force, int* cellStart, int* cellSize, ...) { int i = blockIdx.x * blockDim.x + threadIdx.x; // particle index float3 fi = {0,0,0}; // iterate neighbor cells (27 cells) for each cell in neighbors(cell_of(i)): for j in particles_in(cell): if (i == j) continue; float3 delta = pos[j] - pos[i]; float dist = length(delta); float overlap = radii[i] + radii[j] - dist; if (overlap > 0): fi += contact_force(delta, dist, overlap, vel[i], vel[j]); atomicAdd(&force[i*3], fi.x); // write result atomicAdd(&force[i*3+1], fi.y); atomicAdd(&force[i*3+2], fi.z); }

8. Industrial Applications

DEM is used wherever granular materials play a key role:

Mining & mineral processing: Ball mill grinding (ore + steel balls), conveyor design, silo flow and discharge rates, funnel flow vs mass flow regimes. DEM predicts clogging — where an arch of grains blocks a hopper — critical for preventing costly production stoppages.
Pharmaceutical: Tablet coating in rotating drums, powder blending uniformity, capsule filling. Granule breakage models predict dissolution rate differences. FDA increasingly requires DEM-based process validation.
Geotechnical engineering: Landslide run-out prediction, debris flow, earthquake-induced liquefaction, retaining wall design, underground cave-in mechanics.
Food processing: Cereal, rice, coffee bean handling — predicting segregation by size and density, conveyor wear, bulk storage optimisation.
Additive manufacturing: Metal powder bed fusion (SLM, EBM) — DEM predicts powder bed density and uniformity, which determines laser energy coupling and final part porosity.
Planetary science: Regolith mechanics on asteroids (OSIRIS-REx mission used DEM models to plan the touchdown and sample collection on Bennu), crater formation, ring dynamics in Saturn's rings.

Commercial codes: EDEM (Altair), Rocky DEM (ESSS), PFC (Itasca), LIGGGHTS (open-source LAMMPS-based). Open-source MercuryDPM and YADE are widely used in research. All implement the Hertz–Mindlin model and support GPU acceleration. A minimal DEM from scratch in Python (NumPy + Numba) can run ~10 000 particles in real time for educational purposes.