# Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems

Francisco Armero, Eva Petőcz

► To cite this version:

Francisco Armero, Eva Petőcz. Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems. Computer Methods in Applied Mechanics and Engineering, 1998, 158 (3-4), pp.269-300. 10.1016/S0045-7825(97)00256-9. hal-01435615

HAL Id: hal-01435615

https://hal.science/hal-01435615

Submitted on 14 Jan 2017

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
# Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems

Francisco Armero*, Eva Petőcz¹

*Structural Engineering, Mechanics and Materials, Department of Civil and Environmental Engineering, University of California, Berkeley, CA 94720, USA*

## Abstract

This paper presents the formulation of conserving time-stepping algorithms for frictionless dynamic contact of solids. A new class of finite element methods is proposed for the solution of these problems that exhibit the same conservation laws as the underlying continuum dynamical system. The proposed methods are based on a penalty regularization of the constrained contact problem, and lead to full conservation of the total energy of the system (including the regularization penalty potential) during persistent contact, and restoration of the original energy upon release. Both linear and angular momenta are conserved by the scheme. Furthermore, the newly developed methods have the ability to enforce the associated constraints in the velocity besides the impenetrability constraint in the displacements, while preserving the conservation/restoration properties of the final numerical scheme. A modification of these schemes is described that assures positive energy dissipation if desired (even in the highly nonlinear setting of contact/impact problems), leading to contact schemes with high-frequency energy dissipation. Representative numerical simulations are presented illustrating the performance of the proposed numerical schemes.

## 1. Introduction

The accurate modeling of contact interfaces in solids is one of the main difficulties in common engineering applications. Typical examples are crashworthiness analyses and the simulation of metal forming processes. See the contributions in [17] and [7], respectively, for recent accounts of these considerations. The experience accumulated in the past regarding the numerical analysis of contact problems indicates the inherent difficulty of their solution, the cause being not only the highly nonlinear nature of the problem, but also its *unilaterally constrained* character. The lack of robustness of current implicit methods that impose the contact constraint has led in the past to the consideration of explicit schemes for the numerical solution of contact problems. The difficulties in the enforcement of this constraint appear often as oscillations between contact and released states.

Additional difficulties arise when dynamic problems are considered. The limited conditional stability in time of explicit integration schemes appears as a clear drawback. Implicit schemes may be employed to recover better stability properties but, as it is well known, stable numerical schemes for linear problems may lose this property in the nonlinear context, leading to an unstable increase of the energy during the numerical simulations. Characteristic examples are the trapezoidal and mid-point rules, two energy conserving schemes for linear problems that may result in energy increase (and actual blow up of the computation) in nonlinear problems (see e.g. [18] for representative simulations). These drawbacks have led to the consideration of energy-momentum
conserving schemes that do not suffer of this limited (energy) stability properties, as described in [18,6,19], among others. We can anticipate that the presence of the high nonlinearity due to the contact constraint may lead to similar instabilities, as the simulations of Section 4 show. The goal of the research presented in this paper is the formulation of time-stepping algorithms that possess the desired temporal stability properties by controlling the evolution of the energy and that, at the same time, lead to a stable (non-oscillatory) enforcement of the contact constraints.

A complete account on the numerical analysis of contact problems until the late 1980s can be found in [13]. Finite element methods for dynamic contact problems can be found in [12,10,4,5] for explicit integrators, and [3,13] involving implicit integrators for frictionless and frictional problems, to cite just a few references. The recent works presented in [22,15,16] show the current interest in the formulation of more robust implicit algorithms for frictionless contact. See also the results presented in [1]. The robustness of the numerical scheme requires good stability properties in the limit conservative case, without relying in the physical dissipation introduced by frictional effects.

The approach proposed herein makes use of the properties of the continuum dynamical system for the formulation and analysis of new and more robust implicit time-stepping algorithms for contact problems. Assuming no external forces, the total energy, linear and angular momenta of a system of solids in frictionless contact are conserved. These conservation properties are introduced in the newly developed schemes by construction, thus leading not only to a better modeling of the physical system but also to improved numerical properties. The new schemes are second-order accurate and unconditionally (energy) stable even in the fully nonlinear finite strain range, as implied by the conservation of the total energy of the system. The (unilaterally) constrained problem is regularized via a penalty formulation. Both the constraint in the displacements and the associated constraint in the velocities are enforced in this manner. The (positive) energy corresponding to the penalty potential is taken into account in the evolution of the energy, leading to full restoration of the initial energy of the system of solids upon release (i.e. when the regularization potentials are inactive), while the energy never increases beyond its initial value during persistent contact. The total energy of the system (solids plus regularization potentials) is conserved at all times, leading to the unconditional (energy) stability of the numerical schemes. These properties are combined with full conservation of angular and linear momenta. Numerical experiments have shown that these improved stability properties lead to a superior numerical performance when compared to similar traditional schemes (like the second-order mid-point and trapezoidal rules). As noted above, high-order standard numerical schemes usually involve an artificial increase of the energy, which eventually leads to the actual blow-up of the numerical computation.

Fully energy conserving schemes are appropriate for the long-term simulations of the interactions of solids in contact, where the main interest is the accurate resolution of the configuration of the system in the long-term (and thus its energy content). On the other hand, short term simulations are employed for the study of high-velocity frontal impacts (e.g. a rod impacting a rigid wall), requiring then the resolution of a wide frequency spectrum. In fact, weak shocks (discontinuities in the velocity and strain) dominate completely the solution at these time scales. In these conditions, high frequency energy dissipation is a desired feature. We emphasize that the application of standard dissipative schemes developed typically for linear problems do not assure in general a positive energy dissipation in the numerical scheme, the cause again being the highly nonlinear nature of the contact problem. We propose herein a simple modification of the conservative schemes previously developed that accomplishes this feature, and whose dissipative properties can be proven rigorously.

An outline of the rest of the paper is as follows. Section 2 includes a complete description of the problem under consideration. The governing equations are summarized in Section 2.1, with the conservation laws for frictionless contact described in Section 2.2. The finite element implementation considered in this paper is described in Section 2.3. Next, the formulation of the new energy-restoring, momentum-conserving scheme proposed herein is described in detail in Section 3.1, when only the gap constraint in the displacement is enforced. Rigorous proofs of the conservation properties of the proposed methods as well as extensions involving high-frequency dissipation are described in detail. Section 3.2 considers the enforcement of the velocity gap constraint, arriving to a similar class of conserving algorithms. Representative numerical simulations are presented in Section 4 to assess the performance of the proposed methods. Section 5 includes some concluding remarks. Finally, the consistent linearization of the proposed methods is summarized in Appendix A.
## 2. Problem description

We describe in this section the problem and numerical simulation of dynamic contact of elastic bodies. Section 2.1 summarizes the governing equations. Section 2.2 describes the conservation laws associated to this Hamiltonian system. Finally, the finite element implementation of the governing equations is described in Section 2.3.

### 2.1. The governing equations

Consider the motion of two elastic bodies with a reference placement $\Omega^\alpha$ ($\alpha = 1, 2$), characterized by the deformations $\varphi^\alpha : \Omega^\alpha \times [0, T] \to \mathbb{R}^{n_{\dim}}$ ($n_{\dim} = 1, 2$ or 3). The results presented herein extend trivially to multi-body interactions, as well as to self-contact of solids. We identify the material particles of each solid with the reference coordinate $X \in \cup \Omega^\alpha \subset \mathbb{R}^{n_{\dim}}$. Let $x^\alpha := \varphi^\alpha(X, t)$ be the current placement of the material particle $X \in \Omega^\alpha$ of the solid $\alpha$ at time $t \in [0, T]$, for some time interval $T$.

Denote by $P^\alpha$ the nominal stresses (first Piola–Kirchhoff stresses) in each solid. The case of interest corresponds to two hyperelastic solids characterized by respective stored energy functions $W^\alpha(F^\alpha)$, where $F^\alpha = \text{Grad } \varphi^\alpha$, and

$$P^{\alpha} = \frac{\partial W^{\alpha}}{\partial F^{\alpha}}. \qquad (2.1)$$

By the principle of material frame indifference, the stored energy function is invariant under the action of the proper orthogonal group (the rotation group) $SO(n_{\dim})$, that is

$$W^{\alpha}(QF^{\alpha}) = W^{\alpha}(F^{\alpha}) \quad \forall Q \in SO(n_{\dim}). \qquad (2.2)$$

Considering a one-parameter group of rotations $Q(\eta)$ with $Q(0) = \mathbf{1}_{n_{\dim}}$ (the identity in $\mathbb{R}^{n_{\dim}}$), so

$$\left. \frac{d}{d\eta} Q(\eta) \right|_{\eta=0} = \hat{W} \in \text{so}(n_{\dim}), \qquad (2.3)$$

(where $\text{so}(n_{\dim})$ denotes the linear space of skew-symmetric tensors), and taking the derivative of (2.2) with respect to $\eta$ at $\eta = 0$, we obtain the relation

$$\frac{\partial W^{\alpha}}{\partial F^{\alpha}} F^{\alpha \top} : \hat{W} = 0 \quad \forall \hat{W} \in \text{so}(n_{\dim}), \qquad (2.4)$$

implying the symmetry of the Kirchhoff stress tensor

$$\tau^{\alpha} := \frac{\partial W^{\alpha}}{\partial F^{\alpha}} F^{\alpha \top} = \tau^{\alpha \top}. \qquad (2.5)$$

The symmetry relation (2.2), or equivalently, the relation (2.5), leads to the classical conservation law of the total angular momentum as discussed below. Furthermore, a classical argument (see e.g. [23]) leads then to the dependence of the stored energy function on the Green–Lagrange strain tensor $E = \frac{1}{2}(F^T F - \mathbf{1})$, i.e.

$$W^{\alpha}(F^{\alpha}) = \tilde{W}^{\alpha}(E^{\alpha}). \qquad (2.6)$$

As an example, the simulations presented in Section 4 consider the Saint–Venant Kirchhoff model, characterized by

$$W(F) = \tilde{W}(E) = \frac{1}{2} E : CE, \qquad (2.7)$$

where $C$ denotes the material secant tangent.

Let $V^\alpha := \dot{\varphi}^\alpha$ be the material velocity field of the solid $\alpha$, and $\rho^\alpha$ the corresponding reference density. The superimposed dot (·) refers to the (material) derivative with respect to time $t$. We denote by $\gamma^\alpha := \partial\varphi^\alpha(\Omega^\alpha)$ the boundary of the current configuration of solid $\alpha$, with its reference counterpart $\Gamma^\alpha = \varphi^{\alpha-1}(\gamma^\alpha) = \partial\Omega^\alpha$.
Denoting the current and reference boundaries in contact by $\gamma_c^\alpha$ and $\Gamma_c^\alpha := \varphi^{\alpha-1}(\gamma_c^\alpha)$, respectively, the weak form of the balance of linear momentum equations reads

$$ \sum_{\alpha=1}^{2} \left[ \int_{\Omega^\alpha} \rho^\alpha \dot{\mathbf{V}}^\alpha \cdot \delta\varphi^\alpha \, d\Omega + \int_{\Omega^\alpha} \mathbf{P}^\alpha : \mathrm{Grad}(\delta\varphi^\alpha) \, d\Omega \right] = \sum_{\alpha=1}^{2} \left[ \int_{\Omega^\alpha} \rho^\alpha \mathbf{b} \cdot \delta\varphi^\alpha \, d\Omega + \int_{\Gamma_t^\alpha} \bar{\mathbf{t}} \cdot \delta\varphi^\alpha \, d\Gamma \right] + \int_{\Gamma_t^\alpha} \mathbf{t} \cdot [\delta\varphi^1(\mathbf{X}) - \delta\varphi^2(\tilde{\mathbf{Y}}(\mathbf{X}))] \, d\Gamma, \quad (2.8) $$

for all admissible variations $\delta\varphi^\alpha : \Omega^\alpha \to \mathbb{R}^{n_{\mathrm{dim}}} (\alpha = 1, 2)$ such that

$$ \delta\varphi^\alpha|_{\Gamma_u^\alpha} = \mathbf{0}. \quad (2.9) $$

Here, $\Gamma_u^\alpha$ denotes the part of the reference boundary of solid $\alpha$ with imposed displacements, and $\Gamma_t^\alpha$ is the part of the reference boundary with imposed external tractions $\bar{\mathbf{t}}$. The specific body forces are denoted by $\mathbf{b}$. The decomposition

$$ \Gamma^\alpha = \overline{\Gamma_u^\alpha \cup \Gamma_t^\alpha \cup \Gamma_c^\alpha} \quad \text{with} \quad \Gamma_u^\alpha \cap \Gamma_t^\alpha \cap \Gamma_c^\alpha = \emptyset \quad (\alpha = 1, 2), \quad (2.10) $$

is assumed for a well-posed boundary value problem.

The vector $\mathbf{t}$ in the last term of (2.8) denotes the contact nominal traction between the solids along the contact boundary. For frictionless contact, this traction is given in terms of the (nominal) contact pressure $p \ge 0$ as

$$ \mathbf{t} = p\mathbf{v}, \quad (2.11) $$

with $\mathbf{v}$ denoting the unit outward normal to the current contact boundary $\gamma_c^2$. The contact pressure $p$ corresponds to the Lagrange multiplier imposing the unilateral contact constraint

$$ g(\mathbf{X}) := \mathbf{v} \cdot [\varphi^1(\mathbf{X}) - \varphi^2(\tilde{\mathbf{Y}}(\mathbf{X}))] \ge 0, \quad (2.12) $$

