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

HAL
open science

# An explicit pseudo-energy conservative scheme for contact between deformable solids

Nadine Dirani, Laurent Monasse

► To cite this version:

Nadine Dirani, Laurent Monasse. An explicit pseudo-energy conservative scheme for contact between deformable solids. 2022. hal-03879345

HAL Id: hal-03879345

https://hal.science/hal-03879345

Preprint submitted on 30 Nov 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An explicit pseudo-energy conservative scheme for contact between deformable solids

Nadine Dirani¹, Laurent Monasse¹

¹Université Côte d'Azur, Inria, CNRS, LJAD, EPC COFFEE, 06108 Nice, France

email: nadine.dirani@inria.fr, laurent.monasse@inria.fr

November 30, 2022

## Abstract

We extend the work of [21] on elastodynamics to the treatment of contact. To that end, we propose adequate handling of boundary conditions, either through the resolution of local problems on each of the face displacement unknowns or through interpolation from nearest neighbours. Adapting the time-integration strategy adopted in [20], it is possible to conserve both momentum and a pseudo-energy exactly. Numerical results are presented to illustrate the accuracy of contact treatment and the energy conservation of the system.

## 1 Introduction

The development of numerical methods for dynamic contact problems is a central issue for applications ranging from tire dynamics [12, 13], fracture dynamics [7] and fluid-structure interaction with explosions. An overview of contact problems and applications can be found in [26]. A key issue in this context is to be able to accurately and efficiently capture non-penetration while conserving energy, without spurious oscillations. We focus here on fast dynamics with high-velocity contact for applications in crash-tests and safety assessment.

A large number of works have been devoted to the resolution and numerical analysis of contact variational formulations [22, 16]. A review of contact problems can be found in [19] and [3]. Contact problems involve difficulties on the theoretical and numerical sides and implicate non-linear contact conditions on a part of the boundary. The numerical methods can be classified along three axes: the space discretization method employed, the treatment of the contact boundary condition, and the time-integration scheme used. For the space discretization method,
Finite Element methods have widely been used [5, 9, 26]. The contact boundary condition can be either enforced exactly [6, 14, 18, 23, 24], enforced with penalty [1, 13, 4, 5, 3] or enforced with contact conditions in velocity [1, 17]. Exact enforcement generally hinges on the use of a Lagrange multiplier for non-penetration which can be identified as the contact pressure. Penalty enforcement used to be restricted by the use of a penalization parameter, which had to be chosen carefully since a small value would mean a lax enforcement of non-penetration and a large value would add numerical stiffness to the problem. However, recent works using Nitsche's penalization [4, 5, 3] have demonstrated a robust and efficient choice for the penalization parameter. Finally, the enforcement of contact conditions in velocity translates non-penetration on the relative velocities of contact surfaces, which separates the penetration condition on displacement from its enforcement on velocity. This allows for innovative time-integration strategies.

On the time-discretization side, three main strategies have been applied: implicit schemes, explicit schemes and semi-explicit schemes. We refer to [9] for a review and comparison of these strategies for contact. Implicit schemes, while being generally more stable than explicit schemes, can exhibit spurious oscillations in the case of exact contact conditions enforcement. Explicit schemes generally have difficulties enforcing exactly non-penetration while limiting energy variations, except in the case of modified mass strategies [15, 9]. Semi-explicit schemes try to circumvent this by expliciting the general dynamics, except the non-penetration condition. A useful tool to compare schemes is energy conservation [13, 9]: while dissipative schemes remove spurious oscillations, they are unable to correctly account for multiple bounces, as displacement amplitude is reduced at each contact even without physical dissipation in the model, as the numerical dissipation of energy becomes dominant.

In the present work, we restrict ourselves to frictionless contact in plane linear eleasticity. We choose to use a Discrete Element method (DEM) in the form of a Cell-centered Galerkin variational form [21]. Exact contact non-penetration is enforced on the boundary degrees of freedom, which have no mass. For the time discretization, we use an explicit pseudo-energy and momentum conserving scheme for the time integration of Hamiltonian systems [20]. This scheme is a two-step method of order two which ensures a discrete pseudo-energy conservation, relying on an exact quadrature of force time integrals for nonlinear behavior laws. It is also accompanied by an *a posteriori* stability criterion which allows easier management of the stability of the scheme. When contact occurs, discontinuities occur in the forces which render high-order quadratures inefficient. We propose a modification of the time quadrature to account for contact. It should be noted that the proposed scheme is, to the best of our knowledge, the first explicit scheme to exactly conserve an energy with exact enforcement of contact conditions. Note
that most energy-conserving schemes are implicit to a certain degree for non-linear problems [13, 25, 11], with the exception of a very recent result by Bilbao et al. [2].

This paper is organized as follows. Section 2 briefly recalls the equations of linear elastodynamic in small deformations with frictionless contact. Section 3 presents the space discretization of the governing equations with DEM. In this section, we present two different treatments of boundary conditions, schemes A and B, and their extension to multiple solid contacts. We derive the forces used in DEM. Section 4 addresses the time discretization of the elastodynamic problem with contact, extending the results developed in [20]. Numerical results in one and two space dimensions are presented and discussed in Section 5, demonstrating the accuracy, efficiency and robustness of the proposed method.

## 2 Governing equations for elastodynamics

We consider an isotropic linearly elastic material occupying the domain $\Omega \in \mathbb{R}^d$, $d = \{2, 3\}$ during a finite time interval $(0, T)$ undergoing small deformations under external volumetric forces $f$. The linearized strain tensor is given by $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$, where $\mathbf{u}$ is the displacement field. Hooke's law connects the stress tensor $\sigma$ to the strain tensor $\varepsilon$ by $\sigma = \lambda \text{tr} \varepsilon \mathbf{I} + 2\mu\varepsilon$ where $\lambda$ and $\mu$ are the Lamé coefficients which verify $\mu > 0$ et $2\mu + d\lambda > 0$.

The boundary $\partial\Omega$ of $\Omega$ is polygonal ($d = 2$) or polyhedral ($d = 3$). Let $\partial\Omega = \Gamma^D \cup \Gamma^N(t) \cup \Gamma^C(t)$ be a partition of the boundary of $\Omega$ where $\Gamma^D, \Gamma^N(t)$ and $\Gamma^C(t)$ are the Dirichlet, Neumann and unilateral contact boundary conditions on displacement at time $t$ respectively (see Figure 1). Displacement is imposed on $\Gamma^D$, while the normal stress is imposed on $\Gamma^N(t)$:

$$ \mathbf{u} = \mathbf{u}_D \text{ on } (0, T) \times \Gamma^D, \quad \sigma \cdot \mathbf{n} = \mathbf{0} \text{ on } (0, T) \times \Gamma^N(t), \quad (1) $$

where $\mathbf{u}_D$ is a vector-valued function. For simplicity, we suppose that the contact boundary in 2D is a straight line while in 3D, it is supposed to be a polygon. The normal unit outward vector on $\partial\Omega$ is denoted $\mathbf{n}$. The elastodynamic equations in strong form consist in finding for the displacement field $\mathbf{u}: (0, T) \times \Omega \to \mathbb{R}^d$ such that:

$$ \rho \ddot{\mathbf{u}} - \text{div} \, \sigma(\mathbf{u}) = \mathbf{f} \quad \text{in } (0, T) \times \Omega, \quad (2) $$

where $\rho > 0$ is the density of the elastic material and $\ddot{\mathbf{u}}$ is the acceleration. We impose additionally the initial conditions $\mathbf{u}(0, \cdot) = \mathbf{u}_0$, $\dot{\mathbf{u}}(0, \cdot) = \mathbf{v}_0$ in $\Omega$. In the case of quasi-static computations, the acceleration term is removed from the equations.
## 2.1 Contact condition

