Exact Energy and Momentum Conserving Algorithms for General Models in Nonlinear Elasticity

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.


Introduction.
In this paper we develop conserving time-integration schemes for initial boundary value problems in nonlinear hyperelasticity. We consider both compressible and incompressible hyperelastic material models and employ the formalism presented in Gonzalez 1996a,b] to construct schemes which preserve the total energy of the system, along with the total linear and angular momentum whenever these conservation laws are present in the underlying problem. Within the context of incompressible hyperelasticity, this paper extends to the in nite-dimensional case the formalism for constrained systems presented in Gonzalez 1996b].
In Section 2 we present a displacement-traction initial boundary value problem for a general compressible hyperelastic continuum body, and review some of the Hamiltonian structure for this system. Using the methods of the aforementioned references we then present a conserving time-stepping scheme. As we will see, the resulting scheme is applicable to general models in hyperelasticity and inherits the conservations laws of total energy, linear and angular momentum whenever they are present in the underlying problem. Here we note that conserving schemes for general hyperelastic models in elastodynamics go back to Simo & Tarnow 1992a], who designed schemes based on mean value arguments which involved the solution of an extra nonlinear equation in each time step. This class of schemes were subsequently extended to nonlinear shells in Simo & Tarnow 1992b], and to nonlinear rods in Simo, Tarnow & Doblar e 1993]. As we will see, the class of schemes presented herein do not involve any extra variables or equations; in particular, the conservation properties are inherent to the time discretization.
In Section 3 we develop a formulation of incompressible nonlinear hyperelasticity which ts within the framework for conserving integrators developed in Gonzalez 1996b]; in particular, we will consider a di erential-algebraic formulation of the dynamics (see e.g. Brenan, Campbell & Petzold 1989] for an explanation of this terminology). Such formulations follow naturally from the d'Alembert-Lagrange principle for constrained mechanical systems, and for conservative systems subject to holonomic constraints, it is well-known that such formulations may also be derived from a constrained variational principle (see e.g. Arnold 1989]). While many di erent di erential-algebraic formulations can be used to describe the same system, our attention will be focused on a particular formulation that is \Hamiltonian" in an appropriate sense. For purposes of numerical implementation, we introduce a quasi-incompressible formulation in which the pointwise incompressibility condition is relaxed. In particular, the resulting quasi-incompressible formulation may be viewed as a generalization to elastodynamics of the work of Simo & Taylor 1991]. Using the general framework for conserving schemes mentioned above, we then construct a conserving time-stepping scheme for the quasi-incompressible formulation.
In Section 4 we discuss the nite element implementation of the time-stepping schemes presented herein. Since the nite element implementation for the compressible case is fairly standard, we concentrate on the quasi-incompressible case. In particular, we introduce a general mixed nite element formulation and reduce it to a two-eld formulation. We then discuss a solution strategy which can be interpreted as the generalization to the dynamic case of the augmented multiplier methods used within the context of elastostatics, see e.g. Glowinski & Le Tallec 1981,1989] and Simo & Taylor 1991]. By using a nested iteration scheme, we will see that the two-eld formulation may be e ectively reduced to one eld.
Finally, in Section 5, we present some numerical examples of both the compressible and incompressible formulations. We employ hyperelastic models of the Ogden type (see e.g. Ogden 1982Ogden , 1984), and for the quasi-incompressible case we consider models in terms of deviatoric strain measures as employed in Simo & Taylor 1991].

General Compressible Hyperelasticity.
In this section we present a formulation of general compressible hyperelasticity. We rst review some standard terminology for the subject and then state a dispacement-traction initial boundary value problem for a hyperelastic continuum body. We next introduce a weak formulation, interpreted as a dynamic generalization of the principle of virtual work, which is suitable for numerical approximation. Given this weak formulation we discuss its Hamiltonian structure and then present a conserving time-stepping scheme. For more details on the physical and mathematical structure of the eld equations presented in this section we refer the reader to the introductory text of Gurtin 1981], and the more advanced works of Marsden & Hughes 1983], Ogden 1984], and Ciarlet 1988].

