# On the Moreau–Jean scheme with the Frémond impact law: energy conservation and dissipation properties for elastodynamics with contact, impact and friction

Vincent Acary, Nicholas Anton Collins-Craft

## ► To cite this version:

Vincent Acary, Nicholas Anton 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. Journal of Theoretical, Computational and Applied Mechanics, 2025, pp.1–30. 10.46298/jt-cam.13480. hal-04230941v5

HAL Id: hal-04230941

https://inria.hal.science/hal-04230941v5

Submitted on 28 Jun 2025

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
**Identifiers**

DOI 10.46298/jtcam.13480

HAL hal-04230941v5

**History**

Received Apr 26, 2024

Accepted Oct 28, 2024

Published Jun 30, 2025

**Associate Editor**

Alexander POPP

**Reviewers**

Anonymous

Simon EUGSTER

Mike PUSO

Remco LEINE

**Open Review**

HAL hal-04819313

**Supplementary Material**

10.5281/zenodo.15750400.

See also addendum for
additional information on
datasets and software
solutions.

**Licence**

CC BY 4.0

©The Authors

# On the Moreau-Jean scheme with the Frémond impact law: energy conservation and dissipation properties for elastodynamics with contact, impact and friction

Vincent ACARY and

Nicholas Anton COLLINS-CRAFT

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000, Grenoble

The objective of this paper is to propose a time integration scheme for nonsmooth mechanical systems involving one-sided contact, impact and Coulomb friction, that respects the principles of discrete-time energy balance with positive dissipation. To obtain energetic consistency in the continuous time model when an impact occurs, we work with an impact law with friction inspired by the work of M. Frémond which ensures that dissipation is positive, i.e. that the Clausius-Duhem inequality is satisfied for the impulses and the velocity jumps. On this basis, we propose a time integration method based on the Moreau-Jean scheme with a discrete version of the Frémond impact law, and show that this method has correct dissipation properties.

**Keywords:** contact dynamics, impact, computational contact mechanics, numerical time integration, energy conservation, dissipation properties

**Notations** The following notations are used throughout the paper. The Euclidean norm for a vector $x \in \mathbb{R}^n$ is denoted by $\|x\|$. For two vectors $x, y \in \mathbb{R}^n$, the Hadamard product is denoted by $x \circ y$. For a positive definite (respectively positive semi-definite) matrix $M \in \mathbb{R}^{n \times n}$, $\|x\|_M$ denotes the norm (respectively the semi-norm) in the metric defined by $M$. Let $I$ denote a real time interval of any sort. For a function $f: I \to \mathbb{R}^n$ of Bounded Variation (BV), we denote the right-limit function by $f^+(t) = \lim_{s \to t, s > t} f(s)$, and respectively the left-limit by $f^-(t) = \lim_{s \to t, s < t} f(s)$. We denote by $0 = t_0 < t_1 < \dots < t_k < \dots < t_N = T$ a finite partition (or a subdivision) of the time interval $[0, T]$, where $(T > 0)$. The integer $N$ is the number of time intervals in the subdivision. The length of a time step is denoted by $h_k = t_{k+1} - t_k$. For simplicity's sake, the schemes are presented in the following with a time step denoted by $h$ for short. The value of a real right-continuous function $x(t)$ at the time $t_k$, is approximated by $x_k$. In the same way, the notation $x_{k+\theta} = (1-\theta)x_k + \theta x_{k+1}$ is used for $\theta \in [0, 1]$. The notation $dt$ defines the Lebesgue measure on $\mathbb{R}$.

# 1 Introduction and motivations

The objective of this paper is to propose a time integration scheme for nonsmooth mechanical systems involving one-sided contact, impact and isotropic Coulomb friction, that respects the principles of discrete-time energy balance with positive dissipation. For simplicity, we consider systems that are either discrete or spatially discretised (by the finite element method (FEM) for example), and that have dynamics that are linear, but with possible nonlinear constraints that model contact. In order to develop a consistent integration scheme with correct energy and dissipation properties, the continuous-time model must also have these properties. It is known that discrete systems with Coulomb friction and a kinematic impact law can generate energy, in particular, when the direction of sliding velocity changes during impact (Brogliato 1999; Brogliato 2016; Glocker 2013).

In this work, we therefore propose to work with an impact law inspired by the work of M. Frémond (Frémond 1995; Frémond 2001; Frémond 2002; Frémond 2017; Frémond 2007) which
ensures that dissipation is positive, i.e. that the Clausius-Duhem inequality (and therefore the second law of thermodynamics) is satisfied. The main ingredient is the use of the average of the pre- and post-impact velocities as the primary kinematic variable conjugated with the impulse. This duality pairing is natural when an impact occurs. In the model, the Coulomb law is written with this average velocity. In this way, the work of the contact impulses is always non-positive. When the velocity is continuous, the model reverts to the standard one-sided contact model with Coulomb friction. With the Frémond model including three-dimensional Coulomb friction, positive dissipation is ensured.

On this basis, we propose a time integration method based on the Moreau-Jean scheme (Jean and Moreau 1987; Moreau 1988) and show that this method has desirable energy dissipation properties: the energy is conserved for conservative systems and the discrete dissipation approximated by the numerical scheme is always positive under reasonable assumptions.

Much work in the literature has focused on conservative or energy-dissipating time integration schemes for frictionless one-sided contact systems. These studies are mainly interested in this property because it confers stability properties to the schemes, which are often difficult to obtain. A review of the literature can be found in (Acary 2016) for elasto-dynamics and in (Greenhalgh et al. 2013) for linear complementarity systems. On the other hand, there is comparatively less literature on contact with Coulomb friction. This may be due to the inherent difficulty of Coulomb friction for discrete systems that exhibit impulses. Amongst the work that has been conducted in the context of variational discrete integrators, we may note the article of Fetecau et al. (2003), who briefly consider the case of friction in the continuous-time case, but do not consider it in their discrete algorithm (which in any case requires the precise location of impacts in time, rendering it infeasible for applications with many impacts such as granular flows). The autonomous case for rigid body dynamics has been treated by Johnson et al. (2014), who allow a frictional impulse over a time step to be sufficiently large to cause the contact to stick (but no larger), thus avoiding any change of sliding direction. However, in the general case of non-autonomous mechanical systems (those with externally imposed driving forces or displacements), and when there are interactions with elasticity and damping, changes in the sliding direction may simply be imposed at the contact, regardless of any interaction with frictional forces. Capobianco and Eugster (2018) considered such a system, and developed an integration scheme that preserves a principle of virtual action, but this does not in itself imply no generation of energy during sliding, only that any energy that is created will be compensated elsewhere in the action.

Undoubtedly the most important contributions that have been made to the study of the energetic properties of nonsmooth impact and friction models are those of Leine and Wouw (2008); Glocker (2013). In these works, the authors show that in the case that all of the restitution coefficients in the model are equal (or close to equal), or they are all sufficiently small with respect to the condition number of the Delassus matrix or all the elements of the system are decoupled, the system is guaranteed to be *globally dissipative at impact instants*. However, this does not guarantee that the system is *locally dissipative at every point*. In our opinion, it is this second, stricter condition that it is essential to fulfill in order to have a model that is thermodynamically admissible. In light of these previous works, there would seem to be considerable interest in ensuring that numerical schemes that are event-capturing (that is to say that do not require the precise location in time of each impact event) are able to remain locally dissipative, regardless of the nature (elastic, visco-elastic, rigid, etc.) of the contacting bodies, and the presence of any external forcing in the framework of elasto-dynamics. In this article we propose a numerical scheme that fulfills exactly this property.

The structure of the article is as follows. We begin in Section 2 by presenting the equations of motion of a spatially discrete system with one-sided constraints, Coulomb friction and impacts. After a discussion of the energy balance of the system, we remark that the system with Newton's law of impact and Coulomb friction can be non-dissipative. Finally, we propose an impact model inspired by the work of Frémond which guarantees that the system is dissipative. In Section 3, we propose a new scheme that is dissipative in discrete time, with dissipative terms that are numerical approximations of frictional and impact dissipation. In Section 4, we describe how the schemes are implemented in practice in the Siconos software framework. In Section 5, the dissipation properties of the discrete time scheme are illustrated with various examples and
compared to the classical Moreau-Jean scheme. Section 6 concludes the article and outlines some perspectives.

## 2 Nonsmooth mechanical systems with unilateral contact, Coulomb's friction and impact

In this section we consider systems that are continuous in time but discrete in space, and outline the equations of motion and energy balance with Coulomb friction and impact for the standard model and the Frémond impact model.

### 2.1 Smooth equations of motion for linear visco-elasticity and constitutive laws

The equations of motion of a discrete or discretised Lagrangian mechanical system in the linear visco-elastic case are

$$ q(t_0) = q_0, \quad v(t_0) = v_0, \tag{1a} $$

$$ \dot{q}(t) = v(t), \tag{1b} $$

$$ M\dot{v}(t) + Kq(t) + Cv(t) = F(t) + R(t), \tag{1c} $$

where $q(t) \in \mathbb{R}^n$ is the vector of generalised coordinates and $v(t) = \dot{q}(t)$ is the corresponding vector of generalised velocities, the initial conditions are $q_0 \in \mathbb{R}^n$ and $v_0 \in \mathbb{R}^n$, $M \in \mathbb{R}^{n \times n}$ is the symmetric mass matrix that is assumed to be positive definite, $K \in \mathbb{R}^{n \times n}$ is the positive semi-definite stiffness matrix and $C \in \mathbb{R}^{n \times n}$ the positive semi-definite damping matrix, $F(t)$ is the vector of external applied forces and $R(t)$ is the vector of generalised contact forces. The system is now subjected to a finite set of $m$ unilateral contacts that defines the admissible set for the configuration

$$ C = \{q \in \mathbb{R}^n \mid g_N^\alpha(q) \ge 0, \alpha \in [1, m]\} \subset \mathbb{R}^n, \tag{2} $$

where $g_N^\alpha: \mathbb{R}^n \to \mathbb{R}$ is assumed to be a smooth function with non-vanishing gradients. We denote by the function $g_N: \mathbb{R}^n \to \mathbb{R}^m$ the function with $g_N^\alpha$ as components. For the perfect unilateral constraints, the Signorini condition is written as

$$ 0 \le g_N(q(t)) \perp \lambda(t) \ge 0, \tag{3} $$

where the inequalities involving vectors are understood to hold component-wise and the $\perp$ symbol means that $g_N^\top \lambda = 0$. The Lagrange multiplier $\lambda$ is related to the generalised reaction force $R$ by $R(t) = \nabla_q g_N(q(t)) \lambda(t)$. The Signorini condition at the velocity level can be also be defined as

$$ 0 \le u_N(t) \perp \lambda(t) \ge 0, \text{ if } g_N(q(t)) \le 0, \tag{4} $$

where the local relative velocity $u_N(t)$ is defined by $u_N(t) = \dot{q}(q(t)) = \nabla_q^\top g_N(q(t)) v(t)$. Moreau's viability lemma ensures that Equation (4) implies Equation (3) if the initial conditions are admissible.

**Coulomb friction** For Coulomb friction, the definition of the gradients of the constraints is not sufficient and an orthonormal local basis $(n^\alpha, t_1^\alpha, t_2^\alpha)$ at contact composed of the inward normal vector $n^\alpha \in \mathbb{R}^3$ to the set $C$ and the tangential vectors $t_1^\alpha, t_2^\alpha \in \mathbb{R}^3$ is needed. For an individual contact $\alpha$, the relative velocity and the reaction force $r$ at contact are denoted in this local contact frame as

$$ u^\alpha = \begin{pmatrix} u_N^\alpha \\ u_T^\alpha \end{pmatrix}, \quad r^\alpha = \begin{pmatrix} r_N^\alpha \\ r_T^\alpha \end{pmatrix}, \tag{5} $$

where $u_N^\alpha \in \mathbb{R}$, $u_T^\alpha = (u_{T_1}^\alpha, u_{T_2}^\alpha)^\top \in \mathbb{R}^2$, $r_N^\alpha \in \mathbb{R}$ and $r_T^\alpha = (r_{T_1}^\alpha, r_{T_2}^\alpha)^\top \in \mathbb{R}^2$. The local variables are related to the generalised variables by the contact configuration matrix $H^\alpha(q)$ such that

