A Variational Formulation for Deformable Particle Simulations and its Level Set Discrete Element Method Implementation

This paper presents a computationally efficient deformable Discrete Element Method (DEM) that extends classical rigid-particle dynamics through an energetic variational formulation and level set evolution, enabling physically accurate grain-scale deformation for arbitrary particle geometries at a cost comparable to standard rigid DEM.

The Gist

Imagine you are trying to simulate a pile of sand, a stack of rocks, or even a tumor made of cells on a computer. For decades, scientists have used a tool called the Discrete Element Method (DEM) to do this.

Think of traditional DEM like playing with marbles. You can roll them, bounce them, and stack them. The computer knows exactly where each marble is and how they hit each other. But there's a catch: in this simulation, marbles are perfectly rigid. If you push two marbles together, they don't squish; they just stop. If you want to see a marble flatten or bend, the old software can't do it without breaking the simulation or taking a million years to calculate.

This new paper introduces a clever upgrade: Deformable DEM. It's like turning those rigid marbles into stress balls or gummy bears. Now, when they collide, they squish, flatten, and bounce back, just like real matter.

Here is how they did it, broken down into simple concepts:

1. The "Magic Shape-Shifter" (The Variational Formulation)

Usually, to make a computer understand how something bends, you have to break it into millions of tiny pieces (like a 3D puzzle) and calculate the physics for every single piece. This is incredibly slow.

The authors used a "variational formulation." Think of this as a master recipe for energy. Instead of calculating every tiny atom, they asked: "What is the simplest way to describe how this object changes shape?"

They realized that most objects don't twist into weird, random shapes. They usually bend, stretch, or squash in predictable ways. So, instead of tracking millions of points, they track just a few "modes" (like musical notes).

• Analogy: Imagine a guitar string. You don't need to track every atom in the string to know how it vibrates. You just need to know the "fundamental note" (the main vibration) and maybe a few "harmonics."
• The Innovation: They added these "vibration notes" (deformation modes) to the standard rules of motion. Now, a particle has three types of movement:
  1. Translation: Moving left/right.
  2. Rotation: Spinning.
  3. Deformation: Squishing or bending (the new "note").

2. The "Ghost Outline" (Level Set Method)

When a marble squishes, its shape changes. How does the computer know where the edge is now?

They use something called a Level Set.

• Analogy: Imagine a particle is a balloon floating in a foggy room. The "Level Set" is the invisible boundary where the fog turns into clear air.
• The Trick: When the balloon (particle) squishes, the computer doesn't redraw the whole balloon. Instead, it just updates the "foggy outline" based on the "notes" (modes) we mentioned earlier. It's like using a magic marker to trace the new shape instantly, rather than rebuilding the balloon from scratch. This keeps the simulation fast, even when the particles are changing shape.

3. Why This Matters: The "Goldilocks" Solution

Before this, scientists had two bad options:

• Option A (Rigid DEM): Fast, but unrealistic. Like playing with marbles that never dent.
• Option B (Finite Element Method - FEM): Super realistic, but painfully slow. It's like trying to simulate a pile of sand by modeling every single grain as a complex 3D sculpture. You could only simulate a few grains before your computer crashed.

This new method is the "Goldilocks" solution:

• It's fast (almost as fast as the rigid marble simulation).
• It's realistic (particles actually squish and change shape).
• It's scalable (you can simulate thousands of particles, not just a handful).

Real-World Examples from the Paper

The authors tested this on two scenarios:

1. Bending a Stick: They simulated a long, thin grain bending like a ruler. The computer predicted exactly how much force was needed to bend it, matching real-world physics perfectly.
2. Squishing a Ball: They simulated a sphere being crushed from all sides (like a grape in a press). The particle flattened out, creating a flat spot where it touched the walls. This is something rigid marbles can never do, but it's crucial for understanding how soil compacts or how biological tissues (like tumors) behave.

The Bottom Line

This paper gives us a new way to simulate the world. It allows computers to see materials not just as hard, unyielding blocks, but as living, breathing, squishy things that change shape when pushed.

In short: They taught the computer to play with gummy bears instead of marbles, using a clever shortcut that makes the simulation run fast enough to model entire piles of sand, rocks, or even biological tissues without needing a supercomputer.

Technical Summary

1. Problem Statement

Traditional Discrete Element Method (DEM) simulations model particles as rigid bodies, which limits their ability to capture bulk deformation (global shape changes) of individual grains. While classical DEM approximates local contact deformation (e.g., Hertzian overlap), it cannot simulate the elastic or inelastic shape changes of the particles themselves.

Existing approaches to model deformable particles face significant trade-offs:

• FEM-DEM Coupling: Solves internal boundary value problems for each particle using Finite Element Method (FEM). While accurate, it is computationally prohibitive for large systems due to the need for individual meshes and remeshing.
• Cluster/Spring Models: Represent particles as aggregates of smaller rigid units connected by springs. These require many sub-particles to achieve realistic deformation, leading to high computational costs.
• Reduced-Order Models: Use simplified kinematics or compliance matrices but often lack a rigorous connection to fundamental mechanics, are restricted to specific geometries, or fail to handle nonlinearities transparently.

The challenge is to develop a framework that captures grain-scale bulk deformation while retaining the computational efficiency, scalability, and geometric clarity of classical rigid-particle DEM.

2. Methodology

