# Multiple-Impact Modeling in Multibody Systems

Bernard Brogliato

► To cite this version:

Bernard Brogliato. Multiple-Impact Modeling in Multibody Systems. A.F. Vakakis; O. Gendelman. Handbook on Nonlinear Dynamics, Vibrations and Acoustics; Volume 1: Nonlinear Dynamics and Vibrations: Fundamental Concepts and Analytical Methods, World Scientific Publising, pp.1-39, 2025. hal-04673788

HAL Id: hal-04673788

https://inria.hal.science/hal-04673788v1

Submitted on 20 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiple-Impact Modeling in Multibody Systems

Bernard Brogliato*

August 20, 2024

## 1 Multibody Complementarity Systems Dynamics

This chapter deals with the modeling of simultaneous collisions in rigid multibody systems, i.e., mechanical systems made of rigid bodies connected by unilateral constraints. In a Lagrangian formalism these systems may be written generically as follows (for simplicity of the presentation bilateral constraints are not considered):

$$M(q)\ddot{q} + F(q, \dot{q}, t) = \nabla h_{nu}(q)\lambda_{nu} + H_{tu}(q)\lambda_{tu} \quad (1)$$

$$\text{Impact model (restitution law)} \quad (2)$$

$$\text{Friction model (tangential effects)} \quad (3)$$

$$\text{Complementarity conditions: } 0 \le h_{nu}(q) \perp \lambda_{nu} \ge 0 \quad (4)$$

where $q \in C_{conf} \subseteq \mathbb{R}^n$ is the vector of generalized coordinates $q_i$, $1 \le i \le n$, assumed to be independent when all constraints are deactivated, i.e., $h_{nu}(q) > 0$, $C_{conf}$ is the system's configuration space, $\dot{q}$ is the vector of generalized velocities, $M(q) = M(q)^T$ is the mass matrix, always assumed to be at least positive semi-definite (it may be assumed non-singular in most cases), $F(q, \dot{q}, t)$ collects internal forces (including forces deriving from a potential, plus Coriolis and centrifugal forces), as well as external actions on the system such as disturbances or control. We consider $m$ unilateral constraints (i.e., potential contact points), with which are associated gap functions or signed distances: $h_{nu,i} : C_{conf} \mapsto \mathbb{R}$, $1 \le i \le m$. The gradient matrix $\nabla h_{nu}(q) = (\nabla h_{nu,1}(q), \dots, \nabla h_{nu,m}(q)) \in \mathbb{R}^{n \times m}$. The vector of

<footnote>*Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France. bernard.brogliato@inria.fr</footnote>
Lagrange multipliers $\lambda_{nu} \in \mathbb{R}^m$ corresponds to the normal components of the contact forces, and is obtained from the local contact kinematics [51, Chapter 10] [4, Chapter 3]. The complementarity condition (4) is to be understood componentwise (i.e., per contact $i$). It models the fact that for each contact $i$, the normal contact force satisfies ($h_{nu,i}(q) > 0$ implies $\lambda_{nu,i} = 0$), and ($\lambda_{nu,i} > 0$ if and only if $h_{nu,i}(q) = 0$). The modeling of impacts will be presented in section 4. In section 3.2 Coulomb's friction is introduced. The matrix $H_{tu}(q)$ is defined from the local contact kinematics, and relates local tangential velocities $v_{t,i}$ to the generalized velocity $\dot{q}$ as $v_{t,i} = H_{tu,i}^{\top}(q)\dot{q}$, $1 \le i \le m$. Same holds for the normal component of the local velocities $v_{n,i}$, with $v_{n,i} = \nabla h_{nu,i}^{\top}(q)\dot{q}$, $1 \le i \le m$, and we denote $U_n = \nabla h_{nu}^{\top}(q)\dot{q} = (v_{n1}, \dots, v_{nm})^{\top}$, $U_t = H_{tu}^{\top}(q)\dot{q} = (v_{tu,1}^{\top}, \dots, v_{tu,m}^{\top})^{\top} \in \mathbb{R}^{md}$, $d=1$ or $d=2$. In the following CoR is for coefficient of restitution, CoF is for coefficient of friction. A useful notion is that of normal cone to a closed convex nonempty set $\Phi \subset \mathbb{R}^n$: $\mathcal{N}_{\Phi}(x) = \{z \in \mathbb{R}^n | z^{\top}(y-x) \le 0, \forall y \in \Phi\}$. Examples are depicted in Fig. 3 and 4(a). The orthogonal projection of $x \in \mathbb{R}^n$ on $\Phi \subset \mathbb{R}^n$ in the metric defined by $M = M^{\top} \succ 0$ is denoted as $\text{proj}_M[\Phi; x] = \underset{z \in \Phi}{\text{argmin}} \frac{1}{2}(z-x)^{\top}M(z-x)$. The boundary of a closed set $\Phi$ is denoted as $\text{bd}(\Phi)$. For a convex proper lower semicontinuous function $f: \mathbb{R}^n \to \mathbb{R}$, $\partial f$ denotes its subdifferential [92].

**Remark 1** Bilateral constraints of the form $h_{nb}(q) = 0$ are usually present in multibody systems. They are not considered in this work for simplicity of the presentation. There are several ways to cope with bilaterally constrained systems, the literature on their analysis and numerical simulation is hudge.

## 2 Classes of Impact Models and Impact Dynamics

There are several ways to classify impact models within large classes [81]:

(a) Algebraic (zero-order) models which assume instantaneous collisions between perfectly rigid bodies, and relate post- and pre-impact velocities through a restitution rule: $\dot{q}(t^+) = \mathcal{F}(q(t), \dot{q}(t^-))$, CoR, CoF). Positions are constant, velocities undergo a jump. At an impact time $t$, the dynamics (1) are equivalent to the velocity-impulse relation:

$$M(q(t))(\dot{q}(t^+) - \dot{q}(t^-)) = \nabla h_{nu}(q(t))\bar{\lambda}_{nu}(t) + H_{tu}(q(t))\bar{\lambda}_{tu}(t). \quad (5)$$

This can be derived rigorously from Schwartz' distributions [22, section 1.1] or measure theory [79, 75]. To simplify the presentation let us assume that
the acceleration is a measure which satisfies: $d\dot{q} = \ddot{q}(t)dt + ds$ with $s(\cdot)$ a jump function, $\ddot{q}(\cdot)$ the derivative of a continuous function of bounded variation, and at $t$: $ds = (\dot{q}(t^+) - \dot{q}(t^-))\delta_t$, $\delta_t$ the Dirac measure, $dt$ the Lebesgue measure. Thus the multiplier is also a measure denoted as: $\bar{\lambda}_{nu} = \lambda_{nu}(t)dt + d\bar{\lambda}_{nu}$, where $d\bar{\lambda}_{nu} = \bar{\lambda}_{nu}(t)\delta_t$ (same for $\lambda_{tu}$). Then (5) has to be completed with (2) and (3) to calculate $\dot{q}(t^+)$, and the percussions $\bar{\lambda}_{nu}(t)$, $\bar{\lambda}_{tu}(t)$.

(b) First-order dynamics following the Darboux-Keller approach [22, section 4.3.5], where the normal contact force infinitesimal impulse $d\bar{\lambda}_{nu} = \lambda_{nu}(t)dt$ is used as a new time-scale, where $\lambda_{nu}(t)$ is the contact force during the impact (assumed to consist of a compression phase followed by a restitution phase in single impacts). The velocity is the state, positions are assumed constant during the shock. In the simplest case of a particle colliding with an obstacle, $m\frac{d\dot{q}}{dt} = \lambda_{nu}(t)$ is transformed during the collision to

$$m \frac{d\dot{q}}{dt} = \lambda_{nu}(t)dt = d\bar{\lambda}_{nu} \Leftrightarrow \frac{d\dot{q}}{d\bar{\lambda}_{nu}}(\bar{\lambda}_{nu}) = \frac{1}{m}. \quad (6)$$

When integrated over compression and expansion phases, (6) is usually called Routh's incremental collision dynamics. When Coulomb's friction is present, a nonlinear nonsmooth ODE is obtained [22, section 4.3.5.2]. In the Lagrangian framework we get the Lagrangian Darboux-Keller dynamics (an extension of Routh's incremental form):

$$M(q(t_i)) \frac{d\dot{q}}{d\bar{\lambda}_{nu,j}} = \nabla h_{nu}(q(t_i)) \frac{d\bar{\lambda}_{nu}}{d\bar{\lambda}_{nu,j}} + H_{tu}(q(t_i)) \frac{d\bar{\lambda}_{tu}}{d\bar{\lambda}_{nu,j}}, \quad (7)$$

where $t_i$ is the initial time of the collision, and $1 \le j \le m$ has to be suitably chosen (see Group 3, section 4.4). The dynamics (7), which is the counterpart of (5), holds on $[t_i, t_f]$. The shock dynamics (7) has to be completed with (2) and (3) so that it can be integrated on $[t_i, t_f]$. CoRs are used to determine the final time $t_f$.

(c) Second-order dynamics which rely on compliant models of lumped flexibilities and dampings (linear or nonlinear spring-dashpot rheological models), through assemblies of various basic components: elastic, visco-elastic, frictional, fractional elastic, visco-plastic, etc, see section 3.1. The state is made of positions and velocities.

It is worth noting that models of class (c) usually are only piecewise smooth because of detachment conditions. The complementarity conditions (4) stem from basic contact assumptions [22, section 5.4.1], and apply to all three
classes of impacts. Classes (a) and (b) usually neglect nonimpact forces $F(q, \dot{q}, t)$. This allows us to get the velocity-impulse relations (5) and (7). It is possible to combine different approaches in the same system [22, Example 5.14]. On the other hand, shocks between bodies may also belong to:

• 1) very-low-velocity impacts,

• 2) low-velocity impacts,

• 3) high-velocity impacts (collisions such that effects like body vibrations, temperature rise, cannot be neglected),

• 4) hypervelocity impacts (collisions which induce irreversible damage like cracks, craters, failures, high deformation, debris formation, etc).

Ranges of velocities may be given [22, Definition 2.1], though they certainly depend on the colliding bodies (materials, shapes, contact geometries). We may say that the models described in this chapter, belong to classes 1) and 2). One major issue is the modelling of kinetic energy dissipation. Tangential dissipation is usually modelled using extensions of set-valued Coulomb's friction, or with tangential coefficients of restitution (CoF and tangential CoR may be related [22, section 4.3.3.1]).

## 3 Fundamental Properties of Impacts

Classically there are three consistency constraints to be respected at an impact time $t$:

• (i) **Energetic consistency**: an impact does not increase kinetic energy $T(q, \dot{q}) = \frac{1}{2}\dot{q}^\top M(q)\dot{q}$:

$$T(q(t_k), \dot{q}(t^+)) - T(q(t_k), \dot{q}(t^-)) \le 0 \quad (8)$$

$$\iff \frac{1}{2}\bar{\lambda}_{nu}^\top(t)(v_n(t^+) + v_n(t^-)) + \frac{1}{2}\bar{\lambda}_{tu}^\top(t)(v_t(t^+) + v_t(t^-)) \le 0 \quad (9)$$

where the second inequality is known as Thomson and Tait (or Kelvin<fnref n="1" /> and Tait) formula [22, section 4.3.12] [102], which rules the work performed by the impact forces ($t^+$ and $t^-$ in (8) (9) are understood as the end and the beginning of the impact process, hence (8) (9) apply to all three above classes of impacts),

<footnote n="1">¹William Thomson, known as Lord Kelvin.</footnote>
• (ii) **Kinematic consistency:** the post-impact velocity has to be admissible: $U_n(t^+) \triangleq \nabla h_{nu}(q(t))^{\top} \dot{q}(t^+) \geq 0$,

• (iii) **Kinetic consistency:** the impact force impulse has to be non-negative: $\bar{\lambda}_{nu}(t) \geq 0$.

The respect of these three constraints does not imply that an impact model possesses good prediction capabilities, however. In particular, multiple impact models should possess an additional property:

• (iv) **Parameter consistency and admissibility:**

- (a) parameters (CoR, CoF) must be able to span the whole admissible post-impact-velocity subspace, and

- (b) admissibility of parameters means that parameters associated with multiple impact should be valid if some of the constraints are removed. Consider for instance a chain of $n$ aligned beads, and an impact law whose parameters have been estimated. If several subchains are extracted from the longer chain, can the same parameters be used? In other words, can the $p$-impact law parameters be identified from independent experiments with $q$-impacts, $q < p$? Or are they attached to the studied $p$-impact?

Finally, less well-known yet important properties are:

• (v) **Uniqueness of the outcome:** post-impact velocity, impact forces impulses, may not be unique for given data (pre-impact velocity, parameters),

• (vi) **impact termination:** it may not be always obvious to find a criterion which allows designers to compute the end of an impact, especially when some multiple impact models are used.

## 3.1 Normal Dissipation

A plethora of approaches has been proposed so far to take normal dissipation (i.e., the dissipation process associated with local normal displacements) into account, among them we may cite:

• linear viscosity (damping), nonlinear damping as in the well-known Simon-Hunt-Crossley, Kuwabara-Kono, Tsuji-Tanaka-Ishida models or the Lankarani-Nikravesh approach [8, 31][22, section 2.2.2] [108], see [98, Fig. 3] for a survey of typical force-indentation curves, and Fig. 1 (a) (b),
• phenomenological parameters like coefficients of restitution (CoR) for the zero and first-order dynamics models, the most popular ones being kinematic (Newton), kinetic (Poisson) and energetic CoRs, and which can depend on other parameters [22, section 4.2],

• varying or piecewise continuous (PWC) stiffnesses in compression and expansion phases (the bistiffness Crook’s approach [34] [22, section 4.2.1.2], or tristiffness [81], see also various PWC stiffnesses in [22, section 4.2.1.1] [17, 73, 24]) (this models plasticity with residual deformation which rules the energy loss, hence the CoR), see Fig. 1 (c) (e) (f),

• plasticity, elastoplasticity à la Johnson with three or four regimes ([61, 19, 103], see also [22, section 4.2.1.1] and [46]), Masin, Persoz and viscoelasto-plastic assemblies with set-valued frictional elements [112, 14]. In fact these models and those in the previous item, belong to the same class of elasto-plastic dissipation and PWC force/indentation relation, see Fig. 1 (d) (f),

• adhesive effects with Dugdale cohesive model or Lennard-Jones potentials [76] [22, section 4.2.2], body bulk post-impact vibrations (see [53, Chapter III] [61, Chapter 11] [22, section 4.2.4] for references), wear, noise, temperature increase, etc.

Interestingly enough, these various models can be identified through their force/indentation curves, some of which are depicted in Fig. 1. An interesting critical survey is in [109], see also [38].

<figure><img src="image_2.png" /><figcaption>Figure 1: Force/indentation curves. Dissipated energies = dashed areas.</figcaption></figure>
## 3.2 Friction Modeling

Two main classes of friction models exist: Coulomb's friction with set-valued part, and possibly varying CoF ($\mu_i(t, v_{t,i}, v_{n,i})$, e.g., Stribeck effects), or some kind of regularization around $v_{t,i} = 0$ (in 2D: saturation, sigmoid, dead-zone, or dynamic models with internal state to cope with presliding, see comments below). Regularization suppresses the friction cone, which may be quite problematic in many applications. In this chapter only set-valued Coulomb's models are considered. At contact $i$, let us assume that the geometry of the surface of contact is convex smooth enough, so that we can define a tangent plane of contact $S$ and the normal (unit) vector $\mathbf{n}_i$ at $S$. The vector $\mathbf{n}_i$, together with two arbitrary orthonormal vectors $\mathbf{t}_{1,i}, \mathbf{t}_{2,i}$ lying on $S$, define a local basis $\mathcal{B}_i$. We assume that each contact $i$ involves only two parts of the whole mechanical system, and denote by $v_i$ the local (space) relative velocity between the two parts in contact, which can be decomposed in $\mathcal{B}_i$ as a normal part $v_{n,i} = v_i^\top \mathbf{n}_i \in \mathbb{R}$ and a tangential part $v_{t,i} = \{v_i^\top \mathbf{t}_{1,i}, v_i^\top \mathbf{t}_{2,i}\} \in \mathbb{R}^2$ (when the contact $i$ is activated, then $v_{n,i} = 0$). The contact force $\lambda_i$ can be similarly decomposed onto its normal part $\lambda_{n,i} \in \mathbb{R}$ and its tangential part $\lambda_{t,i} \in \mathbb{R}^2$. Let $\mu_i$ be the coefficient of friction at the contact point labelled $i$.

### 3.2.1 Coulomb's Friction: Formalisms

Let us assume that the contact $i$ is active, keeping in mind that for the unilateral contacts (4) holds at all times. Coulomb's friction law, symbolically denoted as $(v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i)$, where $\mathcal{C}(\mathbf{n}_i, \mu_i) \subset \mathbb{R}^3$ (or $\mathbb{R}^2$ in 2D) is the friction cone, $\mu_i$ is the CoF, relates the normal and tangential parts of both the local velocity $v_i$ and the local contact force $\lambda_i = (\lambda_{n,i}, \lambda_{t,i}^\top)$ according to the following formulations [78] [22, §5.3] [4, §3.9.1]:

$$ (v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \qquad (10) $$
$$ \left\{ \begin{array}{ll} \text{either} & \|\lambda_{t,i}\| \le \mu_i |\lambda_{n,i}| \text{ and } v_i = 0 & \text{(sticking mode)} \\ \text{or} & \|\lambda_{t,i}\| = \mu_i |\lambda_{n,i}| \text{ and } v_{n,i} = 0, v_{t,i} \ne 0 & \\ & \text{and } \exists \alpha_i > 0, \lambda_{t,i} = -\alpha_i v_{t,i} & \text{(sliding mode).} \end{array} \right. $$

In the sliding mode, one equivalently has $\lambda_{t,i} = -\mu_i |\lambda_{n,i}| \frac{v_{t,i}}{|v_{t,i}|}$. Let us recall that $(v_i, \lambda_i) \in \mathcal{C}(\mathbf{n}_i, \mu_i) \Rightarrow \lambda_{t,i} \in \mathcal{D}(\lambda_{n,i}, \mu_i)$, where $\mathcal{D}(\lambda_{n,i}, \mu_i) = |\lambda_{n,i}|\mathcal{D}_{\mu_i}$ is the Coulomb-Moreau's disk [78] [22, Section 5.3.2], and $\mathcal{D}_{\mu_i} = \{z \in \mathbb{R}^d | z^\top z \le \mu_i^2\}$, $d = 1$ (2D friction) or $d = 2$ (3D friction). Then Coulomb's
law (in velocity) in (10) is equivalently formulated as the variational inequality: Find $\lambda_{t,i} \in \mathcal{D}(\lambda_{n,i}, \mu_i)$ such that

$$v_{t,i}^{\top}(y - \lambda_{t,i}) \geq 0 \quad \text{for all } y \in \mathcal{D}(\lambda_{n,i}, \mu_i), \qquad (11)$$

which shows that Coulomb's friction follows a maximal dissipation principle. This variational inequality is in turn equivalent to:

$$\begin{aligned} v_{t,i} \in -\mathcal{N}_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(\lambda_{t,i}) &\Leftrightarrow \lambda_{t,i} = \operatorname{proj}[\mathcal{D}(\lambda_{n,i}, \mu_i); \lambda_{t,i} - \rho v_{t,i}], \quad \forall \rho > 0 \\ &\Leftrightarrow \lambda_{t,i} \in \partial\sigma_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(-v_{t,i}) \end{aligned} \qquad (12)$$

where $\sigma_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(\cdot)$ is the support function of the disc $\mathcal{D}(\lambda_{n,i}, \mu_i)$ [92]. We have $\sigma_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(y) = \mu_i|\lambda_{n,i}| ||y||$, hence $\partial\sigma_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(y) = \mu_i|\lambda_{n,i}| \partial||y||$, so that $\partial\sigma_{\mathcal{D}(\lambda_{n,i}, \mu_i)}(-v_{t,i}) = -\mu_i|\lambda_{n,i}| \partial||v_{t,i}||$, with

$$\partial\|v_{t,i}\| = \begin{cases} \frac{v_{t,i}}{\|v_{t,i}\|} & \text{if } v_{t,i} \neq 0 \\ \mathcal{B}_2 & \text{if } v_{t,i} = 0 \end{cases}, \qquad (13)$$

where $\mathcal{B}_2 = \{z \in \mathbb{R}^d | z^\top z \leq 1\}$. Let us end this brief section on Coulomb friction with the De Saxcé formulation [37, 56]. The objective is to apply a correction to the velocity so that the modified velocity belongs to the dual to the friction cone (hence, yielding a cone complementarity formulation).

To that end let us define $\hat{U}_i \triangleq \begin{pmatrix} v_{n,i} + \mu_i ||v_{t,i}|| \\ v_{t,i} \end{pmatrix}$ at contact point $i$. Then Coulomb's friction with unilateral contact at contact $i$ is rewritten equivalently as (see Fig. 3(a)):

$$0 \leq v_{n,i} \perp \lambda_{n,i} \geq 0 \Leftrightarrow \lambda_{n,i} \in -\mathcal{N}_{\mathbb{R}_+}(v_{n,i})$$

and:

$$C^*(\mathbf{n}_i, \mu_i) \ni \hat{U}_i \perp \lambda_i = \begin{pmatrix} \lambda_{n,i} \\ \lambda_{t,i} \end{pmatrix} \in C(\mathbf{n}_i, \mu_i) \Leftrightarrow \hat{U}_i \in -\mathcal{N}_{C(\mathbf{n}_i, \mu_i)}(\lambda_i)$$

Given $\lambda_i \in C(\mathbf{n}_i, \mu_i)$, find $\hat{U}_i$ such that $\hat{U}_i^\top(\lambda_i - \gamma) \leq 0$ for all $\gamma \in C(\mathbf{n}_i, \mu_i)$. (14)

Notice that $v_{n,i} > 0 \Rightarrow v_{n,i} + \mu_i ||v_{t,i}|| > 0 \Rightarrow \lambda_{n,i} = 0 \Rightarrow \lambda_{t,i} = 0$. Thus active constraints which are to be deactivated, logically do not produce any contact force. The set of conditions in (14) is equivalent to Coulomb's friction [22, section 5.3.3] [37, 56]. Notice that the normal cone inclusion does not mean that Coulomb's friction has become associated, since the velocity is modified. The first equivalence stems from [92, Corollary 23.5.4], the second equivalence stems from the variational definition of the normal cone (keeping in mind that the friction cone is a closed convex set).
<figure><img src="image_3.png" /><figcaption>(a) 3-ball system.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) $\frac{k_2}{k_1} = 1$, $v_1 = 1$, $v_3 = 0$, Hertz contact.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) $\frac{k_2}{k_1} = 0.1$, $v_1 = 1$, $v_3 = -1$, Hertz contact.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d) $\frac{k_2}{k_1} = 100$, $v_1 = 1$, $v_3 = 0$, Hertz contact.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(e) $\frac{k_2}{k_1} = 1$, $v_1 = 1$, $v_3 = 0.3$, Hertz contact.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(f) $\frac{k_2}{k_1} = 0.1$, $v_1 = 1$, $v_3 = -1$, linear stiffness.
(g) $\frac{k_2}{k_1} = 2.24$, $v_1 = 1$, $v_3 = -1$, linear stiffness.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

