Soft-Sphere DEM: Principles and Applications

• SSDEM is a simulation method where particles overlap slightly to mimic elastic deformation using spring-dashpot, friction, and cohesive laws.
• It bridges hard-sphere collision models and continuum methods by modeling enduring contacts, sliding, and rotational dynamics in granular materials.
• SSDEM is applied in planetary science, granular flow, and fluid-coupled studies, though it requires careful calibration and small timesteps due to computational demands.

Soft-sphere discrete element method (SSDEM) is an overlap-based form of the discrete element method in which particles are allowed to penetrate slightly during contact and the resulting overlap is interpreted as elastic deformation that generates finite-duration contact forces and torques. In this formulation, particles remain kinematically rigid or near-rigid while contact mechanics are regularized through spring-dashpot, frictional, and sometimes cohesive laws; enduring contact, sliding, rolling, twisting, and contact-history effects can therefore be represented directly. Across the literature, SSDEM appears both as a granular dynamics method in its own right and as a bridge between hard-sphere collision models and continuum descriptions, especially when the research question depends on path-dependent deformation, localized failure, or particle-scale force transmission (Zhang et al., 2023, Kim et al., 28 Feb 2025, Klerk et al., 2021).

1. Conceptual scope and distinction from related methods

SSDEM is defined less by a single constitutive law than by a contact paradigm. The defining pattern is direct: contact is active over finite overlap, contact forces are computed from overlap geometry and relative contact kinematics, and the interaction persists over multiple timesteps rather than being treated as an instantaneous impulse. This is why the method is repeatedly contrasted with hard-sphere DEM (HSDEM), whose collision treatment is instantaneous and therefore poorly suited to enduring, frictional, many-contact mechanics in dense or quasi-static aggregates (Zhang et al., 2020, Zhang et al., 2017).

That distinction has methodological consequences. In rubble-pile asteroid problems, HSDEM could reproduce gross breakup behavior, but it could not represent the persistent, frictional, many-contact mechanics that govern deformation, shear resistance, dilation, reshaping, and reaccumulation during spin-up or tidal encounters (Zhang et al., 2020). In quasi-static creep problems, SSDEM is preferred because it can sustain persistent contacts, tangential friction, and rotational resistance, making it appropriate for granular creep and failure rather than only for collision sequences (Zhang et al., 2017). A related contrast arises with continuum theory: limit analysis, Mohr-Coulomb or Drucker-Prager yield conditions, and static spin-shape constraints can delimit failure domains, but they do not directly give the dynamical path by which failure localizes, fragments separate, or remnant shapes and spins emerge (Hirabayashi et al., 2015, Zhang et al., 2020).

SSDEM therefore occupies an intermediate position. It is more dynamically explicit than static continuum theory and more faithful to granular contact physics than hard-sphere collision models. A plausible implication is that its most distinctive scientific value appears when geometry, contact history, and dissipation jointly determine outcome rather than merely set perturbative corrections to otherwise ballistic motion.

2. Contact mechanics and geometric generalizations

The literature shows substantial diversity in constitutive detail while preserving the same soft-contact structure. In one representative sphere-based formulation for self-gravitating rubble piles, the normal contact force is decomposed into Hertzian elasticity, normal damping, and cohesion,

$$\vec{\bf f}_e= k_n\xi^{3/2}{\bf\hat n}, \qquad \vec{\bf f}_d=-\gamma_n\dot\xi{\bf\hat n}, \qquad \vec{\bf f}_n=\vec{\bf f}_e+\vec{\bf f}_c+\vec{\bf f}_d,$$

with tangential force truncated by the local Coulomb condition

$|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$

and rolling resistance limited by

$M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$

In that implementation, cohesion enters as an explicit tensile contact force rather than as a continuum cohesive strength parameter (Hirabayashi et al., 2015).

In other implementations the normal law is linear rather than Hertzian. The Hookean DEM used in the variational-integrator formulation writes the normal and tangential forces as

$$\mathbf{F}_{n_{ij}}= k_n \delta_{ij}\,\mathbf{n}_{ij} - \gamma_n m_\mathrm{eff}\,\mathbf{v}_{n_{ij}}, \qquad \mathbf{F}_{t_{ij}}= k_t \delta_{ij}\,\mathbf{t}_{ij} - \gamma_t m_\mathrm{eff}\,\mathbf{v}_{t_{ij}},$$

