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

HAL
open science

# Design of energy conserving algorithms for frictionless dynamic contact problems

Tod Laursen, Vikas Chawla

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

## ► To cite this version:

Tod Laursen, Vikas Chawla. Design of energy conserving algorithms for frictionless dynamic contact problems. International Journal for Numerical Methods in Engineering, 1997, 40 (5), pp.863-886. 10.1002/(SICI)1097-0207(19970315)40:53.0.CO;2-V. hal-01435617

HAL Id: hal-01435617

https://hal.science/hal-01435617

Submitted on 14 Jan 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
DESIGN OF ENERGY CONSERVING ALGORITHMS FOR
FRICTIONLESS DYNAMIC CONTACT PROBLEMS

T. A. LAURSEN AND V. CHAWLA

Department of Civil and Environmental Engineering, Box 90287, Duke University, Durham, NC 27708-0287, U.S.A.

This paper proposes a formulation of dynamic contact problems which enables exact algorithmic conservation of linear momentum, angular momentum, and energy in finite element simulations. It is seen that a Lagrange multiplier enforcement of an appropriate contact rate constraint produces these conservation properties. A related method is presented in which a penalty regularization of the aforementioned rate constraint is utilized. This penalty method sacrifices the energy conservation property, but is dissipative under all conditions of changing contact so that the global algorithm remains stable. Notably, it is also shown that augmented Lagrangian iteration utilizing this penalty kernel reproduces the energy conserving (i.e. Lagrange multiplier) solution to any desired degree of accuracy. The result is a robust, stable method even in the context of large deformations, as is shown by some representative numerical examples. In particular, the ability of the formulation to produce accurate results where more traditional integration schemes fail is emphasized by the numerical simulations.

KEY WORDS: energy conservation; dynamics; contact; finite elements

# 1. INTRODUCTION

The contact problem attracts considerable attention from the computational mechanics community, due in large part to its highly non-linear and discontinuous nature. Indeed, engineering analysts charged with solving such problems will attest that merely achieving convergence of non-linear solution schemes can be difficult under many circumstances. These difficulties stem primarily from the fact that contact surfaces are unknown *a priori*, and must be evolved in a manner consistent with the equilibrium conditions of the two bodies and the Kuhn–Tucker conditions governing contact evolution. These physical requirements imply a variational inequality for the overall system, with admissible variations being constrained by the conditions governing the solution.¹,² Many finite element modelling efforts have been addressed most directly to the quasi-static contact problem, where inertial effects are essentially negligible and the contact conditions can be effectively devised and implemented on individual configurations, without much regard to the temporal variation of contact kinematic measures.

This work, by contrast, concerns itself with treatment of the dynamic contact problem. In general the most prevalent technique for treating dynamics in non-linear solid mechanics has been the semidiscrete finite element method. Formulations of this type can be roughly grouped into two classes: explicit schemes, which are favoured for highly transient problems; and implicit schemes,
which are more appropriate for systems dominated by low-frequency behaviour. Implementations
of the dynamic contact problem that have been proposed might also be grouped similarly, with
explicit treatments appearing for example in References 3–6, and implicit schemes being described,
e.g. in References 7 and 8. Roughly speaking, the conceptual approach to contact constraints in
many of these works can be described as follows. A semidiscrete finite element system is developed
exactly as would be done in a problem with no contact, except that a contact force vector is also
included. The contact tractions defining this force vector must be subject to some type of contact
conditions, and the description of these conditions is critical in determining the properties of the
resulting method. A common choice, for example, is to take the contact conditions to be the
same as those governing quasi-static response, and to apply any required numerical integration
techniques (e.g. for frictional response) exactly as would be done in the quasi-static case. An
analogue to this treatment is to be found in many implementations of elastoplasticity, where the
integration algorithm for the constitutive equations at each quadrature point does not in general
depend on the global time integration scheme, or even on whether the problem is static or dynamic.
Stretching this analogy a bit further, we might describe most contact implementations as treating
the contact force vector in the global equilibrium equations essentially as an extremely non-linear
(and non-smooth) internal force vector, with the ‘constitutive law’ for the contact being defined
by the complementarity conditions and the friction law, if present.

This general conceptual framework, while widely applied, appears to have several drawbacks.
In explicit calculations, when using the penalty method for constraint enforcement, one finds that
large penalties cannot be used in a fully explicit contact treatment without affecting the Courant
stability criterion.⁵ In softening the penalties, the accuracy of constraint is sacrificed to some degree, and because central difference methods possess no numerical damping, the noise generated
by the contact treatment will tend to obscure the solution as the calculation proceeds. In fact, one
can show⁶ that the Lagrange multiplier formulation of the fully explicit contact treatment is
singular, calling into question the basic validity of the fully explicit penalty approach. Accordingly,
Carpenter *et al*.⁶ have advocated an *implicit* treatment of the (quasi-static) contact constraints, along
with an otherwise explicit time integration of the momentum equations. This approach possesses
a well-defined Lagrange multiplier formulation which unfortunately couples the equilibrium equations in general. Carpenter and co-workers propose a Gauss-Seidel iteration scheme for constraint
enforcement, while Zhong⁹ uncouples the constraints by using an alternative discretization he refers
to as the ‘defense node’ approach. Both alternatives appear to be reasonably well behaved in comparison with the fully explicit approach, but the very fact that the contact conditions are treated
completely differently than the rest of the system calls into question our basic understanding of
the appropriate dynamic contact constraints.

In the implicit arena the state of the art is somewhat similar. For example, it has long been
recognized that use of the second-order accurate trapezoidal rule with a fully implicit treatment of
the contact constraints produces significant oscillations which can become worse as time steps and
spatial discretizations are refined (see, e.g. Reference 6). Some authors have proposed corrections
to the Newmark updates aimed at correcting such difficulties; examples are a correction based on
wave propagation results for linear elastic materials⁷ and a more recent treatment⁸ in which new
contacts are effectively forced to be persistent. While in some cases effective, the need for such
‘corrections’ suggests some sort of inadequacy in our basic understanding of contact constraints
in a dynamic context.

This paper makes an attempt at such understanding by formulating an implicit time integration
method which is fully *conservative*, i.e. it conserves all momenta and the total system energy for
hyperelastic bodies undergoing perfectly elastic impact events. Since we concentrate here on such
conservative systems, we consider only non-dissipative contact and therefore assume frictionless
response. Our approach is an extension of the work of Simo and Tarnow,¹⁰ where conservative algorithms were proposed and demonstrated for hyperelastic systems without contact constraints. In the current context, we will see that a Lagrange multiplier formulation of a particular *rate* constraint on interfaces is completely consistent with global conservation laws. Since the penalty method can be much simpler to apply than Lagrange multipliers in general, we propose a penalty regularization which is unconditionally dissipative, altering the energy conservation property but retaining the stability we desire. We will also demonstrate an augmented Lagrangian update scheme which uses this penalty kernel to reproduce the energy conserving (Lagrange multiplier) solution. In this manner, an algorithm is produced which is stable, and which can conserve energy to any desired degree of accuracy through the augmented Lagrangian iteration procedure. All development will be done without restricting the amount of motion or deformation that can occur, enabling application of the method in a very broad context.

The plan of the paper is as follows. Section 2 outlines the contact problem and discusses the conservation laws and their implications in a continuum context. Section 3 reviews the conservative time integration method proposed by Simo and Tarnow¹⁰ and discusses its extension to the frictionless contact problem. Imposition and regularization of the contact constraints in the context of the conservative algorithm is discussed in Section 4. In particular, Lagrange multiplier, penalty, and augmented Lagrangian algorithms will be proposed in this section. Section 5 briefly discusses the spatial discretization of the contact problem and some associated issues involved in finite element implementation. Finally, Section 6 presents some numerical examples demonstrating the performance of the method and some implications of the conservative scheme.

## 2. PROBLEM FORMULATION AND CONSERVATION LAWS

In the following, we briefly discuss the governing equations and contact conditions for the system of interest. The interested reader should consult Laursen and Simo¹¹ for more details on the continuum formulation of large-deformation contact problems.

### 2.1. Governing equations