Preliminaries.
Let B denote a continuum body in 3-dimensional Euclidean space I R 3 . We assume that B, in its reference or undeformed state, occupies a closed subset B of I R 3 where B is a bounded path-connected open set with piecewise smooth boundary @B. In particular, we identify material particles of B in its reference state with points X = (X 1 ; X 2 ; X 3 ) 2 B and call B the reference con guration of B.
A deformation of B in I R 3 , relative to the reference con guration B, is a mapping ' : B ! I R 3 that is in some sense smooth, orientation-preserving, and injective in B. If we denote by Q the set of all deformations, then a motion of B in a time interval 0; T] is a curve in Q of the form 0; T] 3 t 7 ! ' t 2 Q. Hence, we regard a motion as a mapping ' : B 0; T] ! I R 3 , where for each t 2 0; T] the mapping ' t = '( ; t) : B ! I R is a deformation of B. Also, for any time t 2 0; T] in a motion we de ne a material velocity eld V t : B ! I R 3 by the expression V t = (@=@ )j =t '( ; ) = _ ' t . Let M 3 denote the vector space of all real 3 3 matrices equipped with the standard (Euclidean) inner product, de ned for any E; where det : M 3 ! I R denotes the determinant function. For any deformation ' t : B ! I R in a motion of B we de ne a deformation gradient eld F t : B ! M 3 + by F t (X) = D' t (X); 8X 2 B (2.5) where D' t (X) denotes the derivative of ' t evaluated at X, and we denote by C t : B ! S 3 PD the Cauchy strain eld de ned by the expression For any t 2 0; T] we denote by S t : B ! S 3 the second Piola-Kirchho stress eld for the continuum body. Note that for a hyperelastic body S t is related to the deformation ' t locally via a (frame-invariant) constitutive relation of the form S t (X) = 2D 2 W e (X; A) where W e : B S 3 PD ! I R is a hyperelastic stored energy function for the body relative to its reference state and D i denotes the derivative with respect to the i th argument. (2.8) where superposed dots denote di erentiation with respect to time and Div is the divergence operator in B relative to the cartesian coordinates (X 1 ; X 2 ; X 3 ). For any eld P : B ! I R 3 3 recall that Div P](X) 2 I R 3 for each X and has coordinates given by where @ i denotes the partial derivative with respect to the coordinate X i .
In the above system b : B (0; T] ! I R 3 is a prescribed body force density per unit reference volume, : B ! I R + is the mass density of the body in its reference con guration, N : ? ! I R 3 is the unit outward normal eld on ? @B, f : ? (0; T] ! I R 3 is a prescribed surface force density per unit reference surface area, g : ? ' (0; T] ! I R 3 is a speci ed motion on ? ' @B,' 0 : B ! I R 3 is a prescribed deformation of B, andV 0 : B ! I R 3 is a prescribed material velocity eld. The partial di erential equation (2.8) 1 expresses the local balance of linear momentum in B (whereas the local balance of angular momentum is implied by the symmetry of the stress eld S t ), (2.8) 2;3 are traction and displacement boundary conditions, respectively, and (2.8) 4;5 are initial data for the motion. Remark 2.1. In the above formulation we have excluded the general case in which the body force density b, the surface force density f and the boundary condition g, depend on the motion '. In addition, we can simplify matters further by making the assumption that b, f and g are independent of time, in which case we say that the body B is subjected to dead loads. While this level of simpli cation is not necessary for the developments of this paper, we assume henceforth that all loads are dead unless mentioned otherwise.
Typically, we are interested in weak solutions of (2.8); in particular, weak solutions satisfying an associated weak formulation interpreted as a dynamic generalization of the principle of virtual work, which we introduce next. v w 8v; w 2 L 2 (? ; I R 3 ): (2.13) Taking L 2 (B; I R 3 ) as the setting for our weak formulation, we de ne the conguration space Q for B as a smooth (linear) manifold in L 2 (B; I R 3 ), namely Q = f 2L 2 (B; I R 3 ) j D 2 L 2 (B; M 3 ); D (X) 2 M 3 + 8X 2 B and (X) = g(X) 8X 2 ? ' g (2.14) where g : ? ' ! I R 3 is the speci ed boundary condition appearing in (2.8) 3 . Remark 2.2. Following standard abuse of terminology, we call elements of Q \congurations" of B when really they are deformations of B relative to the reference con guration B.
Given the above de nition of Q, we now view a motion as a C 1 curve 0; T] 3 t 7 ! ' t 2 Q with associated velocity _ ' viewed as a curve in V Q , i.e. 0; T] 3 t 7 ! _ ' t 2 V Q , where the vector space V Q is de ned as V Q = f 2L 2 (B; I R 3 ) j D 2 L 2 (B; M 3 ) and (X) = 0 8X 2 ? ' g: (2.15) With the above notation our weak formulation of (2.8) may be stated as follows: Here we note that, if the deformation ' t 2 Q is su ciently smooth for each t 2 0; T], then the weak formulation given above is equivalent to the initial boundary value problem given in (2.8).
Ideally, before introducing an approximation scheme for (2.8) or (2.16) one would like to know whether or not these systems admit solutions, and if so, whether or not solutions are unique. That is, one would like to know if the systems are wellposed. In this paper we make no attempt to address these questions (the interested reader is referred to Marsden & Hughes 1983, Chapter 6] and Ciarlet 1988 Chapters 6,7] for partial results). We proceed under the assumption that the above systems do admit solutions in some sense, and given this, we seek to construct numerical approximations to them. Before doing so, however, we rst review some of the underlying Hamiltonian structure for the above systems.