with $\delta_{ij}=d-|\mathbf r_{ij}|$ for equal-diameter spheres (Klerk et al., 2021). DEM-Engine and the related Chrono clump simulator instead use a history-dependent Hertz–Mindlin law in which

$${F}_n = f(\bar{R},\delta_n)\left(k_n {u}_n - \gamma_n \bar{m} {v}_n\right), \qquad {F}_t = f(\bar{R},\delta_n)\left(-k_t {u}_t - \gamma_t \bar{m} {v}_t\right), \qquad \|{F}_t\|\le \mu \|{F}_n\|,$$

with $f(\bar R,\delta_n)=\sqrt{\bar R\delta_n}$, so the Hertzian $\delta_n^{3/2}$ scaling is recovered through the factor $\sqrt{\bar R\delta_n}$ (Zhang et al., 2023, Zhang et al., 2023). TinyDEM follows the same visco-elastic Hertz–Mindlin pattern and extends torque exchange to sliding, rolling, and twisting friction within a compact spherical-particle code (Vetter, 16 Jul 2025).

The geometric notion of “overlap” is also broader than center-to-center sphere penetration. ParticLS generalizes SSDEM to arbitrary convex bodies described by signed distance functions $\varphi$, and defines contact points by solving a constrained optimization problem on the two zero level sets. The penetration bond

$|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$0

then plays the role usually played by scalar sphere overlap, and its norm gives the penetration distance (Davis et al., 2022). In spheropolygon DEM, overlap is defined by vertex-edge distance,

$|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$1

so the soft-contact principle survives even when particle geometry is non-spherical and contact is no longer reducible to line-of-centers kinematics (Jiang et al., 2019). This suggests that SSDEM is best understood as a compliant-contact framework rather than a sphere-specific formula.

Formulation	Contact geometry	Representative law
Rubble-pile sphere SSDEM	Sphere-sphere overlap	Hertzian normal force, damping, Coulomb tangential force, rolling resistance, optional cohesion (Hirabayashi et al., 2015)
Level-set DEM	Boundary-point optimization on signed distance functions	Linear elastic or linear viscoelastic force over penetration bond $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$2 (Davis et al., 2022)
Spheropolygon DEM	Vertex-edge distance in Minkowski-sum geometry	Linear normal/tangential springs with Coulomb limit (Jiang et al., 2019)

3. Time integration, contact detection, and software realizations

Many SSDEM implementations remain explicit because contact stiffness and short contact timescales make force evaluation the dominant cost. In asteroid applications, PKDGRAV integrates gravity and contact forces with a second-order leapfrog scheme (Zhang et al., 2017). ParticLS uses second-order Runge-Kutta for generic state objects and second-order velocity Verlet for rigid bodies, while storing tangential-force history as a state variable rather than hard-coding a first-order update (Davis et al., 2022). TinyDEM uses semi-implicit Euler for translation and rotation, together with quaternion-based orientation updates by the SPIRAL method and a linked-cell contact search over a Cartesian grid of cell size $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$3 (Vetter, 16 Jul 2025).

The timestep restriction remains a central numerical issue. DEM-Engine states explicitly that DEM requires very small timesteps, typically $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$4–$|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$5, because elements are small and stiff (Zhang et al., 2023). The same stiffness–timestep tradeoff appears in TinyDEM, where the default timestep is tied to elastic-wave propagation and must often be reduced further when high friction stiffens the tangential response (Vetter, 16 Jul 2025). In the variational-integrator work, the implicit scheme does not remove the practical requirement to resolve fast contact dynamics; it instead reformulates the DEM update as a discrete stationarity problem and, in the quasi-static limit, reduces to potential-energy minimization (Klerk et al., 2021).

Large-scale SSDEM has driven specialized software architectures. DEM-Engine, a Project Chrono submodule, is explicitly soft-sphere, dual-GPU, and built around a delayed contact detection algorithm that decouples contact detection from force computation using an Active-Contact Set; the solver adaptively tunes the refresh interval so the dynamics thread seldom waits for the kinematics thread (Zhang et al., 2023). A related Chrono-based simulator for clump-shaped grains uses two GPUs concurrently, assigns one GPU to kinematics and one to dynamics, and resolves non-spherical particles as rigid clumps of overlapping spheres while retaining Hertz–Mindlin sphere-level contact (Zhang et al., 2023). ParticLS takes a different route: geometry, contact mechanics, particle state, and time integration are cleanly separated in an object-oriented framework, so new shapes are introduced through new signed-distance geometries and new constitutive laws through ContactForce classes rather than by rewriting the entire solver (Davis et al., 2022).