We consider the open sets $\Omega^{(i)} \subset \mathbb{R}^{n_{sd}}$, $i = 1, 2$, which represent the reference configurations of two bodies expected to contact during a time interval of interest $[0, T]$. For each body, we define a portion of the boundary $\Gamma^{(i)} \subset \partial\Omega^{(i)}$ so that all expected areas of contact are included. Adopting a Lagrangian description of the problem, we designate material points in the contact surfaces as $\mathbf{X} \in \Gamma^{(1)}$ and $\mathbf{Y} \in \Gamma^{(2)}$. Writing the unknown configuration mappings at any time $t \in [0, T]$ as $\mathbf{\Phi}_t^{(i)}$, $i = 1, 2$, we can express the spatial positions of the contact surfaces as $\gamma_t^{(i)} = \mathbf{\Phi}_t^{(i)}(\Gamma^{(i)})$. Typical points $\mathbf{x} \in \gamma_t^{(1)}$ and $\mathbf{y} \in \gamma_t^{(2)}$ are then given by $\mathbf{x} = \mathbf{\Phi}_t^{(1)}(\mathbf{X})$ and $\mathbf{y} = \mathbf{\Phi}_t^{(2)}(\mathbf{Y})$. We assume that the $(n_{sd} - 1)$ dimensional manifolds $\Gamma^{(i)}$ are parametrized by mappings $\mathbf{\Psi}_0^{(i)}$ such that $\Gamma^{(i)} = \mathbf{\Psi}_0^{(i)}(\mathcal{A}^{(i)})$, $i = 1, 2$, where $\mathcal{A}^{(i)}$ is a parent domain for the surface in question and the mappings $\mathbf{\Psi}_0^{(i)}$ are assumed to be sufficiently smooth. In particular, considering $\Gamma^{(2)}$, we denote points in $\mathcal{A}^{(2)}$ as $\xi$ and write $\mathbf{Y} = \mathbf{\Psi}_0^{(2)}(\xi)$ and $\mathbf{y} = \mathbf{\Psi}_t^{(2)}(\xi)$, where $\mathbf{\Psi}_t^{(2)} = \mathbf{\Phi}_t^{(2)} \circ \mathbf{\Psi}_0^{(2)}$.

Considering any point $\mathbf{X} \in \Gamma^{(1)}$, the normal (impenetrability) contact conditions are written in terms of a gap function $g(\mathbf{X}, t)$.‡ It is defined at any time $t$ in terms of a *closest point projection*

<footnote>‡ We use the term *gap function* throughout the manuscript to be consistent with most contact mechanics literature. Actually, because of the sign convention chosen, the term penetration function might be more appropriate since the function *g* is positive when interpenetration occurs and negative otherwise</footnote>
in the spatial configuration:

$$g(\mathbf{X}, t) = \mathrm{sign}(g(\mathbf{X}, t))|g(\mathbf{X}, t)|$$

where

$$|g(\mathbf{X}, t)| = \min_{\mathbf{Y} \in \Gamma^{(2)}} \|\boldsymbol{\varphi}^{(1)}(\mathbf{X}, t) - \boldsymbol{\varphi}^{(2)}(\mathbf{Y}, t)\| \qquad (1)$$

and

$$\mathrm{sign}(g(\mathbf{X}, t)) = \begin{cases} -1 & \text{if } \boldsymbol{\varphi}^{(1)}(\mathbf{X}, t) \text{ is admissible} \\ 1 & \text{otherwise} \end{cases}$$

Impenetrability is enforced by the condition $g(\mathbf{X}, t) \le 0$. The point in $\Gamma^{(2)}$ achieving the minimization in (1) is written as $\mathbf{Y}$, with its counterpart in $\mathcal{A}^2$ denoted as $\xi$. It is important to remember that given a point $\mathbf{X}$, the identification of $\mathbf{Y}$ and $\xi$ will both in general vary with time, so that we will often write $\mathbf{Y}(\mathbf{X}, t)$ and $\xi(\mathbf{X}, t)$. Given these definitions, a basis can be constructed at each contact point by defining

$$\boldsymbol{\tau}_{\alpha} := \boldsymbol{\Psi}_{t, \alpha}^{(2)}(\xi) = \mathbf{F}_{t}^{(2)}(\mathbf{Y})\boldsymbol{\Psi}_{0, \alpha}^{(2)}(\xi), \quad \alpha = 1, \dots, n_{sd} - 1 \qquad (2)$$

It will be convenient in the following to augment this basis with a surface normal $\mathbf{v}$, which points out of body 2. In three dimensions, $\mathbf{v}$ would be defined via

$$\mathbf{v} = \frac{\boldsymbol{\tau}_1 \times \boldsymbol{\tau}_2}{\|\boldsymbol{\tau}_1 \times \boldsymbol{\tau}_2\|} \qquad (3)$$

where it is assumed that the parameterizations are defined so that $\mathbf{v}$ has the proper sense. Although both $\boldsymbol{\tau}_\alpha$ and $\mathbf{v}$ are to be associated with $\mathbf{X}$ and vary with time, we suppress these arguments in the following to reduce notation.

The (Piola) contact traction $\mathbf{t}^{(1)}(\mathbf{X}, t)$ is resolved into normal and tangential parts via

$$\mathbf{t}^{(1)}(\mathbf{X}, t) = t_N(\mathbf{X}, t)\mathbf{v} + \mathbb{P}_v \mathbf{t}^{(1)}(\mathbf{X}, t) \qquad (4)$$

where $\mathbf{v}$ is the outward normal to $\gamma_t^{(2)}$ at $\mathbf{y}$ (thus the inward normal to $\gamma_t^{(1)}$), and $\mathbb{P}_v \mathbf{t}^{(1)}$ is the projection of $\mathbf{t}^{(1)}$ onto the associated tangent plane. The contact pressure $t_N(\mathbf{X}, t)$ should be positive for compressive contact. We assume no frictional tractions in the present discussion; therefore $\mathbb{P}_v \mathbf{t}^{(1)} = 0$. The conditions for normal contact can now be written as

$$\begin{aligned} g(\mathbf{X}, t) &\le 0 \\ t_N(\mathbf{X}, t) &\ge 0 \\ t_N(\mathbf{X}, t)g(\mathbf{X}, t) &= 0 \\ t_N(\mathbf{X}, t)\dot{g}(\mathbf{X}, t) &= 0 \end{aligned} \qquad (5)$$

which must hold for all $\mathbf{X} \in \Gamma^{(1)}$ and for all $t \in \mathcal{I}$. Equations (5)$_{1-3}$ represent the classical Kuhn-Tucker complementarity conditions between gap and pressure. Equation (5)$_4$ is a constraint called the persistency condition, and requires that non-zero traction may only be generated during persistent contact. It will be of particular importance in subsequent developments.

With the contact conditions written, we specify the problem to be solved as follows.
Given the following boundary conditions on body force, traction, and boundary displacement:

$$
\begin{align}
\mathbf{f}_t^{(i)} & : \Omega^{(i)} \times \mathbb{I} \to \mathbb{R}^{n_{sd}}, \quad i = 1, 2 \notag \\
\mathbf{t}_t^{(i)} & : \Gamma_{\sigma}^{(i)} \times \mathbb{I} \to \mathbb{R}^{n_{sd}}, \quad i = 1, 2 \tag{6} \\
\boldsymbol{\phi}_t^{(i)} & : \Gamma_{\phi}^{(i)} \times \mathbb{I} \to \mathbb{R}^{n_{sd}}, \quad i = 1, 2 \notag
\end{align}
$$

Find the motions $\boldsymbol{\phi}_t^{(i)} : \Omega^{(i)} \times \mathbb{I} \to \mathbb{R}^{n_{sd}}$, $i = 1,2$, such that

$$
\begin{align}
\rho_0^{(i)} \frac{\partial^2}{\partial t^2} \boldsymbol{\phi}_t^{(i)} &= \mathrm{DIV} \mathbf{P}_t^{(i)} + \mathbf{f}_t^{(i)} \quad \text{on } \Omega^{(i)} \times \mathbb{I}, \quad i = 1,2 \notag \\
\boldsymbol{\phi}_t^{(i)} &= \boldsymbol{\phi}_t^{(i)} \quad \text{on } \Gamma_{\phi}^{(i)} \times \mathbb{I}, \quad i = 1,2 \notag \\
\mathbf{P}_t^{(i)} \mathbf{n}_0^{(i)} &= \mathbf{t}_t^{(i)} \quad \text{on } \Gamma_{\sigma}^{(i)} \times \mathbb{I}, \quad i = 1,2 \tag{7} \\
\boldsymbol{\phi}_t^{(i)} \Big|_{t=0} &= \mathbf{I} \quad \text{(the identity mapping)} \quad \text{in } \Omega^{(i)} \notag \\
\frac{\partial}{\partial t} \Big|_{t=0} \boldsymbol{\phi}_t^{(i)} &= \mathbf{V}_0^{(i)} \quad \text{in } \Omega^{(i)} \notag
\end{align}
$$

and contact conditions (5) are satisfied on $\Gamma^{(1)}$.

In the above, $\rho_0^{(i)}$ is the reference density, $\mathbf{f}_t^{(i)}$ is the body force, $\boldsymbol{\phi}_t^{(i)}$ is the prescribed boundary displacement, $\mathbf{t}_t^{(i)}$ is the prescribed boundary traction, and $\mathbf{V}_0^{(i)}$ is the initial material velocity. The constitutive relations governing $\mathbf{P}_t$ are at this point left unprescribed. The subregions $\Gamma_{\sigma}^{(i)}$, $\Gamma_{\phi}^{(i)}$, and $\Gamma^{(i)}$ are assumed to be non-intersecting and invariant with time, while satisfying

$$
\Gamma_{\sigma}^{(i)} \cup \Gamma_{\phi}^{(i)} \cup \Gamma^{(i)} = \partial \Omega^{(i)} \quad (8)
$$

2.2. Variational principle