Figure 2: Contact/impact forces histories, 3 aligned balls [3] (dashed lines; balls 1-2, solid lines: balls 2-3, $v_i$ is ball i pre-impact velocity. See also [81, Chapter 6] for similar results.

### 3.2.2 Comments

* What about line/line, plane/plane, etc, contacts? Since there is always a nonzero area in which friction acts, it may be important to take its effects into account. Friction coefficients may also depend on contact pressure, temperature, surface topology [15].

* Rigid body models as (1) with set-valued friction, class (a) impact models,
bilateral and unilateral constraints, suffer from dependency of constraints ("too many constraints" and "too many contact force multipliers" which induce hyperstaticity issues, or due uniquely to friction [22, section 5.5.6]) which hampers calculation of contact force multipliers in a unique way. It is customary to solve such problems with compliant models, thus introducing other issues. Alternative solutions exist like arbitrary distribution of contact forces among the contact points. Redundancy may imply the choice of a particular solver, see the exhaustive numerical analysis in [2].

* Contact forces calculation: this is very important in some design processes in order to calculate the dimensions of mechanisms' elements. A force/identification compliant model has to be used. For class (a) models, how can this be obtained? One solution is to use a contact force model once the pre- and post-impact velocities are known, see, e.g., [29].

* During persistent contact phases of motion, unwanted oscillations due to contact stiffness, usually are suppressed or at least decreased by drastic increase of the viscous friction coefficient (to reach a superdamped regime). But the question that arises with such a procedure is: when a damping coefficient is tuned this way, is it a mechanical (physical) parameter or a numerical parameter (like the time step)?

* Several classical spring-dashpot models (linear, Kuwabara-Kono) outcome contact/impact force with wrong sign during the impact expansion phase. This has no clear physical justifications, as adhesive effects can hardly be modeled this way.

* Coulomb's friction is often modeled at the impulse/velocity level for impacts of class (a), i.e., all forces $\lambda_{n,i}$, $\lambda_{t,i}$ in (10)-(14) are replaced by impulses $\bar{\lambda}_{n,i}$, $\bar{\lambda}_{t,i}$. This cannot model some complex tangential dynamics which may occur during the impact process, e.g., when $v_{t,i}$ reverses its sign during the collision at some impact point. But, adding dependency on both pre and post-impact velocities significantly improves predictability capabilities. This has long been known for single impacts [22, sections 4.3.5, 4.3.6] [100]. For multiple impacts, repeated compression/expansion phases as in Fig. 2 add difficulty.

* Three-dimensional Coulomb's friction is often facetized, i.e., the friction cone is approximated by a facetized cone (a pyramidal convex polyhedral cone). This has advantages (this can be expressed as LCPs which are convenient to solve numerically), but it has serious drawbacks as well: how many facets should be chosen? Should the designer choose an inside or an
outside approximation? How much does it modify the system's motion [4, section 13.3.7]? And, an LCP numerical solver like Lemke's algorithm, has a complexity which increases much faster than linearly than the problem's size [10].

* A fundamental property of Coulomb's friction is that it is *nonassociated*. This means that it is not possible to express it as $\begin{pmatrix} \lambda_{n,i} \\ \lambda_{t,i} \end{pmatrix} \in \partial \varphi_i \begin{pmatrix} v_{n,i} \\ v_{t,i} \end{pmatrix}$ for some proper, convex, lower semicontinuous function $(v_{n,i}, v_{t,i}) \mapsto \varphi_i(v_{n,i}, v_{t,i})$. A proof is reported in [22, section 5.3.3]. Sometimes the friction model is modified so that it becomes associated [9], see Fig. 3(b). This is called a convexification. Then the numerical treatment is greatly simplified. This convexification implies that when in sliding motion on an activated constraint $i$, $\lambda_{n,i} > 0 \Rightarrow v_{n,i} > 0$, which implies deactivation of this constraint. Moreover the resulting frictional model yields zero dissipation since the reaction contact force is always orthogonal to the contact (local) velocity. In spite of these drawbacks, this is acceptable in some applications.

<figure><img src="image_10.png" /><figcaption>Figure 3: Set-valued friction models, $C^O$ and $C^*$ are polar and dual cones.</figcaption></figure>

* Dynamic models of friction (Lugre, Bliman-Sorine, Leuven, etc) model correctly some effects like pre-sliding, but also present severe drawbacks: poor modeling of sticking, many parameters per contact (estimation not always easy, physical meaning sometimes unclear), poor predictability capacities in some application cases [91, 97]. They also often correspond to some kind of regularization around zero tangential velocity, which can yield stiff equa-
tions to be integrated with specific schemes (stiff integrators), implying very
small time steps thus large simulation times [97].

• As alluded to above, it is easy to extend set-valued Coulomb’s model to varying CoF (Stribeck effects or else), or different static and dynamic CoF.

4 Multiple Impacts Phenomena

Multiple impacts are very complex mechanical phenomena. They incorpo-
rate the modeling difficulties of single impacts, and specific features not
present in single impacts are added. This is illustrated on aligned chains of
balls in [81] and in [22, section 6.1], where analyses and many references to
the field can be found.

**4.1 Definition**

Roughly speaking, a multiple impact occurs in a multibody system each time
several impacts occur at the same time (thus, it may also be called a simul-
taneous impact). The following definition is made in a Lagrange dynamical
systems framework, i.e.: $C_{conf} \subseteq \mathbb{R}^n$ denotes the systems's configuration
space, $n$ is the number of degrees of freedom, $q \in C_{conf}$ is the vector of in-
dependent generalized coordinates, $h_{nu,i}: C_{conf} \rightarrow \mathbb{R}_+$, $1 \le i \le m$, are gap
functions, or signed distances, each of which represents the distance between
two points of the system (what is meant by "system" here may include a fixed
obstacle), each $h_{nu,i}(\cdot)$ is at least continuously differentiable. Therefore, the
admissible domain is defined as $\Phi = \{q \in C_{conf} | h_{nu,i}(q) \ge 0, 1 \le i \le m\}$.
It is assumed that $\Phi \ne \emptyset$, and that $q(0) \in \Phi$: the system is initialized inside
the admissible domain. A single impact occurs at time $t$ if $h_{nu,i}(q(t)) > 0$
for all $i \ne j$, $h_{nu,j}(q(t)) = 0$ and $v_{n,j}(t^-) = \nabla h_{nu,j}(q(t))^{\top} \dot{q}(t^-) < 0$ for some
$j \in \{1, \dots, m\}$ (this is a transversality condition). At this stage an addi-
tional assumption has to be made on the vector of normal direction to the
constraint boundary: $\nabla h_{nu,j}(q(t)) \ne 0$ when $h_{nu,j}(q(t)) = 0$ and $q(t) \in \Phi$.
This will always be assumed when one or several gap functions satisfy the
above conditions.

**Definition 1** Let the above setting hold true and $1 \le p \le m$. A p-impact
occurs when a codimension p singularity of bd($\Phi$) is attained by a trajectory.

In other words, a *p*-impact occurs when the intersection of *p* constraints
surfaces is reached. This is illustrated in Fig. 4(a) with *p* = 2. Clearly when
two or more gap functions vanish at $t$ (they are said to become active) with
transversality, a multiple impact occurs and the normal direction to $bd(\Phi)$
at the impact point is not uniquely defined. Then the notion of normal
cone, which is a basic tool of Convex and Nonsmooth Analysis, becomes
quite useful. Some multiple impact models rely on normal cones. Defini-
tion 1 encompasses cases when only some gap functions become active with
transversality, while others were active before the impact and are deactivated
after the impact: this is the case for the classical examples of ball-chains,
or rocking blocks as in Fig. 4(b). The latter are true multiple impacts, be-
cause even previously closed contacts can accumulate impulse and produce
body detachment (and velocity jump for class (a) models). Also, what hap-
pens in a neighborhood of the singularity is often worth investigating and a
nontrivial issue.

<figure><img src="image_11.png" /><figcaption>Figure 4: Multiple impacts.</figcaption></figure>

## 4.2 Typical Application Examples

Multiple impacts are encountered in granular matter (Fig. 4(b)): chains of aligned balls (monodisperse, polydisperse, decorated, dimer, disordered, tapered and antitapered chains are studied for their ability to amplify, or attenuate force transmission, or to capture energy), 2D and 3D granular materials. Impact dampers (also known as impact oscillators, or vibro-impact nonlinear energy sinks) use chains of beads [54]. Rocking block systems (2D or 3D) with an excited base for monument or building stability analysis when submitted to earthquakes, or masonry structures which may be seen as blocks assemblies [4, section 10.2.4], or falling rock trajectory prediction [16, 66] (see Fig. 4(b)): they involve line/line or plane/plane contact. Joint
clearance in multibody systems often involves surface/surface contacts that may yield multiple impacts (see [7, Fig. 5–10]). However such conformal contacts can also be modeled with non-Hertz elasticity [113] and Winkler’s elastic foundation [109]. In Robotics, foot/ground contact in biped robots locomotion is another instance, as well as object manipulation, or snake robots locomotion: all these tasks usually involve several contact/impact points. More recent applications can be found in clothes, textiles, hair, see, e.g., [30, 36, 68]. These systems usually involve large numbers of contacts.

## 4.3 Some Specific Features of Multiple Impacts

Multiple impacts possess specific features which are not met in single impacts [22, Chapter 6]:

* Discontinuities of trajectories $(q, \dot{q})$ with respect to initial data (first noticed in 1952 in [44]), while systems with single impacts usually enjoy continuity [96]. The main source of these discontinuities are the kinetic angles. As a consequence, different initial positions/velocities yield different sequences of binary collisions with different outcomes (see Fig. 5(a)). Kinetic angles are defined as:

$$ \theta_{ij}(q) = \pi - \arccos \frac{\nabla h_{nu,i}(q)^{\top} M^{-1}(q) \nabla h_{nu,j}(q)}{\|\nabla h_{nu,i}(q)\|_{M^{-1}(q)} \|\nabla h_{nu,j}(q)\|_{M^{-1}(q)}}, \quad (15) $$

with $\|x\|_{M^{-1}(q)} = \sqrt{x^{\top} M^{-1}(q) x}$. Kinetic angles are quantities that reflect the couplings between the inertial and the geometrical properties of the system with unilateral constraints. Let $U_n(q, \dot{q}) = (v_{n1}, \dots, v_{n,m})^{\top}$, then from (5): $U_n(t^+) - U_n(t^-) = A_{nu}(q)\bar{\lambda}_{nu}(t) + A_{nutu}(q)\bar{\lambda}_{tu}$, where $A_{nu}(q)$ is the Delassus’ matrix, whose entries reflect the kinetic angles. In particular it is diagonal if and only if $\theta_{ij}(q) = 0$ for all $i \neq j$, i.e., all constraints are mutually orthogonal and the impacts are decoupled. Same for $A_{nutu}(q)$ for the couplings between normal and tangential directions.

* Several outcomes for constant energetic behaviour (see Fig. 6(a)): proves the need for other parameters not present in single impacts (stiffness ratios, impulse ratios, off-diagonal terms in Frémond matrices). Energy dispersion is needed in addition to dissipation (see section 4.4).

* Nonlinear wave propagation [25, 118] is essential to model energy dispersion in some multiple impacts like chains of aligned balls (not to be confused with bulk vibrations). Fig. 5(b) depicts soliton wave (force pulse at the end of the 7-big-bead monodisperse part) and soliton wave trains (force pulses
<figure><img src="image_12.png" /><figcaption>(a) Particle hitting a wedge.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Soliton and soliton wave trains.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Stepped chain.</figcaption></figure>

Figure 5: Discontinuity w.r.t. initial data; wave effects in stepped chains
[81, Fig. 5.19].

between the end of the chain and the wall, with 25 and 50 small beads,
respectively) in a stepped chain as in Fig. 5(c). Elasticity coefficient, mass
ratios and local contact geometry play crucial roles. See [81, section 5.4] for
few particular features of aligned chains. Geometry of the collision (kinetic
angles between constraints) is essential in other kinds of multiple impacts
like rocking blocks.

* Unexpected behaviours may occur: ground-chain contact force in vertically falling aligned chains of balls is almost independent of the number of balls [39, 116] (nonlinear waves travelling through the chain are the reason why). When a 2-ball chain made of a big ball (basket ball) with a small ball (tennis ball) at its top, hits the rigid ground, the big ball almost stops while the small one rebounds with velocity much higher than the pre-impact velocity of the 2-ball chain.

* Existence of accumulation of impacts with $e_n = 0$, e.g., particle that hits a V-shaped obstacle with angle $> \frac{\pi}{2}$ and Moreau's impact law, with gravity; or rocking blocks which accumulate rocking oscillations until rest is attained in finite time [49].

* Spanning the whole post-impact admissible velocities space (consistency (iv) (a)), may require CoRs larger than 1, or negative CoRs, at some contact points [81, section 3.1.2], see section 4.4 and Region (B) in Fig. 6(b).

* Possible existence of repeated impacts during a single multiple impact. In
other words, before all local velocities become admissible and the impact terminates, several compression/expansion phases may occur at several contact points (see Fig. 2(b)–2(g) for an illustration).

• Dependence of parameters on initial data, mass ratios (e.g., in granular matter), kinetic angles (15), due to energy dispersion, distance effects (e.g., changing the dissipation behaviour at contact i may influence all other contacts $i \neq j$).

• The behaviour in the neighborhood of the p-impact singularity can be complex [105, 83] [81, section 3.4, Appendix A].

• Redundancy of constraints (or, hyperstaticity) may induce modeling and numerical difficulties.

• The termination of the p-impact process may not be trivial to compute or guarantee for some models.

• Activated contacts $h_{nu,i}(q(t)) = 0$, $v_{n,i}(t^{-}) = 0$, can be deactivated with $v_{n,i}(t^{+}) > 0$, which is impossible to obtain in single impacts with $v_{n,i}(t^{+}) = \mathcal{F}(\text{CoR}, v_{n,i}(t^{-}))$, where $\mathcal{F}(\text{CoR}, 0) = 0$ (even impacts-without-collision in Painlevé paradoxes do not produce detachment: they convert sliding into sticking [45]). Only distance effects between contacts i and j, $i \neq j$, allow this.

## 4.4 Various Modelling Approaches

Multiple impact models may be classified as follows (the three classes (a), (b), (c) refer to section 2):

**Group 1:** models of class (a) which rely on generalized equations, Convex Analysis (solving variational inequalities), Complementarity Theory (solving LCPs) and Optimization (solving QPs). We may cite:

• Moreau's law [79], which is an extension of Newton's impact law, with a global CoR $e_n$ (same for all contacts). In a generalized frictionless framework it reads as:

$$ \dot{q}(t^{+}) = -e_n \dot{q}(t^{-}) + (1 + e_n) \operatorname{argmin}_{z \in V(q(t))} \frac{1}{2}(z - \dot{q}(t^{-}))^{\top} M(q(t))(z - \dot{q}(t^{-})), \quad (16) $$

where $V(q(t)) = \{z \in \mathbb{R}^n | \nabla h_{nu,i}(q(t))^{\top} z \ge 0, i \in \mathcal{I}(q(t))\}$, $\mathcal{I}(q(t)) = \{i \in \{1, \dots, m\} | h_{nu,i}(q(t)) = 0\}$, is the tangent cone linearization cone to $\Phi$ at $q(t)$. When $e_n \in [0, 1]$, the impact law (16) satisfies the basic consistencies
(i)-(iii) in section 3 as well as (iv) (b), but largely fails with (iv) (a). Indeed one drawback is the lack of energy dispersion modeling, since it maximizes dispersion [81, section 3.2]. It is necessary to express (16) in a more tractable way. Let $U_n(q, \dot{q}) \triangleq \nabla h_{nu}(q)^{\top} \dot{q} = (v_{n,1}, \dots, v_{n,m})^{\top} \in \mathbb{R}^m$, and $U_t(q, \dot{q}) \triangleq H_t(q)^{\top} \dot{q} = (v_{t,1}^{\top}, \dots, v_{t,m}^{\top})^{\top} \in \mathbb{R}^{md}$. Then (16) implies the Linear Complementarity Problem (LCP):

$$ 
\begin{aligned}
0 \le (I_m + \mathcal{E}_{nn})U_n(t^-) + A_{nu}(q(t))\bar{\lambda}_{nu}(t) \perp \bar{\lambda}_{nu}(t) \ge 0 \\
\bar{\lambda}_{nu}(t) = \underset{z \ge 0}{\operatorname{argmin}} \frac{1}{2} z^{\top} A_{nu}(q(t))z + z^{\top} (I_m + \mathcal{E}_{nn})U_n(t^-) \quad (17) \\
\bar{\lambda}_{nu}(t) = \operatorname{proj}_{A_{nu}(q(t))} [\mathbb{R}_+^m; -A_{nu}(q(t))^{-1}(I_m + \mathcal{E}_{nn})U_n(t^-)],
\end{aligned}
$$

with $A_{nu}(q) = \nabla h_{nu}(q)^{\top} M^{-1}(q) \nabla h_{nu}(q) \succeq 0$ the Delassus' matrix, and:

$$ 
\begin{aligned}
U_n(t^+) - U_n(t^-) &= A_{nu}(q(t))\bar{\lambda}_{nu}(t) \\
U_t(t^+) - U_t(t^-) &= A_{tunu}(q(t))\bar{\lambda}_{nu}(t),
\end{aligned}
\quad (18)
$$

with $A_{tunu}(q) \triangleq H_{tu}(q)^{\top} M^{-1}(q) \nabla h_{nu}(q)$. Notice that $\bar{\lambda}_{nu}(t)$ is as in (5). To illustrate let us consider the three-ball chain in Fig. 2(a), with initial velocity $\dot{q}(t^-) = (V_s, 0, 0)^{\top}$. The complete post-impact velocity space is depicted in Fig. 6(a). Moreau's line (O,B) only is accessible using (17) (18). In Fig. 6 $KER = \frac{T(t^+)}{T(t^-)}$ is the dissipation index, and $C_{KE} = \frac{1}{T(t^+)} \sqrt{\frac{1}{n} \sum_{i=1}^n (T_i(t^+) - \bar{T}(t^+))^2}$, $T_i(t^+) = \frac{1}{2} m_i \dot{q}_i(t^+)$, $\bar{T}(t^+) = \frac{1}{n} \sum_{i=1}^n T_i(t^+)$, $n=3$ in our case study, is the dispersion index.

When friction is present, one basically uses the formulations in section 3.2, replacing the forces $\lambda_{n,i}$ and $\lambda_{t,i}$ by their impulses $\bar{\lambda}_{n,i}$ and $\bar{\lambda}_{t,i}$. This gives rise to a problem of the form:

$$ 
\begin{aligned}
\underbrace{\begin{pmatrix} U_n(t^+) - U_n(t^-) \\ U_t(t^+) - U_t(t^-) \end{pmatrix}}_{=U(t^+) - U(t^-)} &= \underbrace{\begin{pmatrix} A_{nu}(q) & A_{nutu}(q) \\ A_{tunu}(q) & A_{tu}(q) \end{pmatrix}}_{\triangleq A_{ntu}(q)} \begin{pmatrix} \bar{\lambda}_{nu}(t) \\ \bar{\lambda}_{tu}(t) \end{pmatrix} \quad (19) \\
\bar{\lambda}_{nu}(t) &= \bar{\lambda}_{tu}(t)
\end{aligned}
$$

$$ C^{*}(\mathbf{n}_{i}, \mu_{i}) \ni \hat{U}_{e,i}(t) \perp \bar{\lambda}_{u,i}(t) = \begin{pmatrix} \bar{\lambda}_{nu,i}(t) \\ \bar{\lambda}_{tu,i}(t) \end{pmatrix} \in C(\mathbf{n}_{i}, \mu_{i}), \quad i \in \mathcal{I}(q(t)) \quad (20) $$

$$ C(\mathbf{n}_i, \mu_i) = \{ \bar{\lambda}_i(t) \in \mathbb{R}^{d+1} \mid \| \bar{\lambda}_{tu,i}(t) \| \le \mu_i \bar{\lambda}_{nu,i}(t) \}, \quad (21) $$

where (14) is used in (20). In (19) we have $A_{nutu}(q) = A_{tunu}(q)^{\top}$, $A_{tu}(q) = H_{tu}(q)^{\top} M^{-1}(q) H_{tu}(q)$, $C^{*}(\mathbf{n}_i, \mu_i)$ is the dual cone to $C(\mathbf{n}_i, \mu_i)$. The modified velocity $\hat{U}_{e,i}(t)$ is from (14): $\hat{U}_{e,i} \triangleq \begin{pmatrix} v_{n,i}(t^+) + e_n v_{n,i}(t^-) + \mu_i \|v_{t,i}\| \\ v_{t,i} \end{pmatrix}$.
<figure><img src="image_15.png" /><figcaption>Figure 6: Outcome domains, $\gamma_i = \dot{q}_{i+1} - \dot{q}_i$, $i = 1, 2$, monodisperse 3-ball chain [81].</figcaption></figure>

A drawback is the lack of energetic consistency (i) of (19)–(21) in generic situations, where (9) may not hold, depending on $A_{ntu}(q)$, $e_n$ and $\mu_i$. Analyses of Moreau’s impact law are proposed in [52, 22], conditions under which it provides accurate enough results are studied in [81, 82]. Moreau’s framework allows for extensions with better impact laws, as shown next.

• Frémond [41, 42, 43, 26]: constrain the impulses $\bar{\lambda}_{nu}(t)$ and $\bar{\lambda}_{tu}(t)$ inside the subdifferential of a *pseudo-potential* function [77], whose argument depends on both $v_{n,i}(t^+)$ and $v_{n,i}(t^-)$, and on an “energy dispersion” matrix taking into account distance effects between impact points. This allows us to satisfy energetic consistency (i), using Thomson and Tait’s inequal-
ity (9). Notice that Moreau's law (16) or (17) (18) stems from $\bar{\lambda}_{nu}(t) \in$
$-\mathcal{N}_{V(q(t))} \left( \frac{\dot{q}(t^+)+e_n\dot{q}(t^-)}{1+e_n} \right)$ and (5). While (20) is equivalent to $\bar{\lambda}_{u,i}(t) \in$
$-\mathcal{N}_{C^*(n_i, \mu_i)}(\hat{U}_{e,i}(t))$. Frémond's idea is to extend this as

$$
\bar{\lambda}_u(t) \in -\partial\Gamma(\dot{q}(t^+), \dot{q}(t^-)),
$$

where $\Gamma(\cdot)$ is parameterized by CoRs and CoFs: $\mu_1, \dots, \mu_m, e_{n,1}, \dots, e_{n,m}$
(possibly tangential CoRs also). Coulomb's friction may be extended as

$$
\bar{\lambda}_{tu,i}(t) \in -\mu_i |\bar{\lambda}_{nu,i}(t)| \partial \Gamma_{t,i} (v_{t,i}(t^+) + v_{t,i}(t^-)),
$$

and

$$
\bar{\lambda}_{nu,i}(t) \in -\partial\Gamma_{n,i}(e_{n,i}, v_{n,i}(t^+), v_{n,i}(t^-)),
$$

using suitable convex proper lower semicontinuous functions $\Gamma(\cdot)$, $\Gamma_{t,i}(\cdot)$ and $\Gamma_{n,i}(\cdot)$. This yields $m$ terms like $-\mu_i|\bar{\lambda}_{nu,i}|(v_{t,i}(t^+) + v_{t,i}(t^-))\partial\Gamma_t(v_{t,i}(t^+) + v_{t,i}(t^-))$ in (9), which are guaranteed to be nonpositive if the graph of $\partial\Gamma_{t,i}(\cdot)$ contains $(0,0)$, since this set-valued mapping is maximal monotone [92]. Inserting both constraints in (5) yields a generalized equation with unknown $\dot{q}(t^+)$:

$$
\begin{align}
M(q)(\dot{q}(t^+) - \dot{q}(t^-)) &= (\nabla h_{nu}(q) H_{tu}(q)) \begin{pmatrix} \bar{\lambda}_{nu}(t) \\ \bar{\lambda}_{tu}(t) \end{pmatrix} \nonumber \\
&\in -(\nabla h_{nu}(q) H_{tu}(q)) \partial \Gamma(\dot{q}(t^+), \dot{q}(t^-), \mu_i, e_{n,i}) \tag{22}
\end{align}
$$

Consider Moreau's inclusion, for a planar particle hitting the ground: $\bar{\lambda}_{nu}(t) \in$
$-\mathcal{N}_{\mathbb{R}_+} \left( \frac{\dot{q}(t^+) + e_n \dot{q}(t^-)}{1+e_n} \right)$ (a particular choice for $\Gamma_n(\cdot) = \Psi_{\mathbb{R}_+}(\cdot)$, the indi-
cator function of $\mathbb{R}_+$ [92]), $m(\dot{q}(t^+) - \dot{q}(t^-)) = \bar{\lambda}_{nu}(t)$ (this is (5)), and
$\bar{\lambda}_{tu}(t) \in -\mu |\bar{\lambda}_{nu}(t)| \operatorname{sgn}(\dot{q}(t^+) + \dot{q}(t^-))$ (a particular choice for $\Gamma_t(\cdot) = |\cdot|$),
yields a well-posed impact. The strength of Frémond's approach is the free-
dom in choosing pseudo-potentials $\Gamma_n(\cdot)$ and $\Gamma_t(\cdot)$, its weakness may be too
much freedom and the difficulty in properly choosing $\Gamma_n(\cdot)$ and $\Gamma_t(\cdot)$. Many
examples are provided in [42, 43, 26]. Examples are in [22, Propositions 4.3,
4.7] and a general framework using De Saxcé's approach (14) is proposed
in [5], extending material in [22, section 4.3.1.2]. Consistency (iv) (a) may
fail for some pseudo-potentials, but adding parameters solves this issue [81,
section 3.3].

• Pfeiffer and Glocker [87]: extend the kinetic Poisson's approach, using a
sequence of two LCPs for the compression and expansion phases *per contact*,
respectively, with the same maximum compression times at all impact points
(see [81, section 3.1.2] for an explanation), and kinetic CoR $e_{p,i}$ at contact $i$. In the frictionless case this gives rise to

$$ \bar{\lambda}_{nu}(t) = \mathcal{E}_{nn}\bar{\lambda}_{nu}^{comp} + \mathrm{proj}_{A_{nu}(q)}[\mathbb{R}_{+}^{m}; -A_{nu}(q)\mathcal{E}_{nn}\bar{\lambda}_{nu}^{comp}], $$

with $\bar{\lambda}_{nu}^{comp} = \mathrm{proj}_{A_{nu}(q)}[\mathbb{R}_{+}^{m}; U_{n}(t^{-})]$ the impulse at the end of the compression phase, $\mathcal{E}_{nn}$ is a diagonal kinetic restitution matrix [50]. This can be seen as an extension of (17). However multiple impacts can involve repeated compression and expansion phases with different maximum compression times, see Fig. 2, not taken into account in this approach. The event-capturing time-stepping scheme in [11], inspired from the Moreau-Jean algorithm uses this impact law. The Pfeiffer-Glocker model allows us to extend Moreau's line, however region (B) in Fig. 6(b) is accessible only if $1 \le e_{p,i} \le 2$, $i = 1, 2$. Such parameters cannot be estimated from off-line pairwise collisions (hence consistency (iv) (b) in section 3 fails). Other pitfalls are pointed out in [60, section 10].

• **Other contributions.** [67] use Moreau's and Frémond's ideas by constraining impulses as $\lambda_{nu,i} \in -\partial\sigma_{\mathbb{R}_-}(v_{n,i}(t^+) + e_{n,i}v_{n,i}(t^-))$, and $\lambda_{tu,i} \in -\partial\sigma_{\mathcal{D}(\bar{\lambda}_{nu,i}, \mu_i)}(v_{t,i}(t^+) + e_{ti}v_{t,i}(t^-))$, $e_{n,i} \in [0, 1]$, $e_{t,i} \in [-1, 1]$ (these are particular choices for $\Gamma_n(\cdot)$ and $\Gamma_t(\cdot)$). A careful analysis of consistencies (especially (i), crucial for stability analysis) is made using the Thomson and Tait formula (9). In [111] an impact law for chains of balls is proposed, which uses the splitting of the pre-impact velocity space into subspaces (cones) defined from the post-impact patterns. With each cone a restitution matrix is associated which imposes an energy dispersion pattern. The three consistencies (i)–(iii) are verified, relying on maximal monotonicity and non-expansiveness. The impact law in [89] minimizes the cost function $\frac{1}{2}(\bar{\lambda}_{nu}^{\top}, \bar{\lambda}_{tu}^{\top})A_{ntu}(q) (\bar{\lambda}_{nu}^{\top}, \bar{\lambda}_{tu}^{\top}) + (\bar{\lambda}_{nu}^{\top}, \bar{\lambda}_{tu}^{\top}) (v_{n}(t^{-}), v_{t}(t^{-}))$, under restitution constraints $\sum_{i=1}^{m} v_{nu,i}(t^{+})m_{i}^{\alpha} \ge -\sum_{i=1}^{m} e_{n,i}v_{nu,i}(t^{-})m_{i}^{\alpha}$, $\alpha$ a parameter which rules energy dispersion and is to be estimated, $m_i$ are local normal inertias. [32, 95] use pseudo-potentials as in [67]. A time-stepping nonsmooth generalized-$\alpha$ scheme is proposed which filters out previously active contacts (at which it is assumed null impulse). It usually needs several iterations of pairwise collisions. Thus it shares similarities with the algorithmic methods of Group 2 below, but the iterations are ruled automatically.

Moreau's impact framework, though limited in its basic formulation, furnishes a convenient and powerful framework for extensions.

It is implicitly assumed above that $A_{nu}(q(t)) \succ 0 \Rightarrow m \le n$. In general
$A_{nu}(q(t)) \succ 0$ only because of constraints redundancy. Then (17) may have several solutions, but due to symmetry of $A_{nu}(q)$, $U_n(t^+)$ is unique (using [33, Theorems 3.8.6, 3.1.7]).

**Group 2:** models of class (a) not relying on the above mathematical tools:
* Kinetic quasi-velocities and restitution matrices [20, 23, 21]: the Lagrange dynamics is transformed via a specific state-space change, which is not a generalized coordinate transformation. This allows us to get a canonical dynamics which splits normal and tangential (to be taken in a generalized sense) dynamics in $C_{conf}$ at the singularity where the m-impact occurs. Then (5) becomes:

$$ \dot{q}_{\text{norm}}(t^+) - \dot{q}_{\text{norm}}(t^-) = \mathbf{n}_q^\top M(q) \mathbf{n}_q \bar{\lambda}_{\text{nu}}(t) + \mathbf{n}_q^\top H_t(q) \bar{\lambda}_{\text{tu}}(t) \quad (23) $$

$$ \dot{q}_{\text{tan}}(t^+) - \dot{q}_{\text{tan}}(t^-) = \mathbf{t}_q^\top H_t(q) \bar{\lambda}_{\text{tu}}(t), \quad (24) $$

where $\mathbf{n}_{q,i} = \frac{M^{-1}(q)\nabla h_{nu,i}(q)}{\sqrt{\nabla h_{nu,i}^\top(q)M^{-1}(q)\nabla h_{nu,i}(q)}}$, $1 \le i \le m$, $(\mathbf{n}_q, \mathbf{t}_q)$ is an orthonormal basis of $\mathbb{R}^n$, $(\dot{q}_{\text{norm}}, \dot{q}_{\text{tan}}) = (\mathbf{n}_q^\top \quad \mathbf{t}_q^\top) M(q) \dot{q}$. One of the interests of these quasi-Lagrange dynamics is that the “tangential” part (24) is decoupled from the “normal” part (23). In a chain of aligned balls, $\dot{q}_{\text{tan}}$ is the center of mass velocity, in a rocking block this is the horizontal velocity. The framework in (23) (24) is suitable to formulate a kinematic impact law as:

$$ \left( \begin{array}{c} \dot{q}_{\text{norm}}(t_k^+) \\ \dot{q}_{\text{tan}}(t_k^+) \end{array} \right) = -\mathcal{E} \left( \begin{array}{c} \dot{q}_{\text{norm}}(t_k^-) \\ \dot{q}_{\text{tan}}(t_k^-) \end{array} \right) \quad (25) $$

where $\mathcal{E} = (\mathcal{E}_{nn} \quad \mathcal{E}_{nt})$ is a generalized $n \times n$ restitution matrix. Off-diagonal terms in $\mathcal{E}_{nn}$ play a particular role (see [23]), they are necessary for consistency (iv) (a). Admissible CoRs are shown in Fig. 7(b) with associated post-impact patterns of 3-ball monodisperse chains. The Thomson and Tait formula (9) extends as:

$$ T_L(t_k) = \frac{1}{2}(\dot{q}_{\text{norm}}(t_k^+) + \dot{q}_{\text{norm}}(t_k^-))^T [\bar{p}_n(t_k) + A_n(q)\mathbf{n}_q^T H_t(q) p_t(t_k)] \\ + \frac{1}{2}(\dot{q}_{\text{tan}}(t_k^+) + \dot{q}_{\text{tan}}(t_k^-))^T A_t(q) \mathbf{t}_q^T H_t(q) p_t(t_k), \quad (26) $$

where $A_n(q) = (\mathbf{n}_q^\top M(q) \mathbf{n}_q)^{-1}$, $A_t(q) = (\mathbf{t}_q^\top M(q) \mathbf{t}_q)^{-1}$. The three basic consistencies (i)-(iii) (see section 3) can be analysed, they result in constraints on the entries of $\mathcal{E}$ (very detailed analyses for planar rocking blocks and a 4-ball chain are made in [23]). Use of the kinetic metric for multiple
impacts was advocated in [58], see also [41, 52, 50] for a similar idea with normal restitution matrices (named Frémond matrices in [52]). As shown in [22, Remark 6.8], this approach can be used also in the framework of Group 1, using an impact law as in (17). Consistency (iv) (b) usually fails.

• Algorithmic methods: binary collisions where single impacts are treated sequentially until all post-impact velocities are kinematically admissible for consistency (ii) (but, in general, impacts are not sequential, see Fig. 2), and the Han-Gilmore algorithm (an enumerative algorithm which tracks impacting points through the chain of active contact points in a multibody system). These algorithms are not guaranteed to terminate (consistency (vi) may fail). Binary collision sequences are not unique in general (essentially as a consequence of discontinuity w.r.t. initial data, see section 4.3, so that outcomes depend on the chosen sequence of impacts), and may yield nonunique impact outcome (see an example in [82]): consistency (v) in section 3 fails. Consistency (iv) (a) fails also, see Fig. 6(c).

<figure><img src="image_16.png" /><figcaption>Figure 7: Outcome domains (a) and (b)), restitution matrix CoRs (c), monodisperse 3-ball chain [81].</figcaption></figure>

