Viscoelastic Soft-Sphere Collision Model

• The viscoelastic soft-sphere collision model is defined by a finite overlap and time-dependent Hertzian contact force combined with velocity-based damping.
• It models deformable collisions using a nonlinear ζ^(3/2) force law with dissipation, yielding restitution coefficients below unity and non-instantaneous collision times.
• The model’s one-parameter similarity simplifies calibration for discrete element methods, enhancing simulation efficiency in granular kinetic and continuum-to-particle analyses.

A viscoelastic soft-sphere collision model represents contact between deformable bodies through a finite overlap variable and a finite-duration normal force, rather than through an instantaneous impulse. In the formulations most directly associated with spherical particles, the elastic backbone is Hertzian, scaling as the overlap to the power $3/2$, while dissipation is introduced through a velocity-dependent term. The resulting collision law produces a restitution coefficient below unity, a nontrivial collision time, and a contact dynamics that can be reduced, in important cases, to a one-parameter dimensionless family. This structure makes the model central to discrete element methods (DEM), adaptive collision modeling, granular kinetic theory, and continuum-to-particle links for rate-dependent granular dissipation (Ray et al., 2015).

1. Governing contact law and collision kinematics

For normal collisions of spherical particles, a standard macroscopic formulation uses the overlap $\zeta(t)$ and writes the relative motion as

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$

Here $m_{\rm p}$ is particle mass, $d$ is a linear damping coefficient, and $k$ is the Hertz stiffness parameter. The nonlinear restoring force is $k\zeta^{3/2}$, while dissipation is represented by a term linear in the relative normal velocity. The dry coefficient of restitution is defined by

$$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$

with $T_{\rm c}$ the collision duration (Ray et al., 2015).

A more microscopic viscoelastic sphere model, used for smooth, frictionless, non-adhesive head-on collisions, expresses the normal force as

$$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$

where $\zeta(t)$0 is the compression, $\zeta(t)$1, $\zeta(t)$2 encodes elastic and geometric parameters, and $\zeta(t)$3 is a dissipative material constant. In this class of models the restitution coefficient is velocity-dependent rather than constant, because the dissipative force depends explicitly on $\zeta(t)$4 (Mueller et al., 2011).

Rigorous continuum-mechanics derivations for low-velocity collisions of viscoelastic convex bodies recover the same structural form: an elastic Hertz contribution plus a dissipative term proportional to $\zeta(t)$5. One expression is

$\zeta(t)$6

with $\zeta(t)$7 the Hertz prefactor and $\zeta(t)$8 a material-dependent coefficient obtained from elastic and viscous constants. A closely related derivation gives

$\zeta(t)$9

This suggests equivalence up to sign and notation conventions for the compression rate $$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$0 (Goldobin et al., 2015).

These formulations are distinct from linear spring-dashpot surrogates, in which the normal contact force is proportional to overlap rather than $$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$1, and from purely elastic soft-sphere models in which the damping term is absent.

2. Dimensional reduction and one-parameter similarity

A notable property of the damped Hertzian soft-sphere model is that, after an appropriate choice of units, the collision dynamics collapse to a single dimensionless parameter. Using the scalings

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$2

the governing equation becomes

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$3

with

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$4

The dimensionless collision time is $$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$5, and the restitution coefficient is

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$6

Therefore the physical problem is reduced to a one-parameter family indexed by $$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$7, with dynamic similitude across all collisions that share the same value of this parameter (Ray et al., 2015).

The parameter combination may be summarized as

$$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$8

It combines material stiffness $$m_{\rm p}\ddot{\zeta} = - d\dot{\zeta} - k\zeta^{3/2}, \qquad \zeta(0)=0,\quad \dot{\zeta}(0)=u_{\rm in}.$$9, damping $m_{\rm p}$0, particle mass $m_{\rm p}$1, and impact velocity $m_{\rm p}$2. The significance of this reduction is practical as well as analytical: restitution and collision time can be tabulated, approximated, or inverted as functions of a single scalar instead of four separate physical inputs (Ray et al., 2015).

A related organization appears in the viscoelastic extension of Hertz theory for central collisions of isotropic homogeneous spheres, where the deviations from the elastic Hertz limit are controlled by

$m_{\rm p}$3

In that setting, $m_{\rm p}$4 recovers the purely elastic limit, while nonzero $m_{\rm p}$5 produces a longer collision time, smaller maximum contact area, and $m_{\rm p}$6. This parallel emphasizes that one-parameter similarity is not peculiar to the linear-dashpot Hertz surrogate but is a recurring feature of weakly dissipative viscoelastic contact models (He et al., 2013).

