# Exact energy and momentum conserving algorithms for general models in nonlinear elasticity

O Gonzalez

► To cite this version:

O Gonzalez. Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, 2000, 190 (13-14), pp.1763-1783. 10.1016/S0045-7825(00)00189-4. hal-01363585

HAL Id: hal-01363585

https://hal.science/hal-01363585

Submitted on 10 Sep 2016

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
# Exact energy and momentum conserving algorithms for general models in nonlinear elasticity

O. Gonzalez

Département de Mathématiques, Ecole Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland

Implicit time integration schemes that inherit the conservation laws of total energy, linear and angular momentum are considered for initial boundary-value problems in finite-deformation elastodynamics. Conserving schemes are constructed for general hyperelastic material models, both compressible and incompressible, and are formulated in a way that is independent of spatial discretization. Three numerical examples for Ogden-type material models, implemented using a finite element discretization in space, are given to illustrate the performance of the proposed schemes. These examples show that, relative to the standard implicit mid-point rule, the conserving schemes exhibit superior numerical stability properties without a compromise in accuracy. © 2000 Elsevier Science B.V. All rights reserved.

Numerical integration; Nonlinear elastodynamics; Incompressible elasticity; Integral preservation

## 1. Introduction

In this article we consider energy and momentum conserving time discretization schemes for general initial boundary-value problems in finite-deformation elastodynamics. Conserving schemes are developed in a weak form for both compressible and incompressible hyperelastic material models, and are independent of any particular spatial discretization. Strong forms of the schemes presented here can be obtained using standard arguments involving integration by parts. For finite-deformation, compressible elastodynamics we consider the weak formulation

$$ 
\left.
\begin{aligned}
\int_{\Omega} \dot{\boldsymbol{\varphi}} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \mathbf{V} \cdot \boldsymbol{\eta} \, d\Omega \\
\int_{\Omega} \rho_0 \dot{\mathbf{V}} \cdot \boldsymbol{\eta} \, d\Omega + \int_{\Omega} \mathbf{S}(\boldsymbol{\varphi}) : \mathbf{F}(\boldsymbol{\varphi})^{\mathrm{T}} \nabla \boldsymbol{\eta} \, d\Omega \\
&= \int_{\Omega} \mathbf{b} \cdot \boldsymbol{\eta} \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \boldsymbol{\eta} \, d\Gamma \\
\int_{\Omega} \mathbf{V}|_{t=0} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \mathbf{V}_0 \cdot \boldsymbol{\eta} \, d\Omega \\
\int_{\Omega} \boldsymbol{\varphi}|_{t=0} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \boldsymbol{\varphi}_0 \cdot \boldsymbol{\eta} \, d\Omega
\end{aligned}
\right\} \quad \forall \boldsymbol{\eta} \in \mathcal{V}, \qquad (1)
$$

and for finite-deformation, incompressible elastodynamics
$$
\left.
\begin{aligned}
\int_{\Omega} \dot{\boldsymbol{\phi}} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \mathbf{V} \cdot \boldsymbol{\eta} \, d\Omega \\
\int_{\Omega} \rho_0 \dot{\mathbf{V}} \cdot \boldsymbol{\eta} \, d\Omega + \int_{\Omega} [\mathbf{S}(\boldsymbol{\varphi}) + 2\lambda DG(\mathbf{C}(\boldsymbol{\varphi}))] : \mathbf{F}(\boldsymbol{\varphi})^{\mathrm{T}} \nabla \boldsymbol{\eta} \, d\Omega \\
&= \int_{\Omega} \mathbf{b} \cdot \boldsymbol{\eta} \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \boldsymbol{\eta} \, d\Gamma \\
0 &= \int_{\Omega} \phi G(\mathbf{C}(\boldsymbol{\varphi})) \, d\Omega \\
\int_{\Omega} \mathbf{V}|_{t=0} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \mathbf{V}_0 \cdot \boldsymbol{\eta} \, d\Omega \\
\int_{\Omega} \boldsymbol{\varphi}|_{t=0} \cdot \boldsymbol{\eta} \, d\Omega &= \int_{\Omega} \boldsymbol{\varphi}_0 \cdot \boldsymbol{\eta} \, d\Omega
\end{aligned}
\right\}
\quad
\left\{
\begin{array}{l}
\forall \boldsymbol{\eta} \in \mathcal{V} \\
\forall \phi \in \mathcal{E}
\end{array}
\right.
\qquad (2)
$$

where $\Omega \subset \mathbb{R}^3$ is an open set that represents the reference state of a continuum body with material points $\mathbf{X} \in \Omega$, $\boldsymbol{\varphi}(\mathbf{X}, t) \in \mathbb{R}^3$ is the deformation field for the body relative to its reference state, $\mathbf{V}(\mathbf{X}, t) \in \mathbb{R}^3$ is the material velocity field, $\mathbf{F}(\mathbf{X}, t) := \nabla\boldsymbol{\varphi}(\mathbf{X}, t) \in \mathbb{R}^{3 \times 3}$ is the deformation gradient field and $\mathbf{S}(\mathbf{X}, t) \in \mathbb{R}^{3 \times 3}$ is the second Piola–Kirchhoff constitutive stress field. Here $\nabla$ denotes the gradient operator relative to the Cartesian coordinates $\mathbf{X}$ in $\Omega$, an overdot denotes differentiation with respect to time $t \ge 0$ and $\mathbb{R}^{3 \times 3}$ is the vector space of all real $3 \times 3$ matrices, equipped with the standard (Euclidean) inner product $\mathbf{A} : \mathbf{B} = A_{ij}B_{ij}$ where summation on repeated indices is implied.

For hyperelastic bodies, the stress field $\mathbf{S}$ is connected to the deformation field $\boldsymbol{\varphi}$ via a local constitutive relation of the form

$$ \mathbf{S}(\boldsymbol{\varphi}) = 2DW(\mathbf{C}(\boldsymbol{\varphi})) \qquad (3) $$

where $W(\mathbf{C}) \in \mathbb{R}$ is a given strain energy density function, $DW(\mathbf{C})$ is the derivative of $W$ with respect to its tensor argument evaluated at $\mathbf{C}$, and $\mathbf{C}(\boldsymbol{\varphi})$ is the right Cauchy–Green stretch field associated with the deformation $\boldsymbol{\varphi}$, namely, $\mathbf{C} = \mathbf{F}^{\mathrm{T}}\mathbf{F}$. The strain energy density function $W$ may depend upon the material point $\mathbf{X}$, but this dependence is not shown for clarity of notation. In the incompressible case, the total second Piola–Kirchhoff stress field is of the form

$$ \mathbf{S} + \det[\mathbf{C}]^{1/2} \lambda \mathbf{C}^{-1} \qquad (4) $$

where $\lambda(\mathbf{X}, t) \in \mathbb{R}$ is a material pressure field that enforces the material incompressibility constraint

$$ G(\mathbf{C}) = \det[\mathbf{C}]^{1/2} - 1 = 0 \quad \text{in } \Omega. \qquad (5) $$

Since $DG(\mathbf{C}) = \frac{1}{2}\det[\mathbf{C}]^{1/2}\mathbf{C}^{-1}$, we note that the total second Piola–Kirchhoff stress field (4) can be written in the form $\mathbf{S} + 2\lambda DG(\mathbf{C})$.

In the compressible case, we seek functions $(\boldsymbol{\varphi}, \mathbf{V}) : [0, T] \to \mathcal{U} \times \mathcal{V}$ satisfying (1) where $\mathcal{U}$ is a set of deformations satisfying essential boundary conditions on a given boundary patch $\Gamma_d$, $\mathcal{V}$ is a vector space of functions satisfying the homogeneous counterpart of the essential boundary conditions, $\rho_0(\mathbf{X}) > 0$ is a specified mass density field, $\mathbf{b}(\mathbf{X}, t) \in \mathbb{R}^3$ is a specified body force field, $\mathbf{f}(\mathbf{X}, t) \in \mathbb{R}^3$ is a specified traction field on $\Gamma_d = \partial\Omega \setminus \Gamma_d$, $\mathbf{V}_0(\mathbf{X})$ is a specified initial material velocity field and $\boldsymbol{\varphi}_0(\mathbf{X})$ is a specified initial deformation. In the incompressible case, we seek functions $(\boldsymbol{\varphi}, \mathbf{V}, \lambda) : [0, T] \to \mathcal{U} \times \mathcal{V} \times \mathcal{E}$ satisfying (2) where $\mathcal{E}$ is a vector space of pressure fields restricted only by appropriate smoothness requirements. Various well-posedness results have been established for both (1) and (2), and associated linear systems. See [5, 7, 13, 24, 27, 31, 36] for well-posedness results pertaining to (1), and [9–12, 22, 26, 37] for results pertaining to (2). For more details on the field equations of finite-deformation elasticity see [6, 23, 31, 33, 39, 46].

Depending on boundary conditions and external loads, the nonlinear system (1) may possess various integrals related to the total linear momentum $\mathbf{L}(\boldsymbol{\varphi}, \mathbf{V})$, angular momentum $\mathbf{J}(\boldsymbol{\varphi}, \mathbf{V})$ and energy $H(\boldsymbol{\varphi}, \mathbf{V})$ of the material body. These integrals are typically lost under discretization in time. In [42,43] it was shown that all members of the Newmark family preserve $\mathbf{L}$, but none preserve $\mathbf{H}$ for arbitrary constitutive laws. Moreover, it was shown that the explicit central difference scheme is the only member of the Newmark family that preserves $\mathbf{J}$. The main result in [42,43] was the construction of implicit schemes for (1) that preserve $\mathbf{L}$, $\mathbf{J}$ and $\mathbf{H}$, but at the cost of introducing a nonlinear scalar equation at each time step.

In this article we consider generalizations of the time-reversible, integral-preserving time discretization schemes of [19,20] for the treatment of (1). We show that implicit schemes which preserve $\mathbf{L}$, $\mathbf{J}$ and $\mathbf{H}$ can be constructed for general hyperelastic material models. Moreover, in contrast to the methods developed in
[42,43], the schemes developed here do not involve extra parameters, equations or projection steps. The appeal of preserving an underlying energy integral can be traced back to [35] and the concept of energy stability. Indeed, the computational experiments presented here and in [20,42,43] suggest that certain implicit schemes which preserve L, J and H possess superior stability properties than other standard implicit schemes, such as the trapezoidal rule which generally preserves only L, and the mid-point rule which generally preserves only L and J. For more details on the symmetry and stability properties of integral-preserving schemes see [19,20], and for other applications see [41,44,45].

Due to the incompressibility condition, the initial boundary-value problem (2) for incompressible bodies can be interpreted as a constrained or differential-algebraic evolution equation [3,34] in a function space. Two important sets associated with the constraint (5) are

$$ \mathcal{L} = \{ \boldsymbol{\varphi} \in \mathcal{U} \mid G(\mathbf{C}(\boldsymbol{\varphi})) = 0 \text{ in } \Omega \}, \quad (6) $$

which we call the intrinsic configuration space, and the set

$$ \mathcal{M} = \{ (\boldsymbol{\varphi}, \mathbf{V}) \in \mathcal{U} \times \mathcal{V} \mid G(\mathbf{C}(\boldsymbol{\varphi})) = 0 \text{ and } DG(\mathbf{C}(\boldsymbol{\varphi})) : [\nabla \mathbf{V}^T \nabla \boldsymbol{\varphi} + \nabla \boldsymbol{\varphi}^T \nabla \mathbf{V}] = 0 \text{ in } \Omega \}, \quad (7) $$

which we call the intrinsic state space for the system. Any sufficiently smooth solution of (2) has the property that $\boldsymbol{\varphi}(t) \in \mathcal{L}$ and $(\boldsymbol{\varphi}, \mathbf{V})(t) \in \mathcal{M}$. Moreover, depending on boundary conditions and external loads, the functionals $L|_{\mathcal{M}}$, $J|_{\mathcal{M}}$ and $H|_{\mathcal{M}}$ may be integrals. The formulation (2) is only one of many possible descriptions of incompressible elastodynamics. One may alternatively consider stabilized differential-algebraic formulations along the lines of [3,15,16], or unconstrained formulations in which the constraints appear as invariants along the lines of [1,8,22,28,30]. However, such alternative formulations will not be addressed here.

In this article we seek to construct time discretization schemes for (2) that preserve the sets $\mathcal{L}$ and $\mathcal{M}$, along with the integrals $L|_{\mathcal{M}}$, $J|_{\mathcal{M}}$ and $H|_{\mathcal{M}}$ whenever they are present in the underlying problem. Due to special difficulties associated with direct spatial discretization of (2), we will instead consider a modified quasi-incompressible formulation along the lines of [38,40], and generalize the time-reversible, integral-preserving time discretization schemes developed in [21] to this modified formulation. We show that time discretization schemes can be constructed for this formulation which preserve $\mathcal{L}$, the total energy functional $H|_{\mathcal{L}}$, and the momentum functionals $J|_{\mathcal{L}}$ and $L|_{\mathcal{L}}$. However, the schemes will not generally preserve the velocity-level constraints defined by $\mathcal{M}$. As in the compressible case, conservation will be achieved automatically.

The presentation is structured as follows. In Section 2 we outline conservation laws and associated integrals for (1) and present a time discretization scheme that inherits a discrete version of these laws. In Section 3 we introduce a quasi-incompressible formulation of (2), outline conservation laws and present a conserving time discretization scheme for this modified formulation. In Section 4 we consider the finite element spatial discretization of the time-discretized systems. In the quasi-incompressible case, we introduce a general mixed finite element discretization, reduce the general discretization to a two-field formulation and discuss an efficient solution strategy. In Section 5 we present three numerical examples involving both compressible and incompressible material models of the Ogden type.

## 2. Conserving algorithms for compressible models

### 2.1. Conservation laws and integrals

To any state $(\boldsymbol{\varphi}, \mathbf{V}) \in \mathcal{U} \times \mathcal{V}$ of the system (1) we associate a total linear momentum $\mathbf{L}(\boldsymbol{\varphi}, \mathbf{V})$, a total angular momentum $\mathbf{J}(\boldsymbol{\varphi}, \mathbf{V})$ and a total energy $H(\boldsymbol{\varphi}, \mathbf{V})$ defined as

$$ \left. \begin{aligned} \mathbf{L}(\boldsymbol{\varphi}, \mathbf{V}) &:= \int_{\Omega} \rho_0 \mathbf{V} \, d\Omega, \\ \mathbf{J}(\boldsymbol{\varphi}, \mathbf{V}) &:= \int_{\Omega} \boldsymbol{\varphi} \times \rho_0 \mathbf{V} \, d\Omega, \\ H(\boldsymbol{\varphi}, \mathbf{V}) &:= \int_{\Omega} \frac{1}{2} \rho_0 |\mathbf{V}|^2 + W(\mathbf{C}(\boldsymbol{\varphi})) \, d\Omega - \int_{\Omega} \mathbf{b} \cdot \boldsymbol{\varphi} \, d\Omega - \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \boldsymbol{\varphi} \, d\Gamma. \end{aligned} \right\} \quad (8) $$
Depending on the essential boundary conditions and the external loads, the weak initial boundary-value problem (1) may possess various integrals or constants of motion related to **L**, **J** and **H**. This fact is summarized in the following proposition, whose proof is straightforward.