**Group 3:** models of class (b) (Routh, Darboux-Keller's dynamics).

• LZB model [70, 71, 72, 117]: this is an extension of the Darboux-Keller approach for single impacts [22, section 4.3.5]. It relies on (7), an essential ingredient being the *distributing rule* for infinitesimal impulses ratios:

$$d\gamma_{ij} \triangleq \frac{d\bar{\lambda}_{\text{nu},i}}{d\bar{\lambda}_{\text{nu},j}} = \frac{(1 + \eta_i)^{\frac{\eta_i}{\eta_i+1}} k_i^{\frac{1}{1+\eta_i}} E_i^{\frac{\eta_i}{1+\eta_i}}}{(1 + \eta_j)^{\frac{\eta_j}{\eta_j+1}} k_j^{\frac{1}{1+\eta_j}} E_j^{\frac{\eta_j}{1+\eta_j}}}, \quad (27)$$
where $\eta_i$ is the elasticity coefficient, $k_i$ is the equivalent stiffness, $E_i$ is the elastic potential energy (which can be calculated from the velocities), at contact $i$. Equations (7) and (27) show that the multiple impact is mainly ruled by stiffnesses and elasticity coefficients ratios (a fact confirmed by the analysis of class (c) models [22, section 6.1.3] [81, Appendix C]). One energetic CoR $e_{ni} \in [0, 1]$ is associated with each impact point $i$, and a bi- or tri-stiffness model [81] (see Fig. 1) is used to realize the CoR with plasticity-like effects, elasticity coefficients are parameters: consistency (iv) (b) is satisfied. The primary impulse $d\bar{\lambda}_{nu,j}$ used in (7) has to be properly chosen, and can change during the impact, due to repeated compression/expansion cycles, which are allowed at the impact points. When friction is present, using (12) we get: $\frac{d\bar{\lambda}_{tu,i}}{d\bar{\lambda}_{nu,j}} \in \partial\sigma_{D(|d\gamma_{ij}|,\mu_i)}(-v_{t,i})$ (notice that since $q$ is constant, $dv_{t,i} = \nabla h_{nu,i}(q)^{\top}d\dot{q}$), with $|d\gamma_{ij}| \stackrel{\Delta}{=} \frac{|d\bar{\lambda}_{nu,i}|}{d\bar{\lambda}_{nu,j}}$. A detailed analysis of the LZB model is proposed in [81, Chapter 5], the outcomes for the 3-ball chain when varying CoRs, $\eta_i$ and $\frac{k_1}{k_2}$ are in Fig. 6(d) and 7(a): consistency (iv) (a) holds, consistency (vi) is not proved. Numerical pseudo-codes are provided in [81] for the frictionless case, see [117] for the case with friction.

• Hurmuzlu et al [27, 101]: use Routh's incremental dynamics (6) integrated over $[t_i, t_f]$, energetic CoRs and the ratios of impact forces impulses $\gamma_{ij} \stackrel{\Delta}{=} \frac{\bar{\lambda}_{nu,i}}{\bar{\lambda}_{nu,j}}$, named impulse correlation ratios (ICR), as constant parameters to take into account distance effects between impact points. The drawback is that ICRs are not constant in general when the system varies (e.g., increasing or decreasing the number of balls in a chain) [27, Table 2] [101, Fig. 8,9,10] [3]. Consistency (iv) (b) is analysed experimentally [27, Table 2]. ICRs can be interpreted as parameters monitoring energy dispersion [81, Fig. 3.10] [27, Fig. 6].

• **Other contributions.** An approach similar to LZB is proposed in [60]. [55] propose to encompass the whole process (including persistent contact and open contacts) in a differential inclusion which captures all pairwise collisions outcomes computed with Routh's method. This is integrated with a unified time-scale similar to the time-freezing approach in [85, 84] (restricted to plastic impacts with $v_{n,i}(t^+) = \nabla h_{nu,i}(q)^{\top}\dot{q}(t^+) = 0$) in which velocities are continuous, and time $t$ is frozen during collisions. The whole process is formulated in an LCP formalism. This allows an improved numerical simulation of the whole process, allowing for accumulations of impact times. Consistency (i) is proved, as well as (vi) (impact termination, which is a significant achievement compared to other Group 3 methods). [28] use Routh
approach for the impact description of 3D rocking blocks, with Coulomb's
friction at the infinitesimal impulse level, and a global energetic CoR for im-
pact termination (using a virtual end-of-compression time) and consistency
(i). Several steps (in a way similar to the algorithmic Group 2 method)
may be needed until consistency (ii) holds. Indeterminate configurations in
plane/plane impact (i.e., more contact forces than degrees of freedom, a sort
of hyperstatic situation) are treated. [57, 74] use an enumerative process to
calculate the impact with friction outcome, with different termination cri-
teria and a test of consistencies to select solutions.

**Group 4:** models of class (c). These consist merely of placing linear, non-linear or piecewise continuous unilateral spring-dashpot systems at all *m* potential contact/impact points, and integrating (1) as a usual ODE. This can in some cases be a reasonable option (as in chains of balls where contacts are deactivated after the impact, see, e.g., [118]; hence no constraint stabilization issue occurs). However the dissipation modeling may not be trivial, see section 3.1, and parameter estimation may be a hard task as well. Also the numerical simulation with classical schemes may require some careful analysis [59, section 2.2]. Coulomb's friction may or may not be regularized, see section 3.2.2. [12, 13] model frictionless impacts with a linear vibration dynamics $A_{nu}^{-1}(q)\dot{U}_n + KQ_n = 0$ (see (17) for notation), where $\dot{Q}_n = U_n$. Piecewise linear force/indentation laws are used to model elasto-plasticity, see Fig. 1 (f). Positions *q* are supposed constant (hence it could be in Group 3, though it doesn't rely on the impulse time-scale). Redundancy of constraints with indeterminate configurations is treated. Simple examples demonstrate the complexity of multiple impacts with various compression and expansion phases at each contact (in a way similar to Fig. 2(b)–2(g)). Concerning finite-element modeling of impacts, let us quote [18]: ... *single degree of freedom constitutive models...enable significantly more efficient simulations of impact events than high fidelity finite element simulations as only quantities such as the contact forces and contact areas as functions of penetration depth are calculated. In contrast to the single degree of freedom employed by these models, high-fidelity finite element models of the same phenomena can require up to millions of degrees of freedom to ensure a con- vergent response*. It is reported in [65] that the ABAQUS simulation of a single impact event of a sphere against a beam, requires between 4 and 12 hours.
## 4.5 Numerical Simulation