3. Formal solution, asymptotics, and explicit inversion

The reduced damped Hertzian model does not admit a closed-form trajectory $m_{\rm p}$7, but it does admit an exact formal series solution in $m_{\rm p}$8. The analysis is performed in phase space by treating velocity as a function of penetration, $m_{\rm p}$9, so that

$d$0

With

$d$1

the model is transformed into coupled initial-value problems for $d$2 and $d$3. The authors then expand

$d$4

with zeroth-order term

$d$5

which is exactly the conservative Hertzian result. Maximum penetration $d$6, restitution, and collision time follow from $d$7, $d$8, and

$d$9

The small-$k$0 assumption corresponds to weak damping; the range $k$1 with $k$2 is reported as typically safe. The first-order correction to collision time is technically delicate because differentiation under the integral sign produces divergent-looking terms near the turning point, so the regularization must be handled before evaluation (Ray et al., 2015).

For implementation, the formal asymptotics are replaced by compact approximations. The restitution coefficient is represented as

$k$3

and the collision time as

$k$4

where

$k$5

The corresponding small-$k$6 series are

$k$7

$k$8

$k$9

The same approximations yield a direct inverse calibration for the adaptive collision model of Kempe and Fröhlich. Given a target $k\zeta^{3/2}$0 and $k\zeta^{3/2}$1, one first defines $k\zeta^{3/2}$2 and computes

$k\zeta^{3/2}$3

taking the positive root. Then

$k\zeta^{3/2}$4

This converts the inverse problem from an iterative root-find into explicit algebraic formulas (Ray et al., 2015).

4. Restitution as a velocity-dependent constitutive output

In viscoelastic soft-sphere mechanics, the coefficient of restitution is not a prescribed constant but an output of the contact law. For viscoelastic spheres, Schwager and Pöschel derived the exact infinite series

$k\zeta^{3/2}$5

with

$k\zeta^{3/2}$6

The coefficients $k\zeta^{3/2}$7 are universal, so all material dependence is absorbed into $k\zeta^{3/2}$8. Although exact, the series converges extremely slowly: about 20 terms are needed to obtain $k\zeta^{3/2}$9 up to quadratic order in $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$0, and any finite truncation diverges for sufficiently large $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$1. To obtain simulation-ready formulas, Padé approximants were introduced. The paper identifies $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$2 as a good compromise of simplicity and accuracy, $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$3 as limited by a pole near $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$4, and $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$5 as nearly perfect over the full velocity range (Mueller et al., 2011).

The dependence of collision observables on dissipation is not uniform. In the viscoelastic extension of Hertz theory, weak-dissipation expansions give

$$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$6

$$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$7

$$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$8

This model explains why Hertz’s elastic theory remains accurate for predicting collision time and maximum contact area even when $$e_{\rm dry}=-\frac{u_{\rm out}}{u_{\rm in}} =-\frac{\dot{\zeta}(T_{\rm c})}{u_{\rm in}},$$9 of the kinetic energy is lost due to viscous dissipation: the collision time and maximum contact area have a very weak dependence on impact velocity, and the combinations

$T_{\rm c}$0

have vanishing linear corrections in $T_{\rm c}$1 (He et al., 2013).

For dilute granular gases of viscoelastic particles, the same velocity dependence of restitution carries over to kinetic theory. There the simplest first-principles viscoelastic model gives

$T_{\rm c}$2

and an effective constant restitution coefficient $T_{\rm c}$3 may be defined as a collision average. The reduced model describes granular temperature accurately, but for large dissipation it fails quantitatively for the detailed velocity-distribution corrections represented by Sonine coefficients (Dubey et al., 2012).

5. DEM calibration, numerical efficiency, and alternative implementations

The immediate computational value of the damped Hertzian model lies in parameter calibration. In the adaptive collision model framework of Kempe and Fröhlich, the goal is to prescribe a desired dry restitution coefficient and collision time, then compute the stiffness and damping for each collision. The explicit inversion described above removes the need for quasi-Newton or Broyden-type iteration. For binary collisions, reported restitution errors are typically around $T_{\rm c}$4 for $T_{\rm c}$5, with a maximum relative error about $T_{\rm c}$6 in the practical range; even at $T_{\rm c}$7, the restitution error is about $T_{\rm c}$8. Collision-time errors stay around $T_{\rm c}$9. In a multiple-particle sedimentation test with 100 falling particles colliding with 195 fixed particles in hexagonal packing, the approximate direct method and the original iterative ACTM give almost indistinguishable bulk behavior and essentially the same energy partitioning, while the total CPU time is reduced substantially; one reported case drops from roughly $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$0 s to $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$1 s, with the coefficient-calculation time dropping from $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$2 seconds to $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$3 seconds (Ray et al., 2015).