$$ u^\alpha = H^\alpha(q) v, \quad R^\alpha = H^{\alpha, \top}(q) r^\alpha. \tag{6} $$
Coulomb friction together with the Signorini condition at the velocity level is given by

$$-\tilde{u}^{\alpha} \in N_{\mathcal{K}^{\alpha}}(r^{\alpha}), \quad \text{if } g_{N}^{\alpha}(q) \le 0, \text{ else } r^{\alpha} = 0, \qquad (7)$$

where $N_C$ is the normal cone to a convex set C given by $N_C(x) = \{s \in \mathbb{R}^n | s^T(y - x) \le 0 \text{ for all } y \in C\}$, $\tilde{u}$ is the modified De Saxcé velocity expressed as

$$\tilde{u}^{\alpha} = u^{\alpha} + \begin{pmatrix} \mu^{\alpha} \|u_{T}^{\alpha}\| \\ 0 \end{pmatrix}, \qquad (8)$$

and $\mathcal{K}^{\alpha}$ is the Coulomb cone $\mathcal{K}^{\alpha} = \{r, \|r_T\| \le \mu^{\alpha} r_N\}$, with $\mu^{\alpha} \ge 0$ the isotropic friction coefficient. Let $\mathcal{I} = [1, m] \in \mathbb{N}$ be the set of indices of constraints. Let us define now the index of constraint at the velocity level by $\mathcal{I}^1 = \{\alpha \in \mathcal{I} | g_N^{\alpha}(q(t)) \le 0\}$. By collecting all the variables for each contact in the set $\mathcal{I}^1$, and implicitly introducing the notation dropping $\alpha$ ($x = [x^{\alpha, \top}, \alpha \in \mathcal{I}^1]^{\top}$), we get for the frictional contacts

$$-\tilde{u} \in N_{\mathcal{K}}(r), \qquad (9)$$

considering that $r^{\alpha} = 0$ for $\alpha \notin \mathcal{I}^1$ and $\mathcal{K}$ is the Cartesian product of the cone $\mathcal{K}^{\alpha}$. Using the dual cone of $\mathcal{K}$ denoted by

$$\mathcal{K}^{\star} = \{u, \mu \|u_T\| \le u_N\}, \qquad (10)$$

Equation (9) can be written as $-r \in N_{\mathcal{K}^{\star}}(\tilde{u})$ or as a cone complementarity condition $\mathcal{K}^{\star} \ni \tilde{u} \perp r \in \mathcal{K}$. In the following, we consider a single coefficient of friction, although this does not affect the generality of the results. The reader may refer to (Facchinei and Pang 2003) for more details on complementarity problems on cones and to (Acary et al. 2018) for an application to contact with friction. To be self-contained, the equivalence with the standard form of the Coulomb friction is given in the following Lemma 1.

**Lemma 1** [Acary et al. (2011)] For one contact point, the Coulomb friction law with the Signorini condition at the velocity level fulfills one of the three following conditions:

$$\text{either} \quad r = 0 \text{ and } u_N \ge 0, \qquad \text{(take-off)} \\ \text{or} \quad u = 0 \text{ and } \|r_T\| \le \mu r_N, \qquad \text{(sticking)} \qquad (11)$$

$$\text{or} \quad \|r_T\| = \mu r_N \text{ and } u_N = 0 \text{ and } u_T \ne 0 \text{ and } \|u_T\| \|r_T\| = -\|r_T\| \|u_T\|, \quad \text{(sliding)}$$

and is equivalent to the complementarity problem

$$\mathcal{K}^{\star} \ni \tilde{u} \perp r \in \mathcal{K}, \quad \text{with} \quad \tilde{u} = u + \begin{pmatrix} \mu \|u_T\| \\ 0 \end{pmatrix}. \qquad (12)$$

**Proof 1** $\Rightarrow$ Conditions (11) imply that $r \in \mathcal{K}$ and $u_N \ge 0$. Then $\tilde{u}_N \ge \mu \|u_T\|$ and we conclude that $\tilde{u} \in \mathcal{K}^*$. If $u = 0$, then $\tilde{u} = 0$ and the relation $\tilde{u}^T r = 0$ holds trivially. The same applies for $r = 0$. In the sliding case, we have $\tilde{u}^T r = (u_N + \mu \|u_T\|)r_N + u_T^T r_T = \|u_T\| \|r_T\| + u_T^T r_T$ since $u_N = 0$ and $\mu r_N = \|r_T\|$. Furthermore $\|u_T\| \|r_T\| = -\|r_T\| \|u_T\|$, the vectors $u_T$ and $r_T$ are colinear and have opposite directions, we have $u_T^T r_T = -\|u_T\| \|r_T\|$ and hence $\tilde{u}^T r = 0$.

$\Leftarrow$ Conversely, if (12) holds, we have three separate cases:

1. $r = 0$. The relation $\tilde{u} \in \mathcal{K}^*$ holds. This implies $u_N \ge 0$. We have the take-off case.

2. $\tilde{u} = 0$. This implies that $u = 0$ and we have $r \in \mathcal{K}$. Hence, we have the sticking case.

3. $\tilde{u} \ne 0$ and $r \ne 0$. Using $\tilde{u} \in \mathcal{K}^*$ and $r \in \mathcal{K}$ we get $\mu \|\tilde{u}_T\| \|r_T\| \le \mu \tilde{u}_N r_N$, as $0 \le \mu \|\tilde{u}_T\| \le \tilde{u}_N$ and $0 \le \|r_T\| \le \mu r_N$ implies that the product of the two lesser members of the inequality must be less than the product of the two greater members of the inequality, and that both products are greater than or equal to zero. The relation $\tilde{u}^T r = 0$ implies $\tilde{u}_N r_N = -\tilde{u}_T^T r_T$ so we get $\mu \|\tilde{u}_T\| \|r_T\| \le -\mu \tilde{u}_T^T r_T$. The Cauchy-Schwarz inequality gives $\mu \|\tilde{u}_T\| \|r_T\| \ge -\mu \tilde{u}_T^T r_T$, so we conclude that

$$\mu \|\tilde{u}_T\| \|r_T\| = -\mu \tilde{u}_T^T r_T. \qquad (13)$$
When $\mu \neq 0$, the vectors $\tilde{u}_T$ and $r_T$ are colinear and have the opposite direction to each other, since the Cauchy-Schwarz inequality holds as an equality. Let us assume that $r \in \text{int}(\mathcal{K})$, then $\|r_T\| < \mu r_N$. Since $u \in \mathcal{K}^*$, we have $\mu\|\tilde{u}_T\|\|r_T\| < \mu r_N \tilde{u}_N = -\mu \tilde{u}_T^\top r_T$, which contradicts (13). Hence, we conclude $r \in \partial\mathcal{K}$ and hence $\|r_T\| = \mu r_N$. Let us assume that $\tilde{u} \in \text{int}(\mathcal{K}^*)$, then $\mu\|\tilde{u}_T\| < \tilde{u}_N$. Since $r \in \mathcal{K}$, $\mu\|\tilde{u}_T\|\|r_T\| < \mu r_N \tilde{u}_N$ contradicts (13). Hence, we conclude $\tilde{u} \in \partial\mathcal{K}^*$. From $\tilde{u} \in \partial\mathcal{K}^*$, we conclude that $u_N = 0$.

**Smooth dynamics with Coulomb friction** For the sake of simplicity, we note by

$$ \tilde{u} = \Phi(u) \tag{14} $$

the De Saxcé change of variable (8) for all contacts $\alpha$. The smooth dynamics with Coulomb friction is then given by

$$ q(t_0) = q_0, \quad v(t_0) = v_0, \tag{15a} $$

$$ \dot{q}(t) = v(t), \tag{15b} $$

$$ M\ddot{v}(t) + Kq(t) + Cv(t) = F_{\text{ext}}(t) + H^\top(q(t))r(t), \tag{15c} $$

$$ u(t) = H(q(t))v(t), \tag{15d} $$

$$ -\Phi(u) \in \mathcal{N}_{\mathcal{K}}(r). \tag{15e} $$

**Remark 1** By substituting (14) in the normal cone inclusion in (15), we obtain for a contact $\alpha$

[Normal/tangent decomposed formulation]

$$ - \begin{pmatrix} u_N^\alpha + \mu \|u_T^\alpha\| \\ u_T^\alpha \end{pmatrix} \in \mathcal{N}_{\mathcal{K}} \begin{pmatrix} r_N^\alpha \\ r_T^\alpha \end{pmatrix}. \tag{16} $$

After some algebraic manipulations, the inclusion (16) is equivalent to the decomposed version of the contact law in the form $-u_N^\alpha \in \mathcal{N}_{\mathbb{R}_+}(r_N)$, for the normal part and $-u_T \in \mathcal{N}_{\mathcal{D}(\mu r_N)}(r_T)$, for the tangential part. Moreau's friction disk is defined as $\mathcal{D}(\mu r_N) = \{r_T \mid \|r_T\| \le \mu r_N\} \subset \mathbb{R}^2$. This decomposed formulation is certainly more common in the historical literature on the subject, see for example (Moreau 1988). However, we advocate the use of the formulation introduced by De Saxcé for several reasons. Firstly, this formulation respects the definition of the Coulomb cone without introducing an artificial decomposition of the cone. In addition, it allows the contact law to be written in the form of a standard variational inequality. This avoids the notion of quasi-variational inequalities, which is difficult to handle from a mathematical and numerical point of view. Indeed, it becomes difficult to use standard mathematical results and algorithms for solving variational inequalities in the quasi-variational case. For example, a variational inequality on a second-order cone is directly equivalent to a complementarity problem on the same cone. For this problem, a rich literature exists, see for example (Facchinei and Pang 2003; Alizadeh and Goldfarb 2003; Andersen et al. 2003; Acary et al. 2024).

## 2.2 Nonsmooth dynamics

If a contact $\alpha$ is closing at time $t_i$ with a negative relative velocity, that is $g_N^\alpha(q(t_i)) = 0$ and $u_N^\alpha(t_i) < 0$, we have an impact at $t_i$. The velocity must jump to satisfy the unilateral constraint immediately after the impact ($g_N^\alpha(q(t_i + \varepsilon)) \ge 0$ for arbitrary small $\varepsilon \ge 0$). The velocity $v$ is usually assumed to be a function of locally bounded variations. The same applies for the relative velocity at contact. With finite dimensional systems, the smooth dynamics (15) are generally insufficient to characterise the solution (that is to say that there are infinitely many solutions after impacts). To close the system at the impact time $t_i$ when $v^+(t_i) \ne v^-(t_i)$, an impact law must be added. In this work, we consider Newton's impact law

$$ u_N^{\alpha,+}(t_i) = -e^\alpha u_N^{\alpha,-}(t_i) \text{ if } g_N^\alpha(q(t_i)) \le 0 \text{ and } u_N^{\alpha,-}(t_i) \le 0, \tag{17} $$

where $e^\alpha \in [0, 1]$ is the Newton coefficient of restitution. To address the case of multiple contacts, the Newton impact law is extended in terms of complementarity by Moreau (1988) as

$$ 0 \le u_N^{\alpha,+}(t_i) + e^\alpha u_N^{\alpha,-}(t_i) \perp p_i^\alpha \ge 0, \text{ if } \alpha \in \bar{\mathcal{I}}^1, \tag{18} $$
where $p_i$ is the reaction impulse at time $t_i$ and $\bar{I}^1$ is the index set defined by

$$ \bar{I}^1 = \{ \alpha \in \mathcal{I} \mid g_N^\alpha(q(t)) \le 0 \text{ and } u_N^{\alpha, -}(t) \le 0 \}. \quad (19) $$

In the following, we will assume that the local variables at contact are collected for all indices belonging to $\bar{I}^1$.

The nonsmooth equations of motion are written in terms of a differential measure $dv$ associated with $v$ and the local impulse measure $di$ as follows

$$ q(t_0) = q_0, \quad v^-(t_0) = v_0, \qquad (20a) $$

$$ \dot{q} = v, \qquad (20b) $$

$$ Mdv + Kqdt + Cvdt = Fdt + H^T(q)di, \qquad (20c) $$

$$ u = H(q)v, \qquad (20d) $$

$$ -\Psi(u^+) \in N_{\mathcal{K}}(di), \qquad (20e) $$

where we have made use of the following function that takes into account the impact law (18):