for the gap $g(\mathbf{X})$ of a particle $\mathbf{X} \in \Gamma^1$. In (2.8) and (2.12), the mapping $\mathbf{Y} = \tilde{\mathbf{Y}}(\mathbf{X}) \in \Gamma^2$ defines the closest-point projection of a material point $\mathbf{X} \in \Gamma^1$ on the contact surface to $\Gamma^2$ at the current configuration of the solids, that is

$$ \tilde{\mathbf{Y}}(\mathbf{X}) = \arg \min_{\mathbf{Y} \in \Gamma^2} \{\|\varphi^1(\mathbf{X}) - \varphi^2(\mathbf{Y})\|\}, \quad (2.13) $$

where $\|\cdot\|$ denotes the usual Euclidean vector norm. Fig. 1 depicts the geometrical construction behind the definition of the gap function $g(\mathbf{X})$ and the normal $\mathbf{v}(\tilde{\mathbf{Y}}(\mathbf{X}))$ in (2.12). In the continuum problem, the current contact boundaries are defined by the condition $g = 0$, that is, $\gamma_c^1 = \gamma_c^2 = \cap_\alpha \gamma_c^\alpha$ in this case, thus allowing the use of $\Gamma_c^1$ or, alternatively, $\Gamma_c^2$ in (2.8) without a preference for the role played by each particular body in the governing equations. In the discrete problem, the contact boundaries are approximated by the discrete enforcement of the gap constraint (2.12), as described in Section 2.3.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Definition of the gap function $g(\mathbf{X})$ and unit normal $\mathbf{v}$ through the closest-point projection mapping at the current configuration of the solids in contact.</figcaption></figure>
We note for future use the relation

$$ \varphi^1(X) - \varphi^2(\tilde{Y}(X)) = g(X)\nu(\tilde{Y}(X)) \qquad (2.14) $$

as a consequence of the definition (2.13) for the closest-point projection. The unilaterally constrained system under consideration is then completely characterized by

$$ p \ge 0, \quad g \ge 0, \quad pg = 0, \qquad (2.15) $$

the so-called Kuhn-Tucker conditions (see [20]).

During persistent contact, the time derivative of (2.12), which now holds as an equality, implies

$$ h := \dot{g} = \nu \cdot (V^1(X) - V^2(\tilde{Y}(X))) = 0, \qquad (2.16) $$

where we have made use of the property

$$ \dot{\nu} \cdot (\varphi^1(X) - \varphi^2(\tilde{Y}(X))) = 0, \qquad (2.17) $$

a consequence of the closest-point projection $Y = \tilde{Y}(X)$. Therefore, the velocity field is constrained by (2.16). After using the Kuhn-Tucker conditions (2.15), the condition (2.16) can be recast as

$$ pg = 0, \qquad (2.18) $$

commonly referred to as the persistency condition (see [20]).

## 2.2. The conservation laws

The system of nonlinear elastodynamics equations described in the previous section is a characteristic example of an infinite dimensional Hamiltonian system (see [21]). The consideration of contact states converts the system in a unilaterally constrained Hamiltonian system of evolution. The presence of symmetries like (2.2) leads to the conservation laws described in this section.

Consider the following standard definition of the total linear momentum

$$ L := \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \rho^{\alpha} V^{\alpha} d\Omega, \qquad (2.19) $$

and the total angular momentum

$$ J := \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \varphi^{\alpha} \times \rho^{\alpha} V^{\alpha} d\Omega, \qquad (2.20) $$

of the system of solids. The symbol $\times$ denotes the cross product of two vectors in $\mathbb{R}^3$ if $n_{\dim} = 3$, and their equivalent reductions (embeddings of $\mathbb{R}^{n_{\dim}} \subset \mathbb{R}^3$) in lower dimensions $n_{\dim} < 3$. Similarly, denote the total energy of the system by

$$ \mathcal{E} := \underbrace{\sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \frac{1}{2} \rho^{\alpha} V^{\alpha} \cdot V^{\alpha} d\Omega}_{\mathcal{H}} + \underbrace{\sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} W^{\alpha} d\Omega}_{\mathcal{W}} = \mathcal{H} + \mathcal{W}, \qquad (2.21) $$

for the total kinetic energy $\mathcal{H}$ and strain energy $\mathcal{W}$.

The case of interest for the analysis presented below corresponds to the homogeneous Neumann problem, characterized by no imposed boundary displacements and no external loading. In this case, the total energy $\mathcal{E}$, linear momentum $L$ and angular momentum $J$ are conserved as summarized in the following proposition.

**PROPOSITION 2.1.** Let $\Gamma_u^\alpha = \emptyset$ ($\alpha = 1, 2$), $\tilde{t} = 0$ and $b = 0$. Then, the linear momentum $L$, the angular momentum $J$, and the total energy $\mathcal{E}$ are constants of motion.

**PROOF.** The proof is based on classical arguments, and is included herein for completeness. The discrete counterpart presented in Section 3 follows closely the same arguments.

(i) Conservation of the linear momentum. Since $\Gamma_u^\alpha = \emptyset$, an admissible variation is obtained by
$$ \delta\varphi^\alpha = \mathbf{a} \quad \text{for } \alpha = 1, 2, \qquad (2.22) $$

with $\mathbf{a} \in \mathbb{R}^{n_{\text{dim}}} \text{ constant. Hence, } \text{Grad}(\delta\varphi^\alpha) \equiv 0 \text{ in this case. Using (2.8) with the admissible variations (2.22) and noting that } \tilde{\mathbf{t}} = \mathbf{0} \text{ and } \mathbf{b} = \mathbf{0} \text{ by assumption, we have}$

$$ \mathbf{a} \cdot \frac{d\mathbf{L}}{dt} = \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \rho^{\alpha} \dot{\mathbf{V}}^{\alpha} \cdot \mathbf{a} \, d\Omega = 0 \quad \forall \mathbf{a} \in \mathbb{R}^{n_{\text{dim}}}. \qquad (2.23) $$

Therefore, $d\mathbf{L}/dt = \mathbf{0}$ or equivalently $\mathbf{L}(t) = \mathbf{L}(0) = \text{constant. The conservation of linear momentum follows then from the invariance of the equations under the variations (2.22), i.e. the action of the linear (additive) group } \mathbb{R}^{n_{\text{dim}}} \text{ (spatial translations).}$

(ii) Conservation of the angular momentum. Similarly, we can consider the admissible variations defined by

$$ \delta\varphi^\alpha(\mathbf{X}) := \mathbf{w} \times \mathbf{x}^\alpha \quad \text{for } \alpha = 1, 2, \qquad (2.24) $$

where $\mathbf{w} \in \mathbb{R}^{n_{\text{dim}}} \text{ constant, and } \mathbf{x}^\alpha = \varphi^\alpha(\mathbf{X})$. Thus, we have

$$ \text{Grad}(\delta\varphi^\alpha) = \hat{\mathbf{W}}\mathbf{F}^\alpha, \qquad (2.25) $$

where $\hat{\mathbf{W}} \in \text{so}(n_{\text{dim}})$ is the skew-symmetric tensor with axial vector $\mathbf{w}$ (i.e. $\hat{\mathbf{W}}\mathbf{a} = \mathbf{w} \times \mathbf{a} \quad \forall \mathbf{a} \in \mathbb{R}^{n_{\text{dim}}}$). Writing (2.8) with the variations defined by (2.24), we obtain after making use of (2.20) and (2.25)

$$ \begin{aligned} \mathbf{w} \cdot \frac{d\mathbf{J}}{dt} &= \mathbf{w} \cdot \left[ \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} [\dot{\mathbf{x}}^{\alpha} \times \rho^{\alpha} \mathbf{V}^{\alpha} + \mathbf{x}^{\alpha} \times \rho^{\alpha} \dot{\mathbf{V}}^{\alpha}] d\Omega \right] \\ &= \mathbf{w} \cdot \left[ \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \mathbf{x}^{\alpha} \times \rho^{\alpha} \dot{\mathbf{V}}^{\alpha} d\Omega \right] = \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \rho^{\alpha} \dot{\mathbf{V}}^{\alpha} \cdot (\mathbf{w} \times \mathbf{x}^{\alpha}) d\Omega \\ &= - \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \mathbf{P}^{\alpha} : \hat{\mathbf{W}} \mathbf{F}^{\alpha} d\Omega + \int_{\Gamma_c^1} \mathbf{t} \cdot [\mathbf{w} \times \varphi^1(\mathbf{X}) - \mathbf{w} \times \varphi^2(\tilde{\mathbf{Y}}(\mathbf{X}))] d\Gamma \\ &= - \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \underbrace{\mathbf{P}^{\alpha} \mathbf{F}^{\alpha \top} : \hat{\mathbf{W}}}_{= 0 \text{ by (2.5)}} d\Omega + \int_{\Gamma_c^1} \underbrace{\mathbf{p} \mathbf{v} \cdot [\mathbf{w} \times (\varphi^1(\mathbf{X}) - \varphi^2(\tilde{\mathbf{Y}}(\mathbf{X}))) d\Gamma]}_{= \mathbf{g} \mathbf{v} \text{ by (2.14)}} \\ &= \int_{\Gamma_c^1} \underbrace{\mathbf{p} \mathbf{g} \mathbf{v} \cdot (\mathbf{w} \times \mathbf{v}) d\Gamma}_{= 0} = 0 \quad \forall \mathbf{w} \in \mathbb{R}^{n_{\text{dim}}}. \end{aligned} \qquad (2.26) $$

Therefore, $d\mathbf{J}/dt = \mathbf{0}$ or equivalently $\mathbf{J}(t) = \mathbf{J}(0) = \text{constant. The conservation of angular momentum follows then from the invariance of the equations under the variations (2.24) (infinitesimal rotations), i.e. the action of the rotation group } \text{SO}(n_{\text{dim}})$.

(iii) Conservation of energy. Finally, the evolution of the total energy is obtained using the weak equation (2.8) with the variations $\delta\varphi^\alpha = \mathbf{V}^\alpha$ and (2.21) as

$$ \begin{aligned} \frac{d\mathcal{E}}{dt} &= \sum_{\alpha=1}^{2} \left[ \int_{\Omega^{\alpha}} \rho^{\alpha} \dot{\mathbf{V}}^{\alpha} \cdot \mathbf{V}^{\alpha} d\Omega + \int_{\Omega^{\alpha}} \frac{\partial \mathbf{W}^{\alpha}}{\partial \mathbf{F}^{\alpha}} : \text{Grad}(\mathbf{V}^{\alpha}) d\Omega \right] \\ &= \int_{\Gamma_c^1} \mathbf{t} \cdot (\mathbf{V}^1(\mathbf{X}) - \mathbf{V}^2(\tilde{\mathbf{Y}}(\mathbf{X}))) d\Gamma = \int_{\Gamma_c^1} \mathbf{p} \dot{\mathbf{g}} d\Gamma = 0, \end{aligned} \qquad (2.27) $$

after using the persistency condition (2.18). Therefore, the total energy is conserved $\mathcal{E}(t) = \mathcal{E}(0)$ for all time $t$. $\square$

The goal of this paper is the design of time-stepping algorithms that possess these conservation laws.
## 2.3. Finite element implementation

The weak equation (2.8) is discretized in space through a standard isoparametric finite element formulation,

$$X = \sum_{A=1}^{n_{\text{node}}^{\alpha}} N^{\text{A}}(\xi) X_{\alpha}^{\text{A}} \in \Omega^{\alpha} \quad \text{and} \quad \varphi^{\alpha}(X) = X + \sum_{A=1}^{n_{\text{node}}^{\alpha}} N^{\text{A}}(\xi) d_{\alpha}^{\text{A}} \in \varphi^{\alpha}(\Omega^{\alpha}), \qquad (2.28)$$

($\alpha = 1, 2$) based on the shape functions $N^A: \square \to \mathbb{R}$ defined in the parent domain $\xi \in \square$ for $A = 1, n_{\text{node}}^{\alpha}$, the number of nodes for solid $\alpha$, with references coordinates $X_{\alpha}^{\text{A}} \in \mathbb{R}^{n_{\text{dim}}^{\alpha}}$. The nodal displacements $d_{\alpha}^{\text{A}} \in \mathbb{R}^{n_{\text{dim}}^{\alpha}}$ ($A = 1, n_{\text{node}}^{\alpha}$) are grouped in $d \in \mathbb{R}^{n_{\text{eq}}}$ where $n_{\text{eq}} = n_{\text{dim}} \times n_{\text{node}}$, with $n_{\text{node}} = \sum_{\alpha} n_{\text{node}}^{\alpha}$ being the total number of nodes.

Following a standard procedure, the above interpolations lead to the semi-discrete system of equations

$$\left. \begin{aligned} \dot{d}(t) &= M^{-1}p(t), \\ \dot{p}(t) &= -f_{\text{int}}(d(t)) + f_c(d(t)) + f_{\text{ext}}, \end{aligned} \right\} \qquad (2.29)$$

where we have introduced the nodal (linear) momenta

$$p := Mv, \quad \text{with} \quad v := \dot{d}(t), \qquad (2.30)$$

as an intermediate variable. Here, $M$ is the mass matrix defined by the standard assembly procedure $M = A_{e=1}^{n_{\text{el}}} M^e$ of the elemental mass matrices $M^e$ ($n_{\text{el}}$ = total number of elements). For an element with $n_{\text{en}}$ nodes, we have

$$M^e = \begin{bmatrix} M_{11} & \mathbf{1}_{n_{\text{dim}}} & \cdots & M_{1n_{\text{en}}} & \mathbf{1}_{n_{\text{dim}}} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ M_{n_{\text{en}}1} & \mathbf{1}_{n_{\text{dim}}} & \cdots & M_{n_{\text{en}}n_{\text{en}}} & \mathbf{1}_{n_{\text{dim}}} \end{bmatrix}, \qquad (2.31)$$