The software ecosystem reflected in these works is correspondingly heterogeneous. PKDGRAV dominates self-gravitating asteroid and planetesimal applications; ParticLS emphasizes arbitrary geometries and shared infrastructure with peridynamics; DEM-Engine and the Chrono clump simulator emphasize GPU throughput and user-defined contact kernels; TinyDEM provides a minimal C++11 reference implementation with OpenMP parallelization (Zhang et al., 2023, Davis et al., 2022, Vetter, 16 Jul 2025).

4. Coupling to continuum, fluid, and multiscale models

A recurring pattern in the literature is to use SSDEM where localized contact mechanics matter, while embedding it in broader multiphysics frameworks for stresses, fluids, or deformable continua. In asteroid surface-shedding studies, SSDEM acts as the dynamical counterpart to an analytical limit-analysis model: the continuum framework provides upper and lower failure bounds, while SSDEM resolves the actual deformation path, local rearrangement, and particle lofting during quasi-static spin-up (Hirabayashi et al., 2015). This division of labor is methodologically important because the continuum model supplies admissible failure domains, whereas SSDEM determines where failure nucleates, whether it is local or global, and how failed material moves.

Fluid-coupled methods adapt the same idea. In bonded-sphere CFD–DEM for sediment transport, irregular grains are constructed from a small number of non-overlapping component spheres, all contact detection and contact-force calculations remain sphere-sphere, and fluid forces are computed and applied on each sphere individually; the total force and torque are then aggregated to the rigid composite grain (Sun et al., 2016). In unresolved LBM–DEM, particle motion follows

$|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$6

so dry-contact collision force, hydrodynamic force, and explicit lubrication correction coexist in one force balance (Rettinger et al., 2017). In two-way coupled SPH–DEM, the DEM side may be deliberately minimal—here a linear spring-dashpot normal contact law—but porosity smoothing, drag closure, and conservative force transfer between phases become the dominant technical issues (Robinson et al., 2013).

A further extension appears in hybrid MPM–SDEM. There, discrete near-rigid objects are represented by spheropolygons and the continuum side by material points; contact is computed through DEM-style overlap and friction, and the resulting force is projected back to the MPM grid as an external nodal force (Jiang et al., 2019). This suggests a broader interpretation of SSDEM: not merely a standalone granular solver, but a reusable contact-mechanics module for multiscale methods where rigid or near-rigid objects interact with continua or fluids.

5. Scientific applications and characteristic findings

The most developed SSDEM application domain in the supplied literature is planetary granular mechanics. In studies of rotational disruption of rubble-pile asteroids, SSDEM is used to test the strong-core/weak-shell hypothesis for surface shedding. A layered, self-gravitating aggregate with stronger core contacts and weaker shell contacts shows that critical spin for irreversible failure rises with core size, that failure begins locally rather than globally, and that larger cores produce smaller lofted components (Hirabayashi et al., 2015). In a related Didymos study, SSDEM with quasi-static spin-up yields three critical spin limits associated with reshaping and surface shedding, internal structural deformation, and shear failure, and shows that stability depends strongly on particle arrangement, size distribution, and density structure rather than on nominal friction parameters alone (Zhang et al., 2017).

Tidal encounter and low-speed merger problems make different use of the same method. In Earth flyby simulations of rubble piles, SSDEM replaces earlier hard-sphere treatments and produces closer disruption thresholds, lower relative fragment speeds, and highly elongated remnants with axis ratios near $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$7, while also showing that random close packing behaves more naturally than hexagonal close packing (Zhang et al., 2020). In Arrokoth merger simulations, SSDEM constrains the final contact to be very slow and grazing for rubble-pile progenitors, with intact-lobe outcomes requiring impact speeds less than $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$8 and impact angles greater than $|\vec{\bf f}_t|\leq\mu |\vec{\bf f}_n|,$9 in the direct-impact scenario, while gentle synchronous inspiral gives the most plausible contact-binary morphology (Marohnic et al., 2020). In pebble-cloud gravitational collapse, the replacement of perfect-merger super-particles by PKDGRAV SSDEM allows super-particles to rest on one another, eliminates the need for inflated radii, produces many binary planetesimals from one cloud, and yields new results on remnant shapes and spins, including approximately 10-hour rotation periods and a broad morphological distribution ranging from spherical and oblate to top-shaped and prolate (Barnes et al., 22 Jul 2025).