$$ \Psi^{\alpha}(u^{\alpha}) = u^{\alpha} + \begin{pmatrix} e^{\alpha} u_{N}^{\alpha, -} + \mu \|u_{T}^{\alpha}\| \\ 0 \end{pmatrix}. \qquad (21) $$

**Interpretation of the contact law (20e) in terms of measures** Equation (20e) requires some explanation based on the definition of a normal cone inclusion with measures. By substituting (21) in the normal cone inclusion (20e), we obtain for a contact $\alpha$

$$ - \begin{pmatrix} u_{N}^{\alpha,+} + e^{\alpha} u_{N}^{\alpha,-} + \mu \|u_{T}^{\alpha,+}\| \\ u_{T}^{\alpha,+} \end{pmatrix} \in N_{\mathcal{K}} \begin{pmatrix} di_{N}^{\alpha} \\ di_{T}^{\alpha} \end{pmatrix}. \qquad (22) $$

Almost everywhere with respect to the Lebesgue measure $dt$, we have

$$ -\Psi(u^+) \in N_{\mathcal{K}}(r), \qquad (23) $$

where $r$ is the density of $di$ with respect to $dt$. Since we have $u^- = u^+ = u$ almost everywhere, Equation (23) can be written for a contact $\alpha$ as

$$ \bar{u}^{\alpha} = \begin{pmatrix} (1 + e^{\alpha})u_{N}^{\alpha} + \mu \|u_{T}^{\alpha}\| \\ u_{T}^{\alpha} \end{pmatrix} \in N_{\mathcal{K}^{\alpha}}(r^{\alpha}), \qquad (24) $$

Three cases are possible following the proof of the Lemma 1:

1. The take-off case: $r^\alpha = 0$. In that case, we have $u_N^\alpha \ge 0$.

2. The sticking case: $\bar{u}^\alpha = 0$. In that case, we have $u_T^\alpha = 0$ and $(1+e^\alpha)u_N^\alpha = 0 = u_N^\alpha$ and $r^\alpha \in \mathcal{K}^\alpha$.

3. The sliding case: $r^\alpha \ne 0$ and $\bar{u}^\alpha \ne 0$. Then we have $r^\alpha \in \partial\mathcal{K}^\alpha$ and $\bar{u}^\alpha \in \partial\mathcal{K}^{\alpha, \star}$. This implies $(1+e)u_N^\alpha = 0 = u_N^\alpha$, $\|r_T\| = \mu r_N$ and $\|u_T\|_{r_T} = -\|r_T\|_{u_T}$.

To conclude, we retrieve almost everywhere the Coulomb friction with the Signorini condition at the velocity level.

For the impulsive part, we have at any time $t_i$,

$$ -\Psi(u^+) \in N_{\mathcal{K}}(p_i), \qquad (25) $$

where $p_i$ is the density of $di$ with respect to the Dirac atom at $t_i$, $\delta_{t_i}$. Once again, we have three possible cases for a contact $\alpha$:

1. The take-off case: $p_i^\alpha = 0$. In this case, we have $u_N^{\alpha,+} + eu_N^{\alpha,-} \ge 0$.

2. The sticking case: $\Psi(u^{\alpha,+}) = 0$. In this case, we have $u_T^{\alpha,+} = 0$ and $u_N^{\alpha,+} + eu_N^{\alpha,-} = 0$ and $p_i^\alpha \in \mathcal{K}^\alpha$. The impact law is satisfied since $u_N^{\alpha,+} = -eu_N^{\alpha,-}$.

3. The sliding case: $p_i^\alpha \ne 0$ and $\Psi(u^{\alpha,+}) \ne 0$. Then we have $p_i^\alpha \in \partial\mathcal{K}^\alpha$ and $\Psi(u^{\alpha,+}) \in \partial\mathcal{K}^{\alpha, \star}$. This implies $u_N^{\alpha,+} + eu_N^{\alpha,-} = 0$, $\|p_{i,T}^\alpha\| = \mu p_{i,N}^\alpha$ and $\|u_T^{\alpha,+}\|_{p_{i,T}^\alpha} = -\|p_{i,T}^\alpha\|_{u_T^{\alpha,+}}$.

To summarise, the impact law is satisfied in the sticking and sliding cases and the Coulomb friction is written in terms of impulses and of the right limit of the relative velocity. We will see in the following that these choices have a consequence for the dissipativity of the model.
## 2.3 Energy balance analysis

A detailed analysis of the energy balance for nonsmooth systems can be found in (Leine and Wouw 2008) and in (Brogliato et al. 2013). The energy balance is usually obtained by multiplying the equation of motion by $v^+ + v^-$. After some algebraic manipulations, see (Acary 2016), we obtain

$$2d\mathcal{E} = d(v^T M v) + 2q^T K dq = 2v^T F dt - 2v^T C v dt + (v^+ + v^-)^T H^T(q) di, \quad (26)$$

where the standard definition of the total mechanical energy of the system is

$$\mathcal{E} = \frac{1}{2} v^T M v + \frac{1}{2} q^T K q, \quad (27)$$

and $dq = v^+(t)dt = v^-(t)dt$.

**Energy balance almost everywhere** Almost everywhere with respect to the Lebesgue measure $dt$, the energy balance (26) is

$$\frac{d}{dt} \mathcal{E}(t) = v^T(t) F(t) - v^T(t) C v(t) + v^T(t) H^T(q(t)) r(t), \quad (28)$$

which is the classical energy balance for a smooth dynamical system subjected to some constraints. The term $v^T F$ is the power of the external forces and $-v^T C v$ is the power of the viscous forces. The last term corresponds to the power of the reaction forces. For a contact $\alpha$, it can be decomposed as

$$v^T H^{\alpha, \top}(q) r^\alpha = u^{\alpha, \top} r^\alpha = u_N^\alpha r_N^\alpha + u_T^{\alpha, \top} r_T^\alpha. \quad (29)$$

As we assumed that the constitutive law (23) is satisfied almost everywhere, an examination of the three possible cases for the contact law yields

$$v^T H^{\alpha, \top}(q) r^\alpha = -\mu^\alpha r_N^\alpha \|u_T^\alpha\|, \quad (30)$$

which is the standard dissipated power by a sliding contact with Coulomb friction. Note that the dissipated power is always non-positive, which is consistent with the laws of thermodynamics.

**Energy balance at any time** At any time $t_i$, the energy balance (26) is

$$\mathcal{E}^+(t_i) - \mathcal{E}^-(t_i) = \frac{1}{2} (v^+(t_i) + v^-(t_i))^T H^T (q(t_i)) p_i. \quad (31)$$

For a contact $\alpha$, it can be evaluated as

$$\frac{1}{2} (v^+(t_i) + v^-(t_i))^T H^{\alpha, \top} (q(t_i)) p_i^\alpha = \frac{1}{2} (u^{\alpha, +}(t_i) + u^{\alpha, -}(t_i))^T p_i^\alpha. \quad (32)$$

With the constitutive law given by (23), the dissipated power in (32) vanishes in the take-off case. In the sticking case, the sliding velocity vanishes and we get

$$\begin{aligned} \frac{1}{2} (v^+(t_i) + v^-(t_i))^T H^{\alpha, \top} (q(t_i)) p_i^\alpha &= \frac{1}{2} (u_N^{\alpha, +}(t_i) + u_N^{\alpha, -}(t_i)) p_{N,i}^\alpha + \frac{1}{2} u_T^{\alpha, -}(t_i) p_{T,i}^\alpha, \\ &= \frac{1}{2} (1 - e) u_N^{\alpha, -}(t_i) p_{N,i}^\alpha + \frac{1}{2} u_T^{\alpha, -}(t_i) p_{T,i}^\alpha. \end{aligned} \quad (33)$$

Since $u_N^{\alpha, -}(t_i) \le 0$ and $e \in [0, 1]$, the first term of the dissipated power by an impact is non-positive, but we cannot conclude a priori on the sign of the second term. In the sliding case, we get

$$\begin{aligned} \frac{1}{2} (v^+(t_i) + v^-(t_i))^T H^{\alpha, \top} (q(t_i)) p_i^\alpha &= \frac{1}{2} (1 - e) u_N^{\alpha, -}(t_i) p_{N,i}^\alpha + \frac{1}{2} (u_T^{\alpha, +}(t_i) + u_T^{\alpha, -}(t_i))^T p_{T,i}^\alpha, \\ &= \frac{1}{2} (1 - e) u_N^{\alpha, -}(t_i) p_{N,i}^\alpha - \frac{1}{2} \mu^\alpha p_{N,i}^\alpha \|u_T^{\alpha, +}(t_i)\| + \frac{1}{2} (u_T^{\alpha, -}(t_i))^T p_{T,i}^\alpha. \end{aligned} \quad (34)$$
The first two terms in Equation (34) are non-positive. For the last terms, we cannot conclude a priori on its sign. For instance, if there is a reversal of the sliding direction during the impact, $u_T^{\alpha,+}(t_i) = -u_T^{\alpha,-}(t_i)$, the power of the frictional impulse may generate energy. In other words, the Coulomb law in terms of impulses together with the Newton impact law may be non-dissipative and may not respect the second law of thermodynamics.

This fact is not new and examples where the second law of thermodynamics is not respected are known. The famous Kane (1984) puzzle, see also (Brogliato 1996; Brogliato 1999; Glocker 2013) for details and other examples, shows that some amount of energy can be generated using this type of constitutive law. There does not seem to be any prospect of developing a consistent numerical scheme that dissipates energy when the underlying law does not possess this property.

## 2.4 Frémond's model of frictional impact

To obtain a constitutive law that is consistent with the laws of thermodynamics and the Newton impact law, we propose to use

$$-\Xi\left(\frac{1}{2}(u^{+} + u^{-})\right) \in N_{\mathcal{K}}(di) \qquad (35)$$

with

$$\Xi^{\alpha}(\bar{u}^{\alpha}) = \left( \frac{\bar{u}_{N}^{\alpha} + \frac{1}{2}(e^{\alpha} - 1)u_{N}^{\alpha,-} + \mu^{\alpha} \|\bar{u}_{T}^{\alpha}\|}{\bar{u}_{T}^{\alpha}} \right) \quad \text{and} \quad \bar{u}^{\alpha} = \frac{1}{2}(u^{\alpha,+} + u^{\alpha,-}). \qquad (36)$$

This law is inspired by the energy balance (32) and the work of Frémond (Frémond 1995; Frémond 2001; Frémond 2002), presented most clearly in (Frémond 2017). We will refer to this constitutive law as the Frémond model in the following. We will now demonstrate that the Frémond model always dissipates energy.

**Lemma 2** The Frémond model of contact, impact and friction given by (35) and (36) satisfies the dissipation inequality, that is

$$\frac{1}{2}(v^{+} + v^{-})^{\top}H^{\top}(q)di \le 0. \qquad (37)$$

**Proof 2** Almost everywhere with respect to the Lebesgue measure dt, we have $u^{+} = u^{-} = u$, and the constitutive law (35) reduces to

$$-\left(\frac{1}{2}(1 + e^{\alpha})u_{N}^{\alpha} + \mu^{\alpha}\|u_{T}^{\alpha}\|\right) \in N_{\mathcal{K}^{\alpha}}(r^{\alpha}). \qquad (38)$$

The three possible cases of Equation (38) coincide with those of Equation (23). In other words, for the smooth dynamics, there is no difference between the law in Equation (35) and the standard Coulomb law with the Signorini condition at the velocity level given by Equation (23). At any time $t_i$, we obtain

$$-\Xi(\frac{1}{2}[u^{+}(t_i) + u^{-}(t_i)]) \in N_{\mathcal{K}}(p_i). \qquad (39)$$

To be explicit, let us examine the three possible cases for a contact $\alpha$:

1. The take-off case: $p_i^{\alpha} = 0$. In this case, we have $\frac{1}{2}(u_N^{\alpha,+} + u_N^{\alpha,-}) + \frac{1}{2}(e^{\alpha} - 1)u_N^{\alpha,-} \ge 0$, that is to say $u_N^{\alpha,+} + eu_N^{\alpha,-} \ge 0$.

