DEM Contact Model Fundamentals

• DEM Contact Model is a computational approach that decomposes contact forces into normal and tangential components, enabling accurate simulation of particle interactions.
• It applies bilinear elasto-plastic adhesive–frictional laws with incremental spring–dashpot models to capture hysteresis, cohesion, and frictional slip.
• Proper scaling of model parameters preserves macroscopic responses, enhancing computational efficiency while maintaining predictive fidelity across particle sizes.

The Discrete Element Method (DEM) contact model is the foundational component of DEM simulations, governing the computation of forces and moments during particle–particle and particle–wall interactions. The choice and calibration of the contact law govern the predictive fidelity, computational efficiency, and mechanical realism of DEM simulations for a wide range of granular, cohesive, and complex particulate systems.

1. Mathematical Structure of DEM Contact Models

DEM contact models generally decompose the total contact force into normal and tangential (frictional) components, potentially augmented by damping (viscoelastic), adhesive (cohesive), rolling, and torsional contributions. The normal force is typically history-independent for purely elastic spheres but becomes path-dependent or hysteretic in adhesive, plastic, or frictional models. Tangential forces, even in the simplest stick–slip formulations, are inherently path- and history-dependent.

The most common formulation for frictional, elasto-plastic, and cohesive contacts is as follows:

$\mathbf{f}_n = (f_{\rm hys} + f_{nd})\mathbf{u}$

where

• $f_{\rm hys}$: hysteretic (elasto-plastic) spring force,
• $f_{nd} = -B_n v_n$: viscous normal dashpot,
• $\mathbf{u}$: outward contact normal.

For the tangential direction: $\mathbf{f}_t = (f_{ts} + f_{td})$ where

• $f_{ts}$: incremental tangential spring ($$f_{ts}(n) = f_{ts}(n-1) + k_t \Delta\boldsymbol{\xi}_t$$),
• $f_{td} = -B_t v_t$: tangential dashpot,
• Coulomb bound: $$|\mathbf{f}_t| \leq \mu(f_{\rm hys} + k_{adh}\delta_p) - f_0$$.

This structure enables the representation of a variety of contact behaviors such as load–unload hysteresis, elasto-plastic transitions, stick–slip friction, and cohesive adhesion.

2. Bilinear Elasto-Plastic Adhesive–Frictional Contact Laws

Thakur et al. (2015) developed a scale-bridging DEM contact law for mesoscopic modeling of cohesive powders undergoing quasi-static loading, featuring a bilinear, path-dependent force–displacement response:

• Loading branch:

  $$f_{hys} = f_0 + k_1 \delta, \quad \text{if } k_2(\delta - \delta_p) \geq k_1 \delta$$

• Unloading/reloading branch:

  $$f_{hys} = f_0 + k_2(\delta - \delta_p), \quad \text{if } -k_{adh}\delta_p < k_2(\delta - \delta_p) < k_1\delta$$

• Adhesive (pull-off) branch:

  $$f_{hys} = -k_{adh}\delta_p, \quad \text{if } k_2(\delta - \delta_p) \leq -k_{adh}\delta_p$$

Here,

• $\delta$: normal overlap,
• $\delta_p = (k_1/k_2)\delta_{max}$: plastic offset,
• $f_0$: constant pull-off (van der Waals or JKR-type) force.

This model enables hysteretic, elastic-plastic response with cohesive pull-off. The stiffness parameters satisfy $k_1 \propto d$, $k_2 \propto d$, $k_{adh} \propto d$ (with particle diameter $d$), and $f_0 \propto d^2$.

The tangential force uses an incremental spring–dashpot model truncated by a modified Coulomb–type bound: $$|\mathbf{f}_t| \leq \mu(f_{\rm hys} + k_{adh} \delta_p) - f_0$$ with $k_t = \frac{2}{7} k_1 \propto d$. This construction enables model-scale invariance and accurate reproduction of quasi-static bulk loading response under both confined and unconfined compression (Thakur et al., 2015).

3. Scaling Laws for DEM Parameter Coarse-Graining

Parameter scaling is required to accelerate large-scale DEM simulations while preserving macroscopic response, particularly when replacing fine-grained powders with larger "coarse" particles. Thakur et al. (2015) established that to preserve bulk stress–strain and porosity evolution in cohesive and cohesionless solids, DEM contact model parameters must scale as:

• $k_n \propto d$: normal spring stiffness,
• $k_t \propto d$: tangential spring stiffness,
• $f_0 \propto d^2$: cohesive force (JKR/van der Waals/pull-off),
• $k_{adh} \propto d$: adhesive-branch stiffness.

This scaling results in identical macroscopic mechanical behavior in uniaxial and confined compression across particle sizes, enabling up to an order-of-magnitude reduction in particle number and computational cost while retaining predictive fidelity of coordination number evolution, bulk modulus, and failure strength (Thakur et al., 2015).

4. Contact-Regime Evolution and Frictional Sliding Criteria

Individual contact histories typically involve transitions between elastic loading, elasto-plastic unloading, cohesive pull-off, and frictional slip:

1. Approach: $\delta = 0$, initial adhesive force $f_0$ acts.
2. Elastic loading: $\delta > 0$, branch-1 with stiffness $k_1$ up to $\delta_{max}$.
3. Reverse motion: Unloading follows branch-2 ($k_2$) from $\delta_p$.
4. Adhesive plateau: Further unloading reaches pull-off at $-k_{adh} \delta_p$.
5. Tangential friction: Tangential forces accumulate versus the Coulomb bound; sliding initiated when this is exceeded.

The frictional sliding threshold incorporates adhesion: sliding occurs at

$$|\mathbf{f}_t| = \mu (f_{hys} + k_{adh} \delta_p ) - f_0$$

(Thakur et al., 2015). This construction, verified in experiments and simulations, enables accurate modeling of stick–slip and failure transitions during loading cycles.

5. Model Calibration, Benchmarking, and Bulk Validation

Calibration of DEM contact models proceeds by matching stress–strain and stress–porosity data from fine-particle systems in confined and unconfined compression, treating these as reference ("ground truth"). Parameter choices in (Thakur et al., 2015) include:

• Friction $\mu = 0.5$, rolling friction $\mu_r = 0.001$,
• Bilinear spring constants $k_1 = 5\!-\!10 \times 10^3$ N/m, $k_2 = 2k_1$, $k_{adh} = 5\!-\!7.5 \times 10^3$ N/m,
• Pull-off force $f_0$ set for target tensile strengths ($100\!-\!200$ Pa).

Unscaled stiffness led to significant underprediction of stiffness, strength, and coordination in up-scaled particles. Implementing the above scaling, the macroscopic stress–strain, porosity, strength at failure, and coordination number curves "collapsed" across particle sizes (within a few percent) in both confined and unconfined tests (Thakur et al., 2015). Wall-clock simulation speedups by factors of 7–10 were demonstrated under equivalent physics.

6. Computational Considerations and Practical Implementation

For practical DEM codes, the above bilinear adhesive–frictional contact law can be implemented using incremental updates of normal overlap and tangential history variables for each contact. The sliding and adhesive branches are detected using branch conditions on the overlap and the tangential displacement, with transitions managed using stored maxima and plastic offsets. Scaled-up simulations require updating all relevant model scalings for new particle diameters and verifying that the deformation regime remains within the range tested for validity.

Parameter choices must also respect the stability and dynamic (inertia) constraints of the DEM integration scheme. For quasi-static applications, the time step is typically set according to the highest stiffness and smallest particle mass,

$\Delta t < 0.2\sqrt{m_{min}/k_{max}}$

where $m_{min}$ and $k_{max}$ are the smallest mass and largest contact stiffness in the system.

7. Significance and Applicability

The DEM contact model, particularly in the bilinear adhesive–frictional (hysteretic) form with empirically justified scaling, allows computation of mechanical properties, structural evolution, and failure in granular, powder, and weakly–cohesive solids at laboratory and industrially relevant scales. Such physically justified parameter scaling is essential for bridging length and time scales in multi-scale DEM frameworks, accelerating simulations without compromising the integrity of mechanical predictions across the quasi-static loading regime (Thakur et al., 2015). This approach is extensible to a variety of particle shapes, cohesive mechanisms, and loading histories within the limitations of the underlying model assumptions.

References

• "Scaling of discrete element model parameters for cohesionless and cohesive solid" (Thakur et al., 2015)