At any time $t \in \mathbb{I}$, one can introduce admissible variations $\boldsymbol{\phi}^{*(i)}$ on each body and construct a variational principle. Following the development in Laursen and Simo,¹¹ which should be consulted for more details, this weak form of the equations can be stated as

Find $\boldsymbol{\phi}_t^{(i)} \in \mathcal{C}_t^{(i)}$, such that for all $\boldsymbol{\phi}^{*(i)} \in \mathcal{V}^{(i)}$, $i = 1,2$:

$$
\begin{align}
G^{(i)}(\boldsymbol{\phi}_t^{(i)}, \boldsymbol{\phi}^{*(i)}) &:= \int_{\Omega^{(i)}} \rho_0^{(i)} \mathbf{A}_t^{(i)} \cdot \boldsymbol{\phi}^{*(i)} \, d\Omega^{(i)} + \int_{\Omega^{(i)}} \mathbf{P}_t^{(i)} : \mathrm{GRAD}[\boldsymbol{\phi}^{*(i)}] \, d\Omega^{(i)} \nonumber \\
&\quad - \int_{\Omega^{(i)}} \mathbf{f}_t^{(i)} \cdot \boldsymbol{\phi}^{*(i)} \, d\Omega^{(i)} - \int_{\Gamma_{\sigma}^{(i)}} \mathbf{t}_t^{(i)} \cdot \boldsymbol{\phi}^{*(i)} \, d\Gamma_{\sigma}^{(i)} \tag{9} \\
&:= \int_{\Gamma^{(i)}} \mathbf{t}^{(i)} \cdot \boldsymbol{\phi}^{*(i)} \, d\Gamma^{(i)} \nonumber
\end{align}
$$

where $\mathbf{A}_t^{(i)}$ is the material acceleration of the body, and arguments $\mathbf{X}$ and $t$ are dropped from $\mathbf{t}^{(i)}(\mathbf{X}, t)$ for convenience of notation.

The time-dependent solution spaces $\mathcal{C}_t^{(i)}$ and time-independent variational space $\mathcal{V}^{(i)}$ are defined such that

$$
\mathcal{C}_t^{(i)} = \{ \boldsymbol{\phi}_t^{(i)} : \Omega^{(i)} \to \mathbb{R}^{n_{sd}} \mid J_t^{(i)} = \det[D\boldsymbol{\phi}_t^{(i)}] > 0 \text{ and } \boldsymbol{\phi}_t^{(i)} = \boldsymbol{\phi}_t^{(i)} \text{ on } \Gamma_{\phi}^{(i)} \} \quad (10)
$$
and

$$ \mathcal{V}^{(i)} = \{ \mathbf{\Phi}^{* (i)} : \Omega^{(i)} \to \mathbb{R}^{n_{sd}} | \mathbf{\Phi}^{* (i)} = 0 \text{ on } \Gamma_{\varphi}^{(i)} \} \quad (11) $$

One may add the virtual work expressions implied by (9) to define a global variational principle

$$ \begin{aligned} G(\varphi_t, \varphi^*) &:= G^{(1)}(\varphi_t^{(1)}, \varphi^{* (1)}) + G^{(2)}(\varphi_t^{(2)}, \varphi^{* (2)}) \\ &= \int_{\Gamma^{(1)}} \mathbf{t}^{(1)} \cdot \varphi^{* (1)} \, d\Gamma^{(1)} + \int_{\Gamma^{(2)}} \mathbf{t}^{(2)} \cdot \varphi^{* (2)} \, d\Gamma^{(2)} \end{aligned} \quad (12) $$

where $\varphi_t$ is understood to be the collection of mappings $\varphi_t^{(i)}$, $i = 1,2$ (similarly for $\varphi^*$). In (12) we also use the notations $G^{(i)}(\cdot, \cdot)$ to indicate the sum of the internal virtual work and those of the applied forces. The right-hand side of (12), representing the contact virtual work, can be expressed as a single integral over $\Gamma^{(1)}$ by requiring the contact forces on either side of the interface to be equal and opposite, which leads to

$$ G(\varphi_t, \varphi^*) + G_c(\varphi_t, \varphi^*) = 0 \quad (13) $$

where

$$ G_c(\varphi_t, \varphi^*) := - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\varphi^{* (1)}(\mathbf{X}) - \varphi^{* (2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)} \quad (14) $$

Equation (14) can be simplified further by considering linearized variations of the kinematic quantities, denoted here by the symbol $\delta(\cdot)$:

$$ \delta\beta(\mathbf{X}; \varphi_t, \varphi^*) := \left. \frac{\mathrm{d}}{\mathrm{d}\varepsilon} \right|_{\varepsilon=0} \beta(\varphi_t(\mathbf{X}) + \varepsilon\varphi^*(\mathbf{X})) \quad (15) $$

Direct calculation will verify that

$$ \delta g(\mathbf{X}) = -\mathbf{v} \cdot [\varphi^{* (1)}(\mathbf{X}) - \varphi^{* (2)}(\mathbf{Y}(\mathbf{X}, t))] \quad (16) $$

Using the fact that the frictional traction is zero, (16) can be substituted into (14) to produce

$$ G_c(\varphi_t, \varphi^*) = \int_{\Gamma^{(1)}} t_N(\mathbf{X}, t) \delta g(\mathbf{X}) \, d\Gamma^{(1)} \quad (17) $$

## 2.3. Conservation laws

Before we begin development of the algorithm, it is instructive to consider the global conservation laws in the context of the problem at hand. Specifically, we wish to verify that the linear momentum $\mathbf{L}_t$, the angular momentum $\mathbf{J}_t$, and the total energy $E_t^{\text{tot}}$ are globally conserved by the formulation we propose. The paper of Simo and Tarnow¹⁰ should be consulted for more detail on the general approach. Using the current notation, we define the total system linear momentum $\mathbf{L}_t$ and total system angular momentum $\mathbf{J}_t$ for any time $t$ as

$$ \mathbf{L}_t := \int_{\Omega^{(1)}} \rho_0^{(1)} \mathbf{V}_t^{(1)} \, d\Omega^{(1)} + \int_{\Omega^{(2)}} \rho_0^{(2)} \mathbf{V}_t^{(2)} \, d\Omega^{(2)} \quad (18) $$

and

$$ \mathbf{J}_t := \int_{\Omega^{(1)}} \rho_0^{(1)} \varphi_t^{(1)} \times \mathbf{V}_t^{(1)} \, d\Omega^{(1)} + \int_{\Omega^{(2)}} \rho_0^{(2)} \varphi_t^{(2)} \times \mathbf{V}_t^{(2)} \, d\Omega^{(2)} \quad (19) $$
The theorem of power expended states the global energy balance in rate form

$$ \frac{dK_t}{dt} + \mathcal{P}_t^{\text{int}} = \mathcal{P}_t^{\text{ext}} + \mathcal{P}_t^{\text{con}} \quad (20) $$

where the total kinetic energy is defined as

$$ K_t = \frac{1}{2} \left\{ \int_{\Omega^{(1)}} \rho_0^{(1)} \|\mathbf{V}_t^{(1)}\|^2 d\Omega^{(1)} + \int_{\Omega^{(2)}} \rho_0^{(2)} \|\mathbf{V}_t^{(2)}\|^2 d\Omega^{(2)} \right\} \quad (21) $$

the stress power is written as

$$ \mathcal{P}_t^{\text{int}} = \int_{\Omega^{(1)}} (D\boldsymbol{\phi}_t^{(1)}\mathbf{S}_t^{(1)}) : \text{GRAD}[\mathbf{V}_t^{(1)}] d\Omega^{(1)} + \int_{\Omega^{(2)}} (D\boldsymbol{\phi}_t^{(2)}\mathbf{S}_t^{(2)}) : \text{GRAD}[\mathbf{V}_t^{(2)}] d\Omega^{(2)} \quad (22) $$

the expended power of the external loading is given by

$$ \mathcal{P}_t^{\text{ext}} = \sum_{i=1}^{2} \left\{ \int_{\Omega^{(i)}} \mathbf{f}_t^{(i)} \cdot \mathbf{V}_t^{(i)} d\Omega^{(i)} + \int_{\Gamma_{\sigma}^{(i)}} \mathbf{t}_t^{(i)} \cdot \mathbf{V}_t^{(i)} d\Gamma^{(i)} \right\} \quad (23) $$

and the total power input of the contact stresses is given by

$$ \mathcal{P}_t^{\text{con}} = \sum_{i=1}^{2} \int_{\Gamma^{(i)}} \mathbf{t}^{(i)}(\mathbf{X}, t) \cdot \mathbf{V}_t^{(i)} d\Gamma^{(i)} \quad (24) $$

In (22), $\mathbf{S}_t^{(i)}$ is the second Piola–Kirchhoff stress tensor and is related to $\mathbf{P}_t^{(i)}$ by the relation $\mathbf{P}_t^{(i)} = D\boldsymbol{\phi}_t^{(i)}\mathbf{S}_t^{(i)}$.

Introducing the notations $\mathcal{D}_t^{\text{int}}$ for the internal dissipation function and $E_t^{\text{int}}$ for the total stored internal energy, we can write the following form of the second law of thermodynamics (i.e. a reduced dissipation inequality):