2. The sticking case: $\Xi(\frac{1}{2}[u^{\alpha,+}(t_i) + u^{\alpha,-}(t_i)]) = 0$. In this case, we have $\frac{1}{2}(u_T^{\alpha,+} + u_T^{\alpha,-}) = 0$, $u_N^{\alpha,+} + eu_N^{\alpha,-} = 0$ and $p_i^{\alpha} \in \mathcal{K}^{\alpha}$. The impact law is satisfied since $u_N^{\alpha,+} = -eu_N^{\alpha,-}$. Note that in the sticking case, we have $u_T^{\alpha,+} = -u_T^{\alpha,-}$, which may take non-zero values. Only the velocity average over the time of impact vanishes.

3. The sliding case: $p_i^{\alpha} \ne 0$ and $\Xi(\frac{1}{2}(u^{\alpha,+}(t_i) + u^{\alpha,-}(t_i))) \ne 0$. Then we have $p_i^{\alpha} \in \partial\mathcal{K}^{\alpha}$ and $\Xi(\frac{1}{2}[u^{\alpha,+}(t_i) + u^{\alpha,-}(t_i)]) \in \partial\mathcal{K}^{\alpha,*}$. This implies $u_N^{\alpha,+} + eu_N^{\alpha,-} = 0$, $\|p_{i,T}^{\alpha}\| = \mu p_{i,N}^{\alpha}$ and $\|\frac{1}{2}(u^{+}(t_i) + u^{-}(t_i))\|p_{i,T}^{\alpha} = -\|p_{i,T}^{\alpha}\|\frac{1}{2}(u^{+}(t_i) + u^{-}(t_i))$.
In the take-off case, the impulse $p_i$ vanishes and hence the dissipated power vanishes.

In the sticking case, we obtain

$$ 
\begin{aligned}
\frac{1}{2}(v^+(t_i) + v^-(t_i))^\top H^{\alpha, \top} (q(t_i)) p_i^\alpha &= \frac{1}{2}(u_N^{\alpha,+}(t_i) + u_N^{\alpha,-}(t_i)) p_{N,i}^\alpha, \\
&= \frac{1}{2}(1 - e)u_N^{\alpha,-}(t_i) p_{N,i}^\alpha \le 0,
\end{aligned} 
\quad (40) $$

since $\frac{1}{2}(u_T^{\alpha,+}(t_i) + u_T^{\alpha,-}(t_i))$ vanishes. Since $u_N^{\alpha,-}(t_i) \le 0$ and $e \in [0, 1]$, the dissipated power by an impact is non-positive.

In the sliding case, we get

$$ 
\begin{aligned}
\frac{1}{2}(v^+(t_i) + v^-(t_i))^\top H^{\alpha, \top} (q(t_i)) p_i^\alpha &= \frac{1}{2}(1 - e)u_N^{\alpha,-}(t_i) p_{N,i}^\alpha + \frac{1}{2}(u_T^{\alpha,+}(t_i) + u_T^{\alpha,-}(t_i))^\top p_{T,i}^\alpha, \\
&= \frac{1}{2}(1 - e)u_N^{\alpha,-}(t_i) p_{N,i}^\alpha - \frac{1}{2}\mu^\alpha p_{N,i}^\alpha \|\frac{1}{2}(u_T^{\alpha,+}(t_i) + u_T^{\alpha,-}(t_i))\| \le 0.
\end{aligned} 
\quad (41) $$

We conclude that the dissipated power related to the reaction forces or impulses is always non-negative. In other words, the model satisfies the laws of thermodynamics at the time of impact. The average of the pre-impact and post-impact velocities $\frac{1}{2}(u^+ + u^-)$ is the natural dual variable of the impulse in the energy balance, or in the principle of virtual power. Choosing a constitutive law that relates this velocity to the impulse though a pseudo-potential of dissipation

$$ -p \in \delta_{\mathcal{K}^*}(\Xi(\frac{1}{2}[u^+ + u^-])) \quad (42) $$

is a way to define a dissipative law (Moreau 1971; Moreau 1973). The convex function $\delta_C$ is the indicator function of a convex set $C$, defined by

$$ \delta_C(x) = \begin{cases} 0 & \text{if } x \in C, \\ \infty & \text{if } x \notin C. \end{cases} \quad (43) $$

**Remark 2** The model proposed by Frémond is not only consistent with the laws of thermodynamics, but has also been validated by experimental tests. In (Cholet 1998; Frémond 2001; Frémond 2002; Frémond 2017), we find experimental results showing that percussions belong to the Coulomb cone. Even more interestingly, these experimental tests also show that the Coulomb relationship between impact and velocity is valid for the averaged velocity $\frac{1}{2}(u^+(t_i) + u^-(t_i))$, but not for the velocity after impact $u^+(t_i)$.

## 2.5 Moreau's impact law with friction

In (Moreau 1988), the impact law (16) is expressed in a more general form including a tangential coefficient of restitution. In our notation and with a variational inequality over the Coulomb cone, it can be formulated as

$$ - \left( \begin{matrix} u_N^{\alpha,+} + e^\alpha u_N^{\alpha,-} + \mu \|u_T^{\alpha,+} + e_T^\alpha u_T^{\alpha,-}\| \\ u_T^{\alpha,+} + e_T^\alpha u_T^{\alpha,-} \end{matrix} \right) \in \mathcal{N}_\mathcal{K} \left( \begin{matrix} d_N^\alpha \\ d_T^\alpha \end{matrix} \right), \quad (44) $$

where $e_T^\alpha$ is the tangential coefficient of restitution for the contact $\alpha$.

A first remark is that the Frémond model presented in the previous section appears to be a special case of the general Moreau model for $e_T = 1$. As we have already pointed out, a number of experiments show that the choice $e_T = 1$ is justified for some impactors. Moreau's primary motivation was to introduce a parameter to reproduce multiple impact phenomena (for example, the rocking block that we treat in Section 5.2).

For Frémond, the motivation behind his model was not related to a tangential restitution coefficient, but instead thermodynamic considerations. The Clausius-Duhem inequality which is closely linked to the positive dissipation postulate, gives rise to a pair of dual variables for impacts: the average of the velocities before and after impact and the percussion. In the pioneering works of Moreau (1970); Moreau (1971); Moreau (1973) and Ziegler (1983), a constitutive law that ensures the dissipation of the model is guaranteed by deriving it from a postulated pseudo-potential of
dissipation between these dual variables. We follow this approach to ensure the dissipation of the model. In this work, we take the pseudo-potential proposed by Frémond without discussing its mechanical relevance in depth, considering only the experiments cited in Remark 2. It is not possible to capture all types of impacts with this model (the model is generally intended for rapid impacts between very stiff bodies), and so it may be necessary to adopt a different pseudo-potential in the event that other phenomena (viscosity, cohesion etc.) play a role in the impact process. It is also possible to choose general laws that use the duality between $\frac{1}{2}(u^+ + u^-)$ and $p$ in the form

$$-p \in \partial\varphi(\frac{1}{2}[u^+ + u^-]) \qquad (45)$$

where $\varphi$ is a lower-semi continuous proper convex function such that $\varphi(0) = 0$ (a dissipation pseudo-potential à la Moreau (1973)). It may be for instance motivated by experimental evidence, or a desire to model more complex material behaviour such as cohesion of the interface (Collins-Craft et al. 2022; Collins-Craft et al. 2024). Note that the Frémond model including Coulomb friction does not derive from a convex potential, due to friction. It remains thermodynamically consistent, but in the general case of non-convex potentials, care must be taken to ensure this property.

Moreau’s impact model in terms of inclusion and measure is related to Brach’s model (Brach 1984; Brach 1998), where a tangential restitution coefficient was also introduced. For an exhaustive review of the various impact models, see (Brogliato 1996; Brogliato 1999; Brogliato 2016). In these monographs, different examples of collisions of two rigid solids are carefully detailed in Chapter 4, showing that for central impacts (as for two spheres for example) the model (44) is dissipative. For other types of impact, such as the impacting stick on a plane we treat in Section 5.1, dissipation is no longer assured.

More recently, the model in (44) has been also studied from the dissipation point of view in (Leine and Wouw 2008; Glocker 2013). The results are given explicitly using the impact dynamics equations that relate velocity jumps to impacts in a condensed form at contact:

$$u^+ = Wp + u^- \qquad (46)$$

where $W = H(q)M^{-1}H^{\top}(q)$ is the Delassus matrix. It should be noted that what we propose in this article does not require any assumptions on the structure of the Delassus matrix. There are two categories of results in the contributions of Glocker (2013); Leine and Wouw (2008). The first class of results concerns the local (isolated) characterisation of dissipation, which in our framework means examining the dissipation at each contact point. For that, Glocker (2013) assumes that the contacts are decoupled to obtain the conditions necessary for positive dissipation. These conditions are $0 \le e \le 1$ and $-1 \le e_T \le 1$. The second type of results concerns the global dissipation of the system, which makes it possible to generate local energy provided that it is compensated for by the dissipation of the rest of the system. Leine and Wouw (2008) first investigated this problem.

In order to have positive global dissipation, two sufficient conditions can be distinguished (see (Leine and Wouw 2008, Section 7.1) for more detailed results): a) all the coefficients of restitution are equal, or sufficiently small (relative to the condition number of the Delassus matrix), or have almost equal values (again, relative to the condition number of the Delassus matrix) and b) the contacts are completely decoupled from each other and from the tangent and normal parts of each contact. These conditions are quite restrictive, however they do not guarantee that the second principle of thermodynamics is valid at every point in the body. Note that we differ with Glocker (2013) who writes: “we do not require the individual contributions $W_i$ to be non-positive as done by so-called energetic impact laws, which ... would in our opinion restrict the energy transfer within the system in a too harsh way.” where $W_i$ is the work of the percussion at contact $i$. While it could certainly be interesting in the future to investigate whether there are other pseudo-potentials of dissipation, and other conditions that can be imposed on the Delassus matrix to guarantee positive dissipation locally with conjugate variables $p_T$ and $\frac{1}{2}(u^+ + u^-)$, this is outside the scope of this paper.

The choice of the Frémond model with $e_T = 1$, which is admittedly not necessarily suitable for every possible impact experiment with any possible material, ensures positive dissipation whatever
the choice of restitution coefficient $e$, friction coefficient $\mu$ and whatever the couplings between the contact points. This is the result of an approach which ensures that the Clausius-Duhem inequality, and hence the second principle of thermodynamics, is respected at all points.

# 3 Moreau-Jean scheme with discrete Frémond's impact law

## 3.1 Flaws in the standard Moreau-Jean scheme with Newton's impact law and Coulomb's friction

In this section we consider systems that are discrete in both space and time, demonstrate the problems with the Moreau-Jean scheme and the standard law, and demonstrate that the Moreau-Jean scheme with the Frémond law avoids these problems.

**Formulation of the Moreau-Jean scheme** The Moreau-Jean scheme (Jean and Moreau 1987; Moreau 1988; Jean 1999) on an interval $(t_k, t_{k+1}]$ of length $h$ is as

$$M(v_{k+1} - v_k) + hKq_{k+\theta} + hCv_{k+\theta} - hF_{k+\theta} = H^{\top}(q_k)p_{k+1}, \quad (47a)$$

$$q_{k+1} = q_k + h v_{k+\theta}, \quad (47b)$$

$$u_{k+1} = H(q_k) v_{k+1}, \quad (47c)$$

$$-\Gamma(u_{k+1}) \in N_{\mathcal{K}}(p_{k+1}) \quad (47d)$$

with $\theta \in [0, 1]$ and where

$$\Gamma^{\alpha}(u_{k+1}^{\alpha}) = u_{k+1}^{\alpha} + \begin{pmatrix} eu_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T,k+1}^{\alpha}\| \\ 0 \\ 0 \end{pmatrix}. \quad (48)$$

The following approximations are considered: $v_{k+1} \approx v(t_{k+1})$, $u_{k+1} \approx u^{+}(t_{k+1})$, and $p_{k+1} \approx di((t_k, t_{k+1}))$. The index set $\bar{\mathcal{I}}^1$ is approximated at each time step by $\bar{\mathcal{I}}^1_k = \{\alpha \in \mathcal{I} \mid g^{\alpha}(q_k) \le 0 \text{ and } u_{N,k}^{\alpha} \le 0\}$. The discrete variables at each contact are collected in the unknown variables $u_{k+1}$ and $p_{k+1}$ for the indices such that $\alpha \in \bar{\mathcal{I}}^1_k$.