Hamiltonian Structure and Conservation Laws.
In this section we review the Hamiltonian structure and conservation laws for (2.16). The treatment given here is brief and the interested reader is referred to Marsden & Hughes 1983] and Simo, Marsden & Krishnaprasad 1988] for more details. Since our goal is to eventually treat the incompressible case, we rst review a variational characterization of motion and then use the associated Euler-Lagrange equations to formulate (2.16) as a Hamiltonian system.
To begin, we assume the body B undergoes a motion in a time interval 0; T] under the in uence of conservative external forces with potential V ext : Q ! I R. In particular, for dead loading with body force density b : B ! I R 3 and surface force density f : ? ! I R 3 we have (2.20) With the above notation introduce a Lagrangian functional L : Q V Q ! I R for the body by the expression (2.21) Then, under the assumption of conservative loading, the weak formulation presented in (2.16) is related to a variational principle which we brie y outline next.
Consider the set C of all possible motions in 0; T] between two prescribed con gurations' 0 ;' T 2 Q, i.e. C = f' : 0; T] ! Q j ' 2C 1 ( 0; T]; Q); ' 0 =' 0 ; and ' T =' T g; (2.22) and note that, for any ' 2 C, the tangent space T ' C is de ned as T ' C = fu : 0; T] ! V Q j u 2C 1 ( 0; T]; V Q ) and u 0 = u T = 0g (2.23) where C 1 (Y; Z) denotes the set of continuously di erentiable mappings from Y into Z. If we introduce the action functional A : C ! I R by then Hamilton's principle states that the actual motion satis es the following variational problem: (V) Find ' 2 C such that DA(') u = 0 for all u 2 T ' C, where DA(') u denotes the directional derivative of A at ' 2 C in the direction u 2 T ' C. Using standard arguments from the calculus of variations, we deduce that any solution of (V) must satisfy the Euler-Lagrange equation in (0; T); 8 2 V Q : (2.25) Now, while the above variational principle (i.e. Hamilton's principle) pertains to the boundary value problem of motion (i.e. motion between xed con gurations), the Euler-Lagrange equation generated by the principle is suitable for the formulation of the initial value problem of motion in which ' 0 and _ ' 0 are speci ed. In particular, one can verify that (2.25) may be written as h ' t ; i + hF t S t ; D i = hb; i + hf; i ? in (0; T); 8 2 V Q ; (2.26) which is the weak equation (2.16) 1 with dead loads b and f. Given that (2.16) 1 is formally equivalent to (2.25) we can now easily put the weak equations (2.16) in Hamiltonian form as follows.
For any xed ' t 2 Q consider the functional L(' t ; ) : V Q ! I R and introduce the conjugate momenta t 2 V Q by the weak expression h t ; i = D 2 L(' t ; _ where (' 0 ;^ 0 ) 2 P is given initial data.
Since our weak initial boundary value problem for a general hyperelastic body is Hamiltonian, we can use the general theory for Hamiltonian systems to establish the following conservation laws: Conservation of the Hamiltonian. The Hamiltonian functional H is conserved along any solution ('; ) : 0; T max ) ! Q V Q in the sense that 1) G I = I R 3 with action I : G I P ! P de ned as I (c; (' t ; t )) = (' t + c; t ) and associated momentum map J I : P ! I R 3 de ned by called the linear momentum map.
2) G II = SO(3) with action II : G II P ! P de ned as II ( ; (' t ; t )) = ( ' t ; t ) and associated momentum map J II : P ! I R 3 de ned by J II (' t ; t ) = Z B ' t t (2.38) called the angular momentum map.
Given that H is invariant under the actions of G I and G II , and that these actions possess momentum maps, we have two more conservation laws for our system. In particular, for the case we are considering, the linear and angular momentum maps are conserved along any solution ('; ) : 0; T max ) ! Q V Q in the sense that d dt J I (' t ; t ) = 0 and d dt J II (' t ; t ) = 0 in (0; T max ): We next construct a time-stepping scheme for (2.35) which inherits the above conservation laws.

Conserving Time-Stepping Scheme.
For concreteness consider the case for which ? ' = ; and ? = @B, and suppose b : B ! I R 3 and f : ? ! I R 3 are the zero mappings. In this case the underlying system possesses three conservation laws which, ideally, we would like to design into a time-stepping scheme.
To begin, let P = Q V Q with points denoted by z t = (' t ; t ) and let DH(z t ) denote the derivative of H at z t . In view the results presented in Gonzalez 1996a], the construction of a conserving time integration scheme for (2.35) depends on the construction of a G I;II -equivariant discrete derivative for the function H : P ! I R.
This object (to be de ned in Section 3) is a two-point approximation to the exact derivative of H. For example, if we denote an equivariant discrete derivative for H by D G H( ; ), then D G H(z a ; z b ) DH ? z a +z b 2 (2.39) for any two points z a ; z b 2 P. Given an equivariant discrete derivative D G H we de ne a set of weak Hamiltonian di erence equations in P = Q V Q analogous to (2.34), namely h' n+1 ? ' n ; i = tD G 2 H(z n ; z n+1 ) h n+1 ? n ; i = ? tD G 1 H(z n ; z n+1 ) ) 8 2 V Q : (2.40) where t > 0 is the time step and the partial discrete derivatives D G 1 H and D G 2 H are de ned by the relation Using techniques which are discussed in detail when we consider the incompressible formulation, we construct the following equivariant discrete derivative where S algo : B ! S 3 , interpreted as an approximation to the second Piola-Kirchho stress eld, is de ned as follows.
For any X 2 B let W X e = W e (X; ) : S 3 PD ! I R where W e is the (frame-invariant) hyperelastic stored energy function presented in (2.7) and de ne DW X e : S 3 PD S 3 PD ! I R by DW X e (A n ; A n+1 ) = DW X e (A n+ 1 2 ) where ( ) n+ 1 2 = 1 2 ( ) n + ( ) n+1 ] and A = A n+1 ? A n . With the above notation the approximate stress eld S algo is de ned as S algo (X) = 2DW X e (C n (X); C n+1 (X) 8 2 V Q : (2.45) The above equations were derived for the case in which ? ' = ;, ? = @B, and b : B ! I R 3 and f : ? ! I R 3 are the zero mappings. As such, they represent a discrete Hamiltonian system with symmetry in the sense of Gonzalez 1996a], and inherit the conservation laws discussed above. Extending the above results to the general displacement-traction problem with dead loads yields the following timestepping scheme for (2.35): Given b : B ! I R 3 , f : ? ! I R 3 , and t > 0 nd a sequence (' n ; n ) N n=0 in Q V Q such that h' n+1 ? ' n ; i = t where' 0 and^ 0 are initial data. As for the underlying system, the di erence equations in (2.46) possess the following discrete conservation laws: Conservation of the Hamiltonian. The Hamiltonian functional H is conserved along any solution sequence (' n ; n ) N n=0 in Q V Q in the sense that H(' n ; n ) = H(' 0 ; 0 ) for 0 n N: To see this, choose = n+1 ? n in (2.46) 1 and choose = ' n+1 ? ' n in (2.46) 2 .
The result then follows by subtracting the two expressions. Conservation of Linear and Angular Momentum. Consider the case for which ? ' = ;, ? = @B, and b : B ! I R 3 and f : ? ! I R 3 are the zero mappings. For this case we have the following additional conservation laws: 1) The linear momentum map J I is conserved along any solution sequence (' n ; n ) N n=0 in Q V Q in the sense that J I (' n ; n ) = J I (' 0 ; 0 ) for 0 n N: To see this, consider any 2 I R 3 and choose (X) = for all X 2 B. The result then follows directly from (2.46) 2 .
2) The angular momentum map J II is conserved along any solution sequence (' n ; n ) N n=0 in Q V Q in the sense that J II (' n ; n ) = J II (' 0 ; 0 ) for 0 n N: (2.49) To see this, consider any 2 I R 3 and choose (X) = ' n+ 1 2 (X) in (2.46) 2 , and choose (X) = n+ 1 2 (X) in (2.46) 1 . The result then follows using standard manipulations and the de nition of the angular momentum map. Remark 2.3. The time-stepping scheme in (2.46) was developed under the assumption that the loads b and f were dead. Note, however, that the general methods presented in Gonzalez 1996a] are not limited to this case. That is, conserving schemes can be constructed for general external loads b and f that are derivable from a potential. In this case we replace the term thb; i + thf; i ? in (2.46) 2 by ? tDV ext (' n ; ' n+1 ) , where DV ext is a discrete derivative for the potential V ext .
We next present a formulation of incompressible hyperelasticity and develop a conserving time-stepping scheme using the same ideas as above.