where $\mathbf{1}_{n_{\text{dim}}}$ is the rank-two identity matrix in $\mathbb{R}^{n_{\text{dim}}}$, and the mass coefficients $M_{AB}$ are given by the usual expression

$$M_{AB} = \int_{\Omega_e^\alpha} \rho^\alpha N^A N^B \, d\Omega. \qquad (2.32)$$

In Section 3.2 we make use of the lumped mass matrix obtained, for instance, by the standard row-sum technique

$$M_{AB} = M_A \delta_{AB} \quad (\text{no sum}), \quad \text{where} \quad M_A := \int_{\Omega_e^\alpha} \rho^\alpha N^A \, d\Omega, \qquad (2.33)$$

for the element $\Omega_e^\alpha$.

The external force $f_{\text{ext}} \in \mathbb{R}^{n_{\text{eq}}}$ corresponds to the contributions of the volumetric external force $b$ and imposed external tractions $\bar{t}$. The internal force vector $f_{\text{int}} \in \mathbb{R}^{n_{\text{eq}}}$ corresponds to the stress-divergence term in the continuum, and is given by the usual expression

$$f_{\text{int}} = \sum_{\alpha=1}^{2} \int_{\Omega_e^\alpha} \mathbf{B}_t^\text{T} \tau_t \, d\Omega, \qquad (2.34)$$

for the linearized strain operator $B_t$, with

$$B_t \delta d = \nabla_t^s \delta u = \text{sym}[\text{Grad}[\delta u] F_t^{-1}]$$

for the displacement field $u := \varphi_t(X) - X$. Here, the subscript $t$ refers to the configuration at time $t$.

The force of contact $f_c$ in (2.29) is obtained with the use of the now standard master/slave data structure (see [10] for details). In this context, $S$ denotes the slave node in contact with a master surface, at a point located in a master surface element defined by nodes $\{M1, M2, \dots\}$; see Fig. 2. Thus, we can assign two or more master nodes (belonging to the same master surface element) to each slave node in contact, thus establishing a contact element. The slave and master surfaces approximate the contact boundaries $\Gamma_c^1$ and $\Gamma_c^2$, respectively. Double
<figure><img src="image_3.png" /><figcaption>Fig. 2. A slave node contacts a master segment consisting of four master nodes.</figcaption></figure>

pass schemes can be easily accommodated to avoid the bias associated to a particular ordering of the two surfaces (see [10]).

The force of contact $f_c$ is then expressed as

$$f_c = \mathbf{A}_{s=1}^{n_{\text{slave}}} \hat{f}_{s,c}, \quad \text{with} \quad \hat{f}_{s,c} = p_s \hat{G}_s, \qquad (2.35)$$

where $\mathbf{A}_{s=1}^{n_{\text{slave}}}$ denotes the assembly over the $n_{\text{slave}}$ slave nodes/master segment pairs, and

$$\hat{G}_s = \begin{bmatrix} \mathbf{v}_s \\ -N^{M1}(\xi_s)\mathbf{v}_s \\ \vdots \\ -N^{M2}(\xi_s)\mathbf{v}_s \\ \vdots \end{bmatrix} \in \mathbb{R}^{(1+n_{\text{master}})^2 \times n_{\text{dim}}}, \qquad (2.36)$$

where $n_{\text{master}}^s$ is the number of master nodes in the master segment in contact with the slave node S. The normal component of the contact force (= nominal contact pressure × nominal contact area) at the slave node S has been denoted by $p_s$ in (2.35) and what follows, with a slight abuse of notation given the symbol $p$ employed in the previous developments for the nominal contact pressure in the continuum.

In (2.36), $N^{MI}(\xi_s)$ denotes the standard shape function of node MI in the master segment at the point of contact $\xi_s$ with normal $\mathbf{v}_s$, obtained by the closest-point projection mapping as in Eq. (2.13); see Fig. 2. The discrete counterpart of (2.14) holds as

$$\mathbf{x}^S - \sum_{l=1}^{n_{\text{master}}^s} N^{MI}(\xi_l) \mathbf{x}^{MI} = g(X^S) \mathbf{v}_s, \qquad (2.37)$$

where $\mathbf{x}^S = \boldsymbol{\phi}^1(X^S)$ and $\mathbf{x}^{MI} = \boldsymbol{\phi}^2(X^{MI})$ are the current positions of the slave and master nodes, respectively. We note that

$$\sum_{l=1}^{n_{\text{master}}^s} N^{MI}(\xi_l) = 1, \qquad (2.38)$$

at any point $\xi$ of the master segment.

For later use, we introduce the notation

$$\hat{\mathbf{d}}_s := \begin{bmatrix} \mathbf{d}^S \\ \mathbf{d}^{M1} \\ \mathbf{d}^{M2} \\ \vdots \end{bmatrix} \quad \text{and} \quad \hat{\mathbf{v}}_s := \begin{bmatrix} \mathbf{v}^S \\ \mathbf{v}^{M1} \\ \mathbf{v}^{M2} \\ \vdots \end{bmatrix}, \qquad (2.39)$$

that is, objects denoted by $(\hat{\cdot})_s$ refer to individual slave nodes/master segment pairs. The simulations presented in Section 4 consider linear master elements consisting of two master nodes (i.e. $n_{\text{master}}^s = 2$). The subscript $\alpha$
referring to the solid under consideration has been omitted in (2.37), (2.39) and the following developments, since it should be clear from the different role played by the slave and master nodes.

### 2.3.1. Temporal discretization of the continuum contributions

We consider a mid-point approximation of (2.29) and (2.30) leading to the discrete equations

$$ \left. \begin{aligned} \frac{1}{\Delta t} (\mathbf{d}_{n+1} - \mathbf{d}_n) &= \mathbf{v}_{n+1/2}, \\ \frac{1}{\Delta t} \mathbf{M}(\mathbf{v}_{n+1} - \mathbf{v}_n) &= -\mathbf{f}_{\text{int}}^{(n+1/2)} + \mathbf{f}_{\text{c}}^{(n+1/2)} + \mathbf{f}_{\text{ext}}^{(n+1/2)}, \end{aligned} \right\} \qquad (2.40) $$

where $\Delta t = t_{n+1} - t_n$ for a given time partition $\cup_n \{t_n, t_{n+1}\} = \{0, t_1, \dots\}$ of the time interval of interest, $\mathbf{d}_n \approx \mathbf{d}(t_n)$, $\mathbf{v}_n \approx \mathbf{v}(t_n)$, and $\mathbf{v}_{n+1/2} = (\mathbf{v}_{n+1} + \mathbf{v}_n)/2$. The momenta

$$ \mathbf{p}_t = \mathbf{M} \mathbf{v}_t \quad \text{for } t \in \cup_n \{t_n, t_{n+1}\}, \qquad (2.41) $$

have been eliminated in (2.40).

The discrete force of contact $f_c^{(n+1/2)}$ is defined in the following section. The vector $f_{\text{int}}^{(n+1/2)}$ in (2.40) corresponds to the time discretization proposed in [18]. It defines a second-order *conserving* approximation of the internal force vector at $t_{n+1/2}$, and is given by (2.34) with $\mathbf{B}_t$ evaluated at the mid-point configuration $\boldsymbol{\varphi}_{n+1/2} := (\boldsymbol{\varphi}_{n+1} + \boldsymbol{\varphi}_n)/2$ as

$$ f_{\text{int}}^{(n+1/2)} = \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \mathbf{B}_{n+1/2}^{\text{T}} \boldsymbol{\tau}^{(n+1/2)} \, \mathrm{d}\Omega, \qquad (2.42) $$

with the discrete Kirchhoff stresses $\boldsymbol{\tau}^{(n+1/2)}$ calculated as

$$ \boldsymbol{\tau}^{(n+1/2)} := \mathbf{F}_{n+1/2} \left( \frac{1}{2} \mathbf{C}(\mathbf{E}_n + \mathbf{E}_{n+1}) \right) \mathbf{F}_{n+1/2}^{\text{T}}, \qquad (2.43) $$

for the Saint-Venant-Kirchhoff model defined by Eq. (2.7). In (2.43), the deformation gradient $\mathbf{F}_{n+1/2} := \text{Grad } \boldsymbol{\varphi}_{n+1/2}$ is computed at the mid-point configuration, and the Green-Lagrange strain tensors $\mathbf{E}_n$ and $\mathbf{E}_{n+1}$ are evaluated at the configurations $n$ and $n+1$, respectively. The case involving a general stored energy function $W(E)$ can be found in [9].

As shown in [18], the following properties hold for the time discrete internal forces (2.42):

(i) Internal linear momentum contributions

$$ \left[ \sum_{A=1}^{n_{\text{node}}} f_{\text{int}}^{(n+1/2)A} \right] \cdot \mathbf{a} = \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \nabla_{n+1/2}(\mathbf{a}) : \boldsymbol{\tau}^{(n+1/2)} \, \mathrm{d}\Omega = 0 \quad \forall \mathbf{a} \in \mathbb{R}^{n_{\text{dim}}} \qquad (2.44) $$

where the vectors $f_{\text{int}}^{A(n+1/2)} \in \mathbb{R}^{n_{\text{dim}}}$ refer to the nodal forces corresponding to (2.42). We conclude that the summation in the left-hand side of (2.44) vanishes.

(ii) Internal angular momentum contributions

$$ \begin{aligned} & \left[ \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_{n+1/2}^{A} \times f_{\text{int}}^{(n+1/2)A} \right] \cdot \mathbf{w} = \sum_{\alpha=1}^{2} \sum_{A=1}^{n_{\text{node}}} \int_{\Omega^{\alpha}} [\mathbf{B}_{n+1/2}^{A}(\mathbf{x}_{n+1/2}^{A} \times \mathbf{w})]^{\text{T}} \boldsymbol{\tau}^{(n+1/2)} \, \mathrm{d}\Omega \\ &= \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \nabla_{n+1/2}(\mathbf{x}_{n+1/2}^{A} \times \mathbf{w}) : \boldsymbol{\tau}^{(n+1/2)} \, \mathrm{d}\Omega \\ &= \sum_{\alpha=1}^{2} \int_{\Omega^{\alpha}} \hat{\mathbf{W}} : \boldsymbol{\tau}^{(n+1/2)} \, \mathrm{d}\Omega = 0 \quad \forall \mathbf{w} \in \mathbb{R}^{n_{\text{dim}}}, \end{aligned} \qquad (2.45) $$

given the symmetry of $\boldsymbol{\tau}^{(n+1/2)}$ in (2.43). We conclude that the summation in the left-hand-side of (2.45) vanishes.

(iii) Internal energy contributions
$$
\begin{align}
\sum_{A=1}^{n_{\text{node}}} [f_{\text{int}}^{(n+1/2)A} \cdot (\mathbf{d}_{n+1}^A - \mathbf{d}_n^A)] &= \sum_{\alpha=1}^{2} \int_{\Omega^\alpha} \nabla_{n+1/2}(\mathbf{u}_{n+1} - \mathbf{u}_n) : \boldsymbol{\tau}^{(n+1/2)} \, d\Omega \nonumber \\
&= \sum_{\alpha=1}^{2} \int_{\Omega^\alpha} [E_{n+1} - E_n] : \frac{1}{2} C[E_{n+1} + E_n] \, d\Omega \nonumber \\
&= \mathcal{W}_{n+1}^h - \mathcal{W}_n^h, \tag{2.46}
\end{align}
$$

where the superscript (·)h refers to the discretized system of solids.

General (non-conserving) time discretizations of the internal force term, involving in particular high-
frequency dissipation, are considered in Section 4 in combination with the contact scheme developed next.

3. Conserving algorithms for frictionless dynamic contact

Our goal is the design of the time-discrete counterpart of the nodal contact forces (2.35), that enforces the
unilateral contact constraint and retains the conservation properties of the final algorithm. In Section 3.1 we
develop a penalty scheme that possesses these properties. An extension is presented in Section 3.2 that imposes
the velocity constraint (2.16).

3.1. An energy-restoring, momentum-conserving scheme

Consider for a typical time interval $[t_n, t_{n+1}]$ the second-order approximation of the gap evolution equation (2.16) given by

$$
g_{s,n+1}^d = g_{s,n}^d + \nu_{s,n+1/2} \cdot [ (\varphi_{n+1}^1(X^S) - \varphi_{n+1}^2(\tilde{Y}_{n+1/2}(X^S))) - (\varphi_n^1(X^S) - \varphi_n^2(\tilde{Y}_{n+1/2}(X^S))) ], \quad (3.1)
$$

involving the unit normal $\mathbf{v}_{s,n+1/2}$ defined by the closest-point projection $\tilde{\mathbf{Y}}_{n+1/2}(\mathbf{X}^S)$ of the slave node $S$ at the configuration $\boldsymbol{\varphi}_{n+1/2}$. The evaluation of the current positions of the contact particle $\tilde{\mathbf{Y}}_{n+1/2}(\mathbf{X}^S)$ at the times $t_n$ and $t_{n+1}$ is to be noted. We refer to the scalar quantity $g_{s,n}^d$ the *dynamic gap* (at $t_n$) in contrast with the real gap $g_{s,n}$ defined by the closest-point projection algorithm given by (2.13) at $t_n$.

Expression (3.1) can be written equivalently in the notation introduced in the previous section as