**Discrete dissipation analysis** Following (Acary 2016, Lemma 5.1), the discrete-time dissipation equality of the Moreau-Jean scheme (47) over a time-step $[t_k, t_{k+1}]$ is given by

$$\Delta\mathcal{E} - W_{k+1}^{\text{ext}} - W_{k+1}^{\text{damping}} = (\frac{1}{2} - \theta) [\|v_{k+1} - v_k\|_M^2 + \|q_{k+1} - q_k\|_K^2] + u_{k+\theta}^{\top} p_{k+1}, \quad (49)$$

and the discrete approximation of the work done by the external forces within the step by

$$W_{k+1}^{\text{ext}} = h v_{k+\theta}^{\top} F_{k+\theta} \approx \int_{t_k}^{t_{k+1}} F v dt, \quad (50)$$

and the discrete approximation of the work done by the damping term by

$$W_{k+1}^{\text{damping}} = -h v_{k+\theta}^{\top} C v_{k+\theta} \approx - \int_{t_k}^{t_{k+1}} v^{\top} C v dt. \quad (51)$$

To ensure that the scheme dissipates energy, the first condition is $\theta \ge \frac{1}{2}$. The second condition is related to the sign of the discrete work of the reaction impulse given by

$$u_{k+\theta}^{\top} p_{k+1} = u_{N,k+\theta}^{\top} p_{N,k+1} + u_{T,k+\theta}^{\top} p_{T,k+1}. \quad (52)$$

For a contact $\alpha$, the normal term $u_{N,k+\theta}^{\alpha,\top} p_{N,k+1}^{\alpha}$ vanishes in the take-off case, and is given by

$$u_{N,k+\theta}^{\alpha,\top} p_{N,k+1}^{\alpha} = [1 - \theta(1 + e^{\alpha})] u_{N,k}^{\alpha} p_{N,k+1}^{\alpha}, \quad (53)$$

in the sliding and sticking cases. With the condition

$$\theta \le \frac{1}{1 + e^{\alpha}} \le 1, \text{ for all } \alpha \in \mathcal{I}. \quad (54)$$

we conclude that $u_{N,k+\theta}^{\alpha,\top} p_{N,k+1}^{\alpha} \le 0$.
**Remark 3** In the case of $e = 1$, we have the very restrictive condition that $\theta = \frac{1}{2}$, which means that in the case of perfectly elastic contacts, we cannot impose artificial dissipation via the numerical scheme. In the case that additional dissipation is required to stabilise the integration, it must be imposed by adding viscous behaviour in the bulk, or by recourse to generalised nonsmooth integration schemes (Brüls et al. 2014).

For the frictional term, $u_{T,k+\theta}^{\alpha,\top}p_{T,k+1}^{\alpha}$, we get

$$u_{T,k+\theta}^{\alpha,\top}p_{T,k+1}^{\alpha} = \theta u_{T,k+1}^{\alpha,\top}p_{T,k+1}^{\alpha} + (1-\theta)u_{T,k}^{\alpha,\top}p_{T,k+1}^{\alpha}. \qquad (55)$$

The first term is non-positive but we cannot draw any conclusions about the sign of the second term. We will demonstrate that on a simple example this term can be positive and that the scheme will generate energy.

This fact is not particularly surprising. If we accept that the scheme is consistent and converges towards the continuous-time model, then we must converge towards a model that can generate energy. In practice, energy can be generated even in the case where we have reduced the time steps to obtain a more accurate solution.

**Remark 4** A possible way to dissipate energy is to cancel the second term in (55) using $\theta = 1$. In this case, the condition (54) will imply that $e^{\alpha} = 0$ for all $\alpha$. This result coincides with the results obtained by Leine and Wouw (2008); Glocker (2013), who demonstrated that if all restitution coefficients are equal, the system will be globally dissipative. Hence, in this special case, the system is also locally dissipative at all points.

## 3.2 Principles of the proposed scheme

The following proposal for a new time-stepping scheme differs in the way the contact law is discretised to obtain contact dissipation that is always positive. The Moreau-Jean based scheme with the Frémond impact law is defined as follows:

$$M(v_{k+1} - v_k) + hKq_{k+\theta} + hCv_{k+\theta} - hF_{k+\theta} = H^{\top}(q_k)p_{k+1}, \qquad (56a)$$

$$q_{k+1} = q_k + h v_{k+\theta}, \quad u_{k+1} = H(q_k) v_{k+1}, \qquad (56b)$$

$$-\Theta(u_{k+\theta}) \in \mathcal{N}_{\mathcal{K}}(p_{k+1}) \qquad (56c)$$

with the function $\Theta$ as

$$\Theta^{\alpha}(u^{\alpha}) = u^{\alpha} + \begin{pmatrix} (\theta(1+e^{\alpha}) - 1)u_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T}^{\alpha}\| \\ 0 \\ 0 \end{pmatrix}. \qquad (57)$$

Let us enumerate what is solved in each of the three possible cases of the contact law:

1. The take-off case: $p_{k+1}^{\alpha} = 0$. In this case, we have

$$u_{N,k+\theta}^{\alpha} + (\theta(1+e^{\alpha}) - 1)u_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T,k+\theta}^{\alpha}\| \geq \mu^{\alpha} \|u_{T,k+\theta}^{\alpha}\|, \qquad (58)$$

which simplifies to

$$u_{N,k+1}^{\alpha} + e^{\alpha} u_{N,k}^{\alpha} \geq 0. \qquad (59)$$

2. The sticking case: $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) = 0$. In this case, we have $u_{T,k+\theta}^{\alpha} = 0$ and

$$u_{N,k+\theta}^{\alpha} + (\theta(1+e^{\alpha}) - 1)u_{N,k}^{\alpha} = 0, \qquad (60)$$

which simplifies to

$$u_{N,k+1}^{\alpha} = -e^{\alpha} u_{N,k}^{\alpha}. \qquad (61)$$
3. The sliding case: $p_{k+\theta}^{\alpha} \neq 0$ and $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) \neq 0$. Then we have $p_{k+\theta}^{\alpha} \in \partial\mathcal{K}^{\alpha}$ and $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) \in \partial\mathcal{K}^{\alpha, \star}$. This implies

$$u_{N,k+\theta}^{\alpha} + (\theta(1+e^{\alpha}) - 1)u_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T,k+\theta}^{\alpha}\| = \mu^{\alpha} \|u_{T,k+\theta}^{\alpha}\|, \quad (62)$$

which simplifies to

$$u_{N,k+\theta}^{\alpha} + (\theta(1+e^{\alpha}) - 1)u_{N,k}^{\alpha} = 0, \quad (63)$$

and then

$$u_{N,k+1}^{\alpha} = -e^{\alpha} u_{N,k}^{\alpha}. \quad (64)$$

For the frictional part, we have

$$\|p_{T,k+1}^{\alpha}\| = \mu p_{N,k+1}^{\alpha} \quad \text{and} \quad \|u_{T,k+\theta}^{\alpha}\| \|p_{T,k+1}^{\alpha}\| = -\|p_{T,k+1}^{\alpha}\| \|u_{T,k+\theta}^{\alpha}\|. \quad (65)$$

To summarise, the Newton impact law is defined in the same way as in the Moreau–Jean time-stepping scheme. For the normal part of the law, the pre-impact velocity is the velocity at the beginning of the time-step $u_{N,k}^{\alpha}$ and the post-impact velocity is the velocity at the end of the time step $u_{N,k+1}^{\alpha}$. However, in the new Moreau–Jean scheme with Frémond impact law time-stepping scheme the Coulomb friction law is written with the average velocity $u_{T,k+\theta}^{\alpha}$, rather than the velocity at the end of the time-step $u_{T,k+1}^{\alpha}$, as is the case in the classical Moreau–Jean scheme.

When $\theta = 1/2$, the function $\Theta^{\alpha}$ simplifies to

$$\Theta^{\alpha}(u^{\alpha}) = \left( u_{N}^{\alpha} + \frac{1}{2}(e^{\alpha} - 1)u_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T}^{\alpha}\| \right) \quad (66)$$

and the discrete contact law is

$$-\Theta(\frac{1}{2}[u_{k+1} + u_k]) \in N_{\mathcal{K}}(p_{k+1}) \quad (67)$$

We can observe that the expression (66) is very similar to the expression of the function $\Xi$ in (36). In practice, if the scheme converges, we can expect to approximate the Frémond impact law. When $\theta = 1$, the function $\Theta^{\alpha}$ is

$$\Theta^{\alpha}(u^{\alpha}) = \left( u_{N}^{\alpha} + e^{\alpha} u_{N,k}^{\alpha} + \mu^{\alpha} \|u_{T}^{\alpha}\| \right) = \Gamma(u^{\alpha}), \quad (68)$$

and the discrete contact law is

$$-\Theta(u_{k+1}) = -\Gamma(u_{k+1}) \in N_{\mathcal{K}}(p_{k+1}). \quad (69)$$

We observe that the scheme is identical to the Moreau–Jean time-stepping scheme for the case of $\theta = 1$.

### 3.3 Discrete dissipation properties and energy balance

Let us now give a result concerning the dissipation of this new Moreau–Jean scheme with the Frémond impact law.

**Proposition 1**
[Energy dissipation]

The new Moreau–Jean scheme with the Frémond impact law dissipates energy in the sense that

$$\mathcal{E}(t_{k+1}) - \mathcal{E}(t_k) \leq W_{k+1}^{\text{ext}} + W_{k+1}^{\text{damping}}, \quad (70)$$

if

$$\frac{1}{2} \leq \theta \leq \frac{1}{1 + \bar{e}} \leq 1. \quad (71)$$

where $\bar{e} = \max e^{\alpha}, \alpha \in \mathcal{I}$.
In other words, provided that Equation (71) is satisfied, the variation of the total mechanical energy of the system is always less than or equal to the energy supplied by the external and damping forces.

**Proof 3** The discrete-time dissipation equality of the new Moreau-Jean scheme (56) over a time-step $[t_k, t_{k+1}]$ is also given by (49). Hence, we have

$$ (\frac{1}{2} - \theta)[\|v_{k+1} - v_k\|_M^2 + \|q_{k+1} - q_k\|_K^2] \le 0, \text{ if and only if } \theta \ge \frac{1}{2}. \quad (72) $$

All that remains to prove is that $u_{k+\theta}^{\alpha, \top} p_{k+1} \le 0$. Let us examine the three possible cases for a contact $\alpha$.

1. Take-off case: $p_{k+1}^{\alpha} = 0$. In that case, we have $u_{k+\theta}^{\alpha, \top} p_{k+1}^{\alpha} = 0$.

2. Sticking case: $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) = 0$. In that case, we have $u_{T,k+\theta}^{\alpha} = 0$ and then $u_{k+\theta}^{\alpha, \top} p_{k+1}^{\alpha} = u_{N,k+\theta}^{\alpha} p_{N,k+1}^{\alpha}$. From (60), we conclude that $u_{N,k+\theta}^{\alpha} \le 0$ since $\theta \le \frac{1}{1+\bar{e}}$ and $u_{N,k}^{\alpha} \le 0$ and then $u_{k+\theta}^{\alpha, \top} p_{k+1}^{\alpha} \le 0$ since $p_{N,k+1}^{\alpha} \ge 0$.

3. Sliding case: $p_{k+\theta}^{\alpha} \ne 0$ and $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) \ne 0$. Then we have $p_{k+\theta}^{\alpha} \in \partial K^{\alpha}$ and $\Theta^{\alpha}(u_{k+\theta}^{\alpha}) \in \partial K^{\alpha, \star}$. From (63), we conclude that $u_{N,k+\theta}^{\alpha} \le 0$ since $\theta \le \frac{1}{1+\bar{e}}$ and $u_{N,k}^{\alpha} \le 0$, and then $u_{N,k+\theta}^{\alpha} p_{N,k+1}^{\alpha} \le 0$. For the frictional part, we have from (65):

$$ u_{T,k+\theta}^{\alpha, \top} p_{T,k+1}^{\alpha} = -\mu^{\alpha} p_{N,k+1}^{\alpha} \|u_{T,k+\theta}^{\alpha}\| \le 0. \quad (73) $$

Given that we demonstrated that the Frémond impact law is dissipative in continuous time for any physically meaningful values of $e$ and $\mu$, it is once again no surprise that an approximation of this law using a Moreau-Jean scheme is also dissipative when $\theta \ge \frac{1}{2}$. From a practical perspective, we are thus able to conduct numerical integrations of complex systems with multiple impacts and friction (such as arise in simulating granular media or finite-element-discretised bodies undergoing contact) with guarantees on the stability of the solver, increasing both the speed and accuracy of the resolution.