Outside planetary science, SSDEM appears as a benchmark and model-selection tool for terrestrial granular flow. In comparative DEM–SPH simulations of granular collapse and rotating drums, the DEM side uses Hertz–Mindlin contact with overlap-based normal and tangential forces, rolling resistance, and restitution-dependent damping; the study concludes that DEM is more versatile for complex flow regimes even though calibration is computationally expensive and non-unique (Kim et al., 28 Feb 2025). In rover-terrain and digital-simulant studies, clump-based SSDEM reproduces repose angle, cone penetration, drawbar pull, and slope–slip behavior, and shows that single-sphere terrain still shears too easily even when the sphere friction coefficient is increased from $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$0 to $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$1, indicating that geometric interlocking cannot always be replaced by friction tuning alone (Zhang et al., 2023). DEM-Engine extends the same granular mechanics to tens or hundreds of millions of sphere primitives on two GPUs, validating contact-law behavior in rolling, stacking, hopper, and impact tests while keeping the stock physics within the damped Hertz–Mindlin family (Zhang et al., 2023).

6. Limitations, calibration issues, and recurrent misconceptions

Several recurrent misconceptions are corrected by the cited literature. SSDEM is not synonymous with a single Hertz–Mindlin law: linear spring-dashpot models, viscoelastic Hertzian laws, cohesive variants, and arbitrary-shape level-set or spheropolygon formulations all appear within the same soft-contact paradigm (Davis et al., 2022, Klerk et al., 2021). Nor does SSDEM imply realistic grain shape by default. Many planetary applications still use spherical or equal-sized particles, with rolling friction, twisting friction, or a shape parameter such as $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$2 introduced as a proxy for irregular grain interlocking (Zhang et al., 2020, Hirabayashi et al., 2015). Clump methods mitigate this, but at increased memory, contact-detection, and calibration cost (Zhang et al., 2023).

Calibration remains a persistent difficulty. In the DEM–SPH comparison study, the uncertain microscopic parameters are $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$3, $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$4, and $M_r^m=\mu_r R_r |{\vec{\bf f}_n}|.$5, and calibration is described as challenging, non-unique, and computationally expensive because these are microscopic rather than directly measured bulk parameters (Kim et al., 28 Feb 2025). The same paper shows that a collapse test can be insensitive to restitution while a rotating drum is sensitive to it, so parameter identifiability is strongly problem-dependent. A common practical implication is that SSDEM calibration must be benchmark-specific rather than assumed transferable across flow regimes.

Resolution and implementation details also matter. The asteroid surface-shedding study uses only 3000 particles, notes that shell thickness and core boundary are relatively coarse at particle scale, and does not report timestep, integrator, or detailed convergence tests, leaving implementation uncertainty for exact reproduction (Hirabayashi et al., 2015). The tidal-disruption study emphasizes that ordered HCP packings introduce orientation-sensitive crystalline artifacts not seen in random close packing, and that increasing particle count can alter frictional resistance in ways opposite to some hard-sphere results because contact mechanics are modeled explicitly (Zhang et al., 2020). The Didymos work goes further and argues that failure mode and mechanism are mainly affected by internal configuration, whereas shear strength depends on both configuration and material parameters (Zhang et al., 2017).

The computational burden is structural rather than incidental. Explicit SSDEM requires small timesteps tied to contact stiffness, and even sophisticated GPU implementations remain dominated by contact search, force evaluation, and history-variable maintenance (Zhang et al., 2023, Vetter, 16 Jul 2025). A plausible summary is that SSDEM is most informative precisely in the regimes where it is most expensive: dense, frictional, path-dependent systems in which enduring contact and localized granular rearrangement cannot be reduced to impulse-based collision rules or continuum yield envelopes alone.