$$
g_{s,n+1}^d = g_{s,n}^d + \hat{\mathbf{G}}_{s,n+1/2}^T [\hat{\mathbf{d}}_{s,n+1} - \hat{\mathbf{d}}_{s,n}], \quad (3.2)
$$

for the corresponding displacements of slaves and master nodes at $t_n$ and $t_{n+1}$. The evolution of the dynamic gap (3.2) is initialized with the real gap $g_{s,n}$ for the last time step before contact. See details in Table 1 and the discussion below for the contact/release logic.

<table><caption>Table 1
Contact/release logic</caption>
  <tbody>
    <tr>
      <td>Let cont<sub>s,n</sub> = contact flag at t<sub>n</sub> (.true. or .false.), and g<sub>s,n+1</sub> the (real) gap at t<sub>n+1</sub> for slave node S. Then, cont<sub>s,n+1</sub> is defined as</td>
    </tr>
    <tr>
      <td>IF (cont<sub>s,n</sub> .or. (g<sub>s,n+1</sub>.le.0)) THEN</td>
    </tr>
    <tr>
      <td>Compute g<sup>d</sup><sub>s,n-1</sub> using (3.1).</td>
    </tr>
    <tr>
      <td>IF (g<sup>d</sup><sub>s,n+1</sub>.le.0) THEN</td>
    </tr>
    <tr>
      <td>cont<sub>s,n+1</sub> = .true.</td>
    </tr>
    <tr>
      <td>ELSE</td>
    </tr>
    <tr>
      <td>cont<sub>s,n-1</sub> = .false.</td>
    </tr>
    <tr>
      <td>IF (g<sub>s,n+1</sub>.ge.0) THEN</td>
    </tr>
    <tr>
      <td>The dynamic gap will be initialized to the current g<sub>s,n-1</sub> when evaluating (3.1) in the next time step.</td>
    </tr>
    <tr>
      <td>ELSE</td>
    </tr>
    <tr>
      <td>The dynamic gap will be initialized to the current g<sup>d</sup><sub>s,n+1</sub> when evaluating (3.1) in the next time step.</td>
    </tr>
    <tr>
      <td>ENDIF</td>
    </tr>
    <tr>
      <td>ENDIF</td>
    </tr>
    <tr>
      <td>ELSE</td>
    </tr>
    <tr>
      <td>cont<sub>s,n+1</sub> = .false.</td>
    </tr>
    <tr>
      <td>ENDIF</td>
    </tr>
  </tbody>
</table>
The difference $g_{s,n+1}^d \neq g_{s,n+1}^d$ (real gap) as employed in traditional treatments of the problem is to be noted. We point out that (3.1) corresponds to a second-order approximation of Eq. (2.16) for the evolution of the real gap $g_s$, and accounts for the (geometric) change of normal during contact. In one-dimensional problems, for instance, both gaps coincide. No loss of accuracy has been observed because of this approximation. We note in this regard that the definition of the gap function in terms of the closest-point projection (2.13) is, from a physical point of view, completely arbitrary.

The normal component of the contact force $p_s$ for the slave node S at the time step $[t_n, t_{n+1}]$ is defined by the penalty regularization of the contact constraint (2.12) given by the difference quotient

$$p_s = \begin{cases} -\frac{U(g_{s,n+1}^d) - U(g_{s,n}^d)}{g_{s,n+1}^d - g_{s,n}^d} & \text{if } g_{s,n+1}^d \neq g_{s,n}^d, \\ -U'\left(\frac{1}{2}(g_{s,n}^d + g_{s,n+1}^d)\right) & \text{if } g_{s,n+1}^d = g_{s,n}^d, \end{cases}$$

where $U(g)$ is a penalty regularization potential of the form, e.g.

$$U(g) := \begin{cases} \frac{1}{2} \kappa_p g^2 & \text{if } g \le 0, \\ 0 & \text{otherwise,} \end{cases} \qquad (3.4)$$

with a (large) penalty parameter $\kappa_p$. We observe that, given the approximation (3.3) of the derivative of the (decreasing) potential (3.4), $p_s \ge 0$ as required by (2.15)$_1$. The force of contact is then given by

$$\hat{f}_{s,c}^{(n+1/2)} = p_s \hat{G}_{s,n+1/2}, \qquad (3.5)$$

with $p_s$ as in (3.3). The evaluation of the normal contributions $\hat{G}_s$ in (3.2) and (3.5) at the mid-point $\varphi_{n+1/2} = (\varphi_n + \varphi_{n+1})/2$ configuration becomes crucial for the conservation of the total angular momenta as shown in the following section. A standard calculation shows that the final numerical scheme is formally second-order accurate in time.

The contact/release logic is summarized in Table 1 and proceeds as follows. The computation of the dynamic gap $g_{s,n+1}^d$ begins when a negative $g_{s,n+1}^d$ is encountered. As noted above, the dynamic gap $g_{s,n}^d$ is initialized with the value of the real gap at the last converged value before initial contact. Contact is detected if $g_{s,n+1}^d$ is negative, as implied by the check in (3.4) following (3.3) with $g^d$. We note that the normal contact force depends on the contact states at $t_{n+1}$ and $t_n$, and vanishes when both states at $t_n$ and $t_{n+1}$ are released states. We observe that $p_s \neq 0$ while releasing (i.e. $\text{cont}_{s,n} = \text{true}$ and $\text{cont}_{s,n+1} = \text{false}$, following the notation in Table 1). It has a positive value given by the contribution $U(g_{s,n}^d)$ at $t_n$. This final 'kick' restores the energy to the system of solids upon release.

Observe also that the same contribution to (3.3), $U(g_{s,n}^d)$, vanishes in the first contact increment. Therefore, the proposed penalty formulation enforces the *gap constraint* at the end of the time step $t_{n+1}$. This situation is to be contrasted with schemes enforcing only the velocity constraint (2.16) (the rate of the gap), thus requiring small time steps to avoid excessive penetrations of the solids, like in the conservative schemes presented in [14,24], as it has come to our attention recently.

### 3.1.1. Properties of the proposed scheme

The consideration of the interpolation functions in the definition of the linear momentum (2.19) leads to the expression

$$L_t^h := \sum_{A.B=1}^{n_{\text{node}}} M_{AB} \boldsymbol{v}_t^B, \qquad (3.6)$$

for its discrete counterpart at $t \in \bigcup_n \{t_n, t_{n+1}\}$, where $\boldsymbol{v}_t^B$ ($B=1, n_{\text{node}}$) denote the nodal velocities. We note that the same expression is reached by the consistent mass$_1$ (2.32) or lumped mass (2.33).
We define the total angular momentum for the discretized system at $t \in \cup_n \{t_n, t_{n+1}\}$ as

$$ J_t^h := \sum_{A,B=1}^{n_{\text{node}}} M_{AB} \mathbf{x}_t^A \times \mathbf{v}_t^B. \qquad (3.7) $$

For the consistent mass matrix, this expression follows from the inclusion of the isoparametric interpolations in (2.20). Similarly, we define the total energy of the discretized solids as

$$ \mathcal{E}_t^h := \mathcal{H}_t^h + \mathcal{W}_t^h, \quad \text{with} \quad \mathcal{H}_t^h := \frac{1}{2} \mathbf{v}_t^T \mathbf{M} \mathbf{v}_t, \quad \text{and} \quad \mathcal{W}_t^h := \sum_{\alpha=1}^{2} \int_{\Omega^\alpha} W^\alpha(\mathbf{F}^\alpha(\mathbf{d}_t)) \, d\Omega, \qquad (3.8) $$

for the mass matrix considered in the numerical simulation. The superscript $(\cdot)^h$ refers to (finite element) discrete quantities.

Noting that by (2.41)

$$ \mathbf{p}_t^A = \sum_{B=1}^{n_{\text{node}}} M_{AB} \mathbf{v}_t^B \quad A = 1, n_{\text{node}}, \qquad (3.9) $$

we can write the equivalent expressions

$$ L_t^h := \sum_{A=1}^{n_{\text{node}}} \mathbf{p}_t^A, \qquad (3.10) $$

for the discrete linear momentum, and

$$ J_t^h := \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_t^A \times \mathbf{p}_t^A, \qquad (3.11) $$

for the discrete angular momentum. The evolution of these quantities in the scheme defined by Eqs. (3.2)–(3.4) is characterized by the following proposition.

**PROPOSITION 3.1.** Let $\Gamma_u^\alpha = \emptyset$ ($\alpha = 1, 2$), and $f_{\text{ext}}^{(n+1/2)} = \mathbf{0}$ for a time increment $[t_n, t_{n+1}]$ i.e. homogeneous Neumann problem in that time interval). Then, the following evolution relations hold

(i) The linear momentum is conserved, i.e.

$$ L_{n+1}^h = L_n^h. \qquad (3.12) $$

(ii) The angular momentum is conserved, i.e.

$$ J_{n+1}^h = J_n^h. \qquad (3.13) $$

(iii) The energy evolves as

$$ \mathcal{E}_{n+1}^h + \mathcal{P}_{n+1}^h = \mathcal{E}_n^h + \mathcal{P}_n^h, \qquad (3.14) $$

where

$$ \mathcal{P}_t^h := \sum_{s=1}^{n_{\text{slave}}} U(g_{s,t}^d) \ge 0 \quad \text{for } t \in \cup_n \{t_n, t_{n+1}\}, \qquad (3.15) $$

with $\mathcal{P}_t^h = 0$ in a released state.

**PROOF.** The proof follows closely the proof of the Proposition 2.1, its continuum counterpart.

(i) Conservation of linear momentum. Adding the nodal components of Eq. (2.40)$_2$, we obtain

$$ \begin{aligned} L_{n+1}^h - L_n^h &= \sum_{A=1}^{n_{\text{node}}} (\mathbf{p}_{n+1}^A - \mathbf{p}_n^A) = \sum_{A,B=1}^{n_{\text{node}}} M_{AB} (\mathbf{v}_{n+1}^B - \mathbf{v}_n^B) \\ &= -\Delta t \underbrace{\sum_{A=1}^{n_{\text{node}}} f_{\text{int}}^{A,(n+1/2)}}_{= 0 \text{ by (2.44)}} + \Delta t \underbrace{\sum_{A=1}^{n_{\text{node}}} f_c^{A,(n+1/2)}}_{= 0 \text{ by (2.44)}} \end{aligned} $$
$$ = \Delta t \sum_{s=1}^{n_{\text{slave}}} p_s \underbrace{\left( - \sum_{l=1}^{n_{\text{master}}} N^{\text{MI}}(\xi_{s,n+1/2}) \right)}_{= 0 \text{ by (2.38)}} \mathbf{v}_{s,n \cdot 1/2} = 0, \qquad (3.16) $$

after using definition (3.5) of the contact force $f_c^{N+1/2}$.

(ii) Conservation of angular momentum. We first note the algebraic identity

$$ \mathbf{x}_{n+1}^{A} \times \mathbf{p}_{n+1}^{A} - \mathbf{x}_{n}^{A} \times \mathbf{p}_{n}^{A} = \mathbf{x}_{n+1/2}^{A} \times (\mathbf{p}_{n+1}^{A} - \mathbf{p}_{n}^{A}) + (\mathbf{x}_{n+1}^{A} - \mathbf{x}_{n}^{A}) \times \mathbf{p}_{n+1/2}^{A}, \qquad (3.17) $$

for $A = 1, n_{\text{node}}$. Eq. (2.40)$_1$ reads in nodal components

$$ \mathbf{x}_{n+1}^{A} - \mathbf{x}_{n}^{A} = \mathbf{d}_{n+1}^{A} - \mathbf{d}_{n}^{A} = \Delta t \mathbf{v}_{n+1/2}^{A} = \Delta t \sum_{B=1}^{n_{\text{node}}} (M^{-1})_{AB} \mathbf{p}_{n+1/2}^{B}, \qquad (3.18) $$

which leads to

$$ \sum_{A=1}^{n_{\text{node}}} (\mathbf{x}_{n+1}^{A} - \mathbf{x}_{n}^{A}) \times \mathbf{p}_{n+1/2}^{A} = \Delta t \sum_{A,B=1}^{n_{\text{node}}} (M^{-1})_{AB} \mathbf{p}_{n+1/2}^{B} \times \mathbf{p}_{n+1/2}^{A} = 0, \qquad (3.19) $$

by the symmetry of the mass coefficients $(M^{-1})_{AB} = (M^{-1})_{BA}$, and the skew-symmetry of the cross product. By Eq. (2.40)$_2$, we also have

$$ \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_{n+1/2}^{A} \times (\mathbf{p}_{n+1}^{A} - \mathbf{p}_{n}^{A}) = -\Delta t \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_{n+1/2}^{A} \times \mathbf{f}_{\text{int}}^{A,(n+1/2)} + \Delta t \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_{n+1/2}^{A} \times \mathbf{f}_{c}^{A,(n+1/2)}. \qquad (3.20) \\ = 0 \text{ by (2.45)} $$

Combining the definition (3.7) and Eq. (3.17) to (3.20), we obtain

$$ \begin{aligned} J_{n+1}^h - J_n^h &= \sum_{A=1}^{n_{\text{node}}} (\mathbf{x}_{n+1}^A \times \mathbf{p}_{n+1}^A - \mathbf{x}_n^A \times \mathbf{p}_n^A) = \Delta t \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_{n+1/2}^A \times \mathbf{f}_c^{A,(n+1/2)} \\ &= \Delta t \sum_{S=1}^{n_{\text{slave}}} p_s \left( \mathbf{x}_{n+1/2}^S - \sum_{A=1}^{n_{\text{master}}} N^{\text{MI}}(\xi_{s,n+1/2}) \mathbf{x}_{n+1/2}^{MI} \right) \times \mathbf{v}_{s,n+1/2} \\ &= \Delta t \sum_{S=1}^{n_{\text{slave}}} p_s g_{s,n+1/2} (\mathbf{v}_{s,n+1/2} \times \mathbf{v}_{s,n+1/2}) = 0, \end{aligned} \qquad (3.21) $$