**Proposition 1.** Let $(\varphi, \mathbf{V})(t)$ be a solution of (1) for given boundary conditions specified in $\mathcal{U}$ and $\mathcal{V}$, and given external loads **b** and **f**. Moreover, let **R** and **T** denote the resultant force and moment due to the external loads, namely,

$$ 
\begin{align} 
\mathbf{R} &= \int_{\Omega} \mathbf{b} \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \, d\Gamma, \\ 
\mathbf{T} &= \int_{\Omega} \varphi \times \mathbf{b} \, d\Omega + \int_{\Gamma_{\sigma}} \varphi \times \mathbf{f} \, d\Gamma. 
\end{align} 
\tag{9} $$

1. If $\eta = \xi \in \mathcal{V}$ for a fixed vector $\xi$, then

$$ \frac{d}{dt}[\xi \cdot \mathbf{L}] = \xi \cdot \mathbf{R}. $$

2. If $\eta = \xi \times \varphi \in \mathcal{V}$ for a fixed vector $\xi$ and any fixed time $t$, then

$$ \frac{d}{dt}[\xi \cdot \mathbf{J}] = \xi \cdot \mathbf{T}. $$

3. If the boundary conditions and external loads are time-independent, then

$$ \frac{d}{dt} H = 0. $$

From the above result we see that the maximum number of integrals arises in the case of a free body, for which there are no essential boundary conditions and no external loads. In particular, we find that $H$ and each component of **L** and **J** are integrals or constants of motion in this case.

## 2.2. Energy-momentum conserving time discretization

We seek to approximate solutions $(\varphi, \mathbf{V}) : [0, T] \to \mathcal{U} \times \mathcal{V}$ of (1) at discrete times $t_n = n \Delta t$ where $n$ is a positive integer and $\Delta t > 0$. Denoting approximations to $\varphi(t_n)$ and $\mathbf{V}(t_n)$ by $\varphi_n$ and $\mathbf{V}_n$, the algorithmic problem is to construct approximations at time level $n+1$ given those at time level $n$.

Our goal is to develop time discretization schemes for which a discrete form of Proposition 1 can be established. We can accomplish this goal by replacing the exact derivative in (3) by an appropriate discrete derivative [19,20]. The final result is an implicit, one-step scheme of the form

$$ 
\left. 
\begin{aligned}
\int_{\Omega} \frac{1}{\Delta t} (\varphi_{n+1} - \varphi_n) \cdot \eta \, d\Omega &= \int_{\Omega} \mathbf{V}_{n+1/2} \cdot \eta \, d\Omega \\
\int_{\Omega} \frac{\rho_0}{\Delta t} (\mathbf{V}_{n+1} - \mathbf{V}_n) \cdot \eta \, d\Omega \\
&= - \int_{\Omega} \tilde{\mathbf{S}}(\varphi_n, \varphi_{n+1}) : \mathbf{F}(\varphi_{n+1/2})^{\mathrm{T}} \nabla \eta \, d\Omega \\
&\quad + \int_{\Omega} \mathbf{b}_{n+1/2} \cdot \eta \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_{n+1/2} \cdot \eta \, d\Gamma
\end{aligned}
\right\} \quad \forall \eta \in \mathcal{V}. 
\tag{10} $$

Here $(\cdot)_{n+1/2} := \frac{1}{2}[(\cdot)_n + (\cdot)_{n+1}]$, $\tilde{\mathbf{S}}(\varphi_n, \varphi_{n+1})$ is an algorithmic stress field defined by

$$ \tilde{\mathbf{S}}(\varphi_n, \varphi_{n+1}) := 2 \, dW(\mathbf{C}_n, \mathbf{C}_{n+1}), \tag{11} $$

and $dW(\mathbf{C}_n, \mathbf{C}_{n+1})$ is a discrete derivative defined by

$$ dW(\mathbf{C}_n, \mathbf{C}_{n+1}) := DW(\mathbf{C}_{n+1/2}) + \left[ \frac{W(\mathbf{C}_{n+1}) - W(\mathbf{C}_n) - DW(\mathbf{C}_{n+1/2}) : \mathbf{M}}{\|\mathbf{M}\|^2} \right] \mathbf{M}, \tag{12} $$

where $\mathbf{M} := \mathbf{C}_{n+1} - \mathbf{C}_n$ and $\|\mathbf{M}\|^2 = \mathbf{M} : \mathbf{M}$.
Just as in the time-continuous problem, the algorithmic equations (10) may possess various integrals or constants of motion related to **L**, **J** and **H** depending on the essential boundary conditions and the external loads. This fact is summarized in the following proposition.

**Proposition 2.** Let $(\varphi_n, \mathbf{V}_n) \in \mathcal{U} \times \mathcal{V}(n \ge 0)$ be a solution sequence of (10) for given boundary conditions, external loads **b** and **f**, and time step $\Delta t > 0$. Moreover, let $\mathbf{R}_{n+1/2}$ and $\mathbf{T}_{n+1/2}$ denote a resultant force and moment defined as

$$ 
\begin{aligned}
\mathbf{R}_{n+1/2} &= \int_{\Omega} \mathbf{b}_{n+1/2} \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_{n+1/2} \, d\Gamma, \\
\mathbf{T}_{n+1/2} &= \int_{\Omega} \varphi_{n+1/2} \times \mathbf{b}_{n+1/2} \, d\Omega + \int_{\Gamma_{\sigma}} \varphi_{n+1/2} \times \mathbf{f}_{n+1/2} \, d\Gamma.
\end{aligned}
\quad (13) $$

1. If $\eta = \xi \in \mathcal{V}$ for a fixed vector $\xi$, then

$$ \frac{1}{\Delta t} [\xi \cdot \mathbf{L}_{n+1} - \xi \cdot \mathbf{L}_n] = \xi \cdot \mathbf{R}_{n+1/2}. $$

2. If $\eta = \xi \times \varphi_n \in \mathcal{V}$ and $\eta = \xi \times \mathbf{V}_n \in \mathcal{V}$ for a fixed vector $\xi$ and any $n$, then