$$ \mathcal{D}_t^{\text{int}} := \mathcal{P}_t^{\text{int}} - \frac{dE_t^{\text{int}}}{dt} \ge 0 \quad \text{for all } t \in \mathbb{R} \quad (25) $$

Combining (20) and (25) yields

$$ \frac{dE_t^{\text{tot}}}{dt} = \mathcal{P}_t^{\text{ext}} + \mathcal{P}_t^{\text{con}} - \mathcal{D}_t^{\text{int}} \quad (26) $$

where $E_t^{\text{tot}}$ is the sum of the kinetic and internal energies of the two bodies in question. In the study considered in this paper, we will make the following assumptions about the system at hand:

1. The bodies are subject to no body forces, so that $\mathbf{f}_t^{(i)} = \mathbf{0}$ on $\Omega^{(i)} \times \mathbb{R}$ for $i = 1, 2$.

2. There are no Dirichlet (displacement) boundary conditions, and the tractions are zero on the Neumann boundaries, so that $\Gamma_{\phi}^{(i)} = \emptyset$ and $\mathbf{t}_t^{(i)} = \mathbf{0}$ on $\Gamma_{\sigma}^{(i)} \times \mathbb{R}$, $i = 1, 2$.

3. There is no internal dissipation in the bodies under consideration, so that $\mathcal{D}_t^{\text{int}} = 0$ for all time $t$.

Under this set of assumptions, we can examine the conservation properties of our system by making appropriate substitutions into the virtual work expression (13). For example, substituting
$\Phi^* = \eta$ into (13), where $\eta$ is a constant arbitrary vector on $\Omega^{(1)} \cup \Omega^{(2)}$, gives

$$
\begin{align}
0 &= G(\Phi_t, \eta) + G_c(\Phi_t, \eta) \nonumber \\
  &= \sum_{i=1}^{2} \int_{\Omega^{(i)}} \rho_0^{(i)} \mathbf{A}_t^{(i)} \cdot \eta \, d\Omega^{(i)} - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot (\eta - \eta) \, d\Gamma^{(1)} \tag{27} \\
  &= \left( \frac{d\mathbf{L}_t}{dt} \right) \cdot \eta \nonumber
\end{align}
$$

which in turn implies that $d\mathbf{L}_t/dt=\mathbf{0}$ ($\eta$ is arbitrary). In a similar manner, considering the variation
$\Phi^* = \eta \times \Phi_t$, where $\eta$ is a fixed vector, gives

$$
\begin{align}
0 &= G(\boldsymbol{\phi}_t, \boldsymbol{\eta} \times \boldsymbol{\phi}_t) + G_c(\boldsymbol{\phi}_t, \boldsymbol{\eta} \times \boldsymbol{\phi}_t) \nonumber \\
  &= \sum_{i=1}^{2} \left\{ \int_{\Omega^{(i)}} (D\boldsymbol{\phi}_t^{(i)}\mathbf{S}_t^{(i)}): (\hat{\boldsymbol{\eta}}D\boldsymbol{\phi}_t^{(i)}) d\Omega^{(i)} + \int_{\Omega^{(i)}} (\boldsymbol{\phi}_t^{(i)} \times \rho_0^{(i)}\mathbf{A}_t^{(i)}) \cdot \boldsymbol{\eta} d\Omega^{(i)} \right\} \nonumber \\
  &\quad - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\boldsymbol{\eta} \times (\boldsymbol{\phi}_t^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t)))] d\Gamma^{(1)} \nonumber \\
  &= \sum_{i=1}^{2} \int_{\Omega^{(i)}} (D\boldsymbol{\phi}_t^{(i)}\mathbf{S}_t^{(i)}D\boldsymbol{\phi}_t^{(i)\mathrm{T}}): \hat{\boldsymbol{\eta}} d\Omega^{(i)} + \left(\frac{d\mathbf{J}_t}{dt}\right) \cdot \boldsymbol{\eta} \nonumber \\
  &\quad + \int_{\Gamma^{(1)}} \boldsymbol{\eta} \cdot [\mathbf{t}^{(1)}(\mathbf{X}, t) \times (\boldsymbol{\phi}_t^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t)))] d\Gamma^{(1)} \tag{28} \\
  &= \frac{d\mathbf{J}_t}{dt} \cdot \boldsymbol{\eta} \nonumber
\end{align}
$$

where $\hat{\eta}$ is the skew-symmetric tensor whose axial vector is $\eta$. In the last step of (28), the first
term disappears because the contraction of a skew symmetric tensor ($\hat{\eta}$) with a symmetric tensor
$(D\boldsymbol{\phi}_t^{(i)}\mathbf{S}_t^{(i)}D\boldsymbol{\phi}_t^{(i)\mathrm{T}})$ is zero, while the last term, the contact contribution, disappears because in the
frictionless case the traction $\mathbf{t}^{(1)}(\mathbf{X}, t)$ is collinear with $(\boldsymbol{\phi}_t^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t)))$. Again, since $\eta$ is
arbitrary, we can conclude that $d\mathbf{J}_t/dt = \mathbf{0}$, so that angular momentum is conserved.

To examine the energy, we consider use of the material velocity field $\mathbf{V}_t$ as the variation.
Accordingly, we compute

$$
\begin{align}
0 &= G(\boldsymbol{\phi}_t, \mathbf{V}_t) + G_c(\boldsymbol{\phi}_t, \mathbf{V}_t) \nonumber \\
  &= \sum_{i=1}^{2} \left\{ \int_{\Omega^{(i)}} (D\boldsymbol{\phi}_t^{(i)}\mathbf{S}_t^{(i)}): \mathrm{GRAD}[\mathbf{V}_t^{(i)}] d\Omega^{(i)} + \int_{\Omega^{(i)}} \rho_0^{(i)}\mathbf{A}_t^{(i)} \cdot \mathbf{V}_t^{(i)} d\Omega^{(i)} \right\} \nonumber \\
  &\quad - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)} \tag{29} \\
  &= \mathcal{P}_t^{\mathrm{int}} + \frac{dK_t}{dt} - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)} \nonumber
\end{align}
$$

In view of our assumptions above, $\mathcal{P}^{\text{ext}} = 0$. Using this fact and comparing (29) with (20) leads to the identification $\mathcal{P}_t^{\text{con}} = \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)}$. By similar reasoning used to calculate $\delta g$ in (16), the material time derivative $\dot{g}(\mathbf{X}, t)$ is given by

$$
\dot{g}(\mathbf{X}, t) = -\mathbf{v} \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] \quad (30)
$$
so that the expression for $\mathcal{P}_t^{\text{con}}$ becomes

$$ \begin{aligned} \mathcal{P}_t^{\text{con}} &= \int_{\Gamma^{(1)}} \mathbf{t}^{(1)}(\mathbf{X}, t) \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)} \\ &= \int_{\Gamma^{(1)}} t_N(\mathbf{X}, t) \mathbf{v} \cdot [\mathbf{V}_t^{(1)}(\mathbf{X}) - \mathbf{V}_t^{(2)}(\mathbf{Y}(\mathbf{X}, t))] d\Gamma^{(1)} \quad (31) \\ &= - \int_{\Gamma^{(1)}} t_N(\mathbf{X}, t) \dot{g}(\mathbf{X}, t) d\Gamma^{(1)} \end{aligned} $$

If we now examine equation (5)₄, the persistency condition for frictionless contact, we are led to conclude

$$ \mathcal{P}_t^{\text{con}} = 0 \quad (32) $$

which, in view of (26) and the fact that $\mathcal{P}_t^{\text{ext}} = \mathcal{P}_t^{\text{int}} = 0$, implies

$$ \frac{dE_t^{\text{tot}}}{dt} = 0 \quad (33) $$

In other words, total energy is conserved for the system at hand as a direct result of persistency condition (5)₄. Another way to say this is that if we wish to ensure that all contacts are perfectly elastic, so that the net contact power input to the system is zero, the persistency condition should be satisfied. This observation will be key to the ensuing algorithmic development, where a counterpart of this condition in the temporally discrete framework must be found.

### 3. CONSERVATIVE DISCRETIZATION SCHEMES

In this section we extend the energy-momentum conserving scheme of Simo and Tarnow¹⁰ to the contact problem posed in the last section. The reader is referred to that reference for extensive details on the approach, which will be only briefly outlined here. In the following, we pay particular attention to the temporal discretization of the contact conditions, with the aim of algorithmically reproducing the conservation properties outlined for the continuous case in the last section.

We subdivide the time interval of interest $\ll$ into intervals $[t_n, t_{n+1}]$, where $n$ is an index on time steps. Given a time step $n$, we sometimes write $\Delta t = t_{n+1} - t_n$, and note that in general $\Delta t$ need not be uniform throughout a problem. Focusing our attention on a typical time interval $[t_n, t_{n+1}]$ we will use the notation $(\cdot)_n$ to mean the algorithmic (i.e. time discrete) approximation to the continuum variable $(\cdot)(t_n)$. The conservation properties we wish to maintain are