**Corollary 1** In the conservative case with friction, or in the conservative case without friction, that is $e^{\alpha} = 1$ for all $\alpha$ and all contacts are in the sticking or take-off modes, the new Moreau-Jean scheme with the Frémond impact law is an energy conserving scheme in the sense that

$$ \mathcal{E}(t_{k+1}) - \mathcal{E}(t_k) = W_{k+1}^{\text{ext}} + W_{k+1}^{\text{damping}}, \quad (74) $$

provided that $\theta = 1/2$.

# 4 Implementation details in Siconos

We will give some details of the implementation of the numerical solver as it is implemented in Siconos (Acary et al. 2019), where it is possible to solve the systems (47) and (56) directly. If the inverse of the mass matrix is simple (i.e. it is a diagonal or block diagonal matrix), we can also decide to reduce the system to the contact variable, as we do in the following. We will ultimately arrive at a second order cone complementarity problem for each scheme. Siconos/numerics implements a variety of solvers for second order cone complementarity problems, specifically for those arising in non-associative plasticity (Acary et al. 2023) or friction problems (Acary et al. 2018) defined by the system

$$ \begin{cases} w = Ar + b, \\ \tilde{w} = \Phi(w), \\ -\tilde{w} \in N_{K}(r), \end{cases} \quad (75) $$

where $w$ and $r$ are the vectors that solve the system, $A$ is a matrix, $b$ is a vector and $\Phi$ is a function that transforms $w$ to the vector $\tilde{w}$ that has the complementarity relationship with $r$. In the following, $\Phi$ is the De Saxcé velocity function. Siconos has a built-in frictional contact problem type that automatically applies the De Saxcé function, which we exploit in our numerical solutions. Among the available solvers for this type of system, we may cite a) projection-based
solvers for variational inequalities, b) iterative solvers such as block projected Gauss-Seidel, c) semi-smooth Newton methods, and d) interior point methods (Acary et al. 2024). The solver finds the values of r and w that resolve the system to a given tolerance, using the stopping criteria $||r - \text{proj}_{\mathcal{K}}(r - \Phi(Ar + b))||_2/(1 + ||b||_2) < \epsilon$ with a user error tolerance $\epsilon$ set to $10^{-10}$ (unless otherwise specified), i.e. we resolve our systems with very small error. The point $\text{proj}_{\mathcal{K}}(r - \Phi(Ar + b))$ is the closest point $x \in \mathcal{K}$ to $r - \Phi(Ar + b)$ in the sense of the Euclidean distance. We will exploit the nonsmooth Gauss-Seidel solver (Jourdan et al. 1998) with an internal solver based on hybrid semi-smooth Newton with a Goldstein-Price line search (Acary et al. 2018) to resolve our systems, as it is a very robust solver.

## 4.1 Implementation of the Moreau-Jean scheme

To arrive at the appropriate cone complementarity system for the classical Moreau-Jean scheme, we expand the first line of (56a) to obtain

$$ (M+h\theta C+h^2\theta^2 K)v_{k+1}=(M-h(1-\theta)C-h^2\theta(1-\theta)K)v_k-hKq_k+hF_{k+\theta}+H^\top(q_k)p_{k+1}. \quad (76) $$

Then, we construct the matrix $\hat{M} = M + h\theta C + h^2\theta^2 K$ and obtain an expression for the relative velocities given by

$$ u_{k+1}=H(q_k)\hat{M}^{-1}[(M-h(1-\theta)C-h^2\theta(1-\theta)K)v_k-hKq_k+hF_{k+\theta}]+H(q_k)\hat{M}^{-1}H^\top(q_k)p_{k+1}. \quad (77) $$

We will denote the discrete Delassus matrix by $W = H(q_k)\hat{M}^{-1}H^\top(q_k)$, and a constant term by

$$ y = H(q_k)\hat{M}^{-1}[(M - h(1 - \theta)C - h^2\theta(1 - \theta)K)v_k - hKq_k + hF_{k+\theta}]. \quad (78) $$

We obtain a second order cone complementarity problem that must be solved at each time step:

$$ \begin{cases} u_{k+1} = Wp_{k+1} + y, \\ -\Gamma(u_{k+1}) \in N_{\mathcal{K}}(p_{k+1}). \end{cases} \quad (79) $$

Now, considering the generic system (75), we set $A = W$, $\Phi = \Gamma$,

$$ b^{\alpha} = y^{\alpha} + \begin{bmatrix} eu_{N,k}^{\alpha} \\ 0 \\ 0 \end{bmatrix} \quad (80) $$

$r$ is the vector of percussions over the time step $p_{k+1}$ and $w$ the vector of De Saécé relative velocities at the end of the time step $u_{k+1}$, respectively. The solutions $r$ and $w$ will be resolved by the solver so that they fulfill Equation (79).

## 4.2 Implementation of the new Moreau-Jean scheme

For the new scheme, we expand the first line of Equation (56) to obtain the same expression as Equation (76). Then, adding $(M + h\theta C + h^2\theta^2 K)v_k$ to each side and making the appropriate rearrangements and simplifications leads to

$$ (M + h\theta C + h^2\theta^2 K)v_{k+\theta} = Mv_k - h\theta Kq_k + h\theta F_{k+\theta} + \theta H^\top(q_k)p_{k+1}. \quad (81) $$

Once again constructing the matrix $\hat{M} = M + h\theta C + h^2\theta^2 K$, inverting it and multiplying by $H(q_k)$, we obtain an expression for the relative velocities:

$$ u_{k+\theta} = H(q_k)\hat{M}^{-1}(Mv_k - h\theta Kq_k + h\theta F_{k+\theta}) + \theta H(q_k)\hat{M}^{-1}H^\top(q_k)p_{k+1}. \quad (82) $$

By denoting the discrete Delassus matrix $W = \theta H(q_k)\hat{M}^{-1}H^\top(q_k)$, and the constant term $y = H(q_k)\hat{M}^{-1}(Mv_k - h\theta Kq_k + h\theta F_{k+\theta})$, the second order cone complementarity problem

$$ \begin{cases} u_{k+\theta} = Wp_{k+1} + y, \\ -\Theta(u_{k+\theta}) \in N_{\mathcal{K}}(p_{k+1}), \end{cases} \quad (83) $$
has to be solved at each time step. We set $A = W$,

$$b^{\alpha} = y^{\alpha} + \begin{bmatrix} (\theta(1 + e^{\alpha}) - 1)u_{N,k}^{\alpha} \\ 0 \\ 0 \end{bmatrix} \qquad (84)$$

and find the vectors $r$ and $w$ (respectively the vectors of percussions over the time step $p_{k+1}$ and the weighted average local relative velocities $u_{k+\theta}$) that provide us with a solution to Equation (83).

# 5 Numerical illustrations

Here, we demonstrate on several different example systems the difference between the results obtained using the classical Moreau-Jean scheme with discrete Moreau's impact law and Coulomb's friction, against those obtained using the new proposed scheme with discrete Moreau's impact law and Coulomb's friction, for systems described in continuous time by Equation (20).

## 5.1 Impacting stick

The first example is a rigid bar impacting a fixed, rigid obstacle. This very simple example shows the errors in calculating dissipation in discrete time using the classical Moreau-Jean scheme. We consider a bar of length $l$ and mass $m$ subjected to gravity. Its centre of mass $G$ is given by the coordinates $(x, y)$ and the angle of the bar $\gamma$ complete the generalised coordinates so that $q = (x, y, \gamma)$, see Figure 1. For simplicity's sake, the contact point $C$ is considered to be on the

<figure><img src="image_2.png" /><figcaption>Figure 1 Geometry of the impacting stick of length $l$ with a centroid at $G$ and making contact with the ground at point $C$, at an angle of $\gamma$.</figcaption></figure>

middle fibre of the bar. The numerical values (in dimensionless units) are as follows: $m = 1$, $l = 1$, $g = 10$, $e = 1.0$, $\mu = 0.01$ and the initial conditions are $q_0 = (\frac{1}{2}l \sin(\pi/4), \frac{1}{2}l \cos(\pi/4) + 0.01, \pi/4)$ and $v_0 = (-0.5, 0.1, 0.1)$. The simulation is performed with a time step $h = 10^{-4}$, $\theta = 0.5$ and a simulation interval $[0, 0.2]$ to observe the first two impacts in detail, using both the classical and new schemes. The results are displayed in Figures 2 and 3. It is observed that the energy dissipated by friction can have the wrong sign for the Moreau-Jean scheme with Newton law at first impact, and therefore generate energy, which is not the case for the proposed Moreau-Jean scheme with the Frémond law. We also observe here that while the normal percussions at impact occur, they do no net work on the system due to the choice of $e = 1$.

## 5.2 A rocking block

We now consider a rigid block of mass $m$, length $b$ and height $l$, subjected to gravity. The centroid of the block is located at $G(x, y)$, and the generalised coordinates are completed by the angle to the horizontal plane $\gamma$ so that $q = (x, y, \gamma)$. It is assumed that the block makes contact with the rigid ground at its two corner points $C_1$ and $C_2$. The geometry is depicted in Figure 4. The numerical values of the system (once again in dimensionless units) are $b = 1$, $l = 1$, $m = 1$, $g = 10$, $e = 1.0$, $\mu = 0.1$. The initial conditions are $q_0 = (0, 0.6, 0)$ and $v_0 = (0, -0.2, 1)$. The simulation is performed with a time step $h = 10^{-4}$, $\theta = 0.5$ and a simulation interval $[0, 1]$, using both the classical and new schemes. The results are displayed in Figure 5 and Figure 6. As can be clearly seen, the classical Moreau-Jean scheme with the Newton law creates energy at the second, third, sixth and eighth impacts, and in fact the dissipated energy can become positive. The global effect of this can be seen in the total energy curve in Figure 5(f) where the total energy has visible
<figure><img src="image_3.png" /><figcaption>**Figure 2** Results of the simulation of the impacting stick, using the standard Moreau–Jean scheme, with $\theta = 0.5$, $h = 10^{-4}$. (a) The generalised coordinates of the system. (b) The generalised velocities of the system. (c) The local relative normal and tangential velocities. (d) The normal and tangential contact impulses. (e) The work done by the reaction impulses. (f) The total energetic quantities of the system.</figcaption></figure>

"kinks" as a sequence of impacts first dissipates then creates energy. By comparison, Figure 6 demonstrates clearly the difference with the system being purely dissipative, and each impact either doing no dissipative work (as $e = 1$, so the normal part of the impact is purely elastic), or doing negative work due to frictional sliding. The total energy curve clearly decreases at each impact, as would be expected. A comparison of the velocities and positions of the system's centre of mass presented in the subfigures (a) and (b) of each figure demonstrate plainly that the creation of energy in the Moreau–Jean scheme with Newton law has substantial effects on the system trajectory.

## 5.3 A sliding block

In this example we consider a linear elastic block that slides in constant contact with a rigid surface, with the geometry as shown in Figure 7. The block is taken to be made of Poly(methyl methacrylate) (PMMA), and to be linearly elastic. We take the linear elasticity parameters from (Berman et al. 2020). We thus have the Young's modulus $E = 5.75 \times 10^3$ MPa, the Poisson's ratio $\nu = 0.358$ and the density $\rho = 1.17 \times 10^{-3}$ g/mm$^3$. We choose a value of the restitution coefficient $e$ of 0.0 and friction coefficient $\mu = 0.5$, and apply these values to all contact nodes.
We discretise the system in space using the finite element method, generating the mesh
<figure><img src="image_4.png" /><figcaption>**Figure 3** Results of the simulation of the impacting stick, using the new Moreau-Jean scheme with the Frémond law, with $\theta = 0.5$, $h = 10^{-4}$. (a) The generalised coordinates of the system. (b) The generalised velocities of the system. (c) The local relative normal and tangential velocities. (d) The normal and tangential contact impulses. (e) The work done by the reaction impulses. (f) The total energetic quantities of the system.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4** A rocking block of width $b$ and height $h$, with a centroid at $G$ and making contact with the ground at the points $C_1$ and $C_2$, at an angle of $\gamma$.</figcaption></figure>