$$ \frac{1}{\Delta t} [\xi \cdot \mathbf{J}_{n+1} - \xi \cdot \mathbf{J}_n] = \xi \cdot \mathbf{T}_{n+1/2}. $$

3. If the boundary conditions and external loads are time-independent, then

$$ \frac{1}{\Delta t} [H_{n+1} - H_n] = 0. $$

**Proof.** The results follow from the specific form of the assumed variations and may be established as follows.

1. Choosing $\eta = \xi$ in $(10)_2$ gives

$$ \int_{\Omega} \frac{\rho_0}{\Delta t} (\mathbf{V}_{n+1} - \mathbf{V}_n) \cdot \xi \, d\Omega = \int_{\Omega} \mathbf{b}_{n+1/2} \cdot \xi \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_{n+1/2} \cdot \xi \, d\Gamma, \quad (14) $$

which is the required result for **L**.

2. To establish the result for **J**, we first note that, if $\xi \times \varphi_n \in \mathcal{V}$ for any $n$, then $\xi \times \varphi_{n+1/2} \in \mathcal{V}$ by linearity of $\mathcal{V}$, and similarly $\xi \times \mathbf{V}_{n+1/2} \in \mathcal{V}$. Next, we choose a variation $\eta = \xi \times \rho_0 \mathbf{V}_{n+1/2}$ in $(10)_1$ to obtain

$$ \int_{\Omega} \frac{1}{\Delta t} \xi \cdot [\varphi_{n+1} - \varphi_n] \times \rho_0 \mathbf{V}_{n+1/2} \, d\Omega = 0 \quad (15) $$

and choose a variation $\eta = \xi \times \varphi_{n+1/2}$ in $(10)_2$ to obtain

$$ \int_{\Omega} \frac{1}{\Delta t} \xi \cdot \varphi_{n+1/2} \times \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega = - \int_{\Omega} \tilde{\mathbf{S}} : \mathbf{F}_{n+1/2}^{\mathrm{T}} \hat{\xi} \mathbf{F}_{n+1/2} \, d\Omega + \xi \cdot \mathbf{T}_{n+1/2}, \quad (16) $$

where $\hat{\xi}$ is the skew tensor associated with $\xi$, i.e., $\xi \times \mathbf{a} = \hat{\xi}\mathbf{a}$ for any vector **a**. Since $\mathbf{F}_{n+1/2}^{\mathrm{T}} \hat{\xi} \mathbf{F}_{n+1/2}$ is a skew tensor and $\tilde{\mathbf{S}}$ is symmetric, the first term on the right-hand side of (16) vanishes and we are left with

$$ \int_{\Omega} \frac{1}{\Delta t} \xi \cdot \varphi_{n+1/2} \times \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega = \xi \cdot \mathbf{T}_{n+1/2}. \quad (17) $$

The desired result then follows by adding (15) and (17), and expanding the vector products.

3. To establish the result for **H**, we first note that $(\varphi_{n+1} - \varphi_n) \in \mathcal{V}$ for time-independent boundary conditions, and that $\mathbf{V}_n \in \mathcal{V}$ for any $n$ by assumption. Next, we use the variation $\eta = \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n]$ in $(10)_1$, and $\eta = [\varphi_{n+1} - \varphi_n]$ in $(10)_2$ to deduce that
$$
\begin{align}
& \int_{\Omega} \frac{1}{\Delta t} [\boldsymbol{\varphi}_{n+1} - \boldsymbol{\varphi}_n] \cdot \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega = \int_{\Omega} \mathbf{V}_{n+1/2} \cdot \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega, \nonumber \\
& \int_{\Omega} \frac{1}{\Delta t} \rho_0 [\mathbf{V}_{n+1} - \mathbf{V}_n] \cdot [\boldsymbol{\varphi}_{n+1} - \boldsymbol{\varphi}_n] \, d\Omega \nonumber \\
& \qquad = - \int_{\Omega} \tilde{\mathbf{S}} : \mathbf{F}_{n+1/2}^{\mathrm{T}} [\mathbf{F}_{n+1} - \mathbf{F}_n] \, d\Omega + \int_{\Omega} \mathbf{b} \cdot [\boldsymbol{\varphi}_{n+1} - \boldsymbol{\varphi}_n] \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \cdot [\boldsymbol{\varphi}_{n+1} - \boldsymbol{\varphi}_n] \, d\Gamma. \tag{18}
\end{align}
$$

Using the fact that $\tilde{\mathbf{S}}$ is symmetric and the definition of $\mathbf{C}$, we have

$$
\tilde{\mathbf{S}} : \mathbf{F}_{n+1/2}^{\mathrm{T}} [\mathbf{F}_{n+1} - \mathbf{F}_n] = \frac{1}{2} \tilde{\mathbf{S}} : [\mathbf{C}_{n+1} - \mathbf{C}_n] \quad (19)
$$

and by definition of $\tilde{\mathbf{S}}$,

$$
\frac{1}{2} \tilde{\mathbf{S}} : [\mathbf{C}_{n+1} - \mathbf{C}_n] = W(\mathbf{C}_{n+1}) - W(\mathbf{C}_n). \quad (20)
$$

Combining (20) with (19), and subtracting (18)$_2$ from (18)$_1$, leads to

$$
H(\boldsymbol{\varphi}_{n+1}, \mathbf{V}_{n+1}) = H(\boldsymbol{\varphi}_n, \mathbf{V}_n), \quad (21)
$$

which is the required result.

Comparing Proposition 2 with Proposition 1, we see that the algorithmic equations (10) inherit discrete
forms of the conservation laws associated with **L**, **J** and **H**. In particular, **H** and each of the components of
**L** and **J** are integrals of the time-discretized system (10) for the case of a free body.

The conservation properties of (10) are automatic in the sense that no extra parameters, equations or
projection steps are involved. The key feature is the use of the discrete derivative (12) in computing the
stress field $\tilde{\mathbf{S}}$. If the stored energy function $W(\mathbf{C})$ is sufficiently differentiable, then one can show that

$$
\mathrm{d}W(\mathbf{C}_n, \mathbf{C}_{n+1}) = DW(\mathbf{C}_{n+1/2}) + O(\|\mathbf{C}_{n+1} - \mathbf{C}_n\|^2). \quad (22)
$$

Hence, the expression for dW(C_n, C_{n+1}) given in (12) is well-defined in the limit ||C_{n+1} - C_n|| → 0. Moreover, for sufficiently regular solutions of the algorithmic equations we have

$$
\tilde{S}(\boldsymbol{\varphi}_n, \boldsymbol{\varphi}_{n+1}) = S(\boldsymbol{\varphi}_{n+1/2}) + O(\Delta t^2), \quad (23)
$$

and it follows that the conserving scheme outlined in (10) has the same formal order of accuracy as the
implicit mid-point rule.

An important feature of (10) is the use of an interpolated stretch field $C_{n+1/2}$, rather than the stretch field of an interpolated configuration $C(\varphi_{n+1/2})$, in evaluating the stress field. The use of $C(\varphi_{n+1/2})$ introduces an artificial coupling between overall rigid rotations and material stretches that is not present in the time-continuous problem. In [20] it was proved that such coupling can lead to numerical instabilities, while the use of interpolated strains like $C_{n+1/2}$ completely bypasses such problems. For further discussion and numerical examples of such phenomena see [42,43].

The scheme outlined in (10) is independent of any particular spatial discretization. In the stated weak form, the scheme is ideally suited for a finite element discretization in space. However, by employing a standard argument involving an integration by parts, one can obtain a strong form of (10), which could then be discretized in space by an appropriate finite difference method for example.

**3. Conserving algorithms for incompressible models**

3.1. *A quasi-incompressible formulation*

A conserving time discretization scheme can be constructed for the initial boundary-value problem (2) as
stated. However, this formulation may be problematic for subsequent spatial discretization. The basic
difficulty lies in the fact that any smooth solution of this system satisfies the pointwise condition $\det[\mathbf{F}] = 1$ in $\Omega$, and this class of deformations cannot be approximated well by standard, low-order finite element spaces. In particular, one can expect the phenomena of locking to appear in direct treatments of (2) (see, e.g., [25] for a summary account within a linearized context). Here we introduce a formulation of (2) along the lines of [38,40] to alleviate such difficulties. The main idea is to introduce a material field $\Theta(\mathbf{X}, t) > 0$ that approximates the Jacobian field $\zeta(\mathbf{X}, t) := \det[\mathbf{F}(\mathbf{X}, t)] > 0$, and then enforce the condition $\Theta = 1$ in $\Omega$. Since the incompressibility condition is enforced on $\Theta$ rather than $\zeta$, the formulation is called quasi-incompressible.

A quasi-incompressible formulation of (2) may be outlined as follows. To any quasi-Jacobian field $\Theta$ and deformation field $\varphi$ we associate a modified deformation gradient field $\bar{\mathbf{F}}$ by the expression

$$ \bar{\mathbf{F}} := \Theta^{1/3} \mathbf{F}^{\text{dev}}, \qquad (24) $$

where $\mathbf{F}^{\text{dev}}$ is the deviatoric part of the actual deformation gradient field $\mathbf{F}$, i.e.,

$$ \mathbf{F}^{\text{dev}} := \zeta^{-1/3} \mathbf{F}, \qquad (25) $$

and we define an associated modified Cauchy-Green stretch field $\bar{\mathbf{C}}$ by

$$ \bar{\mathbf{C}} := \bar{\mathbf{F}}^{\text{T}} \bar{\mathbf{F}}. \qquad (26) $$

Let $\mathbf{W}(\mathbf{C})$ be the original strain energy function introduced in (3) and define a modified strain energy $\bar{\mathbf{W}}(\mathbf{C}, \Theta)$ by

$$ \bar{\mathbf{W}}(\mathbf{C}, \Theta) := \mathbf{W}(\bar{\mathbf{C}}) = \mathbf{W}(\Theta^{2/3} \zeta^{-2/3} \mathbf{C}). \qquad (27) $$

Then a quasi-incompressible formulation of (2) is to find $(\varphi, \mathbf{V}, \lambda, \Theta, p) : [0, T] \to \mathcal{U} \times \mathcal{V} \times \mathcal{E} \times \mathcal{E} \times \mathcal{E}$ such that