General Incompressible Hyperelasticity.
In this section we present a formulation of general incompressible hyperelasticity. Our goal is to develop a formulation for this system which falls within the framework of Gonzalez 1996b] and to use methods presented therein to construct conserving time-stepping schemes.

Hamiltonian Formulation.
Consider again the variational formulation of general compressible hyperelasticity given in Section 2.4. Let Q denote the con guration space for B considered earlier and let C denote the space of possible motions of B in a time interval 0; T] between xed endpoints. For the moment we assume Q and C contain only smooth mappings to justify the operations in this section.
Now, for the incompressible case, we require that any motion ' 2 C of our body B satisfy a pointwise constraint of the form G(D' t (X)) = 0 in B 0; T] for a smooth function G : M 3 + ! I R given by G(M) = det M] ? 1: In particular, the condition G(D' t (X)) = 0 in B 0; T] says that, for any t 2 0; T], the deformation ' t : B ! I R 3 is volume-preserving at each point X 2 B (see e.g. Ogden 1984] or Marsden & Hughes 1983] for more details).
Given the above constraint, the basis for our formulation of incompressible hyperelasticity is the following constrained variational problem: Find a stationary point ' 2 C of the functional A : C ! I R subject to the pointwise condition G(D' t (X)) = 0 for all (X; t) 2 B 0; T], where A : C ! I R is the action functional de ned in (2.24). Since constrained variational problems are often cumbersome to deal with directly we use the method of (Lagrange) multipliers to transform the constrained variational problem into an equivalent unconstrained one as follows.
Let E denote the set of smooth functions : B ! I R, and de ne a space E as the set of smooth mappings : 0; T] ! E. Next, consider an augmented action functionalÂ : C E ! I R de ned bŷ where L denotes the Lagrangian functional de ned in (2.21). Then, under suitable conditions on the constraint, the constrained problem (CV) is formally equivalent to the following unconstrained variational problem: (UV) Find ('; ) 2 C E such that DÂ('; ) (u; v) = 0 for all (u; v) 2 T ' C T E.
To see that the weak equation (3.3) 1 is a correct description of the motion we substitute the de nition of L and G into (3.3) 1 to get h ' t ; i + F t S t + det F t ] t C ?1 t ; D = hb; i + hf; i ? To pass to a formulation with some Hamiltonian structure we consider, for any xed ' t 2 Q and t 2 E, the functionalL(' t ; ; t ) : V Q ! I R and introduce the conjugate momenta t 2 V Q by the weak expression Next, we construct a mapping^ ' t ; t : V Q ! V Q which solves (3.6) for _ ' t in the sense that h _ ' t ; i = h^ ' t ; t ( t ); i = h^ (' t ; t ; t ); i 8 2 V Q : (3.7) In particular, from (3.6) and the de nition of the Lagrangian we have _ ' t =^ (' t ; t ; t ) = ?1 t : (3.8) At this point, using (3.7) and (3.6), we can write the Euler-Lagrange equations (3.5) as h _ ' t ; i = h^ (' t ; t ; t ); i h_ t ; i = D 1L (' t ; v; t ) v=^ (' t ; t ; t ) 9 = ; 8 2 V Q : (3.9) Introducing a functionalĤ : we see that the Euler-Lagrange equations (3.5) lead to the following initial value problem: ( 3.11) where' 0 and^ 0 are initial data. In the above system (3.11) 1;2 are interpreted as weak Hamiltonian evolution equations for (' t ; t ) 2 Q V Q , and expression (3.11) 3 is viewed as an equation which determines t ; in particular, the system may be interpreted as a Hamiltonian di erentialalgebraic equation within the context of Gonzalez 1996b].
Remark 3.1. The above di erential-algebraic formulation is not the only one which may be used to describe the dynamics of the incompressible system. For example, rather than append the constraint G(F t (X)) = 0 to Hamilton's principle, one may append the time-di erentiated constraint DG(F t (X)) : _ F t (X) = 0 to obtain a di erent formulation. Such formulations, obtained by appending a timedi erentiated con guration constraint to Hamilton's principle, are related to the so-called vakonomic approach to constrained systems. For more details on this approach see Kozlov 1983], Arnold 1988 Here we note that, from the numerical point of view, the above formulation is susceptible to di culties. The main problem lies in the fact that any smooth solution of this system satis es the pointwise condition det F t (X)] = 1 in B 0; T]. It is well-known that this class of deformations cannot be approximated well with standard (low-order) nite element spaces; in particular, one can expect di culties such as the phenomena of locking (see e.g. Hughes 1987] for a summary account). We thus introduce another formulation to alleviate these expected di culties.