$$ \mathbf{L}_{n+1} = \mathbf{L}_n, \quad \mathbf{J}_{n+1} = \mathbf{L}_n, \quad E_{n+1}^{\text{tot}} = E_n^{\text{tot}} \quad (34) $$

where all quantities are as defined in the last section.

Simo and Tarnow described their algorithm by making the following definitions based on convex combinations of variables at $n+1$ and variables at $n$:

$$ \mathbf{\Phi}_{n+\alpha} = \alpha \mathbf{\Phi}_{n+1} + (1-\alpha) \mathbf{\Phi}_{n}, \quad \mathbf{V}_{n+\alpha} = \alpha \mathbf{V}_{n+1} + (1-\alpha) \mathbf{V}_{n} \quad (35) $$

where $\alpha \in [0, 1]$ is an algorithmic parameter. Considering first a problem with only one body, and hence no contact constraints, one could summarize their algorithm in the unforced case as
Given all data at n, find $\varphi_{n+1} \in \mathscr{C}_{n+1}$, such that for all $\varphi^* \in \mathcal{V}$

$$0 = G(\varphi_{n+\alpha}, \varphi^*) := \int_{\Omega} \rho_0 \frac{(V_{n+1} - V_n)}{\Delta t} \cdot \varphi^* d\Omega + \int_{\Omega} (D\varphi_{n+\alpha} S) : \text{GRAD}[\varphi^*] d\Omega \quad (36)$$

$$\mathbf{V}_{n+\alpha} = \frac{1}{\Delta t} [\varphi_{n+1} - \varphi_n]$$

This algorithm has the following properties:

(i) Algorithmic conservation of linear momentum (i.e. satisfaction of $(34)_1$) for any $\alpha \in [0, 1]$ and for an arbitrary constitutive relation describing the (symmetric) second Piola–Kirchhoff tensor $S$;

(ii) Algorithmic conservation of angular momentum (i.e. satisfaction of $(34)_2$) for the case where $\alpha = \frac{1}{2}$ and arbitrary prescription of $S$;

(iii) Algorithmic energy conservation (i.e. satisfaction of $(34)_3$) for the case where $\alpha = \frac{1}{2}$, and where $S$ is defined according to a gradient of a generic stored energy function $\hat{e}(C)$:

$$S := \nabla \hat{e}(C_{n+\beta}) + \nabla \hat{e}(C_{n+1-\beta}) \quad (37)$$

$C_{n+\beta}$ is an algorithmic right Cauchy–Green tensor, defined via

$$C_{n+\beta} := \beta C_{n+1} + (1-\beta)C_n \quad (38)$$

and $\beta \in [0, 1]$ is selected to satisfy

$$\hat{e}(C_{n+1}) - \hat{e}(C_n) = \nabla \hat{e}(C_{n+\beta}) \cdot (C_{n+1} - C_n) \quad (39)$$

Such a $\beta$ must always exist as a direct consequence of the mean value theorem.

(iv) Second-order accuracy under the conditions stated in item (iii).

For more details on the proof of these properties the paper by Simo and Tarnow is recommended. In particular, properties (i)–(iii) can be verified by substitution of $\eta$, $\eta \times \varphi_{n+\alpha}$ and $V_{n+1/2}$, respectively, into $(36)_1$, where $\eta$ is again an arbitrary constant vector. These arguments are directly analogous to those given in the last section for the time continuous case.

Returning once more to the case at hand, where contact constraints are active between two bodies, we propose the following algorithm, which extends that in (36) by including $G_c(\varphi_{n+\alpha}, \varphi^*)$:

Given all data at n, find $\varphi_{n+1} \in \mathscr{C}_{n+1}$, such that for all $\varphi^* \in \mathcal{V}$:

$$0 = G(\varphi_{n+\alpha}, \varphi^*) + G_c(\varphi_{n+\alpha}, \varphi^*) \\ \mathbf{V}_{n+\alpha} = \frac{1}{\Delta t} [\varphi_{n+1} - \varphi_n] \quad (40)$$

where

$$G(\varphi_{n+\alpha}, \varphi^*) := \sum_{i=1}^{2} \left\{ \int_{\Omega^{(i)}} \rho_0 \frac{(V_{n+1}^{(i)} - V_n^{(i)})}{\Delta t} \cdot \varphi^* (\varphi^{(i)}) d\Omega^{(i)} + \int_{\Omega^{(i)}} (D\varphi_{n+\alpha}^{(i)} S^{(i)}) : \text{GRAD}[\varphi^* (\varphi^{(i)})] d\Omega^{(i)} \right\} \quad (41)$$

and

$$G_c(\varphi_{n+\alpha}, \varphi^*) := - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)} \cdot [\varphi^* (\mathbf{X}) - \varphi^* (\mathbf{Y}_{n+\alpha}(\mathbf{X}))] d\Gamma^{(1)} \quad (42)$$
$\mathbf{Y}_{n+\alpha}(\mathbf{X})$ is that point of $\Gamma^{(2)}$ which minimizes $\|\boldsymbol{\phi}_{n+\alpha}^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_{n+\alpha}^{(2)}(\mathbf{Y})\|$; i.e. the closest point projection used to define the contact basis is done in the $n+\alpha$ configuration. One can also infer from this definition that force equilibrium on the interface is enforced in the $n+\alpha$ configuration.

It is readily verified that properties (i) and (ii) of the non-contact algorithm are retained by this algorithm. These verifications are directly analogous to those for the continuum case and will not be repeated here. In examination of the energy conservation property, we set $\alpha = \frac{1}{2}$ and examine the equation

$$G(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2}) + G_c(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2}) = 0 \quad (43)$$

If both bodies are described by the constitutive law outlined in equations (37)–(39), one can show by the same arguments used without contact that

$$G(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2}) = \frac{1}{\Delta t} (E_{n+1}^{\text{tot}} - E_{n}^{\text{tot}}) \quad (44)$$

In view of (43), it is clear that the following must be satisfied for energy to be conserved:

$$G_c(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2}) = - \int_{\Gamma^{(1)}} \mathbf{t}^{(1)} \cdot [\mathbf{V}_{n+1/2}^{(1)}(\mathbf{X}) - \mathbf{V}_{n+1/2}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X}))] d\Gamma^{(1)} = 0 \quad (45)$$

Since we consider frictionless response, and have already made the condition that the basis for $\mathbf{t}^{(1)}$ will be defined in the $n+\frac{1}{2}$ configuration, we can simplify $G_c(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2})$ as

$$\begin{aligned} G_c(\boldsymbol{\phi}_{n+1/2}, \mathbf{V}_{n+1/2}) &= \int_{\Gamma^{(1)}} -t_N \mathbf{v}_{n+1/2} \cdot [\mathbf{V}_{n+1/2}^{(1)}(\mathbf{X}) - \mathbf{V}_{n+1/2}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X}))] d\Gamma^{(1)} \\ &= \int_{\Gamma^{(1)}} t_N \tilde{g}_{n+1/2}(\mathbf{X}) d\Gamma^{(1)} \end{aligned} \quad (46)$$

where $\tilde{g}_{n+1/2}(\mathbf{X})$, an algorithmic gap rate, is defined as

$$\begin{aligned} \tilde{g}_{n+1/2}(\mathbf{X}) &:= -\mathbf{v}_{n+1/2} \cdot [\mathbf{V}_{n+1/2}^{(1)}(\mathbf{X}) - \mathbf{V}_{n+1/2}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X}))] \\ &= -\frac{1}{\Delta t} \mathbf{v}_{n+1/2} \cdot [(\boldsymbol{\phi}_{n+1}^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_{n+1}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X}))) \\ &\qquad -(\boldsymbol{\phi}_{n}^{(1)}(\mathbf{X}) - \boldsymbol{\phi}_{n}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X})))] \end{aligned} \quad (47)$$

Finally, we can conclude from (46) that if the following algorithmic persistency condition is satisfied pointwise on $\Gamma^{(1)}$, then global energy will be conserved for the two-body contact system:

$$t_N \tilde{g}_{n+1/2} = 0 \quad (48)$$

We now have the algorithmic counterpart of (5)$_4$ appropriate for energy conservation. Note that the remainder of the contact conditions, i.e. equations (5)$_{1-3}$, are to this point unenforced by our algorithmic formulation. We discuss some alternatives for carrying out this step in the next section.

**Remark.** Examination of equations (47) and (48) shows that the pointwise contact condition necessary to obtain conservation properties only involves quantities associated with $\mathbf{X}$ and $\mathbf{Y}_{n+1/2}(\mathbf{X})$. In particular, no history terms associated with the point $\mathbf{X}$ need be stored and utilized in calculation of $\tilde{g}_{n+1/2}$, with only algorithmic material velocities and the surface normal at the $n+\frac{1}{2}$ configuration being involved. This is particularly advantageous in large deformation, large slip problems, which are characterized by frequent changes in the elements contacted by individual points $\mathbf{X}$. The fact that the conservation properties are unaffected by such events is a crucial feature of the algorithm we propose.
## 4. DISCRETE CONTACT CONSTRAINTS