$$ \left. \begin{aligned} \int_{\Omega} \dot{\varphi} \cdot \eta \, d\Omega &= \int_{\Omega} \mathbf{V} \cdot \eta \, d\Omega \\ \int_{\Omega} \rho_0 \dot{\mathbf{V}} \cdot \eta \, d\Omega + \int_{\Omega} [\mathbf{S}(\varphi, \Theta) + 2pDG(\mathbf{C}(\varphi))] : \mathbf{F}(\varphi)^{\text{T}} \nabla \eta \, d\Omega \\ &= \int_{\Omega} \mathbf{b} \cdot \eta \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \eta \, d\Gamma \\ 0 &= \int_{\Omega} \varphi [D_{\Theta} \bar{\mathbf{W}}(\mathbf{C}(\varphi), \Theta) + \lambda D \bar{\mathbf{G}}(\Theta) - p] \, d\Omega \\ 0 &= \int_{\Omega} \varphi [\zeta(\varphi) - \Theta] \, d\Omega \\ 0 &= \int_{\Omega} \varphi \bar{\mathbf{G}}(\Theta) \, d\Omega \end{aligned} \right\} \begin{cases} \forall \eta \in \mathcal{V} \\ \forall \varphi \in \mathcal{E} \end{cases} \qquad (28) $$

where

$$ \mathbf{S}(\varphi, \Theta) = 2D_c \bar{\mathbf{W}}(\mathbf{C}(\varphi), \Theta) \quad \text{and} \quad \bar{\mathbf{G}}(\Theta) = \Theta - 1. \qquad (29) $$

The relation between (28) and (2) is established in the following result.

**Proposition 3.** If $(\varphi, \mathbf{V}, \lambda, \Theta, p) : [0, T] \to \mathcal{U} \times \mathcal{V} \times \mathcal{E} \times \mathcal{E} \times \mathcal{E}$ is a smooth solution of the quasi-incompressible formulation (28), then $(\varphi, \mathbf{V}, \lambda) : [0, T] \to \mathcal{U} \times \mathcal{V} \times \mathcal{E}$ is a solution of the original incompressible formulation (2).

**Proof.** From (28)$_4$ we see that, if $\zeta$ and $\Theta$ are continuous, then $\zeta = \Theta$. This implies $\bar{\mathbf{F}} = \mathbf{F}$ and $\bar{\mathbf{W}}(\mathbf{C}, \Theta) = \mathbf{W}(\mathbf{C})$, so $\bar{\mathbf{W}}(\mathbf{C}, \Theta)$ is independent of $\Theta$. By continuity of $\lambda$ and $p$ we deduce from (28)$_3$ that $p = \lambda$, and so (28)$_2$ becomes

$$ \int_{\Omega} \rho_0 \dot{\mathbf{V}} \cdot \eta \, d\Omega + \int_{\Omega} [\mathbf{S}(\varphi) + 2\lambda DG(\mathbf{C}(\varphi))] : \mathbf{F}(\varphi)^{\text{T}} \nabla \eta \, d\Omega = \int_{\Omega} \mathbf{b} \cdot \eta \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \eta \, d\Gamma \quad \forall \eta \in \mathcal{V}, \qquad (30) $$

which is just (2)$_2$. Finally, since $\zeta = \Theta$ we obtain $\bar{\mathbf{G}}(\Theta) = \mathbf{G}(\mathbf{C})$, and so (28)$_5$ is equivalent to (2)$_3$, which establishes the result.
## 3.2. Conservation laws and integrals

Just as in the compressible case, we associate a total linear momentum $L(\varphi, V)$ and a total angular momentum $J(\varphi, V)$ to any state $(\varphi, V, \lambda, \Theta, p)$ of the quasi-incompressible system (28). We also associate to each state a total modified energy $\bar{H}(\varphi, V, \lambda, \Theta, p)$ defined as

$$ \bar{H}(\varphi, V, \lambda, \Theta, p) := \int_{\Omega} \frac{1}{2} \rho_0 |V|^2 + \bar{W}(C(\varphi), \Theta) \, d\Omega + \int_{\Omega} p[\zeta(\varphi) - \Theta] + \lambda \bar{G}(\Theta) \, d\Omega \\ - \int_{\Omega} \mathbf{b} \cdot \varphi \, d\Omega - \int_{\Gamma_{\sigma}} \mathbf{f} \cdot \varphi \, d\Gamma. \qquad (31) $$

As before, the weak initial boundary-value problem (28) may possess various integrals or constants of motion related to $L$, $J$ and $\bar{H}$ depending on the essential boundary conditions and external loads. This fact is summarized in the following proposition, whose proof is straightforward.

**Proposition 4.** Let $(\varphi, V, \lambda, \Theta, p)(t)$ be a solution of (28) for given boundary conditions specified in $\mathcal{U}$ and $\mathcal{V}$, and given external loads $\mathbf{b}$ and $\mathbf{f}$. Moreover, let $\mathbf{R}$ and $\mathbf{T}$ denote the resultant force and moment due to the external loads as defined in (9). Then $L$, $J$ and $\bar{H}$ satisfy conservation laws completely analogous to those in Proposition 1.

## 3.3. Energy-momentum conserving time discretization

Here we seek to approximate solutions $(\varphi, V, \lambda, \Theta, p) : [0, T] \to \mathcal{U} \times \mathcal{V} \times \mathcal{E} \times \mathcal{E} \times \mathcal{E}$ of (28) at discrete times $t_n = n\Delta t$ where $n$ is a positive integer and $\Delta t > 0$. Using the same notation as before, our goal is to develop time discretization schemes for which a discrete form of Proposition 4 can be established. As in the compressible case, we can accomplish this goal by replacing certain exact derivatives in (28) by appropriate discrete derivatives along the lines of [21]. The final result is an implicit, one-step scheme of the form

$$ \left. \begin{aligned} & \int_{\Omega} \frac{1}{\Delta t} (\varphi_{n+1} - \varphi_n) \cdot \eta \, d\Omega = \int_{\Omega} \mathbf{V}_{n+1/2} \cdot \eta \, d\Omega \\ & \int_{\Omega} \frac{\rho_0}{\Delta t} (\mathbf{V}_{n+1} - \mathbf{V}_n) \cdot \eta \, d\Omega = - \int_{\Omega} [\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG] : \mathbf{F}_{n+1/2}^{\mathrm{T}} \nabla \eta \, d\Omega \\ & \qquad + \int_{\Omega} \mathbf{b}_{n+1/2} \cdot \eta \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_{n+1/2} \cdot \eta \, d\Gamma \\ & 0 = \int_{\Omega} \phi \left[ \frac{1}{2} (d\bar{W}_{C_n} + d\bar{W}_{C_{n+1}}) + \lambda_{n+1/2} \, d\bar{G} - p_{n+1/2} \right] \, d\Omega \\ & 0 = \int_{\Omega} \phi [\zeta_{n+1} - \Theta_{n+1}] \, d\Omega \\ & 0 = \int_{\Omega} \phi \bar{G}_{n+1} \, d\Omega \end{aligned} \right\} \begin{cases} \forall \eta \in \mathcal{V} \\ \forall \phi \in \mathcal{E} \end{cases} \qquad (32) $$

where

$$ \tilde{\mathbf{S}} = d\bar{W}_{\Theta_n}(\mathbf{C}_n, \mathbf{C}_{n+1}) + d\bar{W}_{\Theta_{n+1}}(\mathbf{C}_n, \mathbf{C}_{n+1}), $$

$$ d\bar{W}_{\Theta}(\mathbf{C}_n, \mathbf{C}_{n+1}) = D_c \bar{W}(\mathbf{C}_{n+1/2}, \Theta) + \left[ \frac{\bar{W}(\mathbf{C}_{n+1}, \Theta) - \bar{W}(\mathbf{C}_n, \Theta) - D_c \bar{W}(\mathbf{C}_{n+1/2}, \Theta) : \mathbf{M}}{\|\mathbf{M}\|^2} \right] \mathbf{M}, $$

$$ d\bar{W}_C(\Theta_n, \Theta_{n+1}) = D_{\Theta} \bar{W}(\mathbf{C}, \Theta_{n+1/2}) + \left[ \frac{\bar{W}(\mathbf{C}, \Theta_{n+1}) - \bar{W}(\mathbf{C}, \Theta_n) - D_{\Theta} \bar{W}(\mathbf{C}, \Theta_{n+1/2}) \cdot \Delta\Theta}{(\Delta\Theta)^2} \right] \Delta\Theta, \quad (33) $$

$$ dG(\mathbf{C}_n, \mathbf{C}_{n+1}) = DG(\mathbf{C}_{n+1/2}) + \left[ \frac{G(\mathbf{C}_{n+1}) - G(\mathbf{C}_n) - DG(\mathbf{C}_{n+1/2}) : \mathbf{M}}{\|\mathbf{M}\|^2} \right] \mathbf{M}, $$

$$ d\bar{G}(\Theta_n, \Theta_{n+1}) = D\bar{G}(\Theta_{n+1/2}) + \left[ \frac{\bar{G}(\Theta_{n+1}) - \bar{G}(\Theta_n) - D\bar{G}(\Theta_{n+1/2}) \cdot \Delta\Theta}{(\Delta\Theta)^2} \right] \Delta\Theta \equiv 1, $$

$$ \mathbf{M} = \mathbf{C}_{n+1} - \mathbf{C}_n \quad \text{and} \quad \Delta\Theta = \Theta_{n+1} - \Theta_n. $$
Just as in the time-continuous problem, the algorithmic equations (32) may possess various integrals or constants of motion related to **L**, **J** and **H** depending on the essential boundary conditions and the external loads. This fact is summarized in the following proposition.

**Proposition 5.** Let $(\varphi_n, \mathbf{V}_n, \lambda_n, \Theta_n, p_n) \in \mathcal{U} \times \mathcal{V} \times \mathcal{E} \times \mathcal{E} \times \mathcal{E} \ (n \ge 0)$ be a solution sequence of (32) for given boundary conditions, external loads **b** and **f**, and time step $\Delta t > 0$. Moreover, let $\mathbf{R}_{n+1/2}$ and $\mathbf{T}_{n+1/2}$ denote a resultant force and moment as defined in (13). Then **L**, **J** and $\bar{\mathbf{H}}$ satisfy discrete conservation laws completely analogous to those in Proposition 2.

**Proof.** The results for **L** and **J** follow as for the compressible case upon noting that the total algorithmic stress field $\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG$ is symmetric. To establish the result for $\bar{\mathbf{H}}$, choose a variation $\eta = \rho_0[\mathbf{V}_{n+1} - \mathbf{V}_n]$ in (32)$_1$ to obtain