after using (2.37).

(iii) Energy evolution. Combining the evolution equations (2.40) with the symmetry of the mass matrix $M$, we can write

$$ \begin{aligned} \mathcal{H}_{n+1}^h - \mathcal{H}_n^h &= \frac{1}{2} \mathbf{v}_{n+1}^T \mathbf{M} \mathbf{v}_{n+1} - \frac{1}{2} \mathbf{v}_n^T \mathbf{M} \mathbf{v}_n = \mathbf{v}_{n+1/2}^T \mathbf{M} (\mathbf{v}_{n+1} - \mathbf{v}_n) \\ &= -(\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_{\text{int}}^{(n+1/2)} + (\mathbf{d}_{n+1} - \mathbf{d}_n)^T \mathbf{f}_c^{(n+1/2)} \\ &= -(\mathcal{W}_{n+1}^h - \mathcal{W}_n^h) + \sum_{s=1}^{n_{\text{slave}}} p_s \hat{\mathbf{G}}_{n+1/2}^T (\hat{\mathbf{d}}_{s,n+1} - \hat{\mathbf{d}}_{s,n}), \end{aligned} \qquad (3.22) $$

where we have used the relation (2.46) for the internal forces. After noting that $\mathcal{E}^h = \mathcal{H}^h + \mathcal{W}^h$, and using the definitions of the dynamic gap (3.2) and the normal contact force (3.3) from the regularization potential $U(g^d)$, we conclude that

$$ \begin{aligned} \mathcal{E}_{n+1}^h - \mathcal{E}_n^h &= \sum_{s=1}^{n_{\text{slave}}} p_s (g_{s,n+1}^d - g_{s,n}^d) \\ &= - \sum_{s=1}^{n_{\text{slave}}} (U(g_{s,n+1}^d) - U(g_{s,n}^d)) = -(\mathcal{P}_{n+1}^h - \mathcal{P}_n^h), \end{aligned} \qquad (3.23) $$
which proves (3.14). We note that $\mathcal{P}_t^h = 0$ in a released state given the definition (3.4) of the regularization potential. $\square$

Proposition 3.1 shows that the time-stepping defined by (3.2) to (3.4) conserves the total linear and angular momentum of the system of solids in a homogeneous Neumann problem, as the original continuum system does. The relation (3.14) indicates that the total energy of the system solids plus the regularization potential is conserved during persistent contact. We note the important role of the definition of the dynamic gap for this property to hold. Furthermore, given the definition of the regularization potential (3.4), we have $U=0$ in a released state, so we conclude that the energy of the system of solids is conserved upon release. We summarize these observations in the following corollary.

**COROLLARY 3.2.** Let $\mathcal{E}_0^h$ denote the initial energy of the system of solids, corresponding to a released state (in the sense that $\mathcal{P}_0^h = 0$). Consider a homogeneous Neumann problem. Then, the energy at any time $0 \ge t_n \in \cup_n \{t_n, t_{n+1}\}$ is such that $\mathcal{E}_n^h = \mathcal{E}_0^h$ for a released state and $\mathcal{E}_n^h \le \mathcal{E}_0^h$ for a contact state.

**PROOF.** The result follows from (3.14) and the fact that $\mathcal{P}^h \ge 0$. $\square$

We note that Corollary 3.2 indicates that the energy of the system of solids will never increase beyond its initial value during the numerical simulation regardless of the size of the time-step $\Delta t$. We conclude the unconditional (energy) stability of the proposed scheme. The numerical simulations presented in Section 4 illustrate these stability properties.

### 3.1.2. A contact scheme with positive energy dissipation

As noted in the Introduction, when short-term simulations are employed for the study of high-velocity impacts, high-frequency energy dissipation may be a desired feature. We describe in this section a simple modification of the conserving contact scheme developed above that incorporates this property.

During persistent contact, the expression (3.3) for the normal component of the contact force reduces to

$$p_s = -\frac{1}{2} \kappa_p (g_{s,n+1}^d + g_{s,n}^d). \qquad (3.24)$$

A contact scheme with (positive) energy dissipation can be easily obtained by replacing (3.24) during persistent contact (i.e. $\text{cont}_{s,n} = \text{.true.}$ and $\text{cont}_{s,n+1} = \text{.true.}$, following the notation in Table 1) by

$$p_s = -\kappa_p (\vartheta g_{s,n+1}^d + (1-\vartheta) g_{s,n}^d), \qquad (3.25)$$

for $\vartheta > 1/2$. The difference scheme (3.3) is maintained during initial contact and release. Expression (3.24) is recovered with $\vartheta = 1/2$ in (3.25). The accuracy of the scheme drops to first order for $\vartheta \ne 1/2$.

With this modification, the balance of energy (3.23) reads

$$ \begin{aligned} \mathcal{E}_{n+1}^h - \mathcal{E}_n^h &= \sum_{s=1}^{n_{\text{slave}}} p_s (g_{s,n+1}^d - g_{s,n}^d) \\ &= - \sum_{s=1}^{n_{\text{slave}}} \kappa_p (\vartheta g_{s,n+1}^d + (1-\vartheta) g_{s,n}^d) (g_{s,n+1}^d - g_{s,n}^d) \\ &= - \sum_{s=1}^{n_{\text{slave}}} \kappa_p \left[ \frac{1}{2} (g_{s,n+1}^d + g_{s,n}^d) + \left(\vartheta - \frac{1}{2}\right) (g_{s,n+1}^d - g_{s,n}^d) \right] (g_{s,n+1}^d - g_{s,n}^d) \\ &= - \sum_{s=1}^{n_{\text{slave}}} \left( \mathcal{P}_{s,n+1}^h - \mathcal{P}_{s,n}^h \right) - \left(\vartheta - \frac{1}{2}\right) \sum_{s=1}^{n_{\text{slave}}} \kappa_p (g_{s,n+1}^d - g_{s,n}^d)^2, \end{aligned} \qquad (3.26) $$

for a time step in persistent contact. We conclude that

$$\mathcal{E}_{n+1}^h + \mathcal{P}_{n+1}^h \le \mathcal{E}_n^h + \mathcal{P}_n^h, \qquad (3.27)$$

if $\vartheta \ge 1/2$. The conservation of linear and angular momentum still holds, since the proof of these properties in Proposition 3.1 does not depend on the actual value of the normal contact force $p_s$.
It is important to emphasize that energy dissipation is not assured for schemes that are dissipative for linear problems (e.g. HHT type schemes). This fact is illustrated in the numerical simulations presented in Section 4. The normal contact force may create positive work on the initial and final release gaps (see [16]). In contrast, the proposed scheme has the proper dissipative properties as required. We note the important role played by the use of the dynamic gap (3.1) in this argument.

## 3.2. Enforcement of the velocity constraint

In situations where an extended time of contact appears, penalty schemes imposing only the gap constraint are known to lead in general to oscillations of the contact forces. These oscillations are also present in traditional schemes, and their origin can be traced in part to the lack of satisfaction of the constraint in the velocities (2.16). As discussed in Section 2.1, the velocity field is constrained by (2.16) during persistent contact. Finite element formulations where this constraint is enforced explicitly can be found in [22], and [15], among others. It is the goal of this section to present a modification of the penalty scheme described in Section 3.1 that accomplishes the imposition of (2.16) while maintaining the appropriate conservation properties.

To this end, we modify (2.41), and write the nodal linear momenta for a typical slave node/master segment pair as

$$ \hat{\mathbf{p}}_{s,t} = [\hat{\mathbf{M}}_{s,L} + m_{s,t}\hat{\mathbf{G}}_{s,t}\hat{\mathbf{G}}_{s,t}^{\mathrm{T}}]\hat{\mathbf{v}}_{s,t} = \hat{\mathbf{M}}_{s,L}\hat{\mathbf{v}}_{s,t} + m_{s,t}h_{s,t}\hat{\mathbf{G}}_{s,t}, \quad (3.28) $$

for $t \in \bigcup_n \{t_n, t_{n+1}\}$, where

$$ h_{s,t} := \hat{\mathbf{G}}_{s,t}^{\mathrm{T}} \hat{\mathbf{v}}_{s,t} = \mathbf{v}_t \cdot \left[ \mathbf{v}_t^S - \sum_{l=1}^{n_{\mathrm{master}}^{s}} N^{\mathrm{MI}}(\xi_{s,t}) \mathbf{v}_t^{\mathrm{MI}} \right], \quad (3.29) $$

the discrete counterpart of (2.16), the normal gap of the velocity. We consider the lumped mass matrix $\hat{\mathbf{M}}_{s,L}$ of the slave and master segment pair, i.e.

$$ \hat{\mathbf{M}}_{s,L} = \left[ \begin{array}{c} M_S \mathbf{1}_{n_{\mathrm{dim}}} \\ M_{M1} \mathbf{1}_{n_{\mathrm{dim}}} \\ \vdots \\ M_{M2} \mathbf{1}_{n_{\mathrm{dim}}} \\ \vdots \\ \vdots \\ M_S \mathbf{1}_{n_{\mathrm{dim}}} \end{array} \right] \in \mathbb{R}^{(1+n_{\mathrm{master}}^{s}) \times (1+n_{\mathrm{master}}^{s})}, \quad (3.30) $$

to simplify the final numerical implementation. In (3.28), $m_{s,t}$ denotes a mass added to the contacting slave and master nodes, which depends on the contact state as follows:

$$ m_{s,t} := \begin{cases} m_p & \text{if } g_{s,t}^d \le 0 \quad \text{or} \quad p_{s,t} > 0, \\ 0 & \text{otherwise,} \end{cases} \quad (3.31) $$

for a large penalty parameter $m_p > 0$. In (3.31), $p_{s,t}$ denotes the normal contact force for a slave node $S$ at time $t$ obtained via (3.3). We note that we consider the penalty mass active when this normal force component is positive, including the time increment when the contact is released. We have observed a better performance of the final numerical scheme with this combination (less oscillatory response of the final contact force, as described in Section 4). As $m_p \to \infty$, the constraint $h_{s,n+1} = 0$ for a typical time interval $[t_n, t_{n+1}]$ in contact is effectively imposed.

A mid-point approximation of Eqs. (2.29) is considered again. This leads, after the elimination of the momenta $\hat{\mathbf{p}}_{n+1}$, to the following contribution of a typical slave node/master segment pair

$$ \left. \begin{aligned} \frac{1}{\Delta t} (\hat{\mathbf{d}}_{s,n+1} - \hat{\mathbf{d}}_{s,n}) &= \hat{\mathbf{v}}_{s,n+1/2} + \frac{1}{2} \hat{\mathbf{M}}_{s,L}^{-1} (m_{s,n+1} h_{s,n+1} \hat{\mathbf{G}}_{s,n+1} + m_{s,n} h_{s,n} \hat{\mathbf{G}}_{s,n}), \\ \frac{1}{\Delta t} \hat{\mathbf{M}}_{s,L} (\hat{\mathbf{v}}_{s,n+1} - \hat{\mathbf{v}}_{s,n}) &= -\hat{\mathbf{f}}_{s,int}^{(n+1/2)} + \hat{\mathbf{f}}_{s,(c,mass)}^{(n+1/2)} + \hat{\mathbf{f}}_{s,ext}^{(n+1/2)}, \end{aligned} \right\} \quad (3.32) $$

where the modified contact force $\hat{\mathbf{f}}_{s,(c,mass)}^{(n+1/2)}$ is given by$_{15}$
$$ \hat{\mathbf{f}}_{s,(\text{c,mass})}^{(n+1/2)} = \hat{\mathbf{f}}_{s,c}^{(n+1/2)} - \hat{\mathbf{i}}_{s}^{(n+1/2)} \qquad (3.33) $$

with $\hat{\mathbf{f}}_{s,c}^{(n+1/2)}$ given by (3.5), and

$$ \hat{\mathbf{i}}_{s}^{(n+1/2)} := \frac{1}{\Delta t} (m_{s,n+1} h_{s,n+1}^{d} \hat{\mathbf{G}}_{s,n+1} - m_{s,n} h_{s,n}^{d} \hat{\mathbf{G}}_{s,n}). \qquad (3.34) $$

Physically, $\hat{\mathbf{i}}_{s}^{(n+1/2)}$ corresponds to the impulse enforcing the velocity constraint (2.16). We denote by $\hat{\mathbf{i}}_{s}^{(n+1/2),A}$ ($A = 1, n_{\text{node}}$) the corresponding nodal components, which vanish for the nodes not in contact.

### 3.2.1. Properties of the proposed scheme

The evolution of the linear momentum (3.6), the angular momentum (3.7), and the energy (3.8) in the scheme defined by Eqs. (3.32) and (3.33) is characterized by the following proposition.

**PROPOSITION 3.3.** Let $\Gamma_u^\alpha = \emptyset$ ($\alpha = 1, 2$), and $\mathbf{f}_{\text{ext}}^{(n+1/2)} = \mathbf{0}$ for a time increment $[t_n, t_{n+1}]$ (i.e. a homogeneous Neumann problem in that interval). Then,

(i) The linear momentum is conserved, i.e.

$$ \mathbf{L}_{n+1}^{h} = \mathbf{L}_{n}^{h}. \qquad (3.35) $$

(ii) The angular momentum is conserved, i.e.