In the following we discuss three alternative formulations of the contact constraints appropriate for the above conservative framework. In view of the role of the normal pressure-gap rate complementarity condition in energy conservation, we propose the following expression of the contact conditions in the continuum setting. In contrast to (5), they emphasize complementarity conditions between the gap rate and the pressure:

IF $g < 0$ THEN

$$t_N = 0 \tag{49}$$

$\dot{g}$ is unconstrained

ELSE IF $g = 0$ THEN

$$\begin{gathered} t_N \ge 0 \\ \dot{g} \le 0 \end{gathered} \tag{50}$$

$$\dot{g}t_N = 0$$

### ENDIF

We note that $g$ is precluded from becoming positive by equation (50)$_2$, which simply states that if $g = 0$ the only change that can occur is a negative one. In the continuum case, one can verify that conditions (49) and (50) place the same physical restrictions as equations (5). However, these new conditions appear to be more readily accommodated into a conservation context, as we demonstrate below.

### 4.1. Lagrange multiplier formulation

In the time discrete setting we now consider, we let $t_N = \lambda_N$, where $\lambda_N$ is the Lagrange multiplier field on $\Gamma^{(1)}$ satisfying:

IF $g_n < 0$ THEN

$$\lambda_N = 0 \tag{51}$$

$\tilde{g}_{n+1/2}$ is unconstrained

ELSE IF $g_n \ge 0$ THEN

$$\begin{gathered} \lambda_N \ge 0 \\ \tilde{g}_{n+1/2} \le 0 \end{gathered} \tag{52}$$

$$\tilde{g}_{n+1/2} \lambda_N = 0$$

### ENDIF

where $\tilde{g}_{n+1/2}$ is as defined in (47) and $g_n$ is defined via

$$g_n := -\mathbf{v}_n \cdot (\mathbf{\Phi}_n^{(1)}(\mathbf{X}) - \mathbf{\Phi}_n^{(2)}(\mathbf{Y}_n(\mathbf{X}))) \tag{53}$$

We note that the impenetrability condition $g \le 0$ is only enforced in the limit as $\Delta t \to 0$; in particular, examination of (51) and (52) will show that computation of a non-zero contact traction
$t_N$ will lag one time step behind the detection of an interpenetration. Conditions (52) will then act to preclude further violation of the constraint by directly enforcing the gap rate to be non-positive. Note that this prescription, when used in conjunction with the global algorithm of the last section, preserves all the conservation properties discussed. In particular, energy is conserved because $\lambda_N \tilde{g}_{n+1/2} = 0$ at all contact points and for all conditions of contact.

## 4.2. A penalty regularization

Equations (51) and (52) have the advantage that they are readily regularized, enabling a penalty solution which in some cases may be more convenient. Consider the following prescription for $t_N$:

$$t_N = H(g_n) \epsilon_N \langle \tilde{g}_{n+1/2} \rangle \quad (54)$$

where $\epsilon_N > 0$ is a penalty parameter (physically speaking, a viscosity) and $\langle \cdot \rangle$ denotes the MacAuley bracket, or positive part of the operand. One notes that in the limit as $\epsilon_N \to \infty$ conditions (51) and (52) are reproduced. However, it is no longer true that $t_N \tilde{g}_{n+1/2} = 0$, so we cannot expect exact energy conservation in this case. This makes sense; in fact, penalty methods rely on the finite, non-zero energy associated with the penalization to enforce the constraint as the penalty is increased. We would like to ensure, however, that use of (54) with the otherwise conservative framework does not increase energy. This is readily verified by using equations (43), (44) and (46):

$$ 
\begin{align}
E_{n+1}^{\text{tot}} - E_n^{\text{tot}} &= -\Delta t \int_{\Gamma^{(1)}} t_N \tilde{g}_{n+1/2} \, d\Gamma^{(1)} \nonumber \\
&= -\Delta t \int_{\Gamma^{(1)}} H(g_n) \epsilon_N \langle \tilde{g}_{n+1/2} \rangle \tilde{g}_{n+1/2} \, d\Gamma^{(1)} \tag{55} \\
&= -\Delta t \int_{\Gamma^{(1)}} H(g_n) \epsilon_N \langle \tilde{g}_{n+1/2} \rangle^2 \, d\Gamma^{(1)} \nonumber \\
&\le 0 \nonumber
\end{align}
$$

which shows that the algorithm is unconditionally dissipative. From an energy standpoint the integrator thus remains stable. The dissipated energy can also be restored via augmented Lagrangian iteration, as discussed below.

## 4.3. Augmented Lagrangian iteration

Consider the following expression for $t_N$, based on an augmented Lagrangian augmentation of the penalization in (54):

$$t_N = H(g_n) \langle \lambda_N^{(k)} + \epsilon_N \tilde{g}_{n+1/2} \rangle \quad (56)$$

where $\lambda_N^{(k)}$ is a fixed iterate for the Lagrange multiplier satisfying (49) and (50). One may therefore consider solving the following iterative sequence of problems within a time step, in which we begin with some initial estimate $\lambda_N^{(0)}$ of the multipliers on $\Gamma^{(1)}$, and proceed with iterations $(k)$ until some tolerance is satisfied:

1. Solve for $\Phi_{n+1}^{(k)}$, with $\lambda_N^{(k)}$ fixed:

$$G(\Phi_{n+1/2}^{(k)}, \Phi^*) + G_c(\Phi_{n+1/2}^{(k)}, \Phi^*; \lambda_N^{(k)}) = 0 \quad (57)$$
where $G(\cdot, \cdot)$ is as defined in (41), $\varphi_{n+1/2}^{(k)}$ is as defined in (35), and $G_c$ is defined via

$$ G_c(\varphi_{n+1/2}^{(k)}, \varphi^*; \lambda_N^{(k)}) = \int_{\Gamma^{(1)}} -H(g_n)\langle \lambda_N^{(k)} + \varepsilon_N \tilde{g}_{n+1/2}^{(k)} \rangle \cdot \mathbf{v}_{n+1/2}^{(k)} \cdot (\varphi^*_{(1)}{}^{(k)}(\mathbf{X}) - \varphi^*_{(2)}{}^{(k)}(\mathbf{Y}_{n+1/2}(\mathbf{X}))) d\Gamma^{(1)} \quad (58) $$

2. Update the multipliers on $\Gamma^{(1)}$:

$$ \lambda_N^{(k+1)} = H(g_n)\langle \lambda_N^{(k)} + \varepsilon_N \tilde{g}_{n+1/2}^{(k)} \rangle \quad (59) $$
$$ k \leftarrow k + 1 $$

3. IF convergence achieved EXIT, ELSE return to 1.

As we will show in the numerical examples, this algorithm can reproduce conditions (51) and (52) arbitrarily closely as iterations ($k$) proceed. Several convergence criteria might be proposed to define step 3; for example, the error in $E_{n+1}^{\text{tot}} - E_n^{\text{tot}}$ might be a reasonable choice for defining a convergence tolerance for these iterations.

## 5. SPATIAL DISCRETIZATION AND IMPLEMENTATION

In this section we briefly present an overview of the spatial discretization process which, when applied to the conserving algorithm of the last section, produces the same conservation properties in the fully discrete setting.

### 5.1. Finite element discretization

In general terms, the finite element discretization is achieved by introducing $\varphi_t^h$ and $\varphi_t^{*h}$, the finite dimensional approximations of $\varphi_t$ and $\varphi_t^*$. These lie in the discrete spaces $C_t^h$ and $\mathcal{V}^h$ such that

$$ \varphi_t^h = \sum_{A=1}^{n_{\text{nod}}} N_A \mathbf{d}_{A_t} \in C_t^h \subset \mathcal{C}_t \quad (60) $$

and

$$ \varphi^*_{t}^{h} = \sum_{A=1}^{n_{\text{nod}}} N_A \mathbf{c}_A \in \mathcal{V}^h \subset \mathcal{V} \quad (61) $$

where $\mathbf{d}_{A_t}$ is the vector-valued nodal value of the configuration mapping, $N_A$ is a finite element shape function with domain $\Omega^{(1)} \cup \Omega^{(2)}$, and $\mathbf{c}_A$ are nodal constants. Substitution of these finite dimensional approximations into the time discrete weak form and enforcing it for arbitrary combinations of $\mathbf{c}_A$ yields the following fully discrete equations to be solved in each time step $[t_n, t_{n+1}]$:

$$ \mathbf{M}_{\Delta t} \mathbf{1} [\mathbf{V}_{n+1} - \mathbf{V}_n] + \mathbf{F}_{n+1/2}^{\text{int}} + \mathbf{F}_{c_{n+1/2}} = \mathbf{0} \quad (62) $$
$$ \mathbf{1}_{\Delta t} [\mathbf{d}_{n+1} - \mathbf{d}_n] = \mathbf{V}_{n+1/2} $$