$$ \int_{\Omega} \frac{\rho_0}{\Delta t} [\varphi_{n+1} - \varphi_n] \cdot [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega = \int_{\Omega} \rho_0 \mathbf{V}_{n+1/2} \cdot [\mathbf{V}_{n+1} - \mathbf{V}_n] \, d\Omega, \quad (34) $$

a variation $\eta = [\varphi_{n+1} - \varphi_n]$ in (32)$_2$ to obtain

$$ \int_{\Omega} \frac{\rho_0}{\Delta t} [\mathbf{V}_{n+1} - \mathbf{V}_n] \cdot [\varphi_{n+1} - \varphi_n] \, d\Omega = - \int_{\Omega} [\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG] : \mathbf{F}_{n+1/2}^{\mathrm{T}} [\mathbf{F}_{n+1} - \mathbf{F}_n] \, d\Omega + \int_{\Omega} \mathbf{b} \cdot [\varphi_{n+1} - \varphi_n] \, d\Omega + \int_{\Gamma_a} \mathbf{f} \cdot [\varphi_{n+1} - \varphi_n] \, d\Gamma, \quad (35) $$

a variation $\phi = \Theta_{n+1} - \Theta_n$ in (32)$_3$, $\phi = p_{n+1} - p_n$ in (32)$_4$ and $\phi = \lambda_{n+1} - \lambda_n$ in (32)$_5$ to obtain

$$ \left. \begin{aligned} 0 &= \int_{\Omega} \left[ \frac{1}{2}(d\bar{\mathbf{W}}_{\mathbf{C}_n} + d\bar{\mathbf{W}}_{\mathbf{C}_{n+1}}) + \lambda_{n+1/2} \, d\bar{\mathbf{G}} - p_{n+1/2} \right] [\Theta_{n+1} - \Theta_n] \, d\Omega, \\ 0 &= \int_{\Omega} [\zeta_{n+1} - \Theta_{n+1}] [p_{n+1} - p_n] \, d\Omega, \\ 0 &= \int_{\Omega} \bar{\mathbf{G}}_{n+1} [\lambda_{n+1} - \lambda_n] \, d\Omega. \end{aligned} \right\} \quad (36) $$

Using the fact that $\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG$ is symmetric and the definition of **C**, we have

$$ [\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG] : \mathbf{F}_{n+1/2}^{\mathrm{T}} [\mathbf{F}_{n+1} - \mathbf{F}_n] = \frac{1}{2} [\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG] : [\mathbf{C}_{n+1} - \mathbf{C}_n], \quad (37) $$

and by definition of $\tilde{\mathbf{S}}$ and $dG$,

$$ \frac{1}{2} [\tilde{\mathbf{S}} + 2p_{n+1/2} \, dG] : [\mathbf{C}_{n+1} - \mathbf{C}_n] = \frac{1}{2} [\bar{\mathbf{W}}(\mathbf{C}_{n+1}, \Theta_n) - \bar{\mathbf{W}}(\mathbf{C}_n, \Theta_n)] + \frac{1}{2} [\bar{\mathbf{W}}(\mathbf{C}_{n+1}, \Theta_{n+1}) - \bar{\mathbf{W}}(\mathbf{C}_n, \Theta_{n+1})] \\ + p_{n+1/2}[G(\mathbf{C}_{n+1}) - G(\mathbf{C}_n)]. \quad (38) $$

From (36)$_1$ and the definitions of $d\bar{\mathbf{W}}_C$ and $d\bar{\mathbf{G}}$ we obtain

$$ \left[ \frac{1}{2} (d\bar{\mathbf{W}}_{\mathbf{C}_n} + d\bar{\mathbf{W}}_{\mathbf{C}_{n+1}}) + \lambda_{n+1/2} \, d\bar{\mathbf{G}} - p_{n+1/2} \right] [\Theta_{n+1} - \Theta_n] \\ = \frac{1}{2} [\bar{\mathbf{W}}(\mathbf{C}_n, \Theta_{n+1}) - \bar{\mathbf{W}}(\mathbf{C}_n, \Theta_n)] + \frac{1}{2} [\bar{\mathbf{W}}(\mathbf{C}_{n+1}, \Theta_{n+1}) - \bar{\mathbf{W}}(\mathbf{C}_{n+1}, \Theta_n)] + \lambda_{n+1/2} [\Theta_{n+1} - \Theta_n] \\ - p_{n+1/2} [\Theta_{n+1} - \Theta_n]. \quad (39) $$

Subtracting (35) from (34), and adding (36) then leads to the result

$$ \begin{aligned} 0 = \bar{\mathbf{H}}_{n+1} - \bar{\mathbf{H}}_n + \frac{1}{2} \int_{\Omega} [\zeta_{n+1} - \Theta_{n+1}] [p_{n+1} - p_n] \, d\Omega - \frac{1}{2} \int_{\Omega} [\zeta_n - \Theta_n] [p_{n+1} - p_n] \, d\Omega \\ + \frac{1}{2} \int_{\Omega} \bar{\mathbf{G}}_{n+1} [\lambda_{n+1} - \lambda_n] \, d\Omega - \frac{1}{2} \int_{\Omega} \bar{\mathbf{G}}_n [\lambda_{n+1} - \lambda_n] \, d\Omega. \end{aligned} \quad (40) $$

Since $(p_{n+1} - p_n) \in \mathcal{E}$ and $(\lambda_{n+1} - \lambda_n) \in \mathcal{E}$ for any $n$, we obtain $\bar{\mathbf{H}}_{n+1} - \bar{\mathbf{H}}_n = 0$, which is the desired result.
The above result shows that the algorithmic equations (32) inherit discrete forms of the conservation laws associated with L, J and $\tilde{H}$. In particular, $\tilde{H}$ and each of the components of L and J are integrals for the case of a free body. As for the compressible case, an important feature of (32) is the use of an interpolated stretch field $C_{n+1/2}$, rather than the stretch field of an interpolated configuration $C(\varphi_{n+1/2})$, in evaluating the stress field.

The original incompressible formulation (2) and the quasi-incompressible formulation (28) may be interpreted as differential-algebraic evolution equations [3,34] in an appropriate function space. Systems of this form are typically not well-posed for arbitrary initial data. In particular, one expects well-posedness only for data satisfying certain compatibility conditions implied by the constraints [26,34]. For the case of (2), the constraint is $G(C(\varphi)) = 0$ in $\Omega$, and a basic compatibility condition on $(\varphi, V)$ is that $(\varphi, V) \in \mathcal{M}$. A condition on $\lambda$ may be deduced by considering the second time-derivative of the constraint $G(C(\varphi)) = 0$. Finally, since smooth solutions of (28) satisfy (2) with $\Theta = \zeta(\varphi)$ and $p = \lambda$, we see that initial conditions on $\Theta$ and $p$ follow from those on $\varphi$ and $\lambda$.

# 4. Finite element spatial discretization

Here we discuss the finite element spatial discretization of the time integration schemes (10) and (32). Since the finite element discretization of the compressible case (10) is straightforward, we concentrate on the quasi-incompressible case (32). In particular, we introduce a general mixed finite element discretization, reduce it to a two-field formulation and discuss an efficient solution strategy. This strategy can be interpreted as the generalization to elastodynamics of the augmented multiplier methods used in [40] for elastostatics. For more details on these methods see [2,17,18].

## 4.1. Mixed finite element formulation

Consider a partition of the domain $\Omega \subset \mathbb{R}^3$ into nonoverlapping subdomains $\Omega^e$ ($e = 1, \dots, m_{\text{elem}}$), each defined by a set of element nodes $\mathbf{X}_a^e \in \bar{\Omega}^e$ ($a = 1, \dots, m_{\text{elem}}$), and for each $e = 1, \dots, m_{\text{elem}}$ consider the sets

$$ v^e = \{a | a = 1, \dots, m_{\text{elem}}\} \quad \text{and} \quad v_g^e = \{a \in v^e | \mathbf{X}_a^e \in \Gamma_d\}. \qquad (41) $$

In accordance with a standard isoparametric discretization, we parametrize each subdomain $\Omega^e$ with a mapping $\chi^e: \square \to \Omega^e$ of the form

$$ \chi^e(\xi) = \sum_{a=1}^{m_{\text{elem}}} N^a(\xi) \mathbf{X}_a^e, $$

where $\square \subset \mathbb{R}^3$ is typically the bi-unit cube and $N^a: \square \to \mathbb{R}$ ($a = 1, \dots, m_{\text{elem}}$) are standard isoparametric shape functions. In particular, the shape functions satisfy the condition $N^a(\xi_b) = \delta_b^a$, where $\xi_b \in \bar{\Omega}$ ($b = 1, \dots, m_{\text{elem}}$) are element nodes for the parent domain $\square$.

Given a discretization of $\Omega$ as described above, we introduce an approximation $\mathcal{U}^h$ of $\mathcal{U}$, and an approximation $\mathcal{V}^h$ of $\mathcal{V}$, defined as

$$ \begin{aligned} \mathcal{U}^h &= \left\{ \varphi^h : \Omega \to \mathbb{R}^3 \middle| \varphi^h \in C^0(\Omega, \mathbb{R}^3), \varphi^h|_{\Omega^e} \circ \chi^e = \sum_{a \in v^e \setminus v_g^e} N^a \mathbf{d}_a^e + \sum_{a \in v_g^e} N^a \mathbf{g}(\mathbf{X}_a^e), \mathbf{d}_a^e \in \mathbb{R}^3 \right\}, \\ \mathcal{V}^h &= \left\{ \eta^h : \Omega \to \mathbb{R}^3 \middle| \eta^h \in C^0(\Omega, \mathbb{R}^3), \eta^h|_{\Omega^e} \circ \chi^e = \sum_{a \in v^e \setminus v_g^e} N^a \mathbf{c}_a^e, \mathbf{c}_a^e \in \mathbb{R}^3 \right\}. \end{aligned} \qquad (42) $$

where $\mathbf{g}: \Gamma_d \to \mathbb{R}^3$ is a specified deformation on $\Gamma_d$. Similarly, we introduce an approximation $\mathcal{E}^h$ of $\mathcal{E}$, namely,
$$ \mathcal{E}^h = \left\{ \phi^h : \Omega \to \mathbb{R} \mid \phi^h|_{\Omega^e} \circ \chi^e = \sum_{i=1}^{m_{\text{dis}}} \hat{N}^i \alpha_i^e, \alpha_i^e \in \mathbb{R} \right\}, \quad (43) $$

where $\hat{N}^i : \square \to \mathbb{R}$ ($i = 1, \dots, m_{\text{dis}}$) are smooth interpolation functions. Note that the functions in $\mathcal{E}^h$ are smooth within element domains, but are allowed to be discontinuous across element boundaries.

A mixed finite element discretization of the time integration scheme (32) is obtained by restriction to the finite-dimensional spaces $\mathcal{U}^h$, $\mathcal{V}^h$ and $\mathcal{E}^h$. Standard arguments [25] can then be used to reduce this weak formulation to a system of nonlinear algebraic equations for the nodal variables. Finally, we note that the finite element formulation of (32) possesses conservation laws analogous to those of (32) itself.

## 4.2. Reduction to a two-field formulation

Here we eliminate the fields $V_{n+1}^h$, $\Theta_{n+1}^h$, and $p_{n+1}^h$ from the finite element formulation of (32). The result will be a two-field formulation in which we need only solve for $\varphi_{n+1}^h$ and $\lambda_{n+1}^h$ at each step of the algorithm.