A Quasi-Incompressible Formulation.
In this section we reformulate (3.11) to yield a system more suitable for numerical implementations. Following Simo, Taylor & Pister 1985] and Simo & Taylor 1991], the main idea is to introduce a eld t : B ! I R + which in some sense approximates the jacobian eld J t := det D' t ] : B ! I R + , and then enforce the condition t (X) = 1 in B 0; T]. Roughly speaking, since the incompressibility condition is enforced on t rather than J t , the formulation is called quasi-incompressible.
To begin, let Q denote the con guration space for B de ned in (2.14) and consider the spaces E = L 2 (B; I R) and E + = L 2 (B; I R + ). Now, for any t 2 E + and ' t 2 Q de ne a modi ed deformation gradient eld F t : B ! M 3 + by the expression F t (X) = t (X) 1=3 F dev t (X); (3.12) where F dev t : B ! M 3 + is the deviatoric part of the actual deformation gradient eld F t , i.e. F dev t (X) = J t (X) ?1=3 F t (X); (3.13) and de ne an associated modi ed Cauchy strain eld C t : B ! S 3 PD by C t (X) = F T t (X) F t (X): (3.14) Let W : B M 3 + ! I R be the stored energy function de ned in (2.19) and de ne a modi ed internal potential energy V int : Q E + ! I R by V int (' t ; t ) = Z B W(X; F t (X)): (3.15) As before, denote by V ext : Q ! I R the potential of the external loads. Let H : Q V Q E + E E ! I R denote a modi ed Hamiltonian functional of the form H(' t ; t ; t ; p t ; t ) = 1 2 t ; ?1 t + V int (' t ; t ) + V ext (' t ) + hp t ; (J t ? t )i + t ; G( t ) ; (3.16) where G : I R + ! I R is given by G( ) = ? 1; (3.17) and consider the following weak initial value problem: ( 3.18) where' 0 and^ 0 are initial data. The relation between (3.18) and (3.11) is established in the following Proposition 3.1. If for each t 2 0; T] the elds (' t ; t ) 2 Q E + are su ciently smooth, then the formulation in (3.18) is equivalent to that given in (3.11).
Proof. Using (3.16) we arrive at the following directional derivatives of H: ( 3.27) which is precisely the system in (3.11).
From now on we will work exclusively with the quasi-incompressible formulation given in (3.18), which, as we will see later, is better suited for numerical implementations than its counterpart in (3.11). We next show that the quasi-incompressible formulation possesses conservation laws analogous to those of the general compressible formulation in (2.35).

Conservation Laws.
In this section we show that the quasi-incompressible formulation of incompressible elastodynamics given in (3.18) possesses conservation laws analogous to those of the general compressible system. We state these results in a sequence of propositions.
We next construct a time-stepping scheme for the quasi-incompressible formulation that inherits the above conservation laws.

Conserving Time-Stepping Scheme.
For concreteness consider the case for which ? ' = ; and ? = @B, and suppose b : B ! I R 3 and f : ? ! I R 3 are the zero mappings. In this case the underlying system possesses three conservation laws which, ideally, we would like to preserve under time discretization. Since the quasi-incompressible formulation outlined in (3.18) de nes a (weak) Hamiltonian di erential-algebraic system in the sense of Gonzalez 1996b], we use the results of that paper to construct a conserving time-stepping scheme.
To begin, let P = Q V Q E + E E with points denoted by t = (' t ; t ; t ; p t ; t ) and let V P = V Q V Q E E E. The construction of a conserving time-stepping scheme for (3.18) depends on the construction of a G I;II -equivariant discrete derivative for the modi ed Hamiltonian functional H : P ! I R.
For the problem we are considering the modi ed Hamiltonian H is invariant under the actions of G I and G II on Q V Q and the underlying system preserves the momentum maps J I and J II . The following results show that these conservation laws are preserved by (3.39).
The particular case of dead loads is summarized in the next section.

Numerical Implementation.
In this section we discuss the nite element implementation of the time-stepping scheme (3.67) for quasi-incompressible elastodynamics. We rst discuss a mixednite element formulation for our scheme and then reduce it to a two-eld formulation for e cient implementation. We then discuss a solution strategy which can be interpreted as the generalization to the dynamic case of the augmented multiplier methods used in Simo & Taylor 1991]. For more details on these methods we refer the reader to Bertsekas 1982], and to Glowinski & Le Tallec 1981,1989] for applications to nite elasticity.