$$ \mathbf{J}_{n+1}^{h} = \mathbf{J}_{n}^{h}. \qquad (3.36) $$

(iii) The energy evolves as

$$ \mathcal{E}_{n+1}^{h} + \mathcal{P}_{n+1}^{h} + \mathcal{M}_{n+1}^{h} = \mathcal{E}_{n}^{h} + \mathcal{P}_{n}^{h} + \mathcal{M}_{n}^{h}, \qquad (3.37) $$

where $\mathcal{P}_t^h \ge 0$ is defined in (3.15), and

$$ \mathcal{M}_t^h := \sum_{s=1}^{n_{\text{slave}}} m_{s,t} h_{s,t}^2 \left[ 1 + \frac{1}{2} m_{s,t} \hat{\mathbf{G}}_{s,t}^{\text{T}} \hat{\mathbf{M}}_{s,L}^{-1} \hat{\mathbf{G}}_{s,L} \right] \ge 0, \qquad (3.38) $$

for $t \in \cup_n \{t_n, t_{n+1}\}$.

**PROOF.** We first observe that the equivalent expressions (3.10) and (3.11) for the linear and angular momenta in terms of the nodal momenta $\mathbf{p}^A$ ($A = 1, n_{\text{node}}$) still hold for the modified momenta (3.28). Indeed, we have for the linear momentum at any time $t \in \cup_n \{t_n, t_{n+1}\}$

$$ \begin{aligned} \mathbf{L}_t^h &= \sum_{A,B=1}^{n_{\text{node}}} M_{AB} \mathbf{v}_t^B \\ &= \sum_{A=1}^{n_{\text{node}}} \mathbf{p}_t^A - \sum_{s=1}^{n_{\text{slave}}} m_{s,t} h_{s,t} \underbrace{\left(1 - \sum_{I=1}^{n_{\text{master}}} N^{MI}(\boldsymbol{\xi}_{s,t})\right)}_{= 0 \text{ by (2.38)}} \mathbf{v}_t \\ &= \sum_{A=1}^{n_{\text{node}}} \mathbf{p}_t^A, \end{aligned} \qquad (3.39) $$

and for the angular momentum

$$ \begin{aligned} \mathbf{J}_t^h &= \sum_{A,B=1}^{n_{\text{node}}} M_{AB} \mathbf{x}_t^A \times \mathbf{v}_t^B \\ &= \sum_{A=1}^{n_{\text{node}}} \mathbf{x}_t^A \times \mathbf{p}_t^A - \sum_{s=1}^{n_{\text{slave}}} m_{s,t} h_{s,t} \underbrace{\left( \mathbf{x}^S - \sum_{I=1}^{n_{\text{master}}} N^{MI}(\boldsymbol{\xi}_{s,t}) \mathbf{x}_t^{MI} \right)}_{= g_{s,t} \mathbf{v}_t \text{ by (2.37)}} \times \mathbf{v}_t \\ &= \sum_{A=1}^{n_{\text{node}}} \mathbf{x}^A \times \mathbf{p}_t^A - \sum_{s=1}^{n_{\text{slave}}} m_{s,t} h_{s,t} g_{s,t} \underbrace{(\mathbf{v}_t \times \mathbf{v}_t)}_{= 0} = \sum_{A=1}^{n_{\text{node}}} \mathbf{x}^A \times \mathbf{p}_t^A. \end{aligned} \qquad (3.40) $$
The conservation of linear and angular momentum by the scheme follows then easily by rewriting Eqs. (3.32) in terms of the modified momenta p given by (3.28)

$$ \left. \begin{aligned} \frac{1}{\Delta t} (\hat{\mathbf{d}}_{s,n+1} - \hat{\mathbf{d}}_{s,n}) &= \hat{\mathbf{M}}_{s,L}^{-1} \hat{\mathbf{p}}_{s,n+1/2}, \\ \frac{1}{\Delta t} (\hat{\mathbf{p}}_{s,n+1} - \hat{\mathbf{p}}_{s,n}) &= -\hat{\mathbf{f}}_{s,int}^{(n+1/2)} + \hat{\mathbf{f}}_{s,c}^{(n+1/2)}, \end{aligned} \right\} \quad (3.41) $$

for the homogeneous Neumann problem under consideration. After noting that Eqs. (3.32) are the same as the original Eqs. (2.40) in terms of the momenta p, the equivalences (3.39) and (3.40) imply the conservation properties

$$ L_{n+1}^h = L_n^h \quad \text{and} \quad J_{n+1}^h = J_n^h, \quad (3.42) $$

by the results (3.12) and (3.13) of Proposition 3.1 (whose proof has been developed in terms of the momenta p).

Similarly, using again the result (3.14) of Proposition 3.1, we can write for the scheme defined by (3.41) and (3.28) the following relation

$$ \tilde{\mathcal{H}}_{n+1}^h + \mathcal{W}_{n+1}^h + \mathcal{P}_{n+1}^h = \tilde{\mathcal{H}}_{n+1}^h + \mathcal{W}_{n+1}^h + \mathcal{P}_{n+1}^h, \quad (3.43) $$

where

$$ \tilde{\mathcal{H}}_t^h := \frac{1}{2} \mathbf{p}_t^\mathrm{T} \mathbf{M}_L^{-1} \mathbf{p}_t, \quad (3.44) $$

for $t \in \cup_n \{Q_{t_n}, t_{n+1}\}$. With the use of the definition (3.28), we can write

$$ \tilde{\mathcal{H}}_t^h = \mathcal{H}_t^h + \underbrace{\sum_{s=1}^{n_{\text{slave}}} m_{s,t} h_{s,n+1}^2 \left[ 1 + \frac{1}{2} m_{s,t} \hat{\mathbf{G}}_{s,t}^\mathrm{T} \hat{\mathbf{M}}_{s,L}^{-1} \hat{\mathbf{G}}_{s,t} \right]}_{:= \mathcal{M}_t^h}, \quad (3.45) $$

which combined with (3.43) results in (3.37). □

We observe that an extra contribution appears in this case in the energy balance corresponding to a kinetic energy contribution associated to the mass penalty introduced in the formulation. Given the energy balance (3.37) and the fact $m_s = 0$ after full release as defined by (3.31), we conclude that the total energy of the system is restored upon release. We can say that, during persistent contact, part of the energy is stored in the spring-like and the mass-like penalty regularization potentials. In fact, Corollary 3.2 still holds in this case resulting in the no increase of energy beyond its initial value during the numerical simulation and the desired nonlinear energy stability of the proposed method.

## REMARKS 3.1

(1) An augmented Lagrangian scheme for the velocity constraint can be introduced easily by adding to (3.28) a Lagrange multiplier field of the form

$$ \hat{\mathbf{p}}_{s,n+1} = \hat{\mathbf{M}}_{s,L} \hat{\mathbf{v}}_{s,n+1} + (m_{s,n+1} h_{s,n+1} + \lambda_{s,n+1}) \hat{\mathbf{G}}_{s,n+1}. \quad (3.46) $$

The Lagrange multiplier $\lambda_{s,n+1}$ is obtained by the update

$$ \lambda_{s,n+1}^{(k+1)} = \lambda_{s,n+1}^{(k)} + m_{s,n+1} h_{s,n+1}, \quad (3.47) $$

in the iteration (k) of an iteration procedure nested with the solution of the equations of motion, accomplishing the satisfaction of $h_{s,n+1} = 0$ with finite values of the mass penalty $m_{s,n+1}$. See [8,20] among others, for details on augmented Lagrangian methods.

(2) The mass penalty scheme described in this section can be combined with the energy dissipative scheme proposed in Section 3.1.2.
## 4. Representative numerical simulations

The goal of this section is to evaluate the performance of the newly proposed numerical schemes in several representative numerical simulations. To this end, we consider in Section 4.1 the impact of a linear elastic rod on a rigid wall, and the impact of two nonlinear elastic cylinders in Section 4.2.

### 4.1. Impact of a rod on a rigid wall

The purpose of this simulation is to show the important role that an energy restoring contact algorithm plays in the overall stability of the numerical scheme. As noted in Section 3.1, numerical schemes that are (unconditionally) dissipative for linear problems, and consequently (unconditionally) stable, do not possess this property in general nonlinear settings. As an example, we consider the well-known dissipative HHT schemes (or α-method as sometimes called), and show that the energy increases due to contact if the numerical scheme is not used with an adequate contact algorithm.

To this end, we consider a one-dimensional model of a rod impacting a rigid wall using a combination of different continuum and contact algorithms. The problem is sketched in Fig. 3. Linear elasticity is assumed for the one-dimensional continuum, so that the only nonlinearity arises from the contact conditions. We consider general discretizations in time of the continuum to accommodate dissipative schemes. In this setting, the three parameter family of HHT algorithms (see [11])

$$Ma_{n+1} + K[\alpha d_{n+1} + (1-\alpha)d_n] = f_{c,n+\alpha}, \quad (4.1)$$

$$d_{n+1} = d_n + \Delta t \mathbf{v}_n + \frac{1}{2} \Delta t^2 [2\beta a_{n+1} + (1-2\beta)a_n], \quad (4.2)$$

$$\mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t [\gamma a_{n+1} + (1 - \gamma)a_n], \quad (4.3)$$

is considered, where K denotes the usual stiffness matrix of linear elasticity. We note that Eq. (4.1) has been written in the form presented in [19], which differs from the original presentation of the α-method in [11] (the α parameter in (4.1) corresponds to $1 + \alpha$ of [11]).

We consider the following schemes:

(1) Trapezoidal rule: $\alpha = 1.0$, $\beta = 0.25$ and $\gamma = 0.5$.

(2) Midpoint rule: $\alpha = 0.5$, $\beta = 0.5$ and $\gamma = 1.0$.

(3) HHT: $\alpha = 0.51$, $\beta = 0.555025$ and $\gamma = 0.99$.

All three schemes are combined with a standard penalty scheme for the contact, with the contact constraint imposed at $t_{n-\alpha}$, consistent with (4.1). We consider also:

(4) The new energy restoring contact scheme, with midpoint rule for the continuum (as in Algorithm 2).

(5) The new energy dissipative contact scheme of Section 3.1.2 ($\vartheta = 1.0$), with HHT for the continuum (as in Algorithm 3).

We note that for the linear elastic continuum under consideration the conserving algorithm considered in Section 2.3.1 reduces to the midpoint rule and trapezoidal rule, which would coincide in this linear setting. As is well known, both schemes are conservative for linear problems. Similarly, the HHT Algorithm 3 is energy dissipative in the linear elastic case. However, the nonlinearity of the contact conditions when the simulation starts at a non-zero gap, destroys these conservative and dissipative properties, respectively. In essence, the work done by the contact force on the initial gap is not zero, and without control, leading to an increase of energy (see

<figure><img src="image_4.png" /><figcaption>Fig. 3. Impact of a rod on a rigid wall. Problem definition.</figcaption></figure>
[16]). This situation is to be contrasted with the newly proposed schemes. For Algorithms 4 and 5 the energy will not increase during the simulation, and for Algorithm 4 it will be restored completely upon release.

The rod considered in the simulations has unit length ($L = 1$) and unit cross section area ($A = 1$). The Young's modulus is $E = 1$, and density $\rho = 1$. The initial velocity of the rod before impact is $v_0 = 0.5$. The initial configuration of the rod is located at a distance of $d_0 = 7.5 \cdot 10^{-3}$ from the wall. The exact solution consists of a constant stress front propagating along the rod with the elastic wave speed $c = \sqrt{E}/\rho = 1$. The magnitude of the compressive stress is $\sigma = \rho v_0 c = 0.5$. This front reaches the right end at a time $L/c$ after impact where it is reflected. This reflection results in an unloading front that propagates back along the rod reaching the wall at a time $2L/c$ after impact. At this time the rod is released, that is, at

$$t = \frac{d_0}{v_0} + \frac{2L}{c} = 2.015, \qquad (4.4)$$

accounting for the initial time before impact ($d_0/v_0 = 15 \cdot 10^{-3}$). Therefore, the total force of contact is constant, and given by the value

$$f_c = \sigma A = \rho v_0 c A = 0.5, \qquad (4.5)$$

during the contact interval $15 \cdot 10^{-3} < t < 2.015$.

The rod is discretized with 100 linear finite elements in the numerical simulations presented herein. A Courant condition of CFL = 2 is considered, being therefore outside the range of stability of explicit methods like, e.g. central differences ($\alpha = 1.0, \beta = 0, \gamma = 0.5$). In all the cases, the contact penalty parameter is $\kappa_p = 10^6$, and the mass penalty parameter has the value $m_p = 10^3$ for Algorithms 4 and 5.

Figs. 4 to 5 show the results obtained with these schemes. The gap, velocity gap, contact force and total energy of the rod are plotted versus time. With respect to the standard schemes, Figs. 4 and 5 (left column) depict the results for the trapezoidal and midpoint rule, respectively, showing the severe oscillatory behavior associated with these schemes when trying to enforce the contact constraint. Oscillations between contact and released states lead to a clear unsatisfactory performance of the scheme. Furthermore, these oscillations lead to an increase of energy when they occur due to the associated nonlinearity. The trapezoidal rule, with the contact constraint imposed at $t_{n+1}$, improves the performance, as shown in Fig. 4, but the oscillatory response remains, as it does the non-physical increase of energy. Fig. 6 (left column) shows the results for the HHT. We still observe an initial oscillatory response, as well as an energy increase thus leading to potential instabilities of the scheme. Although the oscillations are eventually damped, this is obtained at the cost of a clear energy lost.