A different implementation strategy appears in fully resolved simulations of particle-laden flows. There the solid-solid contact is modeled by a linear soft-sphere spring-dashpot law,

$$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$4

with

$$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$5

The collision time is prescribed as $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$6, so the contact can be artificially stretched over a chosen number of fluid time steps. Tangential contact is treated with an analogous spring-dashpot law plus a Coulomb friction cap. The authors describe this as a linearized version of Hertzian contact theory, computationally attractive because it reduces the number of time steps required to integrate the collision force accurately, provided the prescribed collision time remains much smaller than the characteristic timescale of particle motion (Costa et al., 2015).

The phrase “soft-sphere model” can also denote a purely elastic harmonic approximation rather than a dissipative viscoelastic law. In a stochastic collision model for highly rarefied gases under weak gravity, the near-contact two-body potential is expanded around an equilibrium separation $$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$7 to obtain

$$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$8

The paper notes the more general DEM form

$$F=\min\!\left(0,\,-\rho\,\xi^{3/2}-\frac{3}{2}A\rho\sqrt{\xi}\,\dot{\xi}\right),$$9

but sets $\zeta(t)$00, so the collisions are fully elastic with $\zeta(t)$01. Collisions are further randomized through a control parameter $\zeta(t)$02, yielding an effective pairwise collision probability $\zeta(t)$03 (Tsuzuki, 2024).

6. Validity limits, rigorous corrections, and broader extensions

The damped Hertzian one-parameter model is intended for normal collisions of spheres, weak to moderate damping, restitution not too small, and practical DEM or particle-transport settings in which a computationally cheap calibration of $\zeta(t)$04 and $\zeta(t)$05 is needed. The continuum-mechanics derivations impose an overlapping but more specific regime: low impact velocity, no plastic deformation or fragmentation, homogeneous viscoelastic continua, arbitrary convex smooth shapes, and a small ratio of material relaxation time to collision duration. In this regime, inertia is neglected to leading order while viscous effects are retained perturbatively (Goldobin et al., 2015).

A central correction supplied by the rigorous perturbative theory is that the first-order elastic response induced by viscous deformation must be kept together with the viscous stress. The quasi-static approximation, which uses the viscous stress evaluated on the static Hertz field but drops the associated first-order elastic correction, is incomplete because the omitted elastic correction is of the same order as the viscous contribution. For two different materials, the quasi-static treatment can even violate Newton’s third law by producing unequal first-order stresses on the two bodies. The rigorous solution resolves this inconsistency and also removes the unphysical vanishing of dissipation that the quasi-static formula predicts in limiting cases such as $\zeta(t)$06 (Brilliantov et al., 2014).

Beyond the quasi-static regime, internal elastic modes can qualitatively alter restitution. Numerical simulations of normal head-on collisions of isothermal viscoelastic spheres show that, for small solid viscosity, the restitution coefficient oscillates as a function of impact speed. The oscillation arises from resonance between the contact duration and the eigen-frequencies of the sphere; it disappears when the solid viscosity is sufficiently strong. Under thermal activation, the restitution coefficient can exceed unity if the impact speed is nearly equal to or slower than the thermal speed (Murakami et al., 2013).

A different regime is the viscous-dominant collision of soft polymeric materials and viscoelastic droplets. There the normal deformation obeys

$\zeta(t)$07

which, after scaling, becomes

$\zeta(t)$08

For the parameter choice adopted in that analysis, this reduces to

$\zeta(t)$09

The forward perturbation solution diverges near the turning point, so the restitution coefficient is obtained by matching forward and inverse collision expansions. In that viscous-dominant regime, the velocity scaling differs from the classic elastic-dominant $\zeta(t)$10 law associated with Brilliantov-type results (Kim, 2013).

At the continuum scale, restitution-derived viscoelastic damping has been used to relate a particle-scale restitution coefficient to a continuum bulk viscosity. In a Kelvin–Voigt-type granular continuum model, the viscous stress is written as

$\zeta(t)$11

and the bulk viscosity is estimated from restitution through

$\zeta(t)$12

That framework explicitly separates viscoelastic damping, which governs wave propagation and collision-like oscillations, from $\zeta(t)$13-based plastic flow, which governs dense frictional deformation. A plausible implication is that the soft-sphere collision model now functions not only as a DEM contact rule but also as a constitutive bridge between particle-scale dissipation and macroscopic rate-dependent granular behavior (Chandra et al., 23 Apr 2026).