where $\mathbf{M}$ is the global mass matrix, $\mathbf{d}$ is the global solution vector, $\mathbf{F}_c$ is the contact force vector, and $\mathbf{F}^{\text{int}}$ is the internal force vector. We retain the discrete counterparts of $\varphi_{n+1/2}$ and $\mathbf{V}_{n+1/2}$ used previously, e.g.:

$$ \mathbf{d}_{n+1/2} = \frac{1}{2} [\mathbf{d}_{n+1} + \mathbf{d}_n] \quad (63) $$
$$ \mathbf{V}_{n+1/2} = \frac{1}{2} [\mathbf{V}_{n+1} + \mathbf{V}_n] $$
and remark that the nodal force $\mathbf{F}_{n+1/2}^{\text{int}}$ associated with node A is given by

$$ \mathbf{F}_{n+1/2}^{\text{int}} = \int_{\Omega^{(1)} \cup \Omega^{(2)}} D\boldsymbol{\phi}_{n+1/2}^{h} \mathbf{S}^{h} \text{GRAD}[N_A] \, d\Omega \quad (64) $$

where $\mathbf{S}^h$ is as defined in (37)–(39).

In giving a general expression for $\mathbf{F}_c$ we summarize the main results, referring the reader to Laursen and Simo¹¹ for more details on the general procedure. The integral $G_c$ can be approximated via

$$ G_c(\boldsymbol{\phi}_{n+1/2}^h, \boldsymbol{\phi}^h) \approx - \sum_{l=1}^{n_{\text{int}}} W^l j^l [t_N \mathbf{v}_{n+1/2} \cdot (\boldsymbol{\phi}^{(1)}(\mathbf{X}) - \boldsymbol{\phi}^{(2)}(\mathbf{Y}_{n+1/2}(\mathbf{X})))]^l \quad (65) $$

where $n_{\text{int}}$ is the total number of contact quadrature points on $\Gamma^{(1)}$, $W^l$ is a weight of integration for quadrature point $l$, and $j^l$ is the jacobian resulting from the local to global transformation used to describe $\Gamma^{(1)}$. $\mathbf{F}_c$ can then be expressed as an assembly of individual quadrature point contributions as follows:

$$ \mathbf{F}_{c_{n+1/2}} = \bigwedge_{l=1}^{n_{\text{int}}} W^l j^l \mathbf{f}_{c_{n+1/2}}^l \quad (66) $$

where $\mathbf{A}$ is the standard finite element assembly operator. In the case of interest, frictionless contact, $\mathbf{f}_{c_{n+1/2}}^l$ takes the following form if nodal quadrature is used in (65):

$$ \mathbf{f}_{c_{n+1/2}}^l = -t_N^l \mathbf{n}^l \quad (67) $$

where $\mathbf{n}^l = [\mathbf{v}_{n+1/2}^l, -N_1(\xi_{n+1/2}^l)\mathbf{v}_{n+1/2}^l, \dots, -N_{\text{nel}}(\xi_{n+1/2}^l)\mathbf{v}_{n+1/2}^l]^T$ with $N_a$, $a = 1, \dots, \text{nel}$ representing the shape functions interpolating the element surface containing the projection $\mathbf{Y}_{n+1/2}(\mathbf{X}^l)$. The contact pressure $t_N$ can be described by any of the representations given in the last section. Calculation of the contact stiffness requires exact linearization of the contact force vector. We omit this calculation here, and refer the interested reader to Laursen and Simo¹¹ for elaboration on the general procedure.

## 6. NUMERICAL EXAMPLES

In presenting some results obtained with the proposed algorithm, we give comparisons with results obtained via two other prevalent strategies: the Newmark method¹² and the Hilber–Hughes–Taylor, or HHT method.¹³ As frequently implemented for contact problems, these integrators could be summarized as follows:

Newmark

$$ \begin{align} & \mathbf{M}\mathbf{a}_{n+1} + \mathbf{F}^{\text{int}}(\mathbf{d}_{n+1}) + \mathbf{F}_c(\mathbf{d}_{n+1}) = \mathbf{F}_{n+1}^{\text{ext}} \notag \\ & \mathbf{d}_{n+1} = \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2} [(1 - 2\beta)\mathbf{a}_n + 2\beta \mathbf{a}_{n+1}] \tag{68} \\ & \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta t [(1 - \gamma)\mathbf{a}_n + \gamma \mathbf{a}_{n+1}] \notag \end{align} $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Longitudinal impact between two identical slender bars</figcaption></figure>

HHT

$$
\begin{align}
\mathbf{M}\mathbf{a}_{n+1} + \mathbf{F}^{\text{int}}(\mathbf{d}_{n+\alpha}) + \mathbf{F}_c(\mathbf{d}_{n+\alpha}) &= \mathbf{F}_{n+\alpha}^{\text{ext}} \nonumber \\
\mathbf{d}_{n+\alpha} &= \alpha\mathbf{d}_{n+1} + (1-\alpha)\mathbf{d}_n \nonumber \\
\mathbf{d}_{n+1} &= \mathbf{d}_n + \Delta t \mathbf{v}_n + \frac{\Delta t^2}{2} [(1-2\beta)\mathbf{a}_n + 2\beta\mathbf{a}_{n+1}] \tag{69} \\
\mathbf{v}_{n+1} &= \mathbf{v}_n + \Delta t [(1-\gamma)\mathbf{a}_n + \gamma\mathbf{a}_{n+1}] \nonumber
\end{align}
$$

where $\alpha$, $\beta$ and $\gamma$ represent algorithmic parameters throughout. As discussed in the introduction, the contact constraints are commonly imposed in both strategies without explicit consideration of rate conditions, so that $F_c$ is assembled from tractions $t_N$ satisfying

$$ g \le 0, \quad t_N \ge 0, \quad t_N g = 0 \tag{70} $$

where $g$ in equations (70) is evaluated from the $n+1$ configuration in the case of Newmark and from the $n+\alpha$ configuration in the case of HHT.

## 6.1. Impact between identical bars

In the first example we consider the axial impact of two identical straight elastic bars, which provides a simple yet illustrative demonstration of the algorithm's performance. In the initial state the two bars are collinear but out of contact as shown in Figure 1. Bar A is initially moving at a uniform velocity of 1 unit while bar B is at rest initially. The properties are: density $\rho = 1$ unit, area of cross section $a_c = 1$ unit, length $l = 10$ unit, Young's modulus $E = 1$ unit and Poisson's ratio $\nu = 0$. The problem is driven only by the initial conditions with different amounts of initial separation considered between the bars.

The solution for displacements and velocities of the impacting ends of the two bars based on physical observations is plotted in Figure 2 for the case of small strains. The numerical solution is obtained using small strain bilinear elastic elements, with 100 elements in each bar and with a time step of $\Delta t = 0.1$ unit. We consider four algorithms for the temporal integration: the proposed conservative scheme, the trapezoidal rule (second-order accurate Newmark with $\beta = 0.25$, $\gamma = 0.5$), Newmark with maximum high-frequency dissipation ($\beta = 0.49$, $\gamma = 0.9$) and HHT with optimal high-frequency damping ($\beta = 0.3025$, $\gamma = 0.6$, $\alpha = 0.9$). Solutions for these cases are shown in Figures 3–6. Unless otherwise stated, the method of augmented Lagrangians is used with iterations until the change in the multiplier value is less than 1 per cent.

Examination of the results shows that dissipative Newmark (Figure 4) and HHT (Figure 5) produce oscillations shortly after the initial contact event, characterized by the two bars coming into and out of contact. The high-frequency dissipation of these two integrators eventually damps
<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Displacement vs. time, (b) velocity vs. time at the impact ends of the bars (expected solution)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. (a) Displacement vs. time, (b) velocity vs. time at node *b* and *d* obtained using the trapezoidal rule with penalty method ($\epsilon_N = 10^8$)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 4. (a) Displacement vs. time and (b) velocity vs. time for nodes *b* and *d* obtained using the dissipative Newmark method</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. (a) Displacement vs. time and (b) velocity vs. time for nodes *b* and *d* obtained using the HHT method</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6. (a) Displacement vs. time and (b) velocity vs. time for nodes *b* and *d* obtained using the conservative scheme</figcaption></figure>

out this behaviour, so that a reasonable solution is produced. The trapezoidal rule possesses no such dissipation, so that these numerically induced contact oscillations persist throughout the solution (Figure 3). Notably, the conservative algorithm possesses no dissipation either, and yet does not excite these spurious modes due to the conservative treatment of the contact conditions (Figure 6). In fact, one can show that up to the treatment of the contact conditions, the trapezoidal rule and the conservative algorithm are identical for small-strain elasticity problems. We therefore conclude that the key issue in this problem is the stable treatment of the contact constraints. The slight interpenetration of the bar ends evident in Figure 6 also tells us that insistence on complete impenetrability is unnecessary for accurate prediction of energy and momentum transfer across the interface.