<figure><img src="image_3.png" /><figcaption>Figure 1: Elastic body $\Omega$ in contact with a rigid foundation (left) or in self-contact (right).</figcaption></figure>

Assume that the contact boundary condition $\Gamma^C(t)$ is partitioned in $\Gamma_1^C(t) \cup \Gamma_2^C(t)$ such that the solid boundary points on $\Gamma_1^C(t)$ are in contact with points on $\Gamma_2^C(t)$ (see Figure 1). Let $\mathbf{x}_1^0 \in \Gamma_1^C(t)$ and $\mathbf{x}_2^0 \in \Gamma_2^C(t)$, and set $\mathbf{n}_1$ and $\mathbf{n}_2$ the outward normals to $\Omega$ at $\mathbf{x}_1^0$ and $\mathbf{x}_2^0$ respectively. We use the simplifying assumption that the boundary elements in contact have opposite normals: $\mathbf{n}_2(\mathbf{x}_2^0) = -\mathbf{n}_1(\mathbf{x}_1^0)$.

Set the normal penetration $u_n(\mathbf{x}_1^0, \mathbf{x}_2^0)$ of the material points at $\mathbf{x}_1^0$ and $\mathbf{x}_2^0$:

$$u_n(\mathbf{x}_1^0, \mathbf{x}_2^0) = (\mathbf{x}_2^0 + \mathbf{u}(t, \mathbf{x}_2^0) - \mathbf{x}_1^0 - \mathbf{u}(t, \mathbf{x}_1^0)) \cdot \mathbf{n}_1. \quad (3)$$

Newton's third law gives that $\sigma(t, \mathbf{x}_1^0) \cdot \mathbf{n}_1 = -\sigma(t, \mathbf{x}_2^0) \cdot \mathbf{n}_2$. In the case of frictionless contact, setting the normal contact force $\sigma_{nn}(\mathbf{x}_1^0)$ and the tangential contact force $\sigma_t(\mathbf{x}_1^0)$ as

$$\sigma_{nn}(\mathbf{x}_1^0) = \mathbf{n}_1 \cdot \sigma(t, \mathbf{x}_1^0) \cdot \mathbf{n}_1, \quad \sigma_t(\mathbf{x}_1^0) = \sigma(t, \mathbf{x}_1^0) \cdot \mathbf{n}_1 - \sigma_{nn}\mathbf{n}_1,$$

it is sufficient to check that contact does not result in traction nor friction, i.e. $\sigma_{nn}(\mathbf{x}_1^0) \le 0$ and $\sigma_t(\mathbf{x}_1^0) = 0$.

The unilateral complementarity contact conditions without friction on $\Gamma^C(t)$ at points $(\mathbf{x}_1^0, \mathbf{x}_2^0)$ write

$$u_n(\mathbf{x}_1^0, \mathbf{x}_2^0) \ge 0, \quad \sigma_{nn}(\mathbf{x}_1^0) \le 0, \quad u_n(\mathbf{x}_1^0, \mathbf{x}_2^0)\sigma_{nn}(\mathbf{x}_1^0) = 0, \quad \sigma_t(\mathbf{x}_1^0) = 0. \quad (4)$$
## 2.2 Mechanical energy

The total energy of the solid $\mathcal{E}(t)$ at time $t$ is defined as the sum of the kinetic energy $\mathcal{E}_{kin}(t)$ and the elastic deformation energy $\mathcal{E}_{elas}(t)$

$$ \mathcal{E}(t) = \frac{1}{2} \underbrace{\int_{\Omega} \rho \left| \frac{\partial \mathbf{u}}{\partial t} \right|^2}_{\mathcal{E}_{kin}(t)} + \underbrace{\frac{1}{2} \int_{\Omega} \varepsilon(\mathbf{u}) : \sigma(\mathbf{u})}_{\mathcal{E}_{elas}(t)} \quad (5) $$

In the case of homogeneous Dirichlet boundary conditions and homogeneous Neumann boundary conditions, the total mechanical energy is conserved.

# 3 Discrete element method space discretization

## 3.1 Degrees of freedom and reconstruction operator on facets

Let $\mathcal{T}^h$ be a family of triangulations of the domain $\Omega$. We assume that $\Omega$ is a polygon and we also assume that the mesh is compatible with the partition of the boundary $\partial\Omega$ into the Dirichlet and Neumann parts.

<figure><img src="image_4.png" /><figcaption>Figure 2: Vector-valued degrees of freedom for the displacement in $\Omega$.</figcaption></figure>