The elimination of $V_{n+1}^h$ is straightforward in view of (32)$_1$. In particular, using this equation together with the definition of $V_{n+1/2}^h$ leads to

$$ V_{n+1}^h = \frac{2}{\Delta t} (\varphi_{n+1}^h - \varphi_n^h) - V_n^h. \quad (44) $$

The elimination of $\Theta_{n+1}^h$ makes crucial use of the structure of $\mathcal{E}^h$. Assuming the interpolation functions $\hat{N}^i : \square \to \mathbb{R}$ ($i = 1, \dots, m_{\text{dis}}$) are nontrivial, we introduce for each $e = 1, \dots, m_{\text{elem}}$ a symmetric, positive-definite matrix $\Xi_e \in \mathbb{R}^{m_{\text{dis}} \times m_{\text{dis}}}$ with components

$$ (\Xi_e)^{ij} = \int_{\square} \hat{N}^i \hat{N}^j \det[\nabla \chi^e] d\square. \quad (45) $$

We denote the components of the inverse matrix by $(\Xi_e^{-1})_{ij}$. Since $\mathcal{E}^h$ is spanned by functions allowed to be discontinuous across element boundaries, we can use standard arguments and reduce the finite element version of (32)$_4$ to

$$ \int_{\Omega^e} \Theta_{n+1}^h |_{\Omega^e} \phi^h|_{\Omega^e} d\Omega^e = \int_{\Omega^e} \zeta_{n+1}^h |_{\Omega^e} \phi^h|_{\Omega^e} d\Omega^e \quad \forall \phi^h|_{\Omega_e} (e = 1, \dots, m_{\text{elem}}), \quad (46) $$

where $\zeta_{n+1}^h = \det[F_{n+1}^h]$. For any $e = 1, \dots, m_{\text{elem}}$ we then arrive at the equation

$$ \int_{\square} (\Theta_{n+1}^h|_{\Omega^e} \circ \chi^e) \hat{N}^i \det[\nabla \chi^e] d\square = \int_{\square} (\zeta_{n+1}^h|_{\Omega^e} \circ \chi^e) \hat{N}^i \det[\nabla \chi^e] d\square, \quad (47) $$

which holds for $i = 1, \dots, m_{\text{dis}}$. To solve this equation for $\Theta_{n+1}^h|_{\Omega^e}$ we consider an interpolation of the form

$$ \Theta_{n+1}^h|_{\Omega^e} \circ \chi^e = \sum_{j=1}^{m_{\text{dis}}} \hat{N}^j \Theta_{j,n+1}^e \quad (48) $$

which, in view of (47), leads to the relation

$$ \sum_{j=1}^{m_{\text{dis}}} (\Xi_e)^{ij} \Theta_{j,n+1}^e = \int_{\square} (\zeta_{n+1}^h|_{\Omega^e} \circ \chi^e) \hat{N}^i \det[\nabla \chi^e] d\square \quad (i = 1, \dots, m_{\text{dis}}). \quad (49) $$

Solving this equation for $\Theta_{j,n+1}^e$ ($j = 1, \dots, m_{\text{dis}}$) and substituting the result into (48) then yields

$$ \Theta_{n+1}^h|_{\Omega^e} \circ \chi^e = \sum_{i,j=1}^{m_{\text{dis}}} \hat{N}^j (\Xi_e^{-1})_{ji} \int_{\square} (\zeta_{n+1}^h|_{\Omega^e} \circ \chi^e) \hat{N}^i \det[\nabla \chi^e] d\square, \quad (50) $$

which expresses $\Theta_{n+1}^h|_{\Omega^e}$ as a function of $\varphi_{n+1}^h|_{\Omega^e}$.
The elimination of $p_{n+1}^h$ proceeds along similar lines as the elimination of $\Theta_{n+1}^h$. Working with the finite element version of (32)$_3$, we find that the restriction of $p_{n+1}^h$ to any element $\Omega^e$ is given by

$$ p_{n+1}^h|_{\Omega^e} \circ \chi^e = -p_n^h|_{\Omega^e} \circ \chi^e + 2 \sum_{i,j=1}^{m_{dis}} \left\{ \hat{N}^j(\Xi_e^{-1})_{ji} \int_{\square} \left[ \frac{1}{2} (d\bar{W}_{C_n^h} + d\bar{W}_{C_{n+1}^h}) + \lambda_{n+1/2}^h d\bar{G} \right] \Big|_{\Omega^e} \circ \chi^e \hat{N}^i \det[\nabla\chi^e] d\square \right\}. \quad (51) $$

Our two-field, mixed finite element discretization of (32) can be stated as follows. Given $(\varphi_n^h, \mathbf{V}_n^h, \lambda_n^h, \Theta_n^h, p_n^h)$ find $(\varphi_{n+1}^h, \lambda_{n+1}^h) \in \mathcal{W}^h \times \mathcal{E}^h$ such that

$$ \left. \begin{aligned} & \int_{\Omega} \frac{2p_0}{\Delta t^2} (\varphi_{n+1}^h - \varphi_n^h) \cdot \eta^h \, d\Omega \\ &= \int_{\Omega} \frac{2p_0}{\Delta t} \mathbf{V}_n^h \cdot \eta^h \, d\Omega - \int_{\Omega} [\tilde{\mathbf{S}}^h + 2p_{n+1/2}^h \, d\mathbf{G}^h] : (\mathbf{F}_{n+1/2}^h)^T \nabla \eta^h \, d\Omega \\ & \quad + \int_{\Omega} \mathbf{b}_{n+1/2}^h \cdot \eta^h \, d\Omega + \int_{\Gamma_{\sigma}} \mathbf{f}_{n+1/2}^h \cdot \eta^h \, d\Gamma \\ & 0 = \int_{\Omega} \phi^h \bar{G}_{n+1}^h \, d\Omega \end{aligned} \right\} \begin{cases} \forall \eta^h \in \mathcal{V}^{-h} \\ \forall \phi^h \in \mathcal{E}^h \end{cases} \quad (52) $$

where $\Theta_{n+1}^h$ and $p_{n+1}^h$ are considered as functions of $\varphi_{n+1}^h$ and $\lambda_{n+1}^h$ via expressions (50) and (51), respectively, and $\mathbf{V}_{n+1}^h$ is computed according to (44).

The above algorithm is completed by an appropriate choice of interpolation functions for $\mathcal{W}^h$, $\mathcal{V}^{-h}$ and $\mathcal{E}^h$. Ideally, one would choose functions for which convergence of solutions of (52) to those of (32) could be established. This issue, which is generally a delicate matter [4,17], will not be addressed here. For purposes of illustration, we note that piecewise linear interpolation functions could serve as a basis for $\mathcal{W}^h$ and $\mathcal{V}^{-h}$. Moreover, the constant interpolation function $\hat{N}(\xi) \equiv 1$ ($m_{dis} = 1$) could serve as a basis for $\mathcal{E}^h$. In this case we would have

$$ \Theta_{n+1}^h|_{\Omega^e} = \frac{1}{\mathrm{vol}(\Omega^e)} \int_{\Omega^e} \zeta(\varphi_{n+1}^h) \, d\Omega^e, \quad (53) $$

and the restriction of $\Theta_{n+1}^h$ to any element would be constant and equal to the average Jacobian value over $\Omega^e$. Thus, the incompressibility condition $\bar{G}(\Theta_{n+1}^h) = \Theta_{n+1}^h - 1 = 0$ in $\Omega^e$ ($e = 1, \dots, m_{elem}$) would be a condition on the average Jacobian for each element.

### 4.3. Solution strategy

Following [38] we consider frame-invariant stored energy functions $W(C)$ of the separable form

$$ W(C) = W^{\mathrm{dev}}(C^{\mathrm{dev}}) + U(\det[C]^{1/2}), \quad (54) $$

where $W^{\mathrm{dev}}(C^{\mathrm{dev}})$ is a deviatoric stored energy function and $U(\det[C]^{1/2})$ is a dilatational stored energy function which is convex in its real argument. Here $C^{\mathrm{dev}}$ is the deviatoric part of $C$ defined by $C^{\mathrm{dev}} = \det[C]^{-1/3}C$, in particular $\det[C^{\mathrm{dev}}] = 1$. In view of (54), the modified stored energy function $\bar{W}(C, \Theta)$ defined in (27) becomes

$$ \bar{W}(C, \Theta) = W^{\mathrm{dev}}(C^{\mathrm{dev}}) + U(\Theta). \quad (55) $$

For our quasi-incompressible formulation we can assume without loss of generality that the function $U$ is of the form

$$ U(\Theta) = K\gamma(\Theta). \quad (56) $$

Here $K > 0$ is a constant and $\gamma(\Theta)$ is a convex, nonnegative function that satisfies $\gamma(\Theta) = 0$ if and only if $\Theta = 1$, e.g.,
$$ \gamma(\Theta) = \frac{1}{2}(\Theta^2 - 1) - \ln(\Theta). \qquad (57) $$

This assumption can be made since $\Theta$ is constrained by the condition $\Theta = 1$. Thus, we interpret $U(\Theta)$ as a penalty function in the quasi-incompressible formulation.

With the separable stored energy function (55), we may interpret (52) as an augmented Lagrangian formulation in the sense of optimization theory [2,29]. In particular, the constraint (52)$_2$ is enforced by means of a penalty parameter $K$ and a multiplier $\lambda_{n+1}^h$. There are various efficient solution strategies for such systems [2,14,29]. One example is the classic Uzawa algorithm, which may be summarized as follows using $k$ as an iteration index:

1. Let $(\phi_n^h, \mathbf{V}_n^h, \lambda_n^h, \Theta_n^h, p_n^h)$ be given and set $\lambda_{n+1}^{h,k}|_{k=0} = \lambda_n^h$.

2. For fixed $\lambda_{n+1}^{h,k}$, solve (52)$_1$ for $\phi_{n+1}^h$ and call this solution $\phi_{n+1}^{h,k}$.

3. If (52)$_2$ is not satisfied within a specified tolerance, then update the multiplier field using the rule

$$ \int_{\Omega} \lambda_{n+1}^{h,k+1} \phi^h d\Omega = \int_{\Omega} \lambda_{n+1}^{h,k} \phi^h d\Omega + \int_{\Omega} K \bar{G}(\Theta_{n+1}^{h,k}) \phi^h d\Omega \quad \forall \phi^h \in \mathcal{E}^h, \qquad (58) $$

and repeat Step 2.

Note that the above algorithm effectively reduces the two-field formulation to one field. Hence, Step 2 can be carried out using standard finite element codes suitable for nonlinear problems.

## 5. Numerical examples