Further insight can be gained by looking at Figure 7, which shows the total system energy for the four cases. Clearly, use of constraints (70) with the trapezoidal rule results in a net energy gain with each new impact event. This energy remains in the system and continues to grow until the impacts stop, after which the system energy is constant. The dissipative integrators damp out the contact oscillations as discussed above, but lose system energy in doing so. The conservative algorithm obtains results no more noisy than those provided by the dissipative integrators, while conserving system energy. Importantly, the conservative algorithm in no way forces the contact to be persistent, but instead allows the conditions to evolve as driven by the momentum equations. It would therefore appear that formulations requiring such persistency of new contacts, as in Taylor and Papadopoulos,⁸ may pose needless constraints which could in some problems be non-physical.

Finally, all integrators display oscillations after the separation of bars. This results from the numerical dispersion of the solution near discontinuities causing the higher-frequency oscillations to travel at lower speeds. A common way to deal with such numerical dispersion is to damp out
these oscillations with artificial viscosity (see for example References 14 and 15). An interesting issue we intend to pursue in future work is the incorporation of conserving contact treatments into dissipative global schemes, creating algorithms with high-frequency dissipation which integrate the contact conditions in a stable manner.

<figure><img src="image_10.png" /><figcaption>Figure 7. Total energy vs. time in the bar impact problem for (a) trapezoidal rule, (b) conservative scheme, dissipative Newmark and HHT method</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Ring configurations at intervals of 5 time units for (a) trapezoidal rule; (b) dissipative Newmark; and (c) conservative scheme</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 8. Continued</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9. Total energy vs. time for the conservative scheme, the trapezoidal rule and dissipative Newmark; ring problem</figcaption></figure>

## 6.2. Impact of a ring against rigid surface

This problem is motivated by a system originally discussed in Wriggers *et al.*¹⁶ An elastic ring is thrown with an initial velocity of 2 units at a 45° angle to a flat rigid surface, as depicted in Figure 8. The material properties and dimensions of the ring are as follows: Young's modulus $E = 10^2$ units, density $ρ = 0.01$ unit, Poisson's ratio $v = 0.0001$, outer radius $r_o = 10$ units and inner radius $r_i = 9$ units. Bilinear large strain elements are used to discretize the ring, with the mesh consisting of 64 elements. Results are obtained for the trapezoidal rule, the dissipative Newmark method and the proposed conservative strategy, using a time step $Δt = 0.2$ in all cases. The ring configurations at various time intervals and the total energy are plotted in Figures 8 and 9, respectively.

The configuration plot shows large deformations during and after impact and an unstable blow-up of the simulation for the trapezoidal rule. The dissipative Newmark method provides reasonable
<figure><img src="image_14.png" /></figure>

(a)

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

(b)

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

(c)

<caption>Figure 10. Striker positions at various times for (a) the trapezoidal rule, (b) dissipative Newmark and (c) the conservative scheme</caption>

results for small times, but damps out the low-mode deformation of the ring for larger times, preserving only the rigid-body motion with any accuracy. The conservative scheme is seen to be clearly superior to either of the above two approaches. Examination of Figure 9, a plot of the system energy for all three simulations, reveals clearly the energy instability in the trapezoidal rule simulation, the loss of energy for the dissipative scheme, and the exact conservation for the current method.

### 6.3. The carrom problem

This problem involves perfectly elastic and frictionless impact between rigid bodies. The system can be visualized as a 2-D version of a pool table, where a plastic disk called the striker is
<figure><img src="image_17.png" /><figcaption>Figure 11. Total energy vs. time for the conservative scheme, the trapezoidal rule and dissipative Newmark; carrom problem</figcaption></figure>

manually struck with the aim of hitting other disks and causing them to go into holes in the four corners of the carrom board. The model simulates the motion of the striker in the absence of other disks as it hits the sides of the carrom board assuming no friction (see Figure 10). The sides of the board are simulated with four elements and the striker is placed within the closed space and subjected to an initial velocity. The material properties (contrived so that deformations will be negligible) and dimensions are as follows: striker—Young's modulus $E = 10^2$ units, Poisson's ratio $v = 0.0001$, density $\rho = 0.1$ unit and radius $r = 1$ unit; carrom board—Young's modulus $E = 10^4$ units, Poisson's ratio $v = 0.0001$, density $\rho = 10^2$ units, length of inner side $l_i = 13$ units and length of outer side $l_o = 15$ units. A time step of $\Delta t = 4$ units was utilized in each of three simulations, performed using the trapezoidal rule, dissipative Newmark and the conservative scheme.

The initial condition causes the striker to hit the lower side in the middle and at an angle of 45°, with a velocity of 0.1 units. The results in Figure 10 show that in case of the trapezoidal rule the rebound is too sharp due to gain of energy; with an eventual blow-up of the solution. The dissipative Newmark result has the opposite effect, with the angle of rebound being too flat and the system energy being reduced. Finally, the conservative scheme correctly predicts all rebounds to be at 45° after the initial impact, giving the expected diamond-shaped trajectory of the striker within the carrom board. The plot of total system energy in Figure 11 is consistent with these observations, and again verifies the energy conservation property.

## 7. SUMMARY AND CONCLUSIONS

This work extends the idea of algorithmic energy and momentum conservation to encompass systems featuring mechanical contacts. Importantly, these conservation properties hold for all changing conditions of contact: new contact, release and persistent. The result is an algorithm which is stable in the complete absence of any dissipation. Although the framework we have developed is
by definition only applicable to conservative systems featuring frictionless contacts, the work has
revealed that careful consideration of algorithmic contact work input is crucial to the construction
of stable integration procedures.

It is fair to observe that most physical systems of interest are not conservative, and that few
impact events are reasonably modelled as perfectly elastic. Effective algorithms are needed for
problems featuring physical dissipation arising from inelasticity and/or interface friction. Further-
more, discretization of continuous systems introduces non-physical modes whose energy content
should ordinarily be damped out in transient simulations. When viewed in this light, perhaps the
most important product of this work is its identification of algorithmic persistency conditions that
enable stable numerical integration through conditions of changing contact. We intend to extend
this work so that interface dissipation, consistent with the frictional model used, can be accurately
produced by the algorithm in each time step of a transient simulation. Future work will also explore
the use of such energetically consistent contact formulations in conjunction with otherwise dissip-
tive global integrators (such as HHT). In this manner both physically and numerically dissipative
equations could be integrated without spurious energy input from the contact conditions.

ACKNOWLEDGEMENTS

T. Laursen's work on this project was supported in part by the National Science Foundation under grant number MSS-9308486, while V. Chawla's work was supported by Sandia National Laboratories under contract number AN-9563 with Duke University. This support is gratefully acknowledged.

REFERENCES

1. G. Duvaut and J. L. Lions, *Les Inéquations en Mecanique et en Phisique*, Dunod, Paris, 1972.

2. N. Kikuchi and J. T. Oden, *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods*, SIAM, Philadelphia, 1988.

3. J. O. Hallquist, G. L. Goudreau and D. J. Benson, 'Sliding interfaces with contact-impact in large-scale Lagrangian computations', *Comp. Methods Appl. Mech. Eng.*, **51**, 107-137 (1985).

4. D. J. Benson and J. O. Hallquist, 'A single surface contact algorithm for the post-buckling analysis of shell structures', *Comp. Methods Appl. Mech. Eng.*, **78**, 141-163 (1990).

5. T. Belytschko and M. O. Neal, 'Contact-impact by the pinball algorithm with penalty and Lagrangian methods', *Int. j. numer. methods eng.*, **31**, 547-572 (1991).

6. N. J. Carpenter, R. L. Taylor and M. G. Katona, 'Lagrange constraints for transient finite element surface contact', *Int. j. numer. methods eng.*, **32**, 103-128 (1991).

7. T. J. R. Hughes, R. L. Taylor, J. L. Sackman, A. Curnier and W. Kanoknukulchai, 'A finite element method for a class of contact-impact problems', *Comp. Methods Appl. Mech. Eng.*, **8**, 249-276 (1976).

8. R. L. Taylor and P. Papadopoulos, 'On a finite element method for dynamic contact/impact problems', *Int. j. numer. methods eng.*, **36**, 2123-2140 (1993).

9. Z.-H. Zhong, *Finite Element Procedures for Contact-Impact Problems*, Oxford Science Publications, Oxford, 1993.

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

11. T. A. Laursen and J. C. Simo, 'A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems', *Int. j. numer. methods eng.*, **36**, 3451-3485 (1993).

12. N. M. Newmark, 'A method of computation for structural dynamics', *J. Eng. Mech. Div.*, **ASCE**, **85**, 67-94 (1959).

13. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, 'Improved numerical dissipation for time integration algorithms in structural dynamics', *Earthquake Eng. Struct. Dyn.*, **5**, 283-292 (1977).

14. J. Von Neumann and R. D. Richtmyer, 'A method for the numerical calculation of hydrodynamic shocks', *J. Appl. Phys.*, **21**, 232 (1950).

15. D. J. Benson, 'A new two-dimensional flux-limited shock viscosity for impact calculations', *Comp. Methods Appl. Mech. Eng.*, **93**, 39-95 (1991).

16. P. Wriggers, T. Vu Van and E. Stein, 'Finite element formulation of large deformation impact-contact problems with friction', *Comput. Struct.*, **37**, 319-331 (1990).