It is impossible to separate the models from the discrete-time numerical schemes with which they may be implemented. Each of the above classes possesses peculiar features from this point of view. Group 1 models lend themselves very well to event-capturing time-stepping schemes, see [79, 4, 5, 99], and can allow for very large number of events [1]. This is due to the fact that backward Euler event-capturing time-stepping methods yield discrete-time systems which are very close to the impact dynamics. They cast the computation of the post-impact velocity in a single generalized equation, for which efficient numerical solvers are available. Academic open-source software packages are developed SICONOS<fnref n="2" /> [4, Chapter 14] [1], LMGC90³, CHRONO<fnref n="4" />, SO-BOGUS<fnref n="5" />. Group 3 models usually require a separated time-integration of the impact process itself, which may in turn require advanced numerical schemes, especially when set-valued friction is used (like the LZB approach). Groups 2 and 3 can be used in event-driven algorithms (switching back and forth between impulse and time domains), which seriously limits the number of events and contacts but provides high accuracy in case of few, separated events. In this respect the approach in [55] may represent a significant progress. Group 4 models usually are integrated with classical adaptive high order methods available in MATLAB or other commercial toolboxes (usually adaptive Runge-Kutta integrators). This aptitude and apparent ease of simulation (avoiding the development of specific time-discretization methods) explain their success in several commercial software packages. But, very large stiffnesses may be needed [90], implying very small time-steps. Moreover they require accurate event detection. This limits the number of events. They usually require many parameters whose mechanical meaning may be unclear, and may be quite uneasy to estimate [97]. Finally they still are nonsmooth so that time-integration may require some care. In summary, Group 4 approach is often numerically inefficient.