In this section we illustrate the performance of the time integration schemes (10) and (32). In the compressible case (10) we employ a straightforward Galerkin finite element discretization, and in the quasi-incompressible case (32) we employ the mixed formulation described above.

The example problems we consider are the tumbling block, the stretched plate, and the inverted spherical cap. In all cases, we employ hyperelastic material models of Ogden type [32,33], namely,

$$ W(\mathbf{C}) = \sum_{A=1}^{3} w(\lambda_{A}), \qquad (59) $$

where

$$ w(\lambda_A) = \sum_{m=1}^{k} \left[ \frac{\mu_m}{\alpha_m} (\lambda_A^{\alpha_m} - 1) - \mu_m \ln(\lambda_A) \right]. \qquad (60) $$

Here $\mathbf{C}$ is the right Cauchy-Green stretch tensor, $\lambda_A^2$ ($A = 1, 2, 3$) are the eigenvalues of $\mathbf{C}$, $\mu_m$ and $\alpha_m$ ($m = 1, \dots, k$) are parameters, and $k$ is the number of terms in the Ogden model. Models of this form are particularly useful for rubber-like materials [32,33,47]. For our compressible calculations we use $k = 3$ and the parameters,

$$ \left. \begin{array}{l} \mu_1 = 0.690 \times 10^6 \text{ N/m}^2, \quad \alpha_1 = 1.3, \\ \mu_2 = 0.010 \times 10^6 \text{ N/m}^2, \quad \alpha_2 = 4.0, \\ \mu_3 = -0.012 \times 10^6 \text{ N/m}^2, \quad \alpha_3 = -2.0, \end{array} \right\} \qquad (61) $$

which are taken from [33, p. 498].

An extension of the above model to incompressible materials is given by

$$ W(\mathbf{C}) = W^{\text{dev}}(\mathbf{C}^{\text{dev}}) = \sum_{A=1}^{3} w^{\text{dev}}(\lambda_{A}^{\text{dev}}) \qquad (62) $$
together with the constraint $\zeta = \lambda_1\lambda_2\lambda_3 = 1$, where $w^{\text{dev}}$ is of the form

$$w^{\text{dev}}(\lambda_A^{\text{dev}}) = \sum_{m=1}^{k} \frac{\mu_m}{\alpha_m} \left[ (\lambda_A^{\text{dev}})^{\alpha_m} - 1 \right]. \qquad (63)$$

Here $C^{\text{dev}}$ is the deviatoric part of $C$ and $(\lambda_A^{\text{dev}})^2$ ($A = 1, 2, 3$) are the eigenvalues of $C^{\text{dev}}$. For this model the modified stored energy function $\bar{W}(C, \Theta)$ becomes

$$\bar{W}(C, \Theta) = W^{\text{dev}}(C^{\text{dev}}) + U(\Theta), \qquad (64)$$

where $U$ is interpreted as a penalty function, which we take as

$$U(\Theta) = K\gamma(\Theta). \qquad (65)$$

Here $K > 0$ is a penalty parameter and $\gamma$ is the function given in (57). For our quasi-incompressible calculations we use the same constants as in (61).

For both the compressible and quasi-incompressible cases we take $m_{\text{en}} = 8$ and use standard trilinear interpolation functions as the basis for the spaces $\mathcal{W}^h$ and $\mathcal{V}^{-h}$. For the quasi-incompressible case, we take $m_{\text{dis}} = 1$ and use a constant interpolation function as the basis for the space $\mathcal{E}^h$. Thus, the functions in $\mathcal{W}^h$ and $\mathcal{V}^{-h}$ are piecewise linear and continuous, while those in $\mathcal{E}^h$ are piecewise constant and possibly discontinuous across element boundaries.

### 5.1. Tumbling elastic block

Consider a cube of dimension $l = 0.02$ m as shown in Fig. 1. The cube is composed of a homogeneous, elastic material of Ogden type with parameters as given in (61) and density $\rho_0 = 1000$ kg/m³. We suppose the cube is initially at rest and subject to tractions $\mathbf{f}_1$ and $\mathbf{f}_2$ that vanish after a short period of time. The traction $\mathbf{f}_1$ is a uniform force distribution on the face $X_3 = -l/2$ defined by $\mathbf{f}_1 = p(t)(0, -f/4, -f/8)$ and the traction $\mathbf{f}_2$ is a uniform force distribution on the face $X_3 = l/2$ defined by $\mathbf{f}_2 = p(t)(0, f, f/2)$. Here $f = 32$ N/cm² and $p(t)$ is given by

$$p(t) = \begin{cases} t, & 0 \le t \le 0.005 \text{ s}, \\ 0, & t > 0.005 \text{ s}. \end{cases} \qquad (66)$$

We compute solutions for both the compressible and quasi-incompressible formulations using a uniform spatial discretization of the cube into 27 elements. For the compressible formulation we consider two time integration schemes: the implicit mid-point rule and the energy-momentum preserving scheme (10). For the quasi-incompressible formulation we use the conserving scheme (32) in the reduced form (52). We employ the Uzawa iteration strategy outlined above and iterate until $\max_e |\bar{G}(\Theta_n^e)| < 5 \times 10^{-10}$ for each time step.

Fig. 2 shows snapshots, at time increments of 0.005 s, of the cube motion in the time interval $[0, 0.1]$ s. The data in the figure corresponds to the compressible motion, however, the quasi-incompressible motion is visually identical. For clarity, we show the snapshots in the $X_2-X_3$ plane only.

Figs. 3 and 4 show total energy and angular momentum time histories, computed using two different time steps, for the mid-point rule and conserving scheme applied to the compressible formulation. Note

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic of elastic cube showing applied loads.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Cube motion (left-to-right) in $X_2-X_3$ plane at time increments of 0.005 s.</figcaption></figure>

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

<figure><img src="image_5.png" /><figcaption>Fig. 3. Energy and angular momentum of compressible cube for mid-point rule.</figcaption></figure>

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

<figure><img src="image_7.png" /><figcaption>Fig. 4. Energy and angular momentum of compressible cube for conserving scheme.</figcaption></figure>

that the mid-point rule experiences numerical blow-up for the time step $\Delta t = 0.005$ s in the time interval of the computation.

Fig. 5 shows an accuracy comparison between the conserving scheme and the mid-point rule for the compressible formulation. For this comparison we considered the above motion in the time interval $[0.0100, 0.0125]$ s. For each scheme, we supplied the same initial conditions at $t_0 = 0.0100$ s and computed the solution up to a time $t_1 = 0.0125$ s. Denoting the solution at time $t_1$ for a given scheme and time step by $\varphi_{t_1,\Delta t}$, we define the $L_2$-error in displacements as

$$L_2\text{-error} = \left( \int_{\Omega} \left| \varphi_{t_1, \Delta t} - \varphi_{t_1, \text{conv}} \right|^2 \right)^{1/2}, \qquad (67)$$
<figure><img src="image_8.png" /><figcaption>Fig. 5. L₂-error in displacements of compressible cube versus time step.</figcaption></figure>

where $\varphi_{t_1, \text{conv}}$ is the solution at $t_1$ computed using a time step of $\Delta t = 2.5 \times 10^{-7}$ s. As shown in Fig. 5, the conserving scheme and the mid-point rule have nearly identical accuracy properties for this model problem.

Fig. 6 shows total energy and angular momentum time histories, computed using a time step $\Delta t = 0.005$ s, for the conserving scheme applied to the quasi-incompressible formulation.

## 5.2. Stretched elastic plate

Consider a thick rectangular plate of thickness $w = 0.002$ m, length $l = 0.01$ m and height $h = 0.01$ m as shown in Fig. 7. The plate is composed of a homogeneous, incompressible, elastic material of Ogden type with parameters as given in (61) and density $\rho_0 = 1000$ kg/m³. In this example, we suppose the plate is subject to a zero displacement boundary condition along the shaded region shown in the figure, is initially at rest and is subject to a traction $\mathbf{f}$ that vanishes after a short period of time. The traction $\mathbf{f}$ is applied to the face $X_3 = -h/2$, is uniform across the thickness, and varies in the $X_2$-direction according to

$$ \mathbf{f}(X_2, t) = \begin{cases} p(t)(f/2, 0, -f), & -\frac{l}{2} \le X_2 < 0, \\ p(t)(0, 0, -f), & 0 \le X_2 \le \frac{l}{2}, \end{cases} \qquad (68) $$

where $f = 5$ kN/cm² and $p(t)$ is given by

$$ p(t) = \begin{cases} t, & 0 \le t \le 0.004 \text{ s}, \\ 0, & t > 0.004 \text{ s}. \end{cases} \qquad (69) $$

The zero displacement boundary condition is imposed along the surface defined by $\{X_3 = h/2\} \cap \{-1/8 \le X_2 \le l/2\}$.

<figure><img src="image_9.png" /><figcaption>Fig. 6. Energy and angular momentum of quasi-incompressible cube for conserving scheme with $\Delta t = 0.005$ s.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Schematic of elastic plate showing applied load and zero displacement boundary condition.</figcaption></figure>

In the spatial discretization of the plate we use 2 elements through the thickness and a uniform discretization in the $X_2-X_3$ plane consisting of 25 elements, for a total of 50 elements. Fig. 8 shows snapshots of the plate motion in the interval $[0, 0.0045 \text{ s}]$. This motion was computed with the quasi-incompressible conserving scheme with a step size $\Delta t = 0.0001 \text{ s}$. We employed the Uzawa iteration strategy and iterated until $\max_e |\bar{G}(\Theta_n^e)| < 5 \times 10^{-10}$ for each time step.

Fig. 9 shows total energy and angular momentum time histories. For this system, the total energy is an integral after the loads are removed at $t = 0.0040 \text{ s}$, but the total angular momentum is not an integral due to the displacement boundary conditions. These features are clearly captured by the conserving scheme.

### 5.3. Inversion of a spherical cap

Consider a thick spherical cap with dimensions $r = 0.1 \text{ m}$, $\beta = 0.008 \text{ m}$ and $\alpha = \pi/4$ as shown in Fig. 10. The cap is composed of a homogeneous, compressible, elastic material of Ogden type with parameters as given in (61) and density $\rho_0 = 100 \text{ kg/m}^3$. We suppose the cap is initially at rest and subject to a system of equilibrated force distributions $f_1$ and $f_2$ that vanish after a short period of time. The force distribution $f_1$ is applied to the convex face of the cap, acts in the negative $X_1$-direction, is uniform in a disk centered about the $X_1$-axis and has a resultant given by $p(t)(-f, 0, 0)$ where $f = 160 \text{ kN}$ and

$$p(t) = \begin{cases} t, & 0 \le t \le 0.001 \text{ s}, \\ 0, & t > 0.001 \text{ s}. \end{cases} \qquad (70)$$