with Gmsh 4.11.1 (Geuzaine and Remacle 2009) and setting a characteristic length of 1 mm throughout, then reading it in to Python 3.12.1 with meshio 5.3.5 (Schlömer 2024). We use linear triangle (T3) elements with a single Gauss point and a consistent mass matrix, choose the plane stress simplifying assumption, and work in the linear finite element framework (that is to say that the finite element matrices are calculated once with respect to the undeformed
<figure><img src="image_6.png" /><figcaption>**Figure 5** Results of the simulation of the rocking block, using the standard Moreau–Jean scheme, with $\theta = 0.5$, $h = 10^{-4}$. (a) The generalised coordinates of the system. (b) The generalised velocities of the system. (c) The local relative normal and tangential velocities. (d) The normal and tangential contact percussions. (e) The work done by the contact impulses. (f) The total energetic quantities of the system.</figcaption></figure>

system configuration (shown in Figure 8(a)) and remain constant for the entire simulation). The generalised coordinates correspond to the displacements in the *x* and *y* directions numbered node-by-node, measured relative to the nodes’ positions in the undeformed configuration, so that $q = (x_0, y_0, x_1, y_1, \dots, x_N, y_N)$. The system is first subjected to an applied vertical displacement $d = -0.005$ mm along the top boundary and the initial displacements are resolved elastically with the no slip and no interpenetration conditions enforced at the contact boundary. During the simulation, the block is subjected to its own self-weight (with gravity set to be $-10^{-3}$ mm/ms²), the vertical displacement on the top boundary is held constant and an oscillating traction $\sigma_s$ on the left boundary given by $\sigma_s = 2 \operatorname{sgn}[\sin(4\pi t)]$. We simulate the system for 1 ms, with time steps of $10^{-4}$ ms, so as to ensure a very fine resolution of the system. In this case, as an external force is applied to the system, we use Equation (50) to calculate the work input into the system. When comparing Figure 9 and Figure 10, we can see very clearly that the new scheme completely avoids generating positive tangential work, aligning with the analytical proofs above. However, this is not the case for the original Moreau–Jean scheme. This may seem surprising in light of the results obtained in (Glocker 2013), where setting all $e = 0$ should result in an energetically consistent system. The resolution of this seeming paradox is that the results of Glocker hold only in continuous time, and the discretisation with the Moreau–Jean scheme does not preserve the
<figure><img src="image_7.png" /><figcaption>**Figure 6** Results of the simulation of the rocking, using the new Moreau-Jean scheme with the Frémond law, with $\theta = 0.5$, $h = 10^{-4}$. (a) The generalised coordinates of the system. (b) The generalised velocities of the system. (c) The local relative normal and tangential velocities. (d) The normal and tangential contact percussions. (e) The work done by the contact impulses. (f) The total energetic quantities of the system.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7** A sliding block under self-weight due to a gravity and an applied sliding force $F_s$. The height $h$ is 25 mm and the breadth $b$ is 40 mm. The thickness is set to 15 mm.</figcaption></figure>

energetic invariants of the system. Nevertheless, in this case the rest of the system is dissipative enough to overcome the energy created by the Moreau-Jean scheme, meaning that the global system remains stable and the sum of the total work and dissipation tracks the total energy closely in both cases.

It is also worth remarking on the seeming appearance of Schallamach-type waves (Schallamach 1971) (namely true Schallamach waves, separation pulses and slip pulses) along the bottom surface of the block, visible in Figure 8(b). Previous studies (Barquins 1985; Viswanathan and
<figure><img src="image_9.png" /><figcaption>**Figure 8** Mesh of the sliding block. (a) The block mesh in the undeformed configuration. (b) The block mesh in the deformed configuration at $t = 0.3$, using the Moreau–Jean scheme with Frémond law, with a 100× warp factor applied to the mesh to exaggerate the displacements, which are also shown via the colour bar. Both figures are generated by ParaView (Ahrens et al. 2005).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9** Results of the simulation using the classical Moreau–Jean method with $\theta = 0.5$ and $h = 10^{-4}$ ms. (a) The individual energetic components of the simulation. (b) The sum of the work input to the system with the dissipation, compared to the total energy of the system. (c) The maximum value of the positive part of the tangential work $w_T$ over all contact nodes at each time step. (d) The minimum value of the negative part of the tangential work $w_T$ over all contact nodes at each time step.</figcaption></figure>

Chandrasekar 2022) have emphasised the role played by adhesion, viscoelastic and temperature effects in this phenomenon. While we have not included it in the simulations in this paper, our
<figure><img src="image_11.png" /><figcaption>**Figure 10** Results of the simulation using the new Moreau-Jean method with the Frémond law and with $\theta = 0.5$ and $h = 10^{-4}$ ms. (a) The individual energetic components of the simulation. (b) The sum of the work input to the system with the dissipation, compared to the total energy of the system. c) The maximum value of the positive part of the tangential work $w_T$ over all contact nodes at each time step. (d) The minimum value of the negative part of the tangential work $w_T$ over all contact nodes at each time step.</figcaption></figure>

method includes the possibility of linear viscosity in the damping matrix C, and setting $e = 0$ (a perfectly inelastic contact), as we have in these simulations, we can at least treat the limiting case of zero adhesive force. Further, experimental evidence suggests that propagating frictional ruptures have a complex spatial structure (Berman et al. 2020), suggesting the importance of rate-dependent friction. As this article is dedicated to the description of a theoretical and numerical model, we will not explore this subject any further, beyond noting that having a consistently dissipative numerical method like the new scheme would seem to be a precondition to a robust numerical study of these phenomena, which have more typically been explored theoretically and experimentally.

## 5.4 Impact on a masonry structure

In this example we consider a more practical example of a rockfall impacting on a masonry wall, with the geometry shown in Figure 11. The wall has ten layers of blocks, and is buttressed by six buttresses, composed of the same type of blocks. The wall rests on flat ground, and in front of the wall is a slope at 45°. Each block has a width of 1 m, a depth of 2 m and a height of 1 m. A rock of
approximately 5 m diameter is created with an irregular polyhedral shape. The blocks and the rock are taken to be rigid bodies, subject to gravity of 9.81 m/s². All the bodies are taken to have a density of 2300 kg/m³. The Newton coefficient of restitution for all normal impacts is taken as e = 0.2 and the coefficient of friction as μ = 0.6.

We do not detail here the formulation and the time-integration of three-dimensional rigid solids with finite rotations. The choice that has been made in Siconos is to use unit quaternions to parameterise rotations for computational efficiency. The Euler equations are written in the body frame in order to obtain a constant inertia operator. As a consequence, we can no longer say $\dot{q}(t) = v(t)$, but instead require a transformation depending on the generalised coordinate to be applied i.e. $\dot{q}(t) = T(q)v(t)$, where $T(q)$ is a matrix that acts to transform the time derivative of the quaternions into conventional angular components in the body-fixed frame. For more details on the formulation in Siconos, we refer the reader to Akhadkar et al. (2018). The integration is based on the mid-point method, with gyroscopic forces implicitly taken into account. The integration of the quaternion is done with a Lie group method (Brüls and Cardona 2010). A global Newton method solves these equations to an accuracy of $10^{-12}$. These methods, which are rather expensive to compute with respect to explicit schemes, are stable, conserve the unit quaternions and enable accurate integration that takes account of the very different rotational inertias along the axes (Holm et al. 1984; Hairer et al. 2003).

To obtain sufficiently accurate results, the simulation is run for 4 s, starting at $t = 0$, with time steps of $h = 5 \times 10^{-4}$ and $\theta = 0.5$. The rock is created at $t = 0.01$, its initial position is at the coordinates (30, 0, 35) m and it has an initial translational velocity of (-10, 0, -35) m/s and an angular velocity of (0.5, 10, 0.1) s⁻¹. In the first instance, we highlight that we are pushing the system outside the range for which the proofs above are valid, i.e. rather than a linear Lagrangian system we have a nonlinear Newton–Euler system. However, in practical terms, this does not prevent the model being implemented as a reasonable approximation. In this case, we track the total work done by the external and gyroscopic forces of the system over each increment of time, given by

$$W_{k+1} = v_{k+\theta}^{\top} (F_{k+\theta}^{\text{ext}} + F_{k+\theta}^{\text{gyr}}) \approx \int_{t_k}^{t_{k+1}} \frac{1}{2} (v^{+}(t) + v^{-}(t))^{\top} (F^{\text{ext}}(t) + F^{\text{gyr}}(t)) dt. \quad (85)$$

The system is resolved with Siconos' three dimensional nonsmooth Gauss–Seidel frictional solver, with the tolerance set to $\epsilon = 10^{-4}$, and the contact detection resolved with Siconos' internal version of Bullet (Coumans and Bai 2021). For a full list of the various numerical tolerances implemented in the system, we refer readers to the simulation file (Acary n.d.). The results of the simulation are output at every time step. In contrast to the previous examples, we see that in this case, the proposed scheme does on occasion generate positive tangential work. This is a direct consequence of numerical error related to the relatively large time step and solver tolerance chosen. With a tighter numerical tolerance, the positive part of the tangential work goes towards zero, and in any case the magnitude of the positive part remains three orders of magnitude less than that seen in the classical Moreau–Jean method. We can compare more closely the total dissipation of each method in Figure 14: We can see that where we have clear corresponding impacts in the two systems, the proposed method clearly dissipates more energy, as a consequence of having smaller positive tangential work terms. As a result, the system has less kinetic energy, and so there are fewer impacts with noticeable dissipation than in the corresponding Moreau–Jean system. In particular, at about 3.1 seconds into the simulation, the increased energy in the Moreau–Jean system causes a number of impacts to occur without the new scheme having obviously corresponding impacts (however, the solution with the proposed scheme does have a small accumulation of impacts near the end without an obvious analogue in the Moreau–Jean scheme). Ultimately, the two systems end up having dissipated similar amounts of energy, although the difference in magnitude of each impacting event causes the systems to follow noticeably different paths to get there. Given small differences in the dissipation at contact can induce meaningful changes in the propagation mechanism (e.g. sliding or rolling) and thus the propagation distance of a rockfall event (Bourrier and Acary 2022), we can expect that in practice application of the Moreau–Jean–Frémond scheme may result in substantially different
<figure><img src="image_12.png" /><figcaption>**Figure 11** Impact of a rock on a masonry structure, at the point of creation of the rock, just after impact, and at the end of the simulation.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12** Results of the simulation using the classical Moreau–Jean method with $\theta = 0.5$ and $h = 5 \times 10^{-4}$ s. (a) The individual energetic components of the simulation. (b) The sum of the change in energy work of the system with the cumulative dissipation. (c) The maximum value of the positive part of the tangential work $w_T$ over all contact nodes at each time step. (d) The minimum value of the negative part of the tangential work $w_T$ over all contact nodes at each time step.</figcaption></figure>

predictions to that obtained using the classical scheme, although this question must be studied statistically and in the context of an individual site, which is beyond the scope of this paper.

# 6 Conclusion

In this article, we have demonstrated that for linear visco-elasto-dynamics, the Frémond model of frictional contact gives a treatment of Coulomb friction combined with Newton's impact law that is always dissipative in continuous time. By adapting the Moreau–Jean time integration scheme to approximate the Frémond model in discrete time, we are able to develop a numerical method that is also provably always dissipative. We then considered several examples and demonstrated that in cases where the classical Moreau–Jean scheme creates energy, the new scheme demonstrably does not. Even in cases that are outside the regime for which we have formally proved the dissipativity of the model, we still observe substantially better performance using the new scheme. This work demonstrates new results and techniques that are likely to be of significance across a variety of subfields in mechanics, ranging from simulations of granular
<figure><img src="image_14.png" /><figcaption>**Figure 13** Results of the simulation using the proposed method $\theta = 0.5$ and $h = 5 \times 10^{-4}$ s. (a) The individual energetic components of the simulation. (b) The sum of the change in energy work of the system with the cumulative dissipation. (c) The maximum value of the positive part of the tangential work $w_T$ over all contact nodes at each time step. (d) The minimum value of the negative part of the tangential work $w_T$ over all contact nodes at each time step.</figcaption></figure>

media, rockfall simulation to design of mechanical linkages.