## 4.6 Multiple Impacts with Coulomb's Friction

Roughly speaking, all the formalisms in section 3.2 can be adapted to collisions of class (a) by replacing forces by impulses $\lambda_{nu}(t)$ and $\lambda_{tu}(t)$, as done in (19)–(21). Models of class (b) allow us to write Coulomb’s law with infinitesimal force impulses [71, Equ. (64)-(66)] [28], hence do not suffer from

<footnote n="2">²https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos/</footnote>

<footnote n="3">³http://mimetics-engineering.fr/index.php/en/lmgc90-2/</footnote>

<footnote n="4">⁴https://www.projectchrono.org/</footnote>

<footnote n="5">⁵https://gitlab.inria.fr/elan-public-code/so-bogus</footnote>
inconsistencies met with class (a) models which approximate Coulomb's friction at the impulse level. But, class (a) models allow for very robust and reliable numerical integration even for high-dimension systems, while class (b) models require event-driven schemes which are more delicate to tune. Other class (a) approaches are in [99, 11], with cone facetization and a specific representation of Coulomb's friction inspired from [63, 64].

## 4.7 Comparisons between Multiple-Impact Models

Comparing the models' capabilities is important in order to better determine their respective fields of applications. Let us report few comparisons made on chains of aligned balls in [81, 82], for Moreau's model (Group 1), the binary collision approach (Group 2) and the LZB model (Group 3). Fig. 8(a) and 8(b) compare these models from the dispersion capabilities point of view for the chain in Fig. 2(a). Clearly energy dispersion depends strongly on stiffnesses ratios, elasticity coefficient, and on kinetic angles (see (15)). Moreau and binary collisions models provide outcomes similar to LZB, only in certain ranges of elasticity and kinetic angles. Comparisons on longer chains (100-ball chains) are made in [82]. They show that some conclusions drawn for 3 balls, no longer hold for longer chains, because of different wave behaviours. Other comparative works are presented in [28] (3D rocking block, between LZB and their model), [89] (various simple systems, between Pfeiffer-Glocker's model and their model). Theoretical comparisons are also made in [50, 51, 52].

<figure><img src="image_17.png" /><figcaption>(a) Relation $C_{KE} - k_2/k_1$ for different $\eta$, (b) $C_{KE}$ versus $\theta_{12}$ for different $\eta$, $k_1 = k_2$. $m_1 = m_2 = m_3 = 1$.</figcaption></figure>

Figure 8: Dispersion in an elastic 3-ball chain, [81, Chapter 6].
## 4.8 Experimental Validations

Since Mechanics is primarily an experimental science, testing the validity of mathematical models is mandatory. Studies with thorough experimental val- idations are available in (obviously this list is nonexhaustive): [62] (robotics and rocking blocks, class **a**), [81, 80, 115, 114, 69, 117, 110] (aligned chains of balls -tapered, stepped, dimer, disordered, monodisperse, polydisperse-, 2D rocking blocks, 2D and 3D bouncing dimers, disc-ball impacts, class **b**), [7, 104] (multibody system with joint clearance, class **a**), [107, 106] (robotics, class **a**), [86, 6] (chains of balls, Frémond matrices, class **a**), [40, 118] (monodisperse, multi-stepped chains, class **c**), [30] (clothes dy- namics, class **a**), [93] (clamped slender structures, plates subject to fric- tional contact, class **a**), [94] (3D granular matter, class **a**), [88] (aligned chains of beads, class **c**, group 4), [27, 35, 101] (aligned chains, 2D rocking block, dimers, class **b**), [47, 48] (planar rocking block, class **a**). Cur- rently, Hertz’ elasticity is widely used in the field of chains of aligned balls analysis, where the effect of dissipation is shown to be crucial [25], and cap- turing nonlinear wave effects is necessary (class **b** and class **c**) models allow this).

# 5 Conclusions

This chapter introduces multiple impacts phenomena in multibody dynami- cal systems. Multiple impacts possess very specific features, not encountered in single impacts, which make their study an interesting topic. The main modeling approaches proposed to date are summarized and categorized de- pending on some basic assumptions. Extensive comparative works should be tackled on other systems than chain of aligned balls. Benchmarks could be rocking blocks, chains of rotating balls with friction, packed blocks, to cite a few. The number of contacts and frequency of events are crucial test parameters. Stochastic models constitute a largely open topic.

## References

[1] V. Acary, O. Bonnefon, M. Brémond, O. Huber, F. Pérignon, and S. Sinclair. An introduction to SICONOS. Technical report, INRIA Grenoble, University Grenoble Alpes, Grenoble, France, 2019. https://hal.inria.fr/inria-00162911v3/document.
[2] V. Acary, M. Brémond, and O. Huber. On solving contact problems with Coulomb friction: formulations and numerical comparisons. In R. Leine, V. Acary, and O. Brüls, editors, *Advanced Topics in Nonsmooth Dynamics. Transactions of the European Network for Nonsmooth Dynamics*, pages 375–457, 2018. https://inria.hal.science/hal-01630836/document.