<figure><img src="image_11.png" /><figcaption>Fig. 8. Motion of quasi-incompressible elastic plate computed with conserving scheme. The direction of time is left-to-right, top-to-bottom.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 9. Total energy and angular momentum of quasi-incompressible plate for conserving scheme with $\Delta t = 0.0001$ s.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Schematic of elastic cap showing geometry and applied loads.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Snapshots of cap motion computed with conserving scheme. The direction of time is left-to-right, top-to-bottom.</figcaption></figure>

The force distribution $f_2$ is applied along the edge of the cap as shown in the figure, acts in the positive $X_1$-direction and has a resultant given by $p(t)(f, 0, 0)$.

In this example, we compare two time integration schemes: the implicit mid-point rule and the energy-momentum preserving scheme (10). In the spatial discretization of the cap we use 4 elements through the thickness and a quasi-uniform discretization of each constant-radius surface into 128 elements, for a total of 512 elements. Fig. 11 shows snapshots of the cap motion at time increments of 0.0005 s computed with the conserving scheme.
<figure><img src="image_15.png" /><figcaption>Fig. 12. Energy and angular momentum of cap for conserving scheme with $\Delta t = 0.0005$ s.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Energy and angular momentum of cap for mid-point rule with $\Delta t = 0.0005$ s.</figcaption></figure>

Figs. 12 and 13 show total energy and angular momentum time histories computed using a time step $\Delta t = 0.0005$ s. For the conserving scheme we see that the total energy is conserved after the loads are removed at $t = 0.0010$ s. Also, since the force distributions were equilibrated and symmetric with respect to the $X_1$-axis, we see that the total angular momentum of the cap remains zero throughout the simulation. Note that the mid-point rule experiences numerical blow-up for the given time-step size.

## 6. Closing remarks

In this article we have considered energy and momentum conserving time discretization schemes for general initial boundary-value problems in finite-deformation elastodynamics. Conserving schemes were developed in a weak form for both compressible and incompressible hyperelastic material models, implemented using finite element discretizations in space and applied to three example problems. Compared to the implicit mid-point rule, the conserving schemes were seen to have similar accuracy properties for small time steps, and superior stability properties for large time steps.

Strategies for the solution of the fully discrete system were not addressed in this article. However, we remark that all the example problems were computed using a standard Newton iteration scheme with a consistent linearization. This linearization leads to a generally nonsymmetric Jacobian in the case of the conserving schemes presented here, and a symmetric Jacobian in the case of the implicit mid-point rule. Thus, our implementation of the conserving schemes were more computationally intensive than the mid-point rule. However, it may be possible to lessen the computational costs of these conserving schemes by employing more efficient solution strategies.
# Acknowledgements

The author is deeply indebted to the late Professor Juan C. Simo for his guidance during the early development of this article. The author also gratefully acknowledges the financial support of the US National Science Foundation.

# References

[1] U.M. Ascher, H. Chin, S. Reich, Stabilization of DAEs and invariant manifolds, Numerische Mathematik 67 (1994) 131–149.

[2] D.P. Bertsekas, Constraint Optimization and Multiplier Methods, Academic Press, New York, 1982.

[3] K.E. Brenan, S.L. Campell, L.R. Petzold, Numerical Solution of Initial-value Problems in Differential-Algebraic Equations, Elsevier, Amsterdam, 1989.

[4] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Elements Methods, Springer Series in Computational Mathematics, vol. 15, Springer, Berlin, 1991.

[5] C. Chen, W. von Whal, Das Rand-Anfangswertproblem für quasilineare wellengleichungen in sobolevraumen niedriger ordnung, Journal für die reine und angewandte Mathematik 337 (1982) 77–112.

[6] P.G. Ciarlet, Mathematical Elasticity, Vol. I. Three-dimensional Elasticity, Elsevier, Amsterdam, 1988.

[7] C.M. Dafermos, W.J. Hrusa, Energy methods for quasilinear hyperbolic initial boundary-value problems. Applications to elastodynamics, Archive for Rational Mechanics and Analysis 87 (1985) 267–292.

[8] D.J. Dichmann, J.H. Maddocks, An impetus-striction simulation of the dynamics of an elastica, Journal of Nonlinear Science 6 (1996) 271–292.

[9] D.G. Ebin, R.A. Saxton, The initial-value problem for elastodynamics of incompressible bodies, Archive for Rational Mechanics and Analysis 94 (1986) 15–38.

[10] D.G. Ebin, S.R. Simanca, Small deformations of incompressible bodies with free boundaries, Communication on Partial Differential Equations 15 (1990) 1589–1616.

[11] D.G. Ebin, S.R. Simanca, Deformations of incompressable bodies with free boundaries, Archive for Rational Mechanics and Analysis 120 (1992) 61–97.

[12] D.G. Ebin, Global solutions of the equations of elastodynamics of incompressible neo-Hookean materials, Proceedings of the National Academy of Sciences, USA 90 (1993) 3802–3805.

[13] G. Fichera, Existence Theorems in Elasticity, Handbuch der Physik, vol. VIa/2, Springer, Berlin, 1972.

[14] M. Fortin, A. Fortin, A generalization of Uzawa's algorithm for the solution of the Navier-Stokes equations, Communications in Applied Numerical Methods 1 (1985) 205–208.

[15] C. Füher, B. Leimkuhler, Numerical solution of differential-algebraic equations for constrained mechanical motion, Numerische Mathematik 59 (1991) 55–69.

[16] C.W. Gear, B. Leimkuhler, G.K. Gupta, Automatic integration of Euler-Lagrange equations with constraints, Journal of Computational and Applied Mathematics 12/13 (1985) 77–90.

[17] R. Glowinski, P. Le Tallec, Finite elements in nonlinear incompressible elasticity, in: J.T. Oden, G.F. Carey (Eds.), Finite Elements, Vol. V: Special Problems in Solid Mechanics, Prentice-Hall, Englewood Cliffs, NJ, 1981.

[18] R. Glowinski, P. Le Tallec, Augmented lagrangian and operator splitting methods in nonlinear mechanics, SIAM Studies in Applied Mathematics, 1989.

[19] O. Gonzalez, Time integration and discrete Hamiltonian systems, Journal of Nonlinear Science 6 (1996) 449–467.

[20] O. Gonzalez, J.C. Simo, On the stability of symplectic and energy momentum algorithms for nonlinear Hamiltonian systems with symmetry, Computer Methods in Applied Mechanics and Engineering 134 (1996) 197–222.

[21] O. Gonzalez, Mechanical systems subject to holonomic constraints: differential-algebraic formulations and conservative integration, Physica D 132 (1999) 165–174.

[22] O. Gonzalez, J.H. Maddocks, R.L. Pego, Multi-multiplier ambient-space formulations of constrained dynamical systems with an application to elastodynamics, Archive for Rational Mechanics and Analysis (accepted).

[23] M.E. Gurtin, An Introduction to Continuum Mechanics, Academic Press, New York, 1981.

[24] T.J.R. Hughes, T. Kato, J.E. Marsden, Well-posed quasi-linear second-order hyperbolic systems with applications to nonlinear elastodynamics and general relativity, Archive for Rational Mechanics and Analysis 63 (1977) 273–294.

[25] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, New York, 1987.

[26] W.J. Hrusa, M. Renardy, An existence theorem for the dirichlet problem in the elastodynamics of incompressible materials, Archive for Rational Mechanics and Analysis 102 (1988) 95–117.

[27] T. Kato, Quasi-linear equations of evolution with applications to partial differential equations, in: G. DaPrato, G. Geymonat (Eds.), Hyperbolicity, Centro Internazionale Matematico Estivo, II ciclo, Cortona, pp. 125–191.

[28] B.J. Leimkuhler, S. Reich, Symplectic integration of constrained Hamiltonian systems, Mathematics of Computation 63 (1994) 589–605.

[29] D.G. Luenberger, Linear and Nonlinear Programming, Addison-Wesley, Reading, MA, 1984.

[30] J.H. Maddocks, R.L. Pego, An unconstrained Hamiltonian formulation for incompressible fluid flow, Communications in Mathematical Physics 170 (1995) 207–217.
[31] J.E. Marsden, T.J.R. Hughes, Mathematical Foundations of Elasticity, Dover, New York, 1983.

[32] R.W. Ogden, Elastic deformation of rubber like solids, Mechanics of Solids (1982) 499-537.

[33] R.W. Ogden, Non-Linear Elastic Deformations, Ellis Horwood, Chichester, UK, 1984.

[34] W.C. Rheinboldt, Differential-algebraic systems as differential equations on manifolds, Mathematics of Computation 43 (1984) 473-482.

[35] R.D. Richtmyer, K.W. Morton, Difference Methods for Initial-value Problems, Interscience, New York, 1967.

[36] M. Sablé-Tougeron, Existence pour un problème de l'elasto-dynamique Neumann non lineaire en dimension 2, Archive for Rational Mechanics and Analysis 101 (1988) 261-292.

[37] S. Schochet, The incompressible limit in nonlinear elasticity, Communications In Mathematical Physics 102 (1985) 207-215.

[38] J.C. Simo, R.L. Taylor, K.S. Pister, Variational and projection methods for volume constraint in finite deformation elasto-plasticity, Computer Methods in Applied Mechanics and Engineering 51 (1985) 177-208.

[39] J.C. Simo, J.E. Marsden, P.S. Krishnaprasad, The Hamiltonian structure of nonlinear elasticity: the material and convective representation of solids, rods, and plates, Archives for Rational Mechanics and Analysis 104 (1988) 125-183.

[40] J.C. Simo, R.L Taylor, Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms, Computer Methods in Applied Mechanics and Engineering 85 (1991) 273-310.

[41] J.C. Simo, K.K. Wong, Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and angular momentum, International Journal for Numerical Methods in Engineering 31 (1991) 19-52.

[42] J.C. Simo, N. Tarnow, K.K. Wong, Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics, Computer Methods in Applied Mechanics and Engineering 1 (1992) 63-116.

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

[44] J.C. Simo, N. Tarnow, A new energy and momentum conserving algorithm for the non-linear dynamics of shells, International Journal for Numerical methods in Engineering 37 (1994) 2525-2550.

[45] J.C. Simo, N. Tarnow, M. Doblaré, Non-linear dynamics of three-dimensional rods: exact energy and momentum conserving algorithms, International Journal for Numerical Methods in Engineering 38 (1995) 1431-1474.

[46] C.A. Truesdell, W.M. Noll, The Nonlinear Field Theories of Mechanics, Handbuch der Physik, vol. III/3, Springer, Berlin, 1965.

[47] K.C. Valanis, R.F. Landel, The strain-energy function of a hyperelastic material in terms of the extension ratios, Journal of Applied Physics 38 (1967) 2997-3002.