The performance of the standard schemes is to be contrasted with the newly proposed methods. Fig. 5 (right column) shows the results obtained with the energy restoring contact scheme. The good enforcement of both constraints ($g = 0$ and $\dot{g} = 0$) is to be noted. Even though small oscillations are observed, these are not between contact and released states. Persistent contact is maintained during the theoretical contact interval, as reflected in the persistent positive value of the contact force. The energy of the rod is under control during all the simulation, and it is restored upon final release. The total energy in the discrete system (rod and regularization spring) is conserved at all times. Fig. 6 (right column) shows the results for the energy dissipative scheme proposed in Section 3.1.2. As expected, we observe a damping of the oscillations in this problem involving the high-frequency part of the spectrum in the solution. The energy never increases beyond its initial value, avoiding any type of instabilities.

To gain a better understanding of the proposed methods, we have included in Fig. 7 the results obtained with the previously considered energy restoring scheme, without mass penalty ($m_p = 0$), i.e. no enforcement of the velocity constraint. Whereas the gap constraint is enforced equally for both schemes, we note the improvement accomplished in the imposition of the velocity constraint and the contact force. This improvement is to be traced to the impulse (3.34) introduced by the mass penalty in the definition of the contact force, and leads to better resolution of the small-time scales in problems where the contact intervals need to be resolved.

## 4.2. Impact of two cylinders

We now consider the impact of two nonlinear elastic cylinders in plane strain. The cylinders have a diameter of 2.0, and are discretized with displacement bilinear finite elements, as shown in Fig. 8. The Saint-Venant-Kirchhoff material model is assumed for both cylinders with Lamé constants, $\lambda = 2 \cdot 10^4$, $\mu = 1 \cdot 10^4$, and
TRAPEZOIDAL RULE

<figure><img src="image_5.png" /><figcaption>Fig. 4. Impact of a rod on a rigid wall. Results obtained with the trapezoidal rule, Algorithm 1.</figcaption></figure>

density $\rho=1$. These properties lead to the consideration of quasi-rigid cylinders. A penalty parameter of $\kappa_p = 1 \cdot 10^5$ is considered with $m_p = 0$, i.e. no imposition of the velocity constraint. We note that we are interested in the overall response of the system in this case (the long time scales), rather than resolving the different contact intervals in detail. A constant time step of $\Delta t = 0.1$ is considered.

Fig. 8 depicts the results obtained with the proposed new scheme in a simulation involving rigid walls as depicted. The left cylinder is given an initial velocity of {$v_x, v_y$} = {1, -2}, hitting the bottom rigid wall at $t \approx 1.5$. Fig. 9 shows the plots of the total energy of the cylinders (kinetic plus strain energies), the two components of the linear momentum ($L_x$ and $L_y$), and the angular momentum ($J$). The x-direction corresponds to the horizontal direction in the plots of Fig. 8, with y-direction being the perpendicular direction. We have included the results for the newly proposed contact energy-restoring scheme, and a standard midpoint rule contact (non-conserving), both in combination with the conserving scheme considered in Section 2.3.1 for the continuum. Therefore, the energy and momenta will be conserved for both schemes between contact interactions.

We observe that the initial hit of the left cylinder with the bottom wall leads to an increase in the y component of the linear momentum ($L_y$) and a change of the angular momentum, as expected. The increase of $L_y$ corresponds to the total force applied during contact, positive since it is pointing in the positive y-direction. The x-component of the linear momentum is conserved for both schemes, whereas the energy is only conserved (restored) after bouncing by the newly proposed scheme. In fact, we observe a sudden increase of the energy for the midpoint rule contact (to almost four times the original value), which is accompanied with a large change of
MIDPOINT RULE

<figure><img src="image_6.png" /><figcaption>Time</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Time</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Time</figcaption></figure>

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

ENERGY RESTORING

<figure><img src="image_10.png" /><figcaption>Time</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Time</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Time</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 5. Impact of a rod on a rigid wall. Results obtained with the midpoint rule, Algorithm 2 (left column), and the energy restoring, Algorithm 4 (right column).</figcaption></figure>
HHT

<figure><img src="image_14.png" /><figcaption>Gap</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Velocity Gap</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>Contact Force</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>Total Energy</figcaption></figure>

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

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

ENERGY DISSIPATIVE

<caption>Fig. 6. Impact of a rod on a rigid wall. Results obtained with the HHT scheme, Algorithm 3 (left column), and the energy restoring scheme, Algorithm 5 (right column).</caption>
Gap

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

Velocity Gap

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

Contact Force

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

Total Energy

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

Fig. 7. Impact of a rod on a rigid wall. Results obtained with the energy restoring scheme with no mass-penalty, $m_p = 0$. To be compared with the use of the mass penalty ($m_p = 10^3$) in Fig. 5, right column. As observed in this last case, the addition of the mass-penalty impulse enforces the velocity constraint and eliminates the oscillation of the contact force, as observed in this figure, when trying to resolve the contact time interval (short-time scales).

linear momentum in the y-direction (indicating an excessively large force of contact). The computed solutions will then differ afterwards. We note that due to the quasi-rigid character of the solids, the total energy is mostly kinetic energy.

After bouncing from the bottom wall, the left cylinder impacts the right cylinder which is at rest. This happens at $t \approx 2.2$ for the energy-restoring scheme, and earlier, at $t \approx 2.0$, for the midpoint rule contact, due to the excessive energy that the previous impact added to the left cylinder. As expected, no change of momenta (linear or angular) is associated to this impact for both algorithms. The energy, on the other hand, is increased again for the midpoint rule contact, whereas the energy-restoring scheme recovers again the initial energy after the small interval where the contact constraints are imposed. We note the good resolution of the gap constraint $g = 0$.

Next, the right cylinder impacts the right wall close to the upper right corner (at $t \approx 3.8$ for the energy-restoring scheme, and $t \approx 2.6$ for the midpoint rule). The $x$ component of the linear momentum $L_x$ is reduced due to the application of the contact force (pointing to the negative $x$-direction). A larger contact force is observed again for the midpoint rule, compared to the value obtained with the energy-restoring scheme. The total energy doubles in the former.

After bouncing, the right cylinder hits the upper wall (at $t \approx 4.3$ for the energy-restoring scheme, and $t \approx 3.3$ for the midpoint rule). This can be observed by the corresponding decrease (the contact force points in the negative $y$-direction) of the component $L_y$ of the linear momentum for both schemes. The left cylinder hits the
<figure><img src="image_26.png" /><figcaption>Fig. 8. Impact of two (quasi-rigid) cylinders. Deformed configurations at different times obtained with the newly proposed energy restoring scheme. The left cylinder impacts the right cylinder, which is at rest, after bouncing from the bottom rigid wall.</figcaption></figure>

left wall for the midpoint rule at $t \approx 3.7$, due to the excessive velocity that has gained in the previous impacts, leading to the increase of $L_x$ observed for this case. This does not happen for the energy-restoring scheme.

After these interactions, the two cylinders impact each other again in the middle of the domain. This occurs at $t \approx 6.0$ for the energy-restoring scheme and $t \approx 4.2$ for the midpoint rule. This impact cannot be resolved with the midpoint-rule contact scheme. The numerical computation blows up in this case (no convergence is obtained), with an unrealistic high value of the energy. The computation with the energy-restoring scheme can be continued without problems after the impact of both cylinders (no change of momenta, energy conserved
MIDPOINT RULE CONTACT

ENERGY RESTORING

<figure><img src="image_27.png" /><figcaption>Fig. 9. Impact of two (quasi-rigid) cylinders. Results obtained with the energy conserving scheme for the continuum in combination with a midpoint-rule contact (left column) and energy restoring contact scheme (right column).</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 10. Skew impact of two elastic cylinders. Deformed configurations at different times obtained with the newly proposed energy restoring scheme. Observe the large finite strains.</figcaption></figure>

again upon release). After this impact, the left cylinder hits the left wall at $t \approx 8.0$ (with the change of $L_x$ and angular momenta, no increase of energy again).

These results show the improved stability properties of the newly proposed scheme. These properties are achieved by the proper control of the energy during the computation. The correct conservation of energy and momenta has been verified. On the other hand, the artificial increase in energy for a standard contact scheme,
MIDPOINT RULE CONTACT

ENERGY RESTORING

<figure><img src="image_29.png" /><figcaption>Fig. 11. Skew impact of two elastic cylinders. Results obtained with the energy conserving scheme for the continuum in combination with a midpoint-rule contact (left column) and energy restoring contact scheme (right column).</figcaption></figure>
like the midpoint rule, has been shown to lead to numerical instabilities that force the termination of the computation. We point out that physically dissipative effects (like friction) would not stabilize the computations of standard schemes in general, as observed in [2] in the analysis of the stability of staggered algorithms for thermomechanical problems.

The above results considered quasi-rigid cylinders. In order to test the performance of the scheme with large finite elastic strains (and thus significant changes of the normal to the contact surfaces), we consider the same cylinders with Lamé constants $\lambda = 130$, $\mu = 43.33$, and density $\rho = 8.93$. The left cylinder is given an initial velocity $\{v_x, v_y\} = \{-1, 0.1\}$, while the right cylinder is at rest. Fig. 10 shows the impact of the two cylinders for this case. The large finite strains are apparent. Fig. 11 depicts the evolution of the energy, the two components of the linear momentum, and the angular momentum, for both the midpoint rule contact and the new energy-restoring scheme, both in combination with the conserving scheme developed in Section 2.3.1 for the continuum, as before. The non-physical increase of energy for the former is to be contrasted with the no increase and final conservation for the latter. The two schemes conserve all the momenta for this case. A penalty parameter of $k_p = 10^4$ is assumed, leading to a good satisfaction of the unilateral constraint (2.12), as the small energy associated to the regularization potential $U(g)$ indicates (the ripples in the plot of the energy) in Fig. 11. The same conclusions as for the previous simulations involving quasi-rigid cylinders apply to this case.

# 5. Concluding remarks

We have presented the formulation of a new class of implicit time-stepping algorithms for dynamic contact problems. The main characteristic of the proposed methods is the conservation laws that the discrete numerical schemes inherit from the continuum dynamical system by construction. In particular, it has been shown that the energy is under control at all times during the numerical simulation, leading to the proper (energy) stability properties, while efficiently enforcing the contact constraints.

These properties lead to improved performance in comparison with standard numerical techniques currently in use. The simplicity of the implementation of the proposed scheme, a modification of standard penalty formulations, is to be noted. Modifications involving the imposition of the velocity constraint and the introduction of positive high-frequency energy dissipation have been discussed in detail. Several numerical simulations have been presented that show the improved numerical stability properties of the new schemes over standard time-stepping algorithms.

# Acknowledgments

We are indebted to Professor R.L. Taylor for many helpful discussions, and the use of the finite element code FEAP. Financial support for this research was provided by the AFOSR under contracts no. F49620-97-1-0196 with UC Berkeley and no. 2-DJA-826 with Stanford University. This support is gratefully acknowledged. The first author acknowledges also the financial support of the ONR under contract N00014-96-1-0818 with UC Berkeley.

# Appendix A. Consistent linearization of the proposed schemes

In this appendix we develop the consistent linearization of the time stepping algorithms developed in this paper. To this purpose, we derive in Section A.1 the linearized equations of the problem. The contributions of the contact arrays to the tangent stiffness matrix are derived in Section A.2.

## A.1. *The linearized equations*

We consider the discrete equations (2.29) in terms of the nodal momenta $p$. As indicated below, the final implementation is carried out in terms of the nodal velocities $v$. Only the nodal momenta $p$ of the nodes in contact need to be considered for a non-vanishing mass penalty $m_p$. Define the residuals
$$
\left.
\begin{aligned}
\mathbf{R}_d &:= \frac{1}{\Delta t} \mathbf{M}(\mathbf{d}_{n+1} - \mathbf{d}_n) - \mathbf{p}_{n+1/2}, \\
\mathbf{R}_a &:= \mathbf{f}_{\text{ext}}^{(n+1/2)} + \mathbf{f}_c^{(n+1/2)} - \mathbf{f}_{\text{int}}^{(n+1/2)} - \frac{1}{\Delta t} (\mathbf{p}_{n+1} - \mathbf{p}_n).
\end{aligned}
\right\}
\quad (A.1)
$$

Given the nodal values {$d_n, v_n$}, and corresponding $p_n$ (see below), at time $t_n$, a consistent linearization of these equations leads to

$$
\left.
\begin{aligned}
\mathbf{R}_d^{(i)} + \frac{1}{\Delta t} \mathbf{M} \Delta \mathbf{d}_{n+1}^{(i+1)} - \frac{1}{2} \Delta \mathbf{p}_{n+1}^{(i+1)} &= \mathbf{0}, \\
\mathbf{R}_a^{(i)} + \Delta \mathbf{f}_c^{(n+1/2)} - \Delta \mathbf{f}_{\text{int}}^{(n+1/2)} - \frac{1}{\Delta t} \Delta \mathbf{p}_{n+1}^{(i+1)} &= \mathbf{0},
\end{aligned}
\right\}
\quad (A.2)
$$

with

$$
\mathbf{d}_{n+1}^{(i+1)} = \mathbf{d}_{n+1}^{(i)} + \Delta \mathbf{d}_{n+1}^{(i+1)}, \quad \text{and} \quad \mathbf{p}_{n+1}^{(i+1)} = \mathbf{p}_{n+1}^{(i)} + \Delta \mathbf{p}_{n+1}^{(i+1)}, \qquad (A.3)
$$

for the update between iterations (i) and (i + 1) in time step $[t_n, t_{n+1}]$ of a Newton-Raphson scheme for the solution of (A.1). The elimination of $\Delta \mathbf{p}_{n+1}^{(i+1)}$ leads to the final expression