Mixed Finite Element Formulation.
Consider a partition of the domain B I R 3 into m elem 1 nonoverlapping subdomains B e (e = 1; : : :; m elem ), each de ned by a set of m en 1 element nodes denoted by X e a 2 B e (a = 1; : : :; m en ). We assume the elements B e are de ned such that B = e B e , and we introduce, for each e = 1; : : :; m elem , the sets e = fa j a = 1; : : :; m en g e g = fa 2 e j X e a 2 ? ' g: In accordance with a standard isoparametric discretization,  Note that the functions in E h and E h + are smooth within element domains, but are allowed to be discontinuous across element boundaries.
A mixed nite element formulation of our time-stepping scheme is now constructed by simply restating (3.67) with the nite-dimensional spaces Q h , V h Q , E h and E h + . We then use standard arguments to reduce the formulation to a system of nonlinear algebraic equations. (See e.g. Hughes 1987] for an exposition of these standard ideas.) Remark 4.1. Since the nite element formulation of (3.67) has the same structure as (3.67) itself, it is easy to see that the nite element formulation possesses conservation laws analogous to those of (3.67).

Reduction to a Two-Field Formulation.
In this section we eliminate the elds h n+1 , h n+1 and p h n+1 from the nite element formulation of (3.67). The result will be a two-eld nite element formulation in which we need only solve for ' h n+1 and h n+1 at each step of the algorithm.
To begin, we assume for simplicity that the density : B ! I R + is a constant function. Then, in view of (3.67) 1 , we can eliminate h n+1 from the nite element formulation using the expression Solving (4.12) for e j;n+1 (j = 1; : : :; m dis ) and substituting the result into ( That is, the restriction of h n+1 to any element is constant and equal to the average jacobian value over B e . In this case, the incompressibility condition we are enforcing is h n+1 j B e = 1 (e = 1; : : :; m elem ). Hence, the incompressibility condition is being enforced on an averaged jacobian, and not the point-wise jacobian J h n+1 , which helps circumvent numerical di culties such as locking. (See e.g. Hughes 1987] where W dev e : B S 3 PD ! I R is interpreted as a deviatoric stored energy function, U(X; r 2 ) = U(X; r) for some function U : B I R + ! I R which is convex in r for each X 2 B, and A dev 2 S 3 PD is the deviatoric part of A de ned by A dev = det A] ?1=3 A.
In view of (3.60) and (4.17), the modi ed stored energy function W appearing in our formulation becomes W(X; A; r) = W dev e (X; A dev ) + U(X; r): (4.18) For our quasi-incompressible formulation we can assume, without loss of generality, that the function U is independent of X 2 B and of the form where K > 0 is a constant and : I R + ! I R is convex, non-negative and satis es (r) = 0 if and only if r = 1.
In accordance with standard penalty formulations of constrained variational problems, one expects that any motion of a hyperelastic body with stored energy function as described above satis es t (X) ! 1 in B 0; T] as K ! 1. In fact, if we drop the quasi-incompressibility condition (4.15) 2 and set h n = 0 for all n, then the resulting nite element formulation (with stored energy function as described above) could be interpreted as a standard penalty formulation of quasi-incompressible elastodynamics. These types of formulations, however, are typically ill-conditioned (see e.g. Luenberger 1984] or Bertsekas 1982]).
On the other hand, keeping the formulation in (4.15) as it is (i.e. with a multiplier eld h n enforcing the quasi-incompressibility condition), and employing a stored energy function as described above, we can interpret our system as an augmented lagrangian formulation, in the sense of optimization theory. That is, we enforce the quasi-incompressibility condition (4.15) 2 by means of a multiplier eld and a penalty parameter K. Given this interpretation of our formulation, we have at our disposal various e cient solution strategies (see e.g. Luenberger 1984], Bertsekas 1982] and Fortin & Fortin 1985].) One example is the classical Uzawa algorithm, which, using k as an iteration index, may be summarized as follows: 1) Let (' h n ; h n ; h n ; p h n ; h n ) 2 P h be given and set h;k n+1 j k=0 = h n . (4.20) 2) Note that the above algorithm e ectively reduces the two-eld formulation to one eld. Hence, Step 2 can be carried out using standard nite element codes suitable for nonlinear problems.

Numerical Examples.
In this section we illustrate the performance of the time-stepping schemes presented herein by means of three example problems: the tumbling of an elastic block, the stretching of a rectangular plate, and the inversion of a spherical cap. In all cases we employ hyperelastic models of the Ogden type, namely W e (C) =w( 1 ; 2 ; 3 ) = Here C is the right Cauchy strain tensor, 2 A (A = 1; 2; 3) are the eigenvalues of C, m and m (m = 1; : : :; k) are parameters, and k is the number of terms in the Ogden model. As discussed in Ogden 1982Ogden ,1984 and Valanis & Landel 1967] models of the form presented above are particularly useful for rubber-like materials.
An extension of the above model to incompressible materials is given by For the mixed nite element formulation we construct the spaces Q h and V h Q using standard trilinear shape functions. For the spaces E h and E h + we take m dis = 1 and choose constant interpolation functions over the elements.
Example 5.1. Tumbling Elastic Block. In this simple example we consider a cube of dimension l = 0:02m composed of a homogeneous elastic material of Ogden type with parameters as given in (5.3) and density = 1000kg=m 3 . In the simulation, the cube is initially at rest and is subjected to tractions F1 and F2 for a short period of time. The forces are then removed and the cube is allowed to tumble freely in space. Figure 5.1 shows a schematic of the cube in its reference con guration along with the mentioned forces. Speci cally, the traction F1 is a spatially uniform force distribution on the face X3 = ?l=2 de ned by F1 = p(t)(0; ?f=4; ?f=8) and the traction F2 is a spatially uniform force distribution on the face X3 = l=2 de ned by F2 = p(t)(0; f; f=2) where f = 32N=cm 2 and p(t) is given by p(t) = t; 0 t 0:005s 0; t > 0:005s. For this example problem 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-stepping schemes: the standard mid-point rule and the exact energy-momentum preserving scheme in (2.46). For the quasi-incompressible formulation we use the conserving scheme presented in (3.67). For this case, we employ the nested iteration strategy outlined in the previous section and iterate until max e j G( e n )j < 5 10 ?10 for each time step. Figure 5.2 shows snapshots, at time increments of 0:005s, in the motion of the cube during the time interval 0; 0:1s]. The plot was computed using the compressible formulation, however, the incompressible motion is indistinguishable for the given resolution. For clarity, we show the snapshots in the X2-X3 plane only.   Figure 5.5 shows how the conserving scheme compares with the mid-point rule with respect to accuracy for the compressible formulation. For this comparison we considered the above motion in the time interval 0:0100s; 0:0125s]. For each scheme we supplied the same initial conditions at t 0 = 0:0100s and computed the solution out to a time t 1 = 0:0125s. Denoting the solution at time t 1 for a given scheme and time step by ' t 1 ; t , we de ne the L 2 -error in displacements for that scheme and time step by L2 Error = Z B j' t 1 ; t ? ' t 1 ;conv j 2 1=2 where ' t 1 ;conv is the solution at t 1 computed using a time step of t = 2:5 10 ?7 s. As shown in Figure 5.5 the conserving scheme and the mid-point rule have nearly identical accuracy properties, at least for the model problem considered.  Figure 5.6 shows total energy and angular momentum time histories, computed using a time step size of t = 0:005s, for the conserving scheme applied to the quasiincompressible formulation. As mentioned above, we employed the nested iteration strategy outlined in the previous section and iterated until max e j G( e n )j < 5 10 ?10 for each time step. We next consider another example problem using the quasi-incompressible formulation. In contrast to the last example, the deformations in this next problem are rather severe. Example 5.2. Stretched Elastic Plate. In this example we consider a thick rectangular plate of thickness w = 0:002m, length l = 0:01m and height h = 0:01m as shown in Figure 5.7. The plate is composed of a homogeneous elastic material of Ogden type with parameters as given in (5.3) and density = 1000kg=m 3 . In the simulation, the plate is initially at rest and is subjected to a traction F for a short period of time. The force is then removed and the plate is allowed to vibrate freely subject to a zero displacement boundary condition along the shaded region shown in the gure. Speci cally, the traction F is applied to the face X3 = ?h=2, is uniform across the thickness, and varies in the X2-direction according to F(X2; t) = p(t)(f=2; 0; ?f); ? l 2 X2 < 0 p(t)(0; 0; ?f); 0 X2 l 2 (5.8) where f = 5kN=cm 2 and p(t) is given by p(t) = t; 0 t 0:004s 0; t > 0:004s. (5.9) The zero displacement boundary condition is imposed along the surface de ned by fX3 = h=2g \ f?l=8 X2 l=2g.
In the spatial discretization of the plate we use 2 elements through the thickness and a uniform discretization in the X2-X3 plane consisting of 25 elements, for a total of 50 elements. Figure 5.8 shows snapshots at various times in the motion of the plate computed using the conserving scheme for the quasi-incompressible formulation with a time step of t = 0:0001s. We employed the nested iteration strategy outlined in the previous section and iterated until max e j G( e n )j < 5 10 ?10 for each time step. Figure 5.9 shows total energy and angular momentum time histories for the computed motion of the elastic plate. Note that total energy is conserved after the loads are removed at t = 0:0040s, but the angular momentum is not conserved due to the displacement boundary condition. We next consider another example problem using the compressible formulation. In contrast to the elastic cube example, the deformations will be more severe.
Example 5.3. Inversion of a Spherical Cap. In this example we consider a thick spherical cap with dimensions r = 0:1m, = 0:008m and = =4 as shown in Figure 5.10. The cap is composed of a homogeneous elastic material of Ogden type with parameters as given in (5.3) and density = 100kg=m 3 . In the simulation, the cap is initially at rest and is subjected to a system of equilibrated force distributions F1 and F2 for a short period of time. The forces are then removed and the cap is allowed to vibrate freely in space. Speci cally, the force distribution F1 is applied to the convex face of the cap, acts in the negative X1 direction and is uniform in a disk centered about the X1 axis. The force F1 has a resultant given by p(t)(?f; 0; 0) where f = 160kN and p(t) = t; 0 t 0:001s 0; t > 0:001s. (5.10) The force distribution F2 is applied along the edge of the cap as shown in the gure, acts in the positive X1 direction and has a resultant given by p(t)(f; 0; 0). In the spatial discretization of the cap we use 4 elements through the thickness and a somewhat uniform discretization of each surface of constant radius into 128 elements, for a total of 512 elements. Figure 5.11 show snapshots at time increments of 0:0005s in the motion of the cap computed using the conserving scheme.
Figures 5.12 and 5.13 show total energy and angular momentum time histories, computed using a time step of t = 0:0005s, for the conserving scheme and midpoint rule, respectively, applied to the spherical cap problem. For the conserving scheme we see that the total energy is conserved after the loads are removed at t = 0:0010s. Also, since the force distributions were equilibrated and symmetric with  respect to the X1 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. FIGURE 5.11. Snapshots from the back view at various times in the motion of the elastic spherical cap as computed using the conserving scheme for the compressible formulation. 6. Closing Remarks.
In this paper we have developed time-integration schemes for initial boundary value problems in nonlinear elasticity; in particular, conserving time-stepping schemes for general hyperelastic models, both compressible and incompressible. The approach taken herein was based on Gonzalez 1996a,b] where a general framework for conserving integrators for Hamiltonian systems is established. Using standard nite element discretizations in space for the compressible formulation, and a mixed nite element discretization for the incompressible case, we applied our results to three example problems. We saw that the conserving schemes perform relatively well as compared to a standard (symplectic) scheme: the implicit mid-point rule. In particular, we saw that the mid-point rule is prone to numerical blow-up at time steps not too large when compared with the scales of the overall or average motion. In addition, for relatively small time steps, we saw that the conserving schemes have accuracy properties similar to that of the mid-point rule.
While strategies for the actual solution of the fully discrete system were not addressed in this paper, we note in closing that all the examples in this paper were performed using a standard Newton iteration scheme with a consistent linearization. Within this framework note that the conserving schemes presented herein yield unsymmetric linearizations whereas the mid-point rule yields a symmetric linearization. Hence, in the approach we have taken, the conserving schemes were considerably more computationally intensive. However, we note that it is possible to lessen the computational costs of conserving schemes by employing more e cient solution strategies.