[3] V. Acary and B. Brogliato. Concurrent multiple impacts modelling: case study of a 3-ball chain. In K.J. Bathe, editor, *Computational Fluid and Solid Mechanics 2003*, pages 1842–1847, Oxford, 2003. Elsevier Science Ltd.

[4] V. Acary and B. Brogliato. *Numerical Methods for Nonsmooth Dynamical Systems. Applications in Mechanics and Electronics*, volume 35 of LNACM. Springer-Verlag, Berlin, 2008.

[5] V. Acary and N.A. Collins-Craft. On the Moreau-Jean scheme with the Frémond impact law. Energy conservation and dissipation properties for elastodynamics with contact impact and friction. https://inria.hal.science/hal-04230941v2/document, February 2024.

[6] U. Aeberhard, M. Payr, and C. Glocker. Theoretical and experimental treatment of perfect multi-contact-collisions. In *Proc. 3rd Asian Conf. Multibody Dynamics ACMD06*, University of Tokyo, August 2006. Jap. Soc. Mech. Eng. Paper A000602.

[7] N. Akhadkar, V. Acary, and B. Brogliato. Multibody systems with 3D revolute joints with clearances: an industrial case study with an experimental validation. *Multibody System Dynamics*, 42:249–282, 2018.

[8] J. Alves, N. Peixinho, M. Tavares da Silva, P. Flores, and H.M. Lankarani. A comparative study of the viscoelastic constitutive models for frictionless contact interfaces in solids. *Mechanism and Machine Theory*, 85:172–188, 2015.

[9] M. Anitescu. Optimization-based simulation of nonsmooth rigid multi-body dynamics. *Math. Program.*, 105:113–143, 2006.

[10] M. Anitescu and G.D. Hart. A fixed-point iteration approach for multi-body dynamics with contact and friction. *Mathematical Programming, Series B*, 101(1):3–32, 2004.
[11] M. Anitescu, F.A. Potra, and D.E. Stewart. Time-stepping for three-dimensional rigid body dynamics. *Computer Methods in Applied Mechanics and Engineering*, 177(3):183-197, 1999.

[12] A. Barjau, J.A. Batlle, and J.M. Font-Llagunes. Combining vibrational linear-by-part dynamics and kinetic-based decoupling of the dynamics for multiple smooth impacts with redundancy. *Multibody System Dynamics*, 31:497-517, 2014.

[13] A. Barjau, J.A. Batlle, and J.M. Font-Llagunes. Combining vibrational linear-by-part dynamics and kinetic-based decoupling of the dynamics for multiple elastoplastic smooth impact. *Multibody System Dynamics*, 35:233-256, 2015.

[14] J. Bastien and C.H. Lamarque. Persoz' gephyroidal model model described by a maximal monotone differential inclusion. *Arch. Appl. Mechanics*, 78(5):393-407, 2008.

[15] R. Berthold, J. Burgner-Kahrs, M. Wangenheim, and S. Khams. Investigating frictional contact behavior for soft material robot simulations. *Meccanica*, 58:2165-2176, 2023.

[16] F. Bourrier and V. Acary. Predictive capabilities of 2D and 3D block propagation models integrating block shape assessed from field experiments. *Rock Mechanics and Rock Engineering*, 55:591-609, 2022. https://hal.science/hal-03155240/file/main_figend_NCC_new.pdf.

[17] M.R. Brake. An analytical elastic-perfectly plastic contact model. *Int. J. Solids Struct.*, 49:3129-3141, 2012.

[18] M.R. Brake. The role of epistemic uncertainty of contact models in the design and optimization of mechanical systems with aleatoric uncertainty. *Nonlinear Dynamics*, 77(3):899-922, 2014.

[19] M.R. Brake. An analytical elastic-plastic contact model with strain hardening and frictional effects for normal and oblique impacts. *Int. J. Solids Struct.*, 62:104-123, 2015.

[20] B. Brogliato. *Nonsmooth Impact Mechanics. Models, Dynamics and Control*, volume 220 of LNCIS. Springer-Verlag, London, 1st edition, 1996.
[21] B. Brogliato. Kinetic quasi velocities in unilaterally constrained Lagrangian mechanics with impacts and friction. *Multibody Syst. Dynamics*, 32(2):175–216, 2014.

[22] B. Brogliato. *Nonsmooth Mechanics. Models, Dynamics and Control*. Communications and Control Eng. Springer International Publishing Switzerland, 3rd edition, 2016.

[23] B. Brogliato, H. Zhang, and C. Liu. Analysis of a generalized kinematic impact law for multibody-multicontact systems, with application to the planar rocking block and chains of balls. *Multibody System Dynamics*, 27(3):351–382, 2012.

[24] H.A. Burgoyne, J.A. Newman, W.C. Jackson, and C. Daraio. Guided impact mitigation in 2D and 3D granular crystals. *Procedia Engineering*, 103:52–59, 2015. Proceedings of the 2015 Hypervelocity Impact Symposium (HVIS 2015).

[25] R. Carretero-González, D. Khatri, Mason A. Porter, P. G. Kevrekidis, and C. Daraio. Dissipative solitary waves in granular crystals. *Phys. Rev. Lett.*, 102:024102, Jan 2009.

[26] F. Caselli and M. Frémond. Collision of three balls on a plane. *Comput. Mech.*, 43:743–754, 2009.

[27] V. Ceanga and Y. Hurmuzlu. A new look at an old problem: Newton's cradle. *ASME J. Applied Mech.*, 68(4):575–583, 2001.

[28] A. Chatterjee and A. Bowling. Modeling three-dimensional surface-to-surface rigid contact and impact. *Multibody System Dynamics*, 46:1–40, 2019.

[29] A. Chatterjee, H. Ghaednia, A. Bowling, and M. Brake. Estimation of impact forces during multi-point collisions involving small deformations. *Multibody System Dynamics*, 51(1):45–90, 2021.

[30] F. Coltraro, J. Amoros, M. Alberich-Carraminana, and C. Torras. A novel collision model for inextensible textiles and its experimental validation. *Applied Mathematical Modelling*, 128:287–308, 2024.

[31] E. Corral, R.G. Moreno, M.J.G. García, and C. Castejón. Nonlinear phenomena of contact in multibody systems dynamics: a review. *Nonlinear Dynamics*, 104:1269–1295, 2021.
[32] A. Cosimo, F.J. Cavalieri, A. Cardona, and O. Brüls. On the adaptation of local impact laws for multiple impact problems. *Multibody System Dynamics*, 102:1997–2016, 2020.

[33] R.W. Cottle, J.S. Pang, and R.E. Stone. *The Linear Complementarity Problem*. Academic Press, 1992.

[34] A.W. Crook. A study of some impacts between metal bodies by a piezoelectric method. *Proc. Royal. Soc. A. Math. Phys. Eng. Sci.*, 212(1110):377–390, 1952.

[35] C. Cuneyt, M. Gharib, and Y. Hurmuzlu. Solving frictionless rocking block problem with multiple impacts. *Proc. R. Soc. A*, 465:3323–3339, 2009.

[36] G. Daviet. Interactive hair simulation on the GPU using ADMM. In *SIGGRAPH '23: ACM SIGGRAPH 2023 Conference Proceedings*, pages 1–11, July 2023. Article no 24.

[37] G. De Saxcé and Z.-Q. Feng. The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms. *Mathematical and Computer Modelling*, 28(4):225–245, 1998. Recent Advances in Contact Mechanics.

[38] S. Ding, Y. Hu, B. Jian, Y. Zhang, R. Xia, and G. Hu. A review and comparative analysis of normal contact force models for viscoelastic particles. *International Journal of Impact Engineering*, 189:104968–28, 2024.

[39] E. Falcon, C. Laroche, S. Fauve, and C. Coste. Collision of a 1-d column of beads with a wall. *Eur. Phys. J. B*, 5:111–131, 1998.

[40] Y. Feng, W. Kang, D. Ma, and c. Liu. Multiple impacts and multiple-compression process in the dynamics of granular chains. *Journal of Computational and Nonlinear Dynamics*, 14(12):121002, 10 2019.

[41] M. Frémond. Rigid bodies collisions. *Phys. Letters A*, 204(1):33–41, 1995.

[42] M. Frémond. *Non-smooth Thermomechanics*. Springer, 2002.

[43] M. Frémond. *Collisions*. Università di Roma Tor Vergata, 2007. Divertimento di Ingegneria Civile.
[44] D. Gale. An indeterminate problem in classical mechanics. *Am. Math. Mon.*, 59(5):291-295, 1952.

[45] F. Génot and B. Brogliato. New results on Painlevé paradoxes. *European Journal of Mechanics A/Solids*, 18(4):653-677, 1999.

[46] H. Ghaednia, X. Wang, S. Saha, Y. Xu, A. Sharma, and R.L. Jackson. A review of elastic-plastic contact mechanics. *Applied Mechanics Reviews*, 69(6):060804, 11 2017.

[47] A.I. Giouvanidis and E.G. Dimitrakopoulos. Nonsmooth dynamic analysis of sticking impacts in rocking structures. *Bull. Earthquake Eng.*, 15:2273-2304, 2017.

[48] A.I. Giouvanidis and E.G. Dimitrakopoulos. Nonsmooth modelling of impacts in rocking structures with Poisson's law. In M. Papadrakakis and M. Fragiadakis, editors, *Proc. COMPDYN 2017, 6th ECCO-MAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering*, volume 2, pages 2910-2925, Rhodes Island, Greece, June 2017.

[49] A.I. Giouvanidis, E.G. Dimitrakopoulos, and P.B. Lourenço. Chattering: an overlooked peculiarity of rocking motion. *Nonlinear Dynamics*, 109:459-477, 2022.

[50] C. Glocker. On frictionless impact models in rogod-body systems. *Phil. Trans. Roy. Soc. London*, 359(1789):2385-2404, 2001.

[51] C. Glocker. *Set-Valued Force Laws: Dynamics of Non-Smooth Systems*, volume 1 of LNACM. Springer-Verlag, Heidelberg, 2001.

[52] C. Glocker. An introduction to impacts. In J. Haslinger and G. Stavroulakis, editors, *Nonsmooth Mechanics of Solids*, volume 458 of CISM Courses and Lectures, pages 45-101. Springer, 2006.

[53] W. Goldsmith. *Impact: The Theory and Physical Behavior of Colliding Solids*. E. Arnold Publishers, London, 1960.

[54] E. Gourc, G. Michon, S. Seguy, and A. Berlioz. Targeted energy transfer under harmonic forcing with a vibro-impact nonlinear energy sink: analytical and experimental developments. *Journal of Vibration and Acoustics*, 137(3):031008-7, 2015.
[55] M. Halm and M. Posa. Set-valued rigid-body dynamics for simul-
taneous, inelastic, frictional impacts. *The International Journal of Robotics Research*, 2024. https://dair.seas.upenn.edu/assets/pdf/Halm2023.pdf.

[56] M. Hjiaj, G. de Saxcé, and Z. Mróz. A variational inequality-based formulation of the frictional contact law with a non-associated sliding rule. *European Journal of Mechanics - A/Solids*, 21(1):49-59, 2002.

[57] Y. Hurmuzlu and D.B. Marghitu. Rigid body collisions of planar kinematic chains with multiple contact points. *The International Journal of Robotics Research*, 13(1):82–92, 1994.

[58] A.P. Ivanov. Impacts in a system with certain unilateral couplings. *J. of Applied Mathematics and Mechanics*, 51(4):436–442, 1987.

[59] G. James. Traveling fronts in dissipative granular chains and nonlinear lattices. *Nonlinearity*, 34(3):1758, 2021.

[60] Y.B. Jia, M.T. Mason, and M.E. Erdmann. Multiple impacts: A state transition diagram approach. *The Int. J. of Robotics Research*, 32(1):84–114, 2013.

[61] K.L. Johnson. *Contact Mechanics*. Cambridge University Press, Cambridge, 1985.