In future extensions of this work, application of our results to nonlinear Lagrangian systems and to Hamiltonian systems for solids in finite rotations or more generally, to nonlinear systems that derive from potentials, seems quite feasible without difficulty. For the sake of brevity, this is not done in this article. We can also consider large sliding problems that are generally best treated with a mortar method (Puso et al. 2024), and to consider the possibility of other impact laws such as those including a cohesion acting on the surface, or general standard and implicit material models as illustrated for instance in (45).

## References

[SW Mod.] Acary, V., rock_protection_wall.py, 2024, Inria. SWHID: 〈swh:1:cnt:ade5eb54939407f058 dbaaec28ff285b5e71ba86;origin=https://gricad gitlab.univ-grenoble-alpes.fr/nonsmooth /siconos-tutorials;visit=swh:1:snp:a2682a0bb31f68fd424b93a489b7b171cbebo340;anchor =swh:1:rev:a04ed23ceb61f4d4024ff75fe45aff7c65acbcdb;path=/examples/mechanics/Mas
<figure><img src="image_15.png" /><figcaption>**Figure 14** Total dissipation of the system for each method.</figcaption></figure>

onry/RockProtectionWall/rock_protection_wall.py).

Acary, V., O. Bonnefon, M. Brémond, O. Huber, and S. Sinclair (2019). *An Introduction to Siconos*. Grenoble: INRIA. [HAL].

Acary, V., F. Cadoux, C. Lemaréchal, and J. Malick (2011). A formulation of the linear discrete Coulomb friction problem via convex optimization. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik* 91(2):155–175. [DOI], [HAL].

Acary, V. (2016). Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik* 96(5):585–603. [DOI], [ARXIV].

Acary, V., P. Armand, H. Minh Nguyen, and M. Shpakovych (2024). Second-order cone programming for frictional contact mechanics using interior point algorithm. *Optimization Methods and Software* 39(3):634–663. [DOI], [HAL].

Acary, V., F. Bourrier, and B. Viano (2023). Variational approach for nonsmooth elasto-plastic dynamics with contact and impacts. *Computer Methods in Applied Mechanics and Engineering* 414:116156. [DOI], [OA].

Acary, V., M. Brémond, and O. Huber (2018). On solving frictional contact problems: formulations and comparisons of numerical methods. *Advanced Topics in Nonsmooth Dynamics*. Springer. Chap. 10, pp 375–457. [DOI], [HAL].

Ahrens, J., B. Geveci, and C. Law (2005). *36 - ParaView: An End-User Tool for Large-Data Visualization*. *Visualization Handbook*. Ed. by C. Hansen and C. Johnson. Burlington: Elsevier, pp 717–731. [DOI].

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

Alizadeh, F. and D. Goldfarb (2003). Second-order cone programming. *Mathematical Programming* 95(1):3–51. [DOI].

Andersen, E., C. Roos, and T. Terlaky (2003). On implementing a primal-dual interior-point method for conic quadratic optimization. *Mathematical Programming* 95(2):249–277. [DOI].

Barquins, M. (1985). Sliding Friction of Rubber and Schallamach Waves – A Review. *Materials Science and Engineering* 73:45–63. [DOI].

Berman, N., G. Cohen, and J. Fineberg (2020). Dynamics and Properties of the Cohesive Zone in Rapid Fracture and Friction. *Physical Review Letters* 125(12):125503. pmid: 33016754. [DOI].

Bourrier, F. and V. Acary (2022). Predictive capabilities of 2D and 3D block propagation models integrating block shape assessed from field experiments. *Rock Mechanics and Rock Engineering* 55(2):591–609. [DOI], [HAL].
Brach, R. M. (1984). Friction, restitution, and energy loss in planar collisions. *Journal of Applied Mechanics* 51(1):164-170. [DOI].

Brach, R. M. (1998). Formulation of rigid body impact problems using generalized coefficients. *International Journal of Engineering Science* 36(1):61-71. [DOI], [OA].

Brogliato, B. (1996). *Nonsmooth Mechanics*. Ed. by M. Thoma. 1st ed. Vol. 220. Lecture Notes in Control and Information Sciences. Springer. [DOI].

Brogliato, B. (1999). *Nonsmooth Mechanics*. Ed. by M. Thoma and E. Sontag. 2nd ed. Communications and Control Engineering. Springer. [DOI].

Brogliato, B. (2016). *Nonsmooth Mechanics: Models, Dynamics and Control*. Ed. by A. Isidori, J. H. van Schuppen, E. D. Sontag, and M. Krstic. 3rd ed. Communications and Control Engineering. Cham: Springer International Publishing. [DOI].

Brogliato, B., B. Maschke, R. Lozano, and O. Egeland (2013). *Dissipative systems analysis and control: theory and applications*. Springer. [DOI].

Brüls, O., V. Acary, and A. Cardona (2014). Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme. *Computer Methods in Applied Mechanics and Engineering* 281(1):131-161. [DOI], [OA].

Brüls, O. and A. Cardona (2010). On the use of Lie group time integrators in multibody dynamics. *Journal of Computational and Nonlinear Dynamics* 5(3):031002. [DOI], [OA].

Capobianco, G. and S. R. Eugster (2018). Time finite element based Moreau-type integrators. *International Journal for Numerical Methods in Engineering* 114(3):215-231. [DOI].

Cholet, C. (1998). Chocs de solides rigides. French. PhD thesis. France: Université Paris 6.

Collins-Craft, N. A., F. Bourrier, and V. Acary (2022). On the formulation and implementation of extrinsic cohesive zone models with contact. *Computer Methods in Applied Mechanics and Engineering* 400:115545. [DOI], [OA].

Collins-Craft, N. A., F. Bourrier, and V. Acary (2024). *On the formulation and implementation of mixed mode I and mode II extrinsic cohesive zone models with contact and friction*. [HAL].
Preprint.

[SW] Coumans, E. and Y. Bai, *PyBullet, a Python module for physics simulation for games, robotics and machine learning*, 2021.

Facchinei, F. and J.-S. Pang (2003). *Finite-Dimensional Variational Inequalities and Complementarity Problems – Volume I*. Springer Series in Operations Research. Springer. [DOI].

Fetecau, R. C., J. E. Marsden, M. Ortiz, and M. West (2003). Nonsmooth Lagrangian Mechanics and Variational Collision Integrators. *SIAM Journal on Applied Dynamical Systems* 2(3):381-416. [DOI], [OA].

Frémond, M. (2007). *Collisions*. Universita di Roma Tor Vergata - Dip Ing. Civile.

Frémond, M. (1995). Rigid bodies collisions. *Physics Letters A* 204(1):33-41. [DOI], [HAL].

Frémond, M. (2001). Internal constraints in mechanics. *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 359(1789):2309-2326. [DOI], [HAL].

Frémond, M. (2002). *Non-Smooth Thermomechanics*. Springer. [DOI].

Frémond, M. (2017). The Theory: Mechanics. An Example: Collision of a Point and a Plane. *Collisions Engineering: Theory and Applications*. Springer. Chap. 2, pp 5-32. [DOI].

Geuzaine, C. and J.-F. Remacle (2009). Gmsh: A 3-D Finite Element Mesh Generator with Built-in Pre- and Post-Processing Facilities. *International Journal for Numerical Methods in Engineering* 79(11):1309-1331. [DOI], [HAL].

Glocker, C. (2013). Energetic consistency conditions for standard impacts: Part I: Newton-type inequality impact laws and Kane's example. *Multibody System Dynamics* 29(1):77-117. [DOI], [OA].

Greenhalgh, S., V. Acary, and B. Brogliato (2013). On preserving dissipativity properties of linear complementarity dynamical systems with the θ-method. *Numerische Mathematik* 125(4):601-637. [DOI], [HAL].

Hairer, E., C. Lubich, and G. Wanner (2003). Geometric numerical integration illustrated by the Störmer-Verlet Method. *Acta Numerica* 12:399-450. [DOI], [HAL].

Holm, D., J. Marsden, T. Ratiu, and A. Weinstein (1984). Stability of Rigid Body Motion Using the Energy-Casimir Method. *Fluids and Plasmas: Geometry and Dynamics*. Contemporary
Mathematics (v. 28), 15–23. [DOI], [OA].

Jean, M. (1999). The non-smooth contact dynamics method. *Computer Methods in Applied Mechanics and Engineering* 177(3-4). Special issue on computational modeling of contact and friction:235–257. [DOI], [HAL].

Jean, M. and J. J. Moreau (1987). Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. *Unilateral Problems in Structural Analysis*. Ed. by G. Del Pietro and F. Maceri. 2nd ed. Springer, pp 151–196. [DOI], [HAL].

Johnson, G., S. Leyendecker, and M. Ortiz (2014). Discontinuous Variational Time Integrators for Complex Multibody Collisions. *International Journal for Numerical Methods in Engineering* 100(12):871–913. [DOI].

Jourdan, F., P. Alart, and M. Jean (1998). A Gauss-Seidel like Algorithm to Solve Frictional Contact Problems. *Computer Methods in Applied Mechanics and Engineering* 155(1-2):31–47. [DOI], [HAL].

Kane, T. R. (1984). A dynamics puzzle. *Stanford Mechanics Alumni Club Newsletter* 6:1-4.

Leine, R. and N. van de Wouw (2008). *Stability and Convergence of Mechanical Systems with Unilateral Constraints*. Vol. 36. Lecture Notes in Applied and Computational Mechanics. Springer. [DOI].

Moreau, J. J. (1970). Sur les lois de frottement, de plasticité et de viscosité. *Comptes rendus de l'Académie des sciences. Série A - Sciences mathématiques* 271:608-611. [HAL].

Moreau, J. J. (1971). Sur l'évolution d'un système élasto-visco-plastique. *Comptes rendus hebdomadaires des séances de l'Académie des sciences* 273:118-121. [HAL].

Moreau, J. J. (1973). On Unilateral Constraints, Friction and Plasticity. *New Variational Techniques in Mathematical Physics*. Ed. by G. Capriz and G. Stampacchia. Springer, pp 171–322. [DOI], [HAL].

Moreau, J. J. (1988). Unilateral contact and dry friction in finite freedom dynamics. *Nonsmooth Mechanics and Applications*. Ed. by J. J. Moreau and P. D. Panagiotopoulos. CISM, Courses and lectures 302. Springer, pp 1–82. [DOI], [HAL].

Puso, M. A., J. H. Porter, and T. Slavik (2024). Structure preserving and energy dissipative contact approaches for implicit dynamics. *International Journal for Numerical Methods in Engineering* 125(10):e7454. [DOI].

Schallamach, A. (1971). How Does Rubber Slide? *Wear* 17(4):301-312. [DOI].

[SW] Schlömer, N., *meshio: Tools for mesh files*, 2024. swhID: (swh:1:dir:6e1cd14ba43f09c1e729f41f4b6e2ad5fbod6c24;origin=https://github.com/nschloe/meshio;visit=swh:1:snp:9d55efa5b1f7ff0238d8557c5aee1673b55e3d29;anchor=swh:1:rev:b2ee99842e119901349fdeee06b5bf61e01f450a).

Viswanathan, K. and S. Chandrasekar (2022). Fifty Years of Schallamach Waves: From Rubber Friction to Nanoscale Fracture. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 380(2232):20210339. [DOI], [OA].

Ziegler, H. (1983). *An Introduction to Thermomechanics*. Ed. by E. Becker, B. Budiansky, W. T. Koiter, and H. A. Lauwerier. Rev. ed. Series in Applied Mathematics and Mechanics v. 21. Elsevier, pp v-x. [DOI].
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long

as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** VA: Conceptualisation, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, software, supervision, visualisation, writing – original draft, writing – review and editing NACC: Data curation, formal analysis, funding acquisition, investigation, methodology, software, visualisation, writing – original draft, writing – review and editing

**Supplementary Material** The code required to run the simulations described in this paper is available in a GitHub repository archived on Software Heritage, or alternatively on Zenodo. The sole exception to this is the code required to run the rockfall impact code, which is archived separately on Software Heritage as part of the Siconos tutorial repository. The data outputs of the codes are also available on Zenodo.

**Acknowledgements** The authors acknowledge the support of the Marie Skłodowska-Curie Actions program under the Horizon Europe research and innovation framework program (Grant agreement ID 101064805 LEMMA). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or Horizon Europe. Neither the European Union nor the granting authority can be held responsible for them. The authors also thank Franck Bourrier and Ritesh Gupta for their helpful comments on the manuscript.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.