$$
\left[ \frac{1}{2} (\mathbf{K}_{\text{int}}^{(i)} + \mathbf{K}_{c}^{(i)}) + \frac{2}{\Delta t^2} \mathbf{M} \right] \Delta \mathbf{d}_{n+1}^{(i+1)} = \mathbf{R}_{a}^{(i)} - \frac{2}{\Delta t} \mathbf{R}_{d}^{(i)}, \quad (A.4)
$$

where we have introduced the notation

$$
\Delta \mathbf{f}_{\text{int}}^{(n+1/2)} := \mathbf{K}_{\text{int}}^{(i)} \Delta \mathbf{d}_{n+1/2}^{(i+1)} \quad \left( \text{with } \Delta \mathbf{d}_{n+1/2}^{(i+1)} = \frac{1}{2} \Delta \mathbf{d}_{n+1}^{(i+1)} \right), \qquad (A.5)
$$

for the continuum contributions to the tangent stiffness, and

$$
\Delta \mathbf{f}_c^{(n+1/2)} := -\mathbf{K}_c^{(i)} \Delta \mathbf{d}_{n+1/2}^{(i+1)}, \qquad (A.6)
$$

(note the change of sign) for the contribution of the contact arrays. A closed-form expression for the contact stiffness $K_c^{(i)}$ is derived in Section A.2 below.

Once the updated nodal displacements $d_{n+1}^{(i+1)}$ and nodal momenta $p_{n+1}^{(i+1)}$ are known, the nodal velocities $v_{n+1}^{(i+1)}$ are recovered using the definition (3.28) of $p_{n+1}^{(i+1)}$, i.e.

$$
\mathbf{p}_{n+1}^{(i+1)} = \left( \mathbf{M} + \mathbf{A}_{s=1}^{n_{\text{slave}}} m_{s,n+1}^{(i+1)} \hat{\mathbf{G}}_{s,n+1}^{(i+1)} \otimes \hat{\mathbf{G}}_{n+1}^{(i+1)} \right) \mathbf{v}_{n+1}^{(i+1)}, \quad (A.7)
$$

in the general case involving the mass penalty $m_p \neq 0$. For the case of no mass penalty, $m_p = 0$, the dynamic update equation (A.1)$_1$ is linear, leading to $\mathbf{R}_d^{(i+1)} = \mathbf{0}$ and to the standard update

$$
\mathbf{v}_{n+1}^{(i+1)} = \mathbf{v}_{n+1}^{(i)} + \Delta \mathbf{v}_{n+1}^{(i+1)} = \mathbf{v}_{n+1}^{(i)} + \frac{2}{\Delta t} \Delta \mathbf{d}_{n+1}^{(i+1)}, \quad (A.8)
$$

without the need to consider the extra array $p_{n+1}^{(i+1)}$. As noted in Section 3.2, we consider a lumped mass matrix $\mathbf{M} = \mathbf{M}_L$ for the general case where we enforce the velocity gap constraint ($m_p \neq 0$), leading to the standard update (A.8) for the nodes not in contact, and the update (A.7) involving the nodes in contact only.

**REMARK 1.1.** An implementation avoiding the use of nodal momenta $\mathbf{p}$ for the nodes in contact can be easily devised by considering the linearized version of Eq. (A.7). Details are omitted.
## A.2. The contact stiffness

The linearization of the contact force $f_c$, defined by (2.35) as

$$f_c^{(n+1/2)} = \mathbf{A}_{s=1}^{n \text{ slave}} \hat{f}_{s,c}^{(n+1/2)} \quad \text{with} \quad \hat{f}_{s,c}^{(n+1/2)} = p_s \hat{G}_{n+1/2} \quad (A.9)$$

is given by

$$\Delta \hat{f}_{s,c}^{(n+1/2)} = \underbrace{\Delta p_s \hat{G}_{s,n+1/2}}_{\text{material}} + \underbrace{p_s \Delta \hat{G}_{s,n+1/2}}_{\text{geometric}}. \quad (A.10)$$

The material part of the tangent is defined by

$$\Delta p_s = \begin{cases} -\frac{U'(g_{s,n+1}^d) - p_s}{g_{s,n+1}^d - g_{s,n}^d} \Delta g_{s,n+1}^d & \text{if } g_{s,n+1}^d \neq g_{s,n}^d, \\ -U''(g_{s,n+1}^d) \Delta g_{s,n-1}^d & \text{if } g_{s,n+1}^d = g_{s,n}^d, \end{cases} \quad (A.11)$$

and

$$\Delta g_{s,n+1}^d = [2\hat{G}_{s,n+1/2} - c_1 \hat{D}_{s,n+1/2} - c_2 \hat{T}_{s,n+1/2}]^T \Delta \hat{d}_{s,n+1/2}, \quad (A.12)$$

$$c_1 = \frac{\Delta t}{l_s} \hat{T}_{s,n+1/2}^T \hat{v}_{s,n+1/2} + \frac{g_{s,n+1/2}}{l_s} \nu_{n+1/2}^T (\tau_{n+1} - \tau_n), \quad (A.13)$$

$$c_2 = \nu_{n+1/2}^T (\tau_{n+1} - \tau_n), \quad (A.14)$$

in two dimensional problems, with

$$\hat{D}_{n+1/2} = \begin{bmatrix} \mathbf{0} \\ -N_{,\xi}^{M1}(\xi_c) \boldsymbol{\nu}_{n+1/2} \\ -N_{,\xi}^{M2}(\xi_c) \boldsymbol{\nu}_{n+1/2} \\ \vdots \end{bmatrix}, \quad \text{and} \quad \hat{T}_{n+1/2} = \begin{bmatrix} \tau_{n+1/2} \\ -N_{,\xi}^{M1}(\xi_c) \tau_{n+1/2} \\ -N_{,\xi}^{M2}(\xi_c) \tau_{n+1/2} \\ \vdots \end{bmatrix}. \quad (A.15)$$

Here, we employed the notation

$$\Delta \hat{d}_{s,n+1/2} = \begin{bmatrix} \Delta d_{n-1/2}^s \\ \Delta d_{n+1/2}^{M1} \\ \Delta d_{n+1/2}^{M2} \\ \vdots \end{bmatrix}, \quad \text{and} \quad \hat{v}_{s,n+1/2} = \begin{bmatrix} \boldsymbol{v}_{n+1/2}^s \\ \boldsymbol{v}_{n+1/2}^{M1} \\ \boldsymbol{v}_{n+1/2}^{M2} \\ \vdots \end{bmatrix}. \quad (A.16)$$

The expression $g_{s,n+1/2}$ refers to the real gap found through the closest point projection at the configuration at $t_{n+1/2}$, $\tau_{n+1/2}$ is the normalized tangent vector to the master surface at the point of contact (i.e. $\boldsymbol{v}_{n+1/2}^T \boldsymbol{\tau}_{n+1/2} = 0$), and $l_s$ is the length of the surface element of the master surface corresponding to the given slave node $S$.

The geometric part of the tangent arises from the change of normal and contact point in $\hat{G}_{s,n+1/2}$, and is obtained as follows:

$$\Delta \hat{G}_{s,n+1/2} = - \frac{1}{l_s} \left[ \hat{T}_{s,n+1/2} \otimes \hat{D}_{s,n+1/2} + \hat{D}_{s,n+1/2} \otimes \hat{T}_{s,n+1/2} + \frac{g_{s,n+1/2}}{l_s} \hat{D}_{s,n+1/2} \otimes \hat{D}_{s,n+1/2} \right] \Delta \hat{d}_{s,n+1/2}, \quad (A.17)$$

for straight master segments as employed in the numerical simulations of Section 4. The final expression of the contact stiffness is then given by
$$ K_c = \mathbf{A}_{s=1}^{n_{\text{slave}}} \left\{ \left( \frac{U'(g_{s,n+1}^d) - p_s}{g_{s,n+1}^d - g_{s,n}^d} \right) \hat{G}_{s,n+1/2} \otimes [2\hat{G}_{s,n+1/2} - c_1 \hat{D}_{s,n+1/2} - c_2 \hat{T}_{s,n+1/2}] \right. \\ \left. + \frac{p_s}{l_s} \left[ \hat{T}_{s,n+1/2} \otimes \hat{D}_{s,n+1/2} + \hat{D}_{s,n+1/2} \otimes \hat{T}_{s,n+1/2} + \frac{g_{s,n+1/2}}{l_s} \hat{D}_{s,n+1/2} \otimes \hat{D}_{s,n+1/2} \right] \right\} $$

with the difference quotient in the first term replaced by $U''(g_{s,n+1}^d)$ if $g_{s,n+1}^d = g_{s,n}^d$. We note the non-symmetry of the material part as it occurs with its counterpart for the energy-momentum conserving algorithms considered in this paper for the continuum.

## References

[1] F. Armero and E. Petocz, A new class of conserving algorithms for dynamic contact problems, Numerical Methods in Engineering '96, Proc. Second ECCOMAS Conf. on Numerical Methods in Engineering, 9-13 September, 1996, Paris, France (John Wiley & Sons, 1996).

[2] F. Armero and J.C. Simo, A-priori stability estimates and unconditionally stable product formula algorithms for non-linear coupled thermoplasticity, Int. J. Plasticity 9 (1993) 149-182.

[3] K.J. Bathe and A. Chaudhary, A solution method for planar and axisymmetric contact problems, Int. J. Numer. Methods Engrg. 21 (1985) 65-88.

[4] T. Belytschko and M.O. Neal, Contact-impact by the pinball algorithm with penalty and Lagrangian methods, Int. J. Numer. Methods Engrg. 31 (1991) 547-572.

[5] N.J. Carpenter, R.L. Taylor and M.G. Katona, Lagrange constraints for transient finite element surface contact, Int. J. Numer. Methods Engrg. 32 (1991) 103-128.

[6] M. Crisfield and J. Shi, A co-rotational element/time-integration strategy for non-linear dynamics, Int. J. Numer. Methods Engrg. 37 (1994) 1897-1913.

[7] J.A. Desideri, C. Hirsch, P. Le Tallec, E. Onate, M. Pandolfi, J. Periaux and E. Stein, eds., Computational methods in applied sciences '96. Invited Lectures and Special Technological Sessions of the Second ECCOMAS Conference on Numerical Methods in Engineering, 9-13 September, 1996, Paris, France (John Wiley & Sons, 1996).

[8] R. Glowinski and P. Le Tallec, Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics (SIAM Studies in Applied Mathematics, PA, 1989).

[9] O. Gonzalez and J.C. Simo, Exact energy-momentum conserving algorithms for general models in nonlinear elasticity, Comput. Methods Appl. Mech. Engrg. (1995) to appear.

[10] J.O. Hallquist, G.L. Goudreau and D.J. Benson, Sliding interfaces with contact-impact in large-scale Lagrangian computations, Comput. Methods Appl. Mech. Engrg. 51 (1985) 107-137.

[11] H.M. Hilber, T.J.R. Hughes and R.L. Taylor, Improved numerical dissipation for the time integration algorithms in structural dynamics, Earthquake Engrg. Struct. Dyn. 5 (1977) 283-292.

[12] T.J.R. Hughes, R.L. Taylor, J.L. Sackmann, A. Curnier and W. Kanoknukulcahi, A finite element method for a class of contact-impact problems, Comput. Methods Appl. Mech. Engrg. 8 (1976) 249-276.

[13] N. Kikuchi and J.T. Oden, Contact Problems of Elasticity: A Study of Variational Inequalities and Finite Element Methods (SIAM, Philadelphia, PA, 1988).

[14] T. Laursen and V. Chawla, Design of energy conserving algorithms for frictionless dynamic contact problems, Int. J. Numer. Methods Engrg. (1996) in press.

[15] K. Lee, Numerical solution for dynamic contact problems satisfying the velocity and acceleration compatibilities on the contact surface, Comput. Mech. 15 (1994) 189-200.

[16] A. Munjiza, D.R.J. Owen and J. Bicanic, A combined finite-discrete element method in transient dynamics of fracturing solids, Engrg. Comput. 12 (1995) 145-174.

[17] J.D. Reid and K.H. Yang, Crashworthiness and occupant protection in transportation systems, AMD-Vol. 169 (The 1993 ASME Winter Annual Meeting, New Orleans, 1993).

[18] J.C. Simo and N. Tarnow, The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics, ZAMP 43 (1992) 757-793.

[19] J.C. Simo, N. Tarnow and M. Doblaré, Energy and momentum conserving algorithms for the dynamics of nonlinear rods, Int. J. Numer. Methods Engrg. 38 (1995) 1431-1474.

[20] J.C. Simo and T. Laursen, Augmented Lagrangian treatment of contact problems involving friction, Comput. Struct. 42 (1992) 97-116.

[21] J.C. Simo, J.E. Marsden and P.S. Krishnaprasad, The Hamiltonian structure of nonlinear elasticity. The convected representation of solids, rods and plates, Arch. Rat. Mech. Anal. 104 (1988) 125-183.

[22] R.L. Taylor and P. Papadopoulos, On a finite element method for dynamic contact impact problems, Int. J. Numer. Methods Engrg. 36 (1993) 2123-2140.
[23] C. Truesdell and W. Noll, The nonlinear field theories of mechanics, in: S. Fluege II, ed., Handbuch der Physik Bd. III/3 (Springer-Verlag, Berlin, 1965).

[24] T. Wasfy, Modeling contact impact of flexible manipulators with a fixed rigid surface, Proc. IEEE Int. Conf. on Robotics and Automation, Vol. 1 (IEEE, Piscataway, NJ, 1995) 621-626.