[62] M.J. Jongeneel, L. Poort, N. van de Wouw, and A. Saccon. Experimental validation of nonsmooth dynamics simulations for robotic tossing involving friction and impacts. https://hal.science/hal-03974604, February 2023.

[63] A. Klarbring. A mathematical programming approach to three-dimensional contact problems with friction. *Computer Methods in Applied Mechanics and Engineering*, 58(2):175–200, 1986.

[64] A. Klarbring and G. Björkman. A mathematical programming ap-
proach to contact problems with friction and varying contact surface.
*Computers and Structures*, 30(5):1185–1198, 1988.

[65] M. Krack. *Exploiting the Use of Strong Nonlinearity in Dynamics and Acoustics*, volume 613 of *Courses and Lectures*, chapter Systems with Contact Nonlinearities, pages 235–272. CISM International Center for Mechanical Sciences, 2024.
[66] R.I. Leine, A. Schweizer, M. Christen, J. Glover, P. Bartelt, and W. Gerber. Simulation of rockfall trajectories with consideration of rock shape. *Multibody Syst. Dyn.*, 32:241–271, 2014.

[67] R.I. Leine and N. van de Wouw. *Stability and Convergence of Mechanical Systems with Unilateral Constraints*, volume 36 of LNACM. Springer-Verlag, Berlin, 2008.

[68] J. Li, G. Daviet, R. Narain, F. Bertails-Descoubes, M. Overby, G. Brown, and L. Boissieux. An implicit frictional contact solver for adaptive cloth simulation. *ACM Transactions on Graphics*, 37(4):1–15, 2018. article 52.

[69] C. Liu, H. Zhang, Z. Zhao, and B. Brogliato. Impact-contact dynamics in a disc-ball system. *Proc. R. Soc. A*, 469:1–20, 2013. paper 20120741.

[70] C. Liu, Z. Zhao, and B. Brogliato. Frictionless multiple impacts in multibody systems. I. Theoretical framework. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 464(2100):3193–3211, 2008.

[71] C. Liu, Z. Zhao, and B. Brogliato. Variable structure dynamics in a bouncing dimer. RR-6718, INRIA, https://hal.inria.fr/inria-00337482, November 2008.

[72] C. Liu, Z. Zhao, and B. Brogliato. Frictionless multiple impacts in multibody systems. II. numerical algorithm and simulation results. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 465(2101):1–23, 2009.

[73] D. Ma and C. Liu. Contact law and coefficient of restitution in elasto-plastic spheres. *Journal of Applied Mechanics*, 82(12):121006, 2015.

[74] D.B. Marghitu and Y. Hurmuzlu. Three-dimensional rigid-body collisions with multiple contact points. *Journal of Applied Mechanics*, 62(3):725–732, 09 1995.

[75] M.D.P. Monteiro Marques. *Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction*. Birkhäuser Verlag, Basel, CH, 1993.

[76] D. Maugis. *Contact, Adhesion and Rupture of Elastic Solids*. Solid-State Sciences. Springer, Heidelberg, 2000.
[77] J.J. Moreau. La notion de sur-potentiel et les liaisons unilatérales en élastostatique. *C.R. des séances de l'Académie des sciences*, 267:954-957, 1968. https://hal.science/hal-01868127/document.

[78] J.J. Moreau. Application of convex analysis to some problems of dry friction. In H. Zorski, editor, *Trends in Applications of Pure Math. to Mechanics*, volume 2, pages 99-121, London, 1979. Pitman.

[79] J.J. Moreau. Unilateral contact and dry friction in finite freedom dynamics. In *Nonsmooth Mechanics and Applications*, volume 302 of *CISM Courses and Lectures*, pages 1-82, International Center for Mechanical Sciences, Udine, Italy, 1988. Springer Verlag. https://hal.archives-ouvertes.fr/hal-01713847/document.

[80] N.S. Nguyen and B. Brogliato. Shock dynamics in granular chains: numerical simulations and comparison with experimental tests. *Granular Matter*, 14:341-362, 2012.

[81] N.S. Nguyen and B. Brogliato. *Multiple Impacts in Dissipative Granular Chains*, volume 72 of LNACM. Springer-Verlag, Berlin Heidelberg, 2014.

[82] N.S. Nguyen and B. Brogliato. *Advanced Topics in Nonsmooth Dynamics*, chapter Comparisons of multiple-impact laws for multibody systems: Moreau's law, binary impacts, and the LZB approach, pages 1-45. Transactions of the European Network for Nonsmooth Dynamics. Springer Int. Publishing AG, 2018.

[83] F.Z. Nqi. *Etude Numérique de divers Problèmes Dynamiques avec Impact et de leur Propriétés Qualitatives*. PhD thesis, Université Claude Bernard Lyon 1, Lyon, France, 1997.

[84] A. Nurkanovic, S. Albrecht, B. Brogliato, and M. Diehl. The time-freezing reformulation for numerical optimal control of complementarity Lagrangian systems with state jumps. November 2021. https://hal.inria.fr/hal-03427800/document.

[85] A. Nurkanovic, S. Albrecht, B. Brogliato, and M. Diehl. The time-freezing reformulation for numerical optimal control of complementarity Lagrangian systems with state jumps. *Automatica*, 158:111295, 2023.
[86] M. Payr, C. Glocker, and C. Bösch. Experimental treatment of multiple-contact-collisions. In *Proc. ENOC-2005*, pages 450–459, Eindhoven, Netherlands, August 2005. Euromech.

[87] F. Pfeiffer and C. Glocker. *Multibody Dynamics with Unilateral Contacts*. Wiley, New York, 1996.

[88] M.A. Porter, C. Daraio, I. Szelengowicz, E.B. Herbold, and P.G. Kevrekidis. Highly nonlinear solitary waves in heterogeneous periodic granular media. *Physica D: Nonlinear Phenomena*, 238(6):666–676, 2009.

[89] S. Rakshit and A. Chatterjee. Scalar generalization of Newtonian resti- tution for simultaneous impact. *Int. J. of Mech. Sciences*, 103:141–157, 2015.

[90] A. Raoofian, X. Dai, and J. Kövecses. Contact representation in robotic mechanical systems employing reduced models. *IEEE Robotics and Automation Letters*, 9(2):1556–1563, 2023.

[91] G. Rill, T. Schaeffer, and M. Schuderer. Lugre or not Lugre. *Multibody Syst. Dyn.*, 60:191–218, 2024.

[92] R.T. Rockafellar. *Convex Analysis*. Princeton Landmarks in Mathe- matics. Princeton University Press, New Jersey, 1970.

[93] V. Romero, M. Ly, A.H. Rasheed, R. Charrondière, A. Lazarus, S. Neukrich, and F. Bertails-Descoubes. Physical validation of sim- ulators in computer graphics: A new framework dedicated to slender elastic structures and frictional contact. *ACM Transactions on Graph- ics*, 40(4):1–19, 2021. Article 66.

[94] G. Rousseau, T. Metivet, H. Rousseau, G. Daviet, and F. Bertails- Descoubes. Revisiting the role of friction coefficients in granular collapses: confrontation of 3-D non-smooth simulations with exper- iments. *Journal of Fluid Mechanics*, 975:1–39, 2022. Article A14.

[95] E. Sanchez, A. Cosimo, O. Brüls, A. Cardona, and F.J. Cavalieri. Non- smooth numerical solution for Coulomb friction, rolling and spinning resistance of spheres applied to flexible multibody system dynamics. *Multibody System Dynamics*, 59(1):69–103, 2023.
[96] M. Schatzman. Uniqueness and continuous dependence on data for one-dimensional impact problem. *Math. Comput. Modell.*, 28(4-8):1-18, 1998.

[97] M. Schuderer, G. Rill, T. Schaeffer, and C. Schulz. Friction modeling from a practical point of view. *Multibody System Dynamics*, 2024. https://doi.org/10.1007/s11044-024-09978-0.

[98] L. Skrinjar, J. Slavic, and M. Boltezar. A review of continuous contact-force models in multibody dynamics. *International Journal of Mechanical Sciences*, 145:171-187, 2018.

[99] D.E. Stewart and J.C. Trinkle. An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction. *Int. J. Numerical Methods in Engineering*, 39(15):2673-2691, 1996.

[100] W.J. Stronge. Rigid body collision with friction. *Proc. Royal Soc. Lond. A*, 431(1881):169-181, 1990.

[101] A. Tavakoli, M. Gharib, and Y. Hurmuzlu. Collision of two mass baton with massive external surfaces. *Journal of Applied Mechanics*, 79(5):051019, 2012.

[102] W. Thomson and P.G. Tait. *Treatise on Natural Philosophy*. Clarendon Press, Oxford, 1867.

[103] C. Thornton. Coefficient of restitution for colinear collisions of elastic-perfectly plastic spheres. *ASME J. Applied Mechanics*, 64:383-386, 1997.

[104] T. Thümmel. *Experimentelle Mechanismendynamik: Messung, Modellierung, Simulation, Verifikation, Interpretation und Beeinflussung typischer Schwingungsphänomene an einem Mechanismenprüfstand*. PhD thesis, München, Technische Universität München, Habil.-Schr., 2012. https://mediatum.ub.tum.de/doc/1108538/document.pdf.

[105] D.H. Towne and C.R. Hadlock. One-dimensional collisions and Chebyshev polynomials. *American Journal of Physics*, 45(3):255-259, 1977.

[106] A.A. Transeth, R.I. Leine, C. Glocker, and K.Y. Pettersen. 3-D snake robot motion: Nonsmooth modeling, simulations, and experiments. *IEEE Transactions on Robotics*, 24(2):361-376, 2008.
[107] A.A. Transeth, R.I. Leine, C. Glocker, K.Y. Pettersen, and P. Liljebäck. Snake robot obstacle-aided locomotion: Modeling, simulations, and experiments. *IEEE Transactions on Robotics*, 24(1):88-104, 2008.

[108] Y. Tsuji, T. Tanaka, and T. Ishida. Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. *Powder Technology*, 71(3):239-250, 1992.

[109] G. Wang and C. Liu. Further investigation on improved viscoelastic contact force model extended based on Hertz' law in multibody system. *Mechanism and Machine Theory*, 153:103986, 2020.

[110] J. Wang, C. Liu, and Z. Zhao. Nonsmooth dynamics of a 3D rigid body on a vibrating plate. *Multibody System Dynamics*, 32:217-239, 2014.

[111] T. Winandy and R.I. Leine. A maximal monotone impact law for the 3-ball Newton's cradle. *Multibody System Dynamics*, 39:79-94, 2017.

[112] X. Xiong, R. Kikuuwe, and M. Yamamoto. A multiscale friction model described by continuous differential equations. *Tribol. Letters*, 51:513-523, 2013.

[113] Z. Yang, J. Hong, D. Wang, R. Cheng, and Y. Ma. Vibration analysis of rotor systems with bearing clearance using a novel conformal contact model. *Nonlinear Dynamics*, 2024. https://doi.org/10.1007/s11071-024-09489-9.

[114] H. Zhang, B. Brogliato, and C. Liu. Dynamics of planar rocking-blocks with Coulomb friction and unilateral constraints: comparisons between experimental and numerical data. *Multibody System Dynamics*, 32:1-25, 2014.

[115] H. Zhang, C. Liu, Z. Zhao, and B. Brogliato. Energy evolution in complex impacts with friction. *Sci. China Phys. Mech. Astron.*, 56:875-881, 2013.

[116] Z. Zhao, C. Liu, and B. Brogliato. Energy dissipation and dispersion effects in granular media. *Physical Review E*, 78(3), 2008. Paper 031307.

[117] Z. Zhao, C. Liu, and B. Brogliato. Planar dynamics of a rigid body system with frictional impacts. II. Qualitative analysis and numerical
cal simulations. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 465(2107):2267–2292, 2009.

[118] Z. Zhou, D.M. McFarland, and A.F. Vakakis. One-dimensional granular chains as transmitted force attenuators. *Nonlinear Dynamics*, 111:14713–14730, 2023.