The authors propose a Variational Deformable DEM framework implemented within the Level Set DEM (LS-DEM) formalism. The methodology consists of three core pillars:

A. Variational Formulation (Lagrange–d'Alembert Principle)

The authors extend classical rigid-body dynamics by embedding translational, rotational, and deformation degrees of freedom into a unified energetic description.

• State Variables: Each grain is described by position ($X$), rotation ($\Theta$), and a set of generalized modal coordinates ($\varepsilon_\alpha$) representing internal elastic deformations.
• Reduced Kinematics: Particle deformation is approximated as a linear combination of fixed spatial shape functions (modes) $\Phi_\alpha(x)$ with time-dependent amplitudes $\varepsilon_\alpha(t)$:
  $$u_\varepsilon(x, t) = \varepsilon_\alpha(t) \Phi_\alpha(x)$$
• Energetics: The framework derives kinetic and potential energy terms for these deformation modes.
  • Kinetic Energy: Defined via a generalized mass matrix $M_{\alpha\beta}$.
  • Potential Energy: Defined via a generalized stiffness matrix $K_{\alpha\beta}$ (for linear elasticity) or an arbitrary energy functional $U_\varepsilon(\varepsilon_\alpha)$ (to accommodate nonlinearities).
• Governing Equations: Applying the Lagrange–d'Alembert principle yields a system of Ordinary Differential Equations (ODEs) that couple rigid-body motion with modal evolution:
  $$M_\varepsilon \ddot{\varepsilon} + K_\varepsilon \varepsilon = F_\varepsilon$$
  Where $F_\varepsilon$ represents the generalized forces derived from contact interactions.

B. Level Set Implementation (LS-DEM)

The framework utilizes the LS-DEM formalism, where particle geometry is represented by a level set function $\phi(x)$ on a grid and surface nodes.

• Geometry Update: As particles deform, both the surface nodes and the level set field must evolve.
• Level Set Evolution Strategies: The paper evaluates several methods for updating the level set under deformation:
  1. Recomputation: Accurate but expensive.
  2. Advection: Solves a transport equation; requires velocity extrapolation.
  3. Semi-Lagrangian Advection (Selected): Maps current grid points back to the reference configuration using the inverse displacement field ($x_{ref} = x - u(x,t)$) and interpolates the initial level set. This is computationally efficient and avoids solving PDEs at every step.
  4. Precomputed Interpolation: Uses precomputed level sets for specific deformation states (high memory cost).

C. Mode Extraction

The framework supports both analytical modes (closed-form expressions, e.g., beam bending) and numerical modes extracted from high-fidelity FEM simulations or experiments.

• For complex systems, modes are extracted via Proper Orthogonal Decomposition (POD) of displacement fields from unit cell or Representative Volume Element (RVE) simulations.
• Nonlinear contact effects are handled by constructing a nonlinear energy functional or force-displacement lookup table based on the extracted modes.

3. Key Contributions

1. Unified Variational Framework: A rigorous derivation extending rigid-body dynamics to include deformability using the Lagrange–d'Alembert principle, ensuring consistency with thermodynamics and mechanics.
2. Generalizability: Unlike previous reduced-order models, this formulation is not restricted to specific particle shapes (e.g., spheres) and supports arbitrary geometries and topologies via the Level Set method.
3. Computational Efficiency: By using a reduced-order modal description, the method adds only a small number of scalar ODEs per particle. It maintains the same order of computational complexity as rigid DEM, unlike FEM-DEM which scales poorly with particle count.
4. Nonlinearity Handling: The framework naturally accommodates nonlinear material behavior and, crucially, nonlinear contact-driven effects (e.g., evolving contact areas) through arbitrary energy functionals or data-driven force-displacement relations.
5. Robust Implementation: The integration with LS-DEM allows for the simulation of complex particle interactions and evolving geometries without the need for mesh regeneration.

4. Results and Validation

The authors validated the framework through several benchmarks:

• Analytical Bending Mode: A three-point bending test on a slender particle using an analytical Euler-Bernoulli mode shape. The deformable LS-DEM results showed excellent agreement with the analytical generalized stiffness and force-displacement response.
• Numerical Compaction Mode: A single spherical grain compressed by rigid walls. A deformation mode was extracted from a high-fidelity FEM simulation. The deformable LS-DEM successfully replicated the contact flattening and the nonlinear Hertzian-like force-displacement curve, capturing the transition from point contact to flat contact.
• System-Scale Simulation: A 27-particle cubic packing under hydrostatic compression. The system-level force-displacement response matched a semi-analytical prediction derived from the single-particle response, demonstrating the framework's ability to scale to multi-particle assemblies while capturing grain-scale densification.

5. Significance

This work bridges the gap between high-fidelity continuum mechanics (FEM) and efficient discrete methods (DEM).

• Scientific Impact: It provides a general, extensible theoretical foundation for modeling particulate systems where grain deformation is critical, such as geomaterials (soils/rocks), biological tissues (nacre, tumors), and architected metamaterials.
• Engineering Utility: It enables the simulation of large-scale systems (thousands to millions of particles) with grain-scale physics that were previously inaccessible without prohibitive computational costs.
• Future Potential: The variational structure allows for future extensions to large deformations, rotational degeneracy handling, and coupled multiphysics (e.g., poroelasticity, inelasticity), making it a versatile tool for the next generation of granular matter simulations.