Let us introduce $\mathcal{C}$ the set of mesh cells and $\mathcal{F} = \mathcal{F}^i \cup \mathcal{F}^b$ the set of mesh faces where $\mathcal{F}^i$ denotes the set of interior faces of the solid shared by two mesh cells and $\mathcal{F}^b$ the set of boundary faces, which belong only to one mesh cell, sitting on the boundary of $\partial\Omega$. The vector-valued volumetric degrees of freedom (dofs) for a generic displacement field $\mathbf{u}_\mathcal{C} := (\mathbf{u}_c)_{c \in \mathcal{C}} \in \mathbb{R}^{d\#(\mathcal{C})}$ are placed at the barycenter of every mesh cell $c \in \mathcal{C}$, where $\#(\mathcal{S})$ denotes the cardinality of a set $\mathcal{S}$. In addition, we denote by $\mathbf{u}_\mathcal{F} := (\mathbf{u}_f)_{f \in \mathcal{F}} \in \mathbb{R}^{d\#(\mathcal{F})}$ the reconstructed displacement on the
mesh faces (see Figure 2). For all faces $f \in \mathcal{F}$, we denote by $\mathcal{R}$ the reconstruction operator and we set

$$ \mathbf{u}_f := \mathcal{R}(\mathbf{u}_h)|_f = \sum_{c \in \mathcal{C}_f} \mathbf{T}_{cf} \cdot \mathbf{u}_c \qquad (6) $$

where $\mathcal{C}_f$ denotes the subset of cells used to reconstruct $f$ and $\mathbf{T}_{cf}$ is a second-order tensor. For interior faces $f \in \mathcal{F}^i$, the set $\mathcal{C}_f$ is chosen as $d+1$ cells ($c_1, \dots, c_{d+1}$) whose respective centers ($\mathbf{x}_1, \dots, \mathbf{x}_{d+1}$) form a simplex containing the facet barycenter $\mathbf{x}_f$ (see Figure 3). The reconstruction tensor reduces to a scalar value $\mathbf{T}_{cf} = \alpha_{cf} \mathbf{I}$ and the coefficients $\alpha_{cf}$ are chosen as the barycentric coordinates of $\mathbf{x}_f$ in ($\mathbf{x}_1, \dots, \mathbf{x}_{d+1}$):

$$ \mathbf{x}_f = \sum_{i=1}^{d+1} \alpha_{c_i f} \mathbf{x}_i, \quad \sum_{i=1}^{d+1} \alpha_{c_i f} = 1. \qquad (7) $$

The reconstruction operator is then a linear interpolation on the interior points and the coefficients $\alpha_{cf}$ lie in $(0, 1)$. The definition of reconstruction tensors $\mathbf{T}_{cf}$ for boundary faces $f \in \mathcal{F}^b$ is left for Section 3.4.

<figure><img src="image_5.png" /><figcaption>Figure 3: Dofs associated with the interior facet $f$ used for the reconstruction.</figcaption></figure>

## 3.2 Discrete mechanical energy

It is possible to reconstruct a discrete piecewise-constant gradient $\mathcal{G}_c(\mathbf{u}_f) := (\mathcal{G}_c(\mathbf{u}_f))_{c \in \mathcal{C}}$ in each cell from the reconstructed facet displacements. Using a discrete Stokes formula, for all $c \in \mathcal{C}$, we have

$$ \mathcal{G}_c(\mathcal{R}(\mathbf{u}_h)_{\mathcal{F}}) = \mathcal{G}_c(\mathbf{u}_{\mathcal{F}}) = \sum_{f \in \partial \mathcal{C}} \frac{|f|}{|c|} \mathbf{u}_f \otimes \mathbf{n}_{f,c} \qquad (8) $$
where $|f|$ is the surface of face $f$, $|c|$ is the volume of cell c and $n_{f,c}$ is the outward normal to c on f. The piecewise-constant linearized strain tensor $\varepsilon_c$ and the stress tensor $\sigma_c$ in each cell $c \in \mathcal{C}$ are defined such that

$$ \varepsilon_c(\mathbf{u}_F) = \frac{1}{2} (\mathcal{G}_c(\mathbf{u}_F) + \mathcal{G}_c(\mathbf{u}_F)^T), \quad \sigma_c = \mathbb{C} : \varepsilon_c = \lambda \operatorname{tr} \varepsilon_c \mathbf{I} + 2\mu \varepsilon_c \quad (9) $$

Finally, we define a cellwise nonconforming $P^1$ reconstruction $\mathfrak{R}$ defined for all $c \in \mathcal{C}$ by

$$ \mathfrak{R}(\mathbf{u}_h)_c(\mathbf{x}) := \mathbf{u}_c + \mathcal{G}_c(\mathcal{R}(\mathbf{u}_h)) \cdot (\mathbf{x} - \mathbf{x}_c) \quad (10) $$

where $\mathbf{x} \in c$ is an arbitrary point in the cell and $\mathbf{x}_c$ is the position of the barycentre of the cell c.

The discrete elastic deformation energy $\mathcal{E}_p$ is written as an analogue to the continuous $\mathcal{E}_{elas}$ (5):

$$ \mathcal{E}_p(\mathbf{u}_h) = \sum_{c \in \mathcal{C}} \frac{1}{2} |c| \sigma_c(\mathbf{u}_h) : \varepsilon_c(\mathbf{u}_h) \quad (11) $$

The kernel of the discrete gradient operator $\mathcal{G}_c \circ \mathcal{R}$ contains the vectors $\mathbf{u}_h$ which have identical components on each cell, i.e. the evaluation of a constant function on the mesh. This is due to the geometric property on closed cells c:

$$ \sum_{f \in \partial c} |f| n_{f,c} = 0. $$

However, the kernel is not restricted to constant functions in general, which may lead to hourglass instabilities. The aim of stabilization is to eliminate the spurious modes from the kernel $\mathcal{G}_c \circ \mathcal{R}$. It consists in penalizing the jumps between $\mathbb{P}_1$ reconstructions of local adjacent faces. We denote by $[\cdot]_F$ the jump of a quantity through an interior face $F$: for the reconstruction $\mathfrak{R}$, denoting by $c_1$ and $c_2$ the two mesh cells sharing the interior face $f$ such that $f = \partial c_1 \cap \partial c_2$ (see Figure 4), we define

$$ [\mathfrak{R}(\mathbf{u}_h)]_f = \mathfrak{R}(\mathbf{u}_h)_{c_1}(\mathbf{x}_f) - \mathfrak{R}(\mathbf{u}_h)_{c_2}(\mathbf{x}_f). \quad (12) $$

We define the weakly consistent stabilization energy $\mathcal{E}_S$ as follows:

$$ \mathcal{E}_S(\mathbf{u}_h) = \sum_{f \in \mathcal{F}} \frac{\eta}{2h_f} |f| [\mathfrak{R}(\mathbf{u}_h)]_f^2, \quad (13) $$

where $h_F$ is the diameter of the facet $f \in \mathcal{F}$. The dimension of $\eta$ is the same as that of the Lamé coefficients $\lambda$ and $\mu$ and the Young's modulus $E$. The scheme
<figure><img src="image_6.png" /><figcaption>Figure 4: Two mesh cells $c_1$ and $c_2$ sharing an interior face $f$.</figcaption></figure>

converges to the exact solution for any positive $\eta$. A reasonable choice in practice is to take $\eta = E$ [21].

Setting $v_c = \dot{u}_c$ the cell velocity, the full discrete energy $\mathcal{E}_h$ is defined as

$$ \mathcal{E}_h(\mathbf{u}_h) = \sum_{c \in \mathcal{C}} \rho |c| \mathbf{v}_c^2 + \mathcal{E}_P(\mathbf{u}_h) + \mathcal{E}_S(\mathbf{u}_h). \quad (14) $$

### 3.3 Derivation of forces

The energy $\mathcal{E}_h$ of the system allows to define the force $F_c$ exerted on a particle $c$ deriving from the potential energy in (14):

$$ F_c = -\nabla_{\mathbf{u}_c}(\mathcal{E}_P + \mathcal{E}_S) = F_c^P + F_c^S, \quad (15) $$

where $F_c^P$ and $F_c^S$ denote respectively the force contribution of the elastic deformation and of the stabilization. A simple calculation gives that

$$ F_c^P = - \sum_{f \in \mathcal{F}} T_{cf} \cdot \sum_{c \ni f} \sigma_c \cdot n_{f,c} + F_c^G. \quad (16) $$

Note that for a given cell $c$, $T_{cf}$ is null for almost all facets $f$, except those in which cell $c$ contributes to the reconstruction operator $\mathcal{R}_f$. The right-hand side is in fact a weighted sum of stress jumps on faces. The additional force $F_c^G$ occurs in the case of sliding contact interfaces and originates from the fact that certain reconstruction tensors $T_{cf}$ depend on the $\mathbf{u}_h$. This term is explained in Section 3.4.3 and detailed in equation (34).

Similarly, the stabilization force can be computed with long but straightforward
calculations:

$$ 
\begin{aligned}
\boldsymbol{F}_c^S = & - \sum_{f \in \mathcal{F}^i, f \in c} \frac{\eta}{h_f} |f| \delta_{cf} [\mathfrak{R}(\boldsymbol{u}_h)]_f \\
& - \sum_{\tilde{f} \in \mathcal{F}} \eta |\tilde{f}| \boldsymbol{T}_{c\tilde{f}} \cdot \left( \sum_{\tilde{c} \ni \tilde{f}} \sum_{f \in \mathcal{F}^i, f \in \tilde{c}} \delta_{\tilde{c}\tilde{f}} \frac{|f|}{h_f |\tilde{c}|} (\boldsymbol{n}_{\tilde{c}\tilde{f}} \cdot (\boldsymbol{x}_f - \boldsymbol{x}_{\tilde{c}})) [\mathfrak{R}(\boldsymbol{u}_h)]_f \right).
\end{aligned} 
\quad (17) $$

## 3.4 Treatment of boundary conditions

We detail here the expression of reconstruction tensors $T_{cf}$ for boundary facets $f \in \mathcal{F}^b$. We first propose two different reconstructions for the contact condition between two particles: first, scheme A implements the complementary constraints (4) at the discrete level as a Neumann/Dirichlet dichotomy; second, scheme B is a simplified (but less accurate) approach to this dichotomy. We then explain how to extend the contact condition to multiple particles in contact.

### 3.4.1 Scheme A

The homogeneous Neumann boundary condition for the boundary cells without contact writes $\boldsymbol{\sigma} \cdot \boldsymbol{n} = 0$, which, for a boundary face $f \in \mathcal{F}^b$ on a cell $c \in \mathcal{C}$, translates using equations (8) and (9):

$$ 
\begin{aligned}
\boldsymbol{\sigma}_c \boldsymbol{n}_{f,c} &= \lambda \operatorname{tr} \boldsymbol{\varepsilon}_c \boldsymbol{n}_{f,c} + 2\mu \boldsymbol{\varepsilon}_c \cdot \boldsymbol{n}_{f,c} \quad (18) \\
&= \sum_{\tilde{f} \in \partial c} \frac{|\tilde{f}|}{|c|} \left( \lambda(\boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{\tilde{f},c}) \boldsymbol{n}_{f,c} + \mu(\boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{f,c}) \boldsymbol{n}_{\tilde{f},c} + \mu(\boldsymbol{n}_{f,c} \cdot \boldsymbol{n}_{\tilde{f},c}) \boldsymbol{u}_{\tilde{f}} \right) = \mathbf{0}. \quad (19)
\end{aligned}
 $$

Equation (19) is a linear system on the unknown reconstruction $\boldsymbol{u}_f$ using the interior face reconstructions $\boldsymbol{u}_{\tilde{f}}$ defined in Section (3.1) which has a unique solution when cells $c \in \mathcal{C}$ have at most one boundary facet $f \in \mathcal{F}^b$:

$$ 
\begin{aligned}
\boldsymbol{u}_f = & - \left[ \sum_{\substack{\tilde{f} \in \partial c \\ \tilde{f} \neq f}} \frac{|\tilde{f}|}{|f|} \left( \frac{\lambda}{\lambda + 2\mu} (\boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{\tilde{f},c}) + \frac{2\mu}{\lambda + 2\mu} (\boldsymbol{n}_{f,c} \cdot \boldsymbol{n}_{\tilde{f},c}) \boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{f,c} \right) \right] \boldsymbol{n}_{f,c} \\
& - \sum_{\substack{\tilde{f} \in \partial c \\ \tilde{f} \neq f}} \frac{|\tilde{f}|}{|f|} \left( (\boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{f,c})(\boldsymbol{n}_{\tilde{f},c} - \boldsymbol{n}_{\tilde{f},c} \cdot \boldsymbol{n}_{f,c} \boldsymbol{n}_{f,c}) + (\boldsymbol{n}_{f,c} \cdot \boldsymbol{n}_{\tilde{f},c})(\boldsymbol{u}_{\tilde{f}} - \boldsymbol{u}_{\tilde{f}} \cdot \boldsymbol{n}_{f,c} \boldsymbol{n}_{f,c}) \right).
\end{aligned} 
\quad (20) $$
This gives the expression of the reconstruction tensor $T_{\tilde{c}\tilde{f}}$ in terms of the interpolation coefficients $\alpha_{\tilde{c}\tilde{f}}$ on the other faces $\tilde{f} \in \partial c$ for $c \ni f$:

$$ T_{\tilde{c}\tilde{f}} = -\sum_{\substack{\tilde{f} \in \partial c \\ \tilde{f} \neq f}} \frac{|\tilde{f}|}{|f|} \left( \frac{\lambda}{\lambda + 2\mu} (\mathbf{n}_{f,c} \otimes \mathbf{n}_{\tilde{f},c}) + \mathbf{n}_{\tilde{f},c} \otimes \mathbf{n}_{f,c} + (\mathbf{n}_{f,c} \cdot \mathbf{n}_{\tilde{f},c}) \mathbf{I} \right. \\ \left. - \frac{2(\lambda + \mu)}{\lambda + 2\mu} (\mathbf{n}_{f,c} \cdot \mathbf{n}_{\tilde{f},c}) \mathbf{n}_{f,c} \otimes \mathbf{n}_{f,c} \right) \alpha_{\tilde{c}\tilde{f}}. \quad (21) $$

In the case of contact between two faces $f_1$ and $f_2 \in \mathcal{F}^b$ such that $f_1 \in \partial c_1$ and $f_2 \in \partial c_2$, it is necessary to define a common normal for the two faces. Denote $\mathbf{n}_{f_1 f_2}$ the average normal taken as the contact normal:

$$ \mathbf{n}_{f_1 f_2} = \frac{|f_1|\mathbf{n}_{f_1 c_1} - |f_2|\mathbf{n}_{f_2 c_2}}{|f_1| + |f_2|}. $$

Similarly to (3), define the normal penetration $u_{f_1 f_2}$ as

$$ u_{f_1 f_2} = (\mathbf{x}_{f_2} + \mathbf{u}_{f_2} - \mathbf{x}_{f_1} - \mathbf{u}_{f_1}) \cdot \mathbf{n}_{f_1 f_2}. \quad (22) $$

As long as $u_{f_1 f_2} \ge 0$, the faces $f_1$ and $f_2$ are not in contact, and reconstruction (20)–(21) holds. However, when $u_{f_1 f_2} < 0$, we change the reconstruction operator so that

$$ u_{f_1 f_2} = 0, \quad \sigma_{nn,c_1} = \sigma_{nn,c_2}, \quad (23) $$

$$ \sigma_{t,c_1} = \mathbf{0}, \quad \sigma_{t,c_2} = \mathbf{0}, \quad (24) $$

where

$$ \sigma_{nn,c_1} = \mathbf{n}_{f_1 f_2} \cdot \sigma_{c_1} \cdot \mathbf{n}_{f_1 f_2}, \quad \sigma_{t,c_1} = \sigma_{c_1} \cdot \mathbf{n}_{f_1 f_2} - \sigma_{nn,c_1} \mathbf{n}_{f_1 f_2}, $$

$$ \sigma_{nn,c_2} = \mathbf{n}_{f_1 f_2} \cdot \sigma_{c_2} \cdot \mathbf{n}_{f_1 f_2}, \quad \sigma_{t,c_2} = \sigma_{c_2} \cdot \mathbf{n}_{f_1 f_2} - \sigma_{nn,c_2} \mathbf{n}_{f_1 f_2}. $$

Again, equations (23)–(24) constitute a local linear system which has a unique solution $(\mathbf{u}_{f_1}, \mathbf{u}_{f_2})$. The full expression of the solution is straightforward but lengthy and will not be expressed here.

In effect, the boundary condition changes only the reconstruction tensor $T_{cf}$. Remark that $T_{cf}$ does not depend on the cell displacements without contact in (21), but that contact makes it depend on $\mathbf{u}_h$. In addition, note that the two expressions coincide when $u_{f_1 f_2} = 0$: this translates the continuity of the reconstructed face displacements $\mathbf{u}_{f_1}$ and $\mathbf{u}_{f_2}$ with regards to the cell displacements $\mathbf{u}_h$.
### 3.4.2 Scheme B

The resolution of system (23)–(24) can be complex to implement, since it involves a larger stencil for which $T_{cf}$ is non trivial, with the cells $c$ distributed on both sides of the contact zone. In addition, the solution of the homogeneous Neumann condition (21) is valid only when there is only one Neumann facet per boundary cell. In case the mesh does not satisfy this condition, it cannot be guaranteed that there exists a solution to the problem. This leads us to use a simpler albeit less accurate scheme, noted scheme B.

On Neumann facets $f \in \mathcal{F}^b$ without contact, we set the value of the reconstructed displacement $u_f$ to the value of the adjacent cell $c \ni f$: $u_f = u_c$. This allows a very simple expression for the reconstruction tensor: $T_{cf} = I$, and for all $\tilde{c} \neq c$, $T_{\tilde{c}f} = 0$. This corresponds to a rigid body movement for the boundary cells. Obviously, it cannot be expected that the stress tensor is correct in boundary cells, since the stress boundary condition $\sigma \cdot n = 0$ is incorrectly captured.

In the case of contact between two faces $f_1$ and $f_2 \in \mathcal{F}^b$ such that $f_1 \in \partial c_1$ and $f_2 \in \partial c_2$, we define the common normal $n_{f_1 f_2}$ similarly to scheme A. Contact is detected when $u_{f_1 f_2} < 0$, and the reconstruction operator is modified. The slip boundary condition is written as:

$$u_{f_1 f_2} = 0, \quad (I - n_{f_1 f_2} \otimes n_{f_1 f_2}) \cdot (u_{f_1} - u_{c_1}) = 0, \quad (I - n_{f_1 f_2} \otimes n_{f_1 f_2}) \cdot (u_{f_2} - u_{c_2}) = 0. \tag{25}$$

The first equation translates the non-penetration condition, while the last two equations translate the fact that each cell follows an independent tangential rigid body movement. In fact, system (25) is not well-posed: there exist infinitely many solutions of the form, for any $\theta \in \mathbb{R}$,

$$u_{f_1} = u_{c_1} + \theta(u_{c_2} + x_{f_2} - u_{c_1} - x_{f_1}) \cdot n_{f_1 f_2} n_{f_1 f_2} \tag{26}$$

$$u_{f_2} = u_{c_2} - (1 - \theta)(u_{c_2} + x_{f_2} - u_{c_1} - x_{f_1}) \cdot n_{f_1 f_2} n_{f_1 f_2}. \tag{27}$$

Note that these solution verify $(u_{f_2} + x_{f_2} - u_{f_1} - x_{f_1}) \cdot n_{f_1 f_2} = 0$ and (25). A symmetric choice for parameter $\theta$, which gives the same role to faces $f_1$ and $f_2$, is $\theta = \frac{1}{2}$. We opt for this particular choice in the sequel. The reconstruction tensor then has the following nontrivial expressions:

$$T_{c_1 f_1} = I - \frac{1}{2} n_{f_1 f_2} \otimes n_{f_1 f_2}, \quad T_{c_2 f_1} = \frac{1}{2} n_{f_1 f_2} \otimes n_{f_1 f_2}, \tag{28}$$

$$T_{c_2 f_2} = I - \frac{1}{2} n_{f_1 f_2} \otimes n_{f_1 f_2}, \quad T_{c_1 f_2} = \frac{1}{2} n_{f_1 f_2} \otimes n_{f_1 f_2}. \tag{29}$$

### 3.4.3 Multiple particle contact

In general, contact does not occur only between pairs of facets: most facets are in contact with two or more facets, which in turn are in contact with other facets.
Writing the conditions of non-penetration and continuity of normal stress (23) of scheme A leads to a fully coupled problem on all the facets in a contiguous contact zone, which usually has no solution. This is why we will describe multiple particle contact only in the case of scheme B.

For a boundary facet $f \in \mathcal{F}^b$, $f \in \partial c$, denote $\mathcal{C}(f)$ the set of facets $\tilde{f} \in \partial \tilde{c}$ in contact with $f$, that is, facets which verify

$$u_{\tilde{f}\tilde{f}} = (\mathbf{u}_{\tilde{c}} + \mathbf{x}_{\tilde{f}} - \mathbf{u}_c - \mathbf{x}_{\tilde{f}}) \cdot \mathbf{n}_{\tilde{f}\tilde{f}} < 0 \text{ and } |f \cap \tilde{f}| = 0, \quad (30)$$

where $|f \cap \tilde{f}|$ is the area of the intersection of the orthogonal projections of $f$ and $\tilde{f}$ (respectively displaced by $\mathbf{u}_c$ and $\mathbf{u}_{\tilde{c}}$) along their common normal $\mathbf{n}_{\tilde{f}\tilde{f}}$. In order to extend (26) to multiple particles, we propose the following expression for the face displacement reconstruction:

$$\mathbf{u}_f = \mathbf{u}_c + \frac{1}{2} \sum_{\tilde{f} \in \mathcal{C}(f)} \frac{|f \cap \tilde{f}|}{|f|} (\mathbf{u}_{\tilde{c}} + \mathbf{x}_{\tilde{f}} - \mathbf{u}_c - \mathbf{x}_f) \cdot \mathbf{n}_{\tilde{f}\tilde{f}} \mathbf{n}_{f\tilde{f}}. \quad (31)$$

The rationale behind the expression is threefold. First, note that we get the same expression as (26) when the cardinal of $\mathcal{C}(f)$ is 1, as well as in the case when all faces $\tilde{f} \in \mathcal{C}(f)$ have the same value for $u_{\tilde{f}\tilde{f}}$: the expression is not influenced by the number of facets but only by the average of penetration weighted by facet overlap. Second, when facets have negligible overlap $f \cap \tilde{f}$, the contribution of their mutual contact is also negligible. This leads us to the third point, which is crucial for energy conservation: the reconstruction $\mathbf{u}_f$ is continuous with regards to the $(\mathbf{u}_c)_{c \in \mathcal{C}}$. This can be seen by the fact that the only potential discontinuity can occur when the face displacement reconstruction switches from “non-contact” to “contact” status, which is only possible when the normal penetration $u_{\tilde{f}\tilde{f}}$ is null or when the face overlap $|f \cap \tilde{f}|$ is null. In both cases, the switch from “non-contact” reconstruction $\mathbf{u}_f = \mathbf{u}_c$ to “contact” reconstruction (31) is seamless since the additional term is 0.

Remark that reconstruction (31) introduces an additional geometrical nonlinearity in the expression of the energy (14), since the reconstruction tensors $\mathbf{T}_{cf}$ now depend implicitly on the displacement of cells $c$ and $\tilde{c}$ through the term $|f \cap \tilde{f}|$:

$$\mathbf{T}_{cf} = \mathbf{I} - \frac{1}{2} \sum_{\tilde{f} \in \mathcal{C}(f)} \frac{|f \cap \tilde{f}|}{|f|} \mathbf{n}_{\tilde{f}\tilde{f}} \otimes \mathbf{n}_{f\tilde{f}} \quad (32)$$

$$\forall \tilde{f} \in \mathcal{C}(f), \tilde{f} \in \tilde{c}, \quad \mathbf{T}_{\tilde{c}\tilde{f}} = \frac{|f \cap \tilde{f}|}{2|f|} \mathbf{n}_{\tilde{f}\tilde{f}} \otimes \mathbf{n}_{f\tilde{f}}. \quad (33)$$

This gives an additional force $\mathbf{F}_c^G$ in the case of sliding faces in contact, when
$f \in \partial c$ is a boundary facet of $c$:

$$ \mathbf{F}_c^G = -\frac{1}{2} \sum_{\tilde{f} \in \mathcal{C}(f)} |f| u_{f\tilde{f}}(\mathbf{n}_{f,c} \cdot \boldsymbol{\sigma}_c \mathbf{n}_{f\tilde{f}}) \frac{\partial |f \cap \tilde{f}|}{\partial \mathbf{u}_c} \cdot \mathbf{n}_{f\tilde{f}} \\ -\frac{1}{2} \sum_{\substack{\tilde{f}, f \in \mathcal{C}(\tilde{f}) \\ \tilde{c} \ni \tilde{f}}} |f| u_{f\tilde{f}}(\mathbf{n}_{\tilde{f},\tilde{c}} \cdot \boldsymbol{\sigma}_{\tilde{c}} \mathbf{n}_{f\tilde{f}}) \frac{\partial |f \cap \tilde{f}|}{\partial \mathbf{u}_c} \cdot \mathbf{n}_{f\tilde{f}}. \quad (34) $$

The expression of the derivative $\frac{\partial |f \cap \tilde{f}|}{\partial \mathbf{u}_c}$ is lengthy and will not be given here, but can be computed easily, as the position of each face $f$ and $\tilde{f}$ is a linear function of $\mathbf{u}_c$ and $\mathbf{u}_{\tilde{c}}$ respectively, so that the derivative is constant by parts.

## 4 Time discretization

We use the explicit second-order accurate pseudo-energy conservative MEMM scheme developed in [20] for Hamiltonian systems. The main advantages of the scheme is that there is a pseudo-energy (with calculable expression) which is exactly conserved, even for non-linear systems. The scheme writes, for $c \in C$:

$$ \mathbf{u}_c^{n+1} = \mathbf{u}_c^n + \Delta t \mathbf{v}_c^{n+1/2}, \quad (35a) $$

$$ \mathbf{v}_c^{n+3/2} = \mathbf{v}_c^{n-1/2} + \frac{2}{\rho|c|} \int_0^{\Delta t} \mathbf{F}_c(\mathbf{u}_h^n + \tau \mathbf{v}_h^{n+1/2}) \mathrm{d}\tau, \quad (35b) $$

where $\mathbf{v}_c$ denotes the velocity of cell $c$, and we give explicitly the dependency of force $\mathbf{F}_c$ in terms of the current displacements of cells at time $t^n + \tau$.

The exact pseudo-energy conservation is achieved only when the force integral in (35b) is exactly evaluated. A high-order quadrature in time manages to recover exact pseudo-energy conservation when the forces are polynomial, since the movement of particles $\mathbf{u}_h^n + \tau \mathbf{v}_h^{n+1/2}$ is affine in time. Smooth nonlinear forces can also be accurately captured up to machine-error accuracy with relatively low-order quadratures (Gauss-Legendre 2 or 3) [20]. However, in the case of contact, the forces are not smooth since they can be discontinuous (we have managed to keep a continuous energy though, see Section 3.4.3). As a consequence, even high-order quadratures fail to converge fast to machine-error accuracy in pseudo-energy conservation.

We resort to a modification of the MEMM scheme [20] to carry out the exact evaluation of the integral of forces in (35b). Remark that the forces are smooth except when contact events occur, where there are force discontinuities. These contact events can be of two different natures: either pairs of facets switch from
non-penetration to normal penetration ($u_{f\tilde{f}} = 0$), or the facet overlap $|f \cap \tilde{f}|$ is not differentiable, which occurs when a vertex of $f$ crosses an edge of $\tilde{f}$. These events are in finite number $N_c$, and given the displacement $\mathbf{u}_h^n$ and velocity $\mathbf{v}_h^{n+1/2}$, it is possible to compute the set of times $(\tau_k)_{k \in \{1, \dots, N_c\}}$ such that $\mathbf{u}_h^n + \tau_k \mathbf{v}_h^{n+1/2}$ verifies one of the two contact events condition. The penetration condition $u_{f\tilde{f}}$ and the facet positions in $f \cap \tilde{f}$ being affine in terms of the cell displacements $\mathbf{u}_h$, these are linear problems in $\tau$ for each pair of faces $(f, \tilde{f})$ which change contact status during the time-step.

Given the set of contact event times $(\tau_k)_{k \in \{1, \dots, N_c\}}$, and setting $\tau_0 = 0$ and $\tau_{N_c+1} = \Delta t$, it is easy to carry out an exact or quasi-exact quadrature on each interval $(\tau_k, \tau_{k+1})$, for $k \in \{0, \dots, N_c\}$, since the forces are smooth on each interval. In the present work, we consider linear elasticity, so that the forces are affine on each sub-interval, and a mid-point quadrature is sufficient to obtain exact conservation.

Note that the proposed modification is different from doing a sub-time-stepping matching contact events: cell displacement $\mathbf{u}_c$ and velocity $\mathbf{v}_c$ is not updated on each sub-interval $(\tau_k, \tau_{k+1})$ and remains constant through the whole time-step. To the best of our knowledge, sub-time stepping is difficult in practical and has not been carried out. Remark also that it is not necessary to do a global sub-integration for all cells $c \in \mathcal{C}$: each integral is evaluated independently in (35b) and the only cells $c$ which require sub-integration are those for which $T_{cf} \neq 0$, with $f$ a facet involved in a contact event. A further refinement of the implementation consists in detecting contact events for facets and updating the set of contact events $\tau_{c,k}$ for each cell $c$ involved in the displacement reconstruction on the facet.

Additionally, note that the present scheme has the advantage that the external boundary facet degrees of freedom have no mass. Only cell degrees of freedom have mass, so that contact boundary conditions can implement the Signorini complementary condition on reconstructions without modifying mass-endowed cell degrees of freedom. The altered reconstructions only impact the computation of forces. In [15], the authors use the Finite Element method and test the penetration with edge nodes which have a mass. For this reason, the modification of the mass matrix is necessary for the treatment of contact which conserves the energy. It should be noted that the authors of [15] propose to remove mass from the degrees of freedom in contact. This is a key ingredient for the convergence of contact methods with Finite Element methods [8].
# 5 Numerical results

## 5.1 One-dimensional benchmarks: elastic bar rebound

We first compare the results obtained with the present scheme with well-established discretizations of the dynamic Signorini problem, reviewed in [9]. The article considers two benchmarks: benchmark 1 consists in the impact of an elastic bar on a rigid surface and targets the absence of spurious oscillations in the displacement and contact pressure on the contact zone; benchmark 2 consists in multiple rebounds of an elastic bar on a rigid surfaces and aims at evaluating energy conservation.

For benchmark 1, we consider an elastic bar with initial velocity $v_0 = -10$, length $L = 10$, density $\rho = 1$, Young's modulus $E = 900$ and initial position of the beam tip at $h_0 = 5$. The rigid obstacle sits at $x = 0$. The bar has a rigid body motion until it collides with the obstacle, stays in contact while waves travel through it, and bounces back as the waves reach the contact zone. An exact solution for displacement and contact pressure is given in [9]. We follow the same discretization size as [9] with $\Delta x = 0.1$ and use the same CFL condition number $\nu_c = c_0 \frac{\Delta t}{\Delta x} = 0.75$, where $c_0 = \sqrt{\frac{E}{\rho}}$ is the wave velocity. Unless stated otherwise, the time-step is $\Delta t = 0.0025$.

Figure 5 shows the displacement and contact pressure at the tip of the beam for schemes A and B, with a comparison with the exact solution. We observe that the behavior of schemes A and B is extremely similar both on displacement and contact pressure, except the fact that after bouncing off, scheme A accurately gives a zero contact pressure while small pressure oscillations remain for scheme B. Both schemes present mild pressure oscillations and no displacement oscillations after contact. These can be explained by the Gibbs phenomenon near pressure discontinuity. These are much milder than the implicit schemes presented in [9], with the exception of modified mass methods [15]. They compare with explicit or semi-explicit results, with more marked oscillations probably due to the particular space discretization used here instead of classical linear Finite Element methods.

In Figure 6, the proposed contact event detecting scheme is compared with higher-order quadratures of the force integral in (35b). We plot the absolute value of the pseudo-energy variation in logscale against time for the present scheme (with midpoint quadrature on each subinterval), and Gauss-Legendre (GL) quadratures of increasing order: 1 point (midpoint, order 2), 2 points (order 4), 3 points (order 6) and 5 points (order 10). It can be observed that while the present scheme preserves the pseudo-energy to $10^{-14}$ through the simulation, all other quadratures exhibit a jump when the beam lifts off at $t = 0.3$s. In addition, an increase of the order of quadrature does not result in a crucial improvement of energy variation,
<figure><img src="image_7.png" /><figcaption>Figure 5: Impact of an elastic bar, benchmark 1: Displacement (left) and contact pressure (right) at the beam tip.</figcaption></figure>

which stays of the order of $10^{-9}$. It is important to note that the number of force evaluations in the present scheme is roughly equal to that of the midpoint quadrature (except two additional evaluations when the beam comes into contact and lifts off), whereas the Gauss-Legendre quadratures' cost is proportional to the number of quadrature points used. The proposed scheme is therefore as fast a midpoint quadrature while achieving almost exact conservation.

Using now the proposed scheme, we compare the convergence of the numerical solution to the exact solution for schemes A and B, which differ on the reconstruction of boundary facet dofs, in Figure 7. It can be observed that, while scheme A is supposed to have a higher order than scheme B, the convergence curves are comparable. It should be noted that in the present case, the order 0.7 of convergence is suboptimal, due to the nonregularity of the solution. Since we are interested in contact problems for which the solution is in general nonsmooth, the reduction of order introduced by scheme B cannot be considered as redhibitory.

For benchmark 2, we consider the same elastic beam, but change the initial velocity to $v = 0$ and add a gravity acceleration $g_0 = 10$, and simulate the system over time $T = 20$. The beam has multiple bounces with the same amplitude as the mechanical energy is conserved. With the adequate choice of parameters, again, an exact solution for displacement can be found in [9].

Figure 8 shows the numerical results for schemes A and B with a discretization step $\Delta x = 0.1$ and $\Delta t = 0.0025$. We observe that both treatments of the boundary condition give similar results, with a good capture of the first four bounces and a progressively degrading solution as the simulation advances. The amplitude of the bounces diminishes, in the same way as most methods presented in [9], with the exception of the semi-explicit scheme dealing with the impact (scheme 4.7), the implicit Newmark scheme with penalty contact condition (scheme 5.1) and the schemes with modified mass (schemes 7.1 and 7.2). It should however be noted
<figure><img src="image_8.png" /><figcaption>Figure 6: Energy error for the present contact event detecting scheme, midpoint, Gauss-Legendre 2, 3 and 5 points integrators.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Convergence in $L^\infty$ norm in time and in space for schemes A and B.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Bounces of an elastic bar, benchmark 2: Displacement for schemes A and B, $\Delta x = 0.1$ and $\Delta t = 0.0025$.</figcaption></figure>

that while an energy error analysis explains the evolution of bounce amplitude in [9], we have exact conservation of the pseudo-energy to $10^{-13}$. We assume that the discrepancy comes from the space discretization accuracy rather than the time discretization scheme.

To test this hypothesis, we simulate benchmark 2 with 5 successive space discretization sizes $\Delta x = L/N$, with $N = 100, N = 200, N = 400, N = 800$ and $N = 1600$, and fix a common time-step $\Delta t = 0.0003125$ given by the stability condition on the finer mesh. The time-step is small enough for the time discretization error to be dominated by the space discretization error for coarse meshes. Figure 9 shows the results obtained for schemes A and B. We observe that both schemes give very similar results. Even with a very small time-step, the results on the coarser grid $N = 100$ are not significantly improved compared to Figure 8. As the mesh is refined, the numerical results converge to the exact solution, without any time-step refinement. This shows that energy is in fact not dissipated, but stored in local oscillations which decrease the amplitude of the tip bounce. These oscillations are diminished when refining the mesh.

## 5.2 Static Hertz's contact

A disc of radius 20 cm is considered with a contact boundary $\Gamma_c$ restricted to the lower part of the boundary ($y < 20$ cm). On $\Gamma_c$, we apply a homogeneous Neumann condition and a vertical density of volume forces of 2 MN/m³ with no
<figure><img src="image_11.png" /><figcaption>Figure 9: Bounces of an elastic bar, benchmark 2: Displacement for scheme A (left) and scheme B (right) at the beam tip for different space discretization sizes.</figcaption></figure>

Dirichlet condition ($\Gamma^D = \emptyset$). Homogeneous isotropic linear elasticity in plane strain approximation is considered with a Young modulus fixed at $E = 25\,\text{MPa}$ and $P = 0.25$. Hertz' test-case consists in finding the static solution to this problem. We refer to [3] and [4] for reference Finite Element calculations.

To compute the static solution, we use a standard BFGS method to minimize potential energy [10]. We carry out simulations with mesh sizes $h = 1.5 \times 10^{-2}$, $h = 7.5 \times 10^{-3}$, $h = 3.75 \times 10^{-3}$ and $h = 1.875 \times 10^{-3}$. The results are presented in Figure 10, and can be compared with the reference solution in Figure 11. Note that the present results have the same contact area and vertical displacement as the reference solution, which are the observable characteristics for this problem. The convergence of Von Mises stress is obviously less accurate, since stresses are approximated with piecewise constant values on each cell in the present scheme, as opposed to second-order polynomials in the Finite Element method [3, 4]. The accuracy of the method is nonetheless satisfactory.

## 5.3 Contact between two disks

We consider two identical disks of radius 0.1m centered at (0.14, 2.27) m (for the left disk) (resp. at (0.36, 2.27) m (for the right disk)). The contact boundary is restricted to the right part of the boundary of the disk placed on left hand (resp. the left part of the boundary of the disk placed on the right hand). We assume that the two disks are initially moving with a constant velocity $\vec{V} = 0.3\vec{e}_x \, \text{m.s}^{-1}$ at left when $x < 0.25 \, \text{m}$ (resp. $\vec{V} = -0.3\vec{e}_x \, \text{m.s}^{-1}$ at right when $x > 0.25 \, \text{m}$) and without rotation. The contact axis is considered at $x = 0.25 \, \text{m}$. The different parameters used for the model are a Poisson ratio $\nu = 0.2$ and a Young's modulus fixed at $E = 100 \, \text{Pa}$. The solid density is given by $\rho_s = 100 \, \text{kg.m}^{-1}$. The CFL condition is taken as 0.2. The deformations are larger for the linear elasticity model to remain
<figure><img src="image_12.png" /><figcaption>Figure 10: 2D deformed configurations with mesh size $h = 1.5 \times 10^{-2}$ (top left), $h = 7.5 \times 10^{-3}$ (top right), $h = 3.75 \times 10^{-3}$ (bottom left), and $h = 1.875 \times 10^{-3}$ (bottom right).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Reference solution with FEM, $h = 1.5 \times 10^{-3}$ [3, 4], color plot of Von Mises stress.</figcaption></figure>

valid.

In Figure 12, structural deformation including contact and the normal stress distribution are illustrated at different times. The simulation time is $t = 1$ s. The first contact occurs at time $t = 0.04$ s. The deformation of the solid becomes too large with increasing time. At time $t = 0.19$ s, we notice a remarkable deformation of the solid and a concentration of normal stress close to the contact axis. We observe a wave which propagates in the solids. At time $t = 0.25$ s and $t = 0.26$ s, the contact (resp. normal stress and deformation) reaches the maximum. Beyond time $t = 0.27$ s, the solids bounce and go in the opposite direction. Figure 12 at time $t = 0.45$ s and $t = 0.58$ s shows the detachment of solids and their displacement in the opposite direction of their initial velocity. We still observe the propagation of waves in the material with a deformation linked to the compression waves generated by the contact phenomenon and which gradually decreases over time.

in Figure 13, we present the relative energy conservation error, computed as the difference between the initial energy and the energy at the different time-steps for the solid. The conservation originates from the modified time-integration scheme of the solid presented in Section 4 which ensures the exact conservation of energy. The relative conservation error on energy observed is extremely low (as $10^{-15}$ of the energy exchange in the system).
<figure><img src="image_14.png" /><figcaption>Figure 12: Normal stress distribution in the disks with periodic boundary conditions at times 0.07s, 0.13s, 0.19s, 0.22s, 0.45s and 0.58s.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13: Relative conservation error on energy.</figcaption></figure>

## 5.4 Billiard test-case

This test-case intends to demonstrate the efficiency and capability of the present method to manage contact between several deformable structures with multiple rebounds. The position of the billiard table is fixed and corresponds to the domain $(x, y) \in [2, 10] \times [2, 10]$ m with a thickness of 1m. The balls are 7 disks of radius 0.5m, with initial center positions (4, 5), (5, 4), (8, 4), (7, 6), (5, 7), (4, 8) and (7, 8). We assume that the disks are moving with an initial velocity $\vec{V}$ such that $\vec{V} = 0.1\vec{e}_y \cdot m \cdot s^{-1}$ for the disks with initial center position $(x, y)$ with $x > y$, and $\vec{V} = 0.1\vec{e}_x \cdot m \cdot s^{-1}$ otherwise. The solids density and Young's modulus are, respectively, $\rho_s = 100 \text{ kg} \cdot \text{m}^{-3}$ and $E = 100 \text{ Pa}$. The Poisson ratio is taken $\nu = 0.2$. The total simulation time is $T = 50$ s.

At time $t = 11$s, we observe in Figure 14 the first contact between two couples of billiard balls. We observe, for the particles in contact, an increase in normal stress. At time $t = 16$s, the first contact with the billiard table occurs, we observe the propagation of shock waves which generates the deformation of the right part of the table. This deformation increases with time and spreads over a wider area. Note that, beyond that point, multiple contact events occur almost simultaneously between balls, and some balls are even in contact with two other balls, for instance at $t = 30.5$s. The present scheme is still capable of dealing with multi-contact between structures while ensuring exact conservation of energy as shown in Figure
# 6 Conclusion

In this work, we treated the contact between deformable solids by using a Discrete Elements method coupled with an explicit time integration method which preserves the pseudo-energy at contact time. A particular importance in these problems returns to the conservation of mass, momentum and especially the conservation of energy of the system. Taking into account the contact between two deformable bodies generates complexities in numerical modeling and simulation in solid mechanics. This is mainly related to the non-linearity of the boundary conditions on part of the boundary and to the non-regularity of the simulated physical principles. The challenge was to prove the conservation of energy when contact takes place between two deformable solids. We proposed a modification of the time integration scheme proposed in [20] and adapted the scheme to obtain the conservation of the pseudo-energy at the moment of contact between the two solids. The first complementary development would be the transition to three-dimensional contact, which presents many challenges. To follow up on this work, another perspectives is to add plasticity and friction to our scheme and to couple them with the Discrete Element methods.

# Acknowledgements

The authors would like to thank A. Ern (ENPC, EPC SERENA) for fruitful discussions.

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-17-CE39-0013 (project PRECIS).

# References

[1] Ted Belytschko and Mark O Neal. Contact-impact by the pinball algorithm with penalty and Lagrangian methods. *International Journal for Numerical Methods in Engineering*, 31(3):547-572, 1991.

[2] Stefan Bilbao, Michele Duceschi, and Fabiana Zama. Explicit exactly energy-conserving methods for Hamiltonian systems. arXiv preprint, 2022.

[3] Franz Chouly, Mathieu Fabre, Patrick Hild, Rabii Mlika, Jérôme Pousin, and Yves Renard. An overview of recent results on Nitsche's method for contact
<figure><img src="image_16.png" /><figcaption>Figure 14: Deformation and multi-contact between solids with color plot of Von Mises stresses at 11 s, 16 s, 20 s, 24 s, 30.5 s and 33 s.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15: Relative conservation error on energy.</figcaption></figure>

problems. In *Geometrically unfitted finite element methods and applications*, pages 93–141. Springer, 2017.

[4] Franz Chouly, Patrick Hild, and Yves Renard. Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments. *Mathematics of Computation*, 84(293):1089–1112, 2015.

[5] Franz Chouly, Rabii Mlika, and Yves Renard. An unbiased Nitsche's approximation of the frictional contact between two elastic structures. *Numerische Mathematik*, 139(3):593–631, 2018.

[6] Peter Deuflhard, Rolf Krause, and Susanne Ertel. A contact-stabilized Newmark method for dynamical contact problems. *International Journal for Numerical Methods in Engineering*, 73(9):1274–1290, 2008.

[7] David Doyen. *Méthodes numériques pour des problèmes dynamiques de contact et de fissuration*. PhD thesis, Université Paris-Est, Université Paris-Est, 2010.

[8] David Doyen, Alexandre Ern, et al. Convergence of a space semi-discrete modified mass method for the dynamic Signorini problem. *Communications in Mathematical Sciences*, 7(4):1063–1072, 2009.
[9] David Doyen, Alexandre Ern, and Serge Piperno. Time-integration schemes for the Finite Element dynamic Signorini problem. *SIAM J. Sci. Comput.*, 33(1):223–249, 2011.

[10] Roger Fletcher. *Practical methods of optimization*. John Wiley & Sons, 2013.

[11] O. Gonzalez. Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Comput. Methods Appl. Mech. Eng.*, 190:1763–1783, 2000.

[12] Patrice Hauret. *Méthodes numériques pour la dynamique des structures non-linéaires incompressibles à deux échelles*. PhD thesis, Ecole Polytechnique, Ecole Polytechnique, 2004.

[13] Patrice Hauret and Patrick Le Tallec. Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact. *Computer methods in applied mechanics and engineering*, 195(37-40):4890–4916, 2006.

[14] Couro Kane, Eduardo A Repetto, Michael Ortiz, and Jerrold E Marsden. Finite element analysis of nonsmooth contact. *Computer methods in applied mechanics and engineering*, 180(1-2):1–26, 1999.

[15] Houari Boumediene Khenous, Patrick Laborde, and Yves Renard. Mass re-distribution method for finite element contact problems in elastodynamics. *European Journal of Mechanics-A/Solids*, 27(5):918–932, 2008.

[16] Noboru Kikuchi and John Tinsley Oden. *Contact problems in elasticity: a study of variational inequalities and finite element methods*. SIAM, 1988.

[17] TA Laursen and Vikas Chawla. Design of energy conserving algorithms for frictionless dynamic contact problems. *International Journal for Numerical Methods in Engineering*, 40(5):863–886, 1997.

[18] TA Laursen and GR Love. Improved implicit integrators for transient impact problems—geometric admissibility within the conserving framework. *International Journal for Numerical Methods in Engineering*, 53(2):245–274, 2002.

[19] Tod Alan Laursen. *Formulation and treatment of frictional contact problems using finite elements*. PhD thesis, Stanford University, Stanford University, 1992.

[20] Frédéric Marazzato, Alexandre Ern, Christian Mariotti, and Laurent Monasse. An explicit pseudo-energy conserving time-integration scheme for Hamiltonian dynamics. *Computer Methods in Applied Mechanics and Engineering*, 347:906–927, 2019.
[21] Frédéric Marazzato, Alexandre Ern, and Laurent Monasse. A variational discrete element method for quasistatic and dynamic elastoplasticity. *International Journal for Numerical Methods in Engineering*, 121(23):5295-5319, 2020.

[22] Jean Jacques Moreau. Application of convex analysis to some problems of dry friction. In *Trends in Applications of Pure Mathematics to Mechanics*, pages 263-280. Pitman, 1977.

[23] Laetitia Paoli and Michelle Schatzman. A numerical scheme for impact problems i: The one-dimensional case. *SIAM Journal on Numerical Analysis*, 40(2):702-733, 2002.

[24] Laetitia Paoli and Michelle Schatzman. A numerical scheme for impact problems ii: The multidimensional case. *SIAM journal on numerical analysis*, 40(2):734-768, 2002.

[25] J. C. Simo, N. Tarnow, and K. K. Wong. Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. *Comput. Methods Appl. Mech. Eng.*, 100:63-116, 1992.

[26] Peter Wriggers and Giorgio Zavarise. Computational contact mechanics. *Encyclopedia of computational mechanics*, 2004.