tceic.com

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

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

O. Gonzalez
? D? epartement de Math? ematiques, Ecole Polytechnique F? ed? erale de Lausanne, CH-1015 Lausanne, Switzerland

Abstract Implicit time integration schemes that inherit the conservation laws of total energy, linear and angular momentum are considered for initial boundary-value problems in ?nite-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 ?nite 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. Key words: 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 ?nite-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 ?nite-deformation, compressible elastodynamics we consider
Preprint submitted to Elsevier Preprint 18 April 1999

the weak formulation
? ?

˙ · η d? = ?

?

V · η d?
?

˙ · η d? + ρ0 V

S(?) : F(?)T ?η d? =
?

b · η d? +

Γσ

f · η dΓ

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

?η ∈ V,

(1)

? ?

V|t=0 · η d? = ?|t=0 · η d? =

? ?

V0 · η d? ? 0 · η d?

and for ?nite-deformation, incompressible elastodynamics ˙ · η d? = ??
? ? V · η d? ? [S(?)

˙ · η d? + ρ0 V

+ 2λDG(C(?))] : F(?)T ?η =
?

b · η d? +

Γσ

f ·η

0=

?

φG(C(?)) d?

? ? ? ? ? ? ? ? ? d? ? ? ? ? ? ? ?? ? dΓ ? ? ? ?? ? ?η ? ? ? ? ?φ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

∈V ∈E

(2)

? ?

V|t=0 · η d? = ?|t=0 · η d? =

? ?

V0 · η d? ?0 · η d?

where ? ? R3 is an open set that represents the reference state of a continuum body with material points X ∈ ?, ?(X, t) ∈ R3 is the deformation ?eld for the body relative to its reference state, V(X, t) ∈ R3 is the material velocity ?eld, F(X, t) := ??(X, t) ∈ R3×3 is the deformation gradient ?eld and S(X, t) ∈ R3×3 the second Piola-Kirchho? constitutive stress ?eld. Here ? denotes the gradient operator relative to the cartesian coordinates X in ?, an overdot denotes di?erentiation with respect to time t ≥ 0 and R3×3 is the vector space of all real 3 ×3 matrices, equipped with the standard (Euclidean) inner product A : B = Aij Bij where summation on repeated indices is implied. For hyperelastic bodies, the stress ?eld S is connected to the deformation ?eld ? via a local constitutive relation of the form S(?) = 2DW (C(?)) (3)

where W (C) ∈ R is a given strain energy density function, DW (C) is the derivative of W with respect to its tensor argument evaluated at C, and C(?) is the right Cauchy-Green stretch ?eld associated with the deformation ?, namely C = FT F. The strain energy density function W may depend upon the material point X, but this dependence is not shown for clarity of notation. In the incompressible case, the total second Piola-Kirchho? stress ?eld is of 2

the form S + det[C]1/2 λC?1 (4)

where λ(X, t) ∈ R is a material pressure ?eld that enforces the material incompressibility constraint G(C) = det[C]1/2 ? 1 = 0 in ?. (5)

1 Since DG(C) = 2 det[C]1/2 C?1 , we note that the total second Piola-Kirchho? stress ?eld (4) can be written in the form S + 2λDG(C).

In the compressible case, we seek functions (?, V) : [0, T ] → U × V satisfying (1) where U is a set of deformations satisfying essential boundary conditions on a given boundary patch Γd , V is a vector space of functions satisfying the homogeneous counterpart of the essential boundary conditions, ρ0 (X) > 0 is a speci?ed mass density ?eld, b(X, t) ∈ R3 is a speci?ed body force ?eld, f (X, t) ∈ R3 is a speci?ed traction ?eld on Γσ = ? ?\Γd , V0 (X) is a speci?ed initial material velocity ?eld and ?0 (X) a speci?ed initial deformation. In the incompressible case, we seek functions (?, V, λ) : [0, T ] → U × V × E satisfying (2) where E is a vector space of pressure ?elds 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 ?eld equations of ?nite-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 L(?, V), angular momentum J(?, V) and energy H (?, 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 L, but none preserve H for arbitrary constitutive laws. Moreover, it was shown that the explicit central di?erence scheme is the only member of the Newmark family that preserves J. The main result in [42,43] was the construction of implicit schemes for (1) that preserve L, J and 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, integralpreserving time discretization schemes of [19,20] for the treatment of (1). We show that implicit schemes which preserve L, J and 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 [42,43,20] suggest that 3

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 di?erentialalgebraic evolution equation [3,34] in a function space. Two important sets associated with the constraint (5) are Q = {? ∈ U | G(C(?)) = 0 in ?} which we call the intrinsic con?guration space, and the set M = {(?, V) ∈ U × V | G(C(?)) = 0 and DG(C(?)) : [?VT ?? + ??T ?V] = 0 in ?} (7) (6)

which we call the intrinsic state space for the system. Any su?ciently smooth solution of (2) has the property that ?(t) ∈ Q and (?, V)(t) ∈ M. Moreover, depending on boundary conditions and external loads, the functionals L|M , J|M and H |M may be integrals. The formulation (2) is only one of many possible descriptions of incompressible elastodynamics. One may alternatively consider stabilized di?erential-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 Q and M, along with the integrals L|M , J|M and H |M whenever they are present in the underlying problem. Due to special di?culties associated with direct spatial discretization of (2), we will instead consider a modi?ed quasi-incompressible formulation along the lines of [38,40], and generalize the time-reversible, integral-preserving time discretization schemes developed in [21] to this modi?ed formulation. We show that time discretization schemes can be constructed for this formulation which preserve Q, the total energy functional H |Q , and the momentum functionals J|Q and L|Q . However, the schemes will not generally preserve the velocity-level constraints de?ned by 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 quasiincompressible formulation of (2), outline conservation laws and present a conserving time discretization scheme for this modi?ed formulation. In Section 4 4

we consider the ?nite element spatial discretization of the time-discretized systems. In the quasi-incompressible case, we introduce a general mixed ?nite element discretization, reduce the general discretization to a two-?eld formulation and discuss an e?cient 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 (?, V) ∈ U × V of the system (1) we associate a total linear momentum L(?, V), a total angular momentum J(?, V) and a total energy H (?, V) de?ned as L(?, V) := J(?, V) := H (?, V) :=
? ρ0 V d? ?

? × ρ0 V d? + W (C(?)) d?
Γσ

1 2 ? 2 ρ0 |V |

?

?

b · ? d? ?

f ·?

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 (?, V)(t) be a solution of (1) for given boundary conditions speci?ed in U and 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 R= T= b d? + f dΓ (9)
Γσ

? ? ? ? ? ? ? ? d Γ. ?

? ? ? ? ? ? ? ? ? ?

(8)

? ?

Γσ

? × b d? +

? × f d Γ.

(1) If η = ξ ∈ V for a ?xed vector ξ, then d [ξ · L] = ξ · R. dt 5

(2) If η = ξ × ? ∈ V for a ?xed vector ξ and any ?xed time t, then d [ξ · J] = ξ · T. dt (3) If the boundary conditions and external loads are time-independent, then d H = 0. dt 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 ?nd 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 (?, V) : [0, T ] → U × V of (1) at discrete times tn = n?t where n is a positive integer and ?t > 0. Denoting approximations to ?(tn ) and V(tn ) by ?n and Vn , 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 ?nal result is an implicit, one-step scheme of the form
1 ? ?t (?n+1 ρ0 ? ?t (Vn+1

? ?n ) · η d? = ? Vn ) · η d?

? Vn+1/2 · η d?

=?

T ? ? S(?n , ?n+1 ) : F(?n+1/2 ) ?η d?

+

?

bn+1/2 · η d? +

Γσ fn+1/2

·η

? ? ? ? ? ? ? ? dΓ ?

? ? ? ? ? ? ? ? ? ?

?η ∈ V.

(10)

[(·)n + (·)n+1 ], Here (·)n+1/2 := 1 2 de?ned by

? (?n , ?n+1 ) is an algorithmic stress ?eld S (11)

? (?n , ?n+1 ) := 2 dW (Cn , Cn+1 ), S and dW (Cn , Cn+1 ) is a discrete derivative de?ned by dW (Cn , Cn+1 ) := DW (Cn+1/2 ) W (Cn+1 ) ? W (Cn ) ? DW (Cn+1/2 ) : M + M ||M||2

(12)

6

where M := Cn+1 ? Cn and ||M||2 = M : 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 (?n , Vn ) ∈ U × V (n ≥ 0) be a solution sequence of (10) for given boundary conditions, external loads b and f , and time step ?t > 0. Moreover, let Rn+1/2 and Tn+1/2 denote a resultant force and moment de?ned as Rn+1/2 = Tn+1/2 =
? ?

bn+1/2 d? +

Γσ

fn+1/2 dΓ (13)
Γσ

?n+1/2 × bn+1/2 d? +

?n+1/2 × fn+1/2 dΓ.

(1) If η = ξ ∈ V for a ?xed vector ξ, then 1 [ξ · Ln+1 ? ξ · Ln ] = ξ · Rn+1/2 . ?t (2) If η = ξ × ?n ∈ V and η = ξ × Vn ∈ V for a ?xed vector ξ and any n, then 1 [ξ · Jn+1 ? ξ · Jn ] = ξ · Tn+1/2 . ?t (3) If the boundary conditions and external loads are time-independent, then 1 [Hn+1 ? Hn ] = 0. ?t

PROOF. The results follow from the speci?c form of the assumed variations and may be established as follows. (1) Choosing η = ξ in (10)2 gives ρ0 (Vn+1 ? Vn ) · ξ d? = ?t bn+1/2 · ξ d? + fn+1/2 · ξ dΓ, (14)

?

?

Γσ

which is the required result for L. (2) To establish the result for J, we ?rst note that, if ξ × ?n ∈ V for any n, then ξ × ?n+1/2 ∈ V by linearity of V, and similarly ξ × Vn+1/2 ∈ V. Next, we choose a variation η = ξ × ρ0 Vn+1/2 in (10)1 to obtain 1 ξ · [?n+1 ? ?n ] × ρ0 Vn+1/2 d? = 0 ?t 7 (15)

?

and choose a variation η = ξ × ?n+1/2 in (10)2 to obtain 1 ξ · ?n+1/2 ×ρ0 [Vn+1 ? Vn ] d? ?t ? : FT ? =? S n+1/2 ξ Fn+1/2 d? + ξ · Tn+1/2
?

?

(16)

? is the skew tensor associated with ξ , that is, ξ × a = ξ ?a for any where ξ T ? ? vector a. Since Fn+1/2 ξFn+1/2 is a skew tensor and S is symmetric, the ?rst term in the right-hand side of (16) vanishes and we are left with 1 ξ · ?n+1/2 × ρ0 [Vn+1 ? Vn ] d? = ξ · Tn+1/2 . ?t (17)

?

The desired result then follows by adding (15) and (17), and expanding the vector products. (3) To establish the result for H , we ?rst note that (?n+1 ? ?n ) ∈ V for time-independent boundary conditions, and that Vn ∈ V for any n by assumption. Next, we use the variation η = ρ0 [Vn+1 ? Vn ] in (10)1, and η = [?n+1 ? ?n ] in (10)2 to deduce that
1 ? ?t

?n+1 ? ?n · ρ0 [Vn+1 ? Vn ] d? =
?

Vn+1/2 · ρ0 [Vn+1 ? Vn ] d?

1 ? ?t ρ0

[Vn+1 ? Vn ] · ?n+1 ? ?n d?
? ?

=? +

? : FT S n+1/2 [Fn+1 ? Fn ] d? b · ?n+1 ? ?n d? +
Γσ

f · ?n+1 ? ?n

? is symmetric and the de?nition of C, we have Using the fact that S
1? ? : FT S n+1/2 [Fn+1 ? Fn ] = 2 S : [Cn+1 ? Cn ] ,

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? d Γ. ?

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

(18)

(19)

?, and by de?nition of S
1? S 2

: [Cn+1 ? Cn ] = W (Cn+1 ) ? W (Cn ).

(20)

Combining (20) with (19), and subtracting (18)2 from (18)1 , leads to H (?n+1 , Vn+1 ) = H (?n , Vn ), 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 8 (21)

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 ? . If the stored use of the discrete derivative (12) in computing the stress ?eld S energy function W (C) is su?ciently di?erentiable, then one can show that dW (Cn , Cn+1) = DW (Cn+1/2) + O ||Cn+1 ? Cn ||2 . (22)

Hence, the expression for dW (Cn , Cn+1 ) given in (12) is well-de?ned in the limit ||Cn+1 ? Cn || → 0. Moreover, for su?ciently regular solutions of the algorithmic equations we have ? (?n , ?n+1 ) = S(?n+1/2 ) + O(?t2 ), S (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 ?eld Cn+1/2 , rather than the stretch ?eld of an interpolated con?guration C(?n+1/2 ), in evaluating the stress ?eld. The use of C(?n+1/2 ) introduces an arti?cial 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 Cn+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 ?nite 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 ?nite di?erence 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 di?culty lies in the fact that any smooth solution of this system satis?es the pointwise condition det[F] = 1 in ?, and this class of deformations cannot be well-approximated 9

by standard, low-order ?nite element spaces. In particular, one can expect the phenomena of locking to appear in direct treatments of (2) (see, for example, [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 di?culties. The main idea is to introduce a material ?eld Θ(X, t) > 0 which approximates the jacobian ?eld ζ (X, t) := det[F(X, t)] > 0, and then enforce the condition Θ = 1 in ?. Since the incompressibility condition is enforced on Θ rather than ζ , the formulation is called quasi-incompressible. A quasi-incompressible formulation of (2) may be outlined as follows. To any quasi-jacobian ?eld Θ and deformation ?eld ? we associate a modi?ed defor? by the expression mation gradient ?eld F ? := Θ1/3 Fdev , F (24)

where Fdev is the deviatoric part of the actual deformation gradient ?eld F, that is Fdev := ζ ?1/3 F, ? by and we de?ne an associated modi?ed Cauchy-Green stretch ?eld C ? := F ?T F ?. C (26) (25)

Let W (C) be the original strain energy function introduced in (3) and de?ne ? (C, Θ) by a modi?ed strain energy W ? ) = W (Θ2/3 ζ ?2/3 C). ? (C, Θ) := W (C W (27)

Then a quasi-incompressible formulation of (2) is to ?nd (?, V, λ, Θ, p) : [0, T ] → U × V × E × E × E such that ˙ · η d? = ??
? ? V · η d?

˙ · η d? + ρ0 V

0= 0= 0=

? ? ?

? (C(?), Θ) + λD G ? (Θ) ? p] d? φ[DΘ W φ[ζ (?) ? Θ] d? ? (Θ) d? φG

? ? ? ? ? ? ? ? ? T ? [ S ( ? , Θ) + 2 pDG ( C ( ? ))] : F ( ? ) ? η d ? ? ? ? ? ? ? ? ? ? ? ?η ? f · η dΓ = b · η d? +
? Γσ

∈V ∈E

? ? ? ?φ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

(28)

where ? (C(?), Θ) and G ? (Θ) = Θ ? 1. S(?, Θ) = 2DC W The relation between (28) and (2) is established in the following result. 10 (29)

Proposition 3 If (?, V, λ, Θ, p) : [0, T ] → U × V × E × E × E is a smooth solution of the quasi-incompressible formulation (28), then (?, V, λ) : [0, T ] → U × V × E is a solution of the original incompressible formulation (2).

PROOF. From (28)4 we see that, if ζ and Θ are continuous, then ζ = Θ. ? = F and W ? (C, Θ) = W (C), so W ? (C, Θ) is independent of This implies F Θ. By continuity of λ and p we deduce from (28)3 that p = λ, and so (28)2 becomes
?

˙ · η d? + ρ0 V

? [S(?)

+ 2λDG(C(?))] : F(?)T ?η d? =
? b · η d? + Γσ f · η dΓ,

(30) ?η ∈ V

? (Θ) = G(C), and so (28)5 which is just (2)2 . Finally, since ζ = Θ we obtain G 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(?, V) and a total angular momentum J(?, V) to any state (?, V, λ, Θ, p) of the quasi-incompressible system (28). We also associate to each state a total mod? (?, V, λ, Θ, p) de?ned as i?ed energy H ? (?, V, λ, Θ, p) := H
1 2 ? 2 ρ0 |V |

? (C(?), Θ) d? +W

+

As before, the weak initial boundary-value problem (28) may possess various ? depending on the integrals or constants of motion related to L, J and H essential boundary conditions and external loads. This fact is summarized in the following proposition, whose proof is straightforward. Proposition 4 Let (?, V, λ, Θ, p)(t) be a solution of (28) for given boundary conditions speci?ed in U and V, and given external loads b and f . Moreover, let R and T denote the resultant force and moment due to the external loads as ? satisfy conservation laws completely analogous de?ned in (9). Then L, J and H to those in Proposition 1.

? (Θ) d? p[ζ (?) ? Θ] + λG ? ? ? ? ? ? ? b · ? d? ? Γσ f · ? dΓ. ?
?

? ? ? ? ? ? ?

(31)

11

3.3 Energy-Momentum Conserving Time Discretization Here we seek to approximate solutions (?, V, λ, Θ, p) : [0, T ] → U×V×E×E×E of (28) at discrete times tn = n?t where n is a positive integer and ?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 ?nal result is an implicit, one-step scheme of the form
1 ? ?t (?n+1 ρ0 ? ?t (Vn+1

? ?n ) · η d? = ? Vn ) · η d? =

?

Vn+1/2 · η d?

?

? [S

? + 2pn+1/2 dG] : FT +

n+1/2 ?η

d?

? bn+1/2 · η d? +

Γσ fn+1/2 · η dΓ

0= 0= 0= where

? ? ?

1 ? C ) + λn+1/2 dG ? ? pn+1/2 ] ? Cn + dW φ[ 2 ( dW n+1

φ[ζn+1 ? Θn+1 ] d? ? n+1 d? φG

? ? ? ? ?φ ? ? ? d? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ?η

∈V ∈E

(32)

? = dW ? Θn+1 (Cn , Cn+1 ) ? Θn (Cn , Cn+1 ) + dW S ? Θ (Cn , Cn+1 ) = DC W ? (Cn+1/2 , Θ) dW +
? (Cn+1/2 ,Θ):M ? (Cn+1 ,Θ)?W ? (Cn ,Θ)?DC W W ||M||2

M

? C (Θn , Θn+1 ) = DΘ W ? (C, Θn+1/2 ) dW +
? (C,Θn+1/2 )·?Θ ? (C,Θn+1 )?W ? (C,Θn )?DΘ W W (?Θ)2

(33)

dG(Cn , Cn+1 ) = DG(Cn+1/2 ) +
G(Cn+1 )?G(Cn )?DG(Cn+1/2 ):M ||M||2

M

? (Θn , Θn+1 ) = D G ? (Θn+1/2 ) dG +
? (Θn+1 )?G ? (Θn )?D G ? (Θn+1/2 ):?Θ G 2 (?Θ)

≡1

M = Cn+1 ? Cn and ?Θ = Θn+1 ? Θn . Just as in the time-continuous problem, the algorithmic equations (32) may ? dependpossess various integrals or constants of motion related to L, J and H ing on the essential boundary conditions and the external loads. This fact is summarized in the following proposition. 12

Proposition 5 Let (?n , Vn , λn , Θn , pn ) ∈ U × V × E × E × E (n ≥ 0) be a solution sequence of (32) for given boundary conditions, external loads b and f , and time step ?t > 0. Moreover, let Rn+1/2 and Tn+1/2 denote a resultant force ? satisfy discrete conservation and moment as de?ned in (13). Then L, J and H laws completely analogous to those in Proposition 2.

PROOF. The results for L and J follows as for the compressible case upon ? + 2pn+1/2 dG is symmetric. To noting that the total algorithmic stress ?eld S ? establish the result for H , choose a variation η = ρ0 [Vn+1 ? Vn ] in (32)1 to obtain
ρ0 ? ?t [?n+1

? ?n ] · [Vn+1 ? Vn ] d? =

?

ρ0 Vn+1/2 · [Vn+1 ? Vn ] d?,

(34)

a variation η = [?n+1 ? ?n ] in (32)2 to obtain
ρ0 ? ?t [Vn+1

? Vn ] · [?n+1 ? ?n ] d? = ? + 2pn+1/2 dG] : FT +
? n+1/2 [Fn+1 Γσ

?

? [S

? Fn ] d? f · [?n+1 ? ?n ] dΓ,

(35)

b · [?n+1 ? ?n ] d? +

a variation φ = Θn+1 ? Θn in (32)3 , φ = pn+1 ? pn in (32)4 and φ = λn+1 ? λn in (32)5 to obtain 0= 0= 0=
1 ? ? [ 2 ( dWCn ? [ζn+1 ?

? Θn+1 ][pn+1 ? pn ] d?

? n+1 [λn+1 ? λn ] d?. G

? ? C ) + λn+1/2 dG ? ? pn+1/2 ][Θn+1 ? Θn ] d? ? ? + dW ? n+1 ? ? ? ? ? ? ? ?

?

(36)

? + 2pn+1/2 dG is symmetric and the de?nition of C, we Using the fact that S have
1 ? ? + 2pn+1/2 dG] : FT [S n+1/2 [Fn+1 ? Fn ] = 2 [S + 2pn+1/2 dG] : [Cn+1 ? Cn ], (37)

? and dG, and by de?nition of S
1 ? [S 2 1 ? ? (Cn , Θn )] + 2pn+1/2 dG] : [Cn+1 ? Cn ] = 2 [W (Cn+1 , Θn ) ? W

? (Cn+1 , Θn+1 ) ? W ? (Cn , Θn+1 )] [W +1 2 +pn+1/2 [G(Cn+1 ) ? G(Cn )]. (38)

13

? C and dG ? we obtain From (36)1 and the de?nitions of dW ? Cn+1 ) + λn+1/2 dG ? ? pn+1/2 ][Θn+1 ? Θn ] ? Cn + dW ( dW [1 2
1 ? ? (Cn , Θn )] + 1 [W ? (Cn+1 , Θn+1 ) ? W ? (Cn+1 , Θn )] =2 [W (Cn , Θn+1) ? W 2

+λn+1/2 [Θn+1 ? Θn ] ? pn+1/2 [Θn+1 ? Θn ]. (39) Subtracting (35) from (34), and adding (36) then leads to the result ? n+1 ? H ?n 0=H
1 +2 ? [ζn+1

? Θn+1 ][pn+1 ? pn ] d? ? +1 2
?

1 [ζ 2 ? n

? Θn ][pn+1 ? pn ] d?
1 2 ?

(40)

? n+1 [λn+1 ? λn ] d? ? G

? n [λn+1 ? λn ] d?. G

? n+1 ? H ? n = 0, Since (pn+1 ? pn ) ∈ E and (λn+1 ? λn ) ∈ E for any n, we obtain H which is the desired result. The above result shows that the algorithmic equations (32) inherit discrete ? . In particular, forms of the conservation laws associated with L, J and H ? and each of the components of L and J are integrals for the case of a H free body. As for the compressible case, an important feature of (32) is the use of an interpolated stretch ?eld Cn+1/2 , rather than the stretch ?eld of an interpolated con?guration C(?n+1/2 ), in evaluating the stress ?eld. The original incompressible formulation (2) and the quasi-incompressible formulation (28) may be interpreted as di?erential-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(?)) = 0 in ?, and a basic compatibility condition on (?, V) is that (?, V) ∈ M. A condition on λ may be deduced by considering the second time-derivative of the constraint G(C(?)) = 0. Finally, since smooth solutions of (28) satisfy (2) with Θ = ζ (?) and p = λ, we see that initial conditions on Θ and p follow from those on ? and λ.

4

Finite Element Spatial Discretization

Here we discuss ?nite element spatial discretization of the time integration schemes (10) and (32). Since the ?nite element discretization of the compressible case (10) is straightforward, we concentrate on the quasi-incompressible 14

case (32). In particular, we introduce a general mixed ?nite element discretization, reduce it to a two-?eld formulation and discuss an e?cient 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 ? ? R3 into nonoverlapping subdomains ?e ?e (e = 1, . . . , melem ), each de?ned by a set of element nodes Xe a ∈ ? (a = 1, . . . , men ), and for each e = 1, . . . , melem consider the sets
e ν e = {a | a = 1, . . . , men } and νg = {a ∈ ν e | X e a ∈ Γd }.

(41)

In accordance with a standard isoparametric discretization, we parametrize each subdomain ?e with a mapping χe : → ?e of the form
men

χe (ξ ) =
a=1

N a (ξ )X e a,

where ? R3 is typically the bi-unit cube and N a : → R (a = 1, . . . , men ) are standard isoparametric shape functions. In particular, the shape functions a , where ξb ∈ ? (b = 1, . . . , men ) are element satisfy the condition N a (ξb ) = δb nodes for the parent domain . Given a discretization of ? as described above, we introduce an approximation Uh of U, and an approximation Vh of V, de?ned as Uh = {?h : ? → R3 | ?h ∈ C 0 (?, R3 ), ?h |?e ? χe = η h |?e ? χe =
e a∈ν e \νg

N a de a+

e a∈νg

e 3 N a g (X e a ), da ∈ R }

Vh = {η h : ? → R3 | η h ∈ C 0 (?, R3 ),
e a∈ν e \νg

(42)

e 3 N a ce a , ca ∈ R },

where g : Γd → R3 is a speci?ed deformation on Γd . Similarly, we introduce an approximation Eh of E, namely
mdis

E = { φ : ? → R | φ | ?e ? χ =
i=1

h

h

h

e

? i α e , α e ∈ R} N i i

(43)

? i : → R (i = 1, . . . , mdis ) are smooth interpolation functions. Note where N that the functions in Eh are smooth within element domains, but are allowed to be discontinuous across element boundaries. 15

A mixed ?nite element discretization of the time integration scheme (32) is obtained by restriction to the ?nite-dimensional spaces Uh , Vh and Eh . 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 ?nite element formulation of (32) possesses conservation laws analogous to those of (32) itself.

4.2 Reduction to a Two-Field Formulation
h h h Here we eliminate the ?elds Vn +1 , Θn+1 and pn+1 from the ?nite element formulation of (32). The result will be a two-?eld formulation in which we h need only solve for ?h n+1 and λn+1 at each step of the algorithm. h The elimination of Vn +1 is straightforward in view of (32)1 . In particular, using this equation together with the de?nition of Vn+1/2 leads to

h Vn +1 =

2 h h (?h n+1 ? ?n ) ? Vn . ?t

(44)

h The elimination of Θh n+1 makes crucial use of the structure of E . Assuming ?i : the interpolation functions N → R (i = 1, . . . , mdis ) are non-trivial, we introduce for each e = 1, . . . , melem a symmetric, positive-de?nite matrix Ξe ∈ Rmdis ×mdis with components

(Ξe )ij =

? iN ? j det[?χe ] d . N

(45)

1 h We denote the components of the inverse matrix by (Ξ? e )ij . Since E is spanned by functions allowed to be discontinuous across element boundaries, we can use standard arguments and reduce the ?nite element version of (32)4 to e h Θh n+1 |?e φ |?e d?

?e

=

?e

h ζn +1 |?e

φ |

h

?e

d?

e

?φ |

h

(46)
?e

(e = 1, . . . , melem )

h h where ζn +1 = det[Fn+1 ]. For any e = 1, . . . , melem we then arrive at the equation e ?i e (Θh n+1 |?e ? χ )N det[?χ ] d

=

e ?i e h (ζ n +1 |?e ? χ )N det[?χ ] d

(47)

which holds for i = 1, . . . , mdis . To solve this equation for Θh n+1 |?e we consider 16

an interpolation of the form
mdis e Θh n+1 |?e ? χ = j =1

? j Θe N j,n+1

(48)

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

(Ξe )ij Θe j,n+1 =
j =1

e ?i e h (ζ n +1 |?e ? χ )N det[?χ ] d

(i = 1, . . . , mdis ). (49)

Solving this equation for Θe j,n+1 (j = 1, . . . , mdis ) and substituting the result into (48) then yields
mdis e Θh n+1 |?e ? χ = i,j =1 h which expresses Θh n+1 |?e as a function of ?n+1 |?e . h The elimination of ph n+1 proceeds along similar lines as the elimination of Θn+1 . Working with the ?nite element version of (32)3 , we ?nd that the restriction e of ph n+1 to any element ? is given by mdis e e h ph n+1 |?e ? χ = ?pn |?e ? χ + 2 i,j =1 1 2

? j (Ξ?1 )ji N e

e ?i e h (ζ n +1 |?e ? χ )N det[?χ ] d

(50)

? ? ?

? j (Ξ?1 )ji N e ? dG ? χ N det[?χ ] d
e

? C h + dW ? Ch dW + n n+1

λh n+1/2

?i

e

?e

? ? ?

(51) .

Our two-?eld, mixed ?nite element discretization of (32) can be stated as h h h h h h h h follows. Given (?h n , Vn , λn , Θn , pn ) ?nd (?n+1 , λn+1 ) ∈ U × E such that
2ρ0 h ? ?t2 (?n+1

?

?h n)

· η d? =

h

2ρ0 h ? ?t V n

· η d?

h

?

?h ? [S +

h h T h + 2 ph n+1/2 dG ] : (Fn+1/2 ) ?η d? h h ? bn+1/2 · η d? +

0=

h ?h ? φ Gn+1 d?

h h ? ? ?φh ∈ Eh ? Γσ fn+1/2 · η dΓ ? ?

? ? ? ? ? ? ? ? ? ?? ? ? ?η h ?? ? ? ? ?

∈ Vh

(52)

h h h where Θh n+1 and pn+1 are considered as functions of ?n+1 and λn+1 via expresh sions (50) and (51), respectively, and Vn +1 is computed according to (44).

The above algorithm is completed by an appropriate choice of interpolation functions for Uh , Vh and Eh . 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 17

purposes of illustration, we note that piecewise linear interpolation functions could serve as a basis for Uh and Vh . Moreover, the constant interpolation ? (ξ) ≡ 1 (mdis = 1) could serve as a basis for Eh . In this case we function N would have Θh n+1 |?e = 1 vol(?e )
?e e ζ (?h n+1 ) d? ,

(53)

and the restriction of Θh n+1 to any element would be constant and equal to the average jacobian value over ?e . Thus the incompressibility condition h e ? (Θh G n+1 ) = Θn+1 ? 1 = 0 in ? (e = 1, . . . , melem ) 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 dev (Cdev ) + U (det[C]1/2 ) (54)

where W dev (Cdev ) 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 Cdev is the deviatoric part of C de?ned by Cdev = det[C]?1/3 C, in particular ? (C, Θ) det[Cdev ] = 1. In view of (54), the modi?ed stored energy function W de?ned in (27) becomes ? (C, Θ) = W dev (Cdev ) + U (Θ). W (55)

For our quasi-incompressible formulation we can assume without loss of generality that the function U is of the form U (Θ) = Kγ (Θ). (56)

Here K > 0 is a constant and γ (Θ) is a convex, non-negative function that satis?es γ (Θ) = 0 if and only if Θ = 1, for example
1 γ (Θ) = 2 (Θ2 ? 1) ? ln(Θ).

(57)

This assumption can be made since Θ is constrained by the condition Θ = 1. Thus, we interpret U (Θ) 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 18

K and a multiplier λh n+1 . There are various e?cient 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:
h,k h h h h h (1) Let (?h n , Vn , λn , Θn , pn ) be given and set λn+1 |k =0 = λn . h,k h (2) For ?xed λh,k n+1 , solve (52)1 for ?n+1 , and call this solution ?n+1 .

(3) If (52)2 is not satis?ed within a speci?ed tolerance, then update the multiplier ?eld using the rule
+1 h λh,k n+1 φ d? = h λh,k n+1 φ d? + h h h ? (Θh,k KG n+1 )φ d? ?φ ∈ E , (58)

?

?

?

and repeat Step 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.

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 ?nite element discretization, and in the quasi-incompressible (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 the Ogden type [32,33], namely
3

W ( C) =
A=1

w (λ A )

(59)

where
k

w (λ A ) =
m=1

?m αm (λ ? 1) ? ?m ln(λA ) . αm A

(60)

Here C is the right Cauchy-Green stretch 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. Models of this form are particularly useful for rubber-like materials [32,33,47]. For our compressible calculations 19

we use k = 3 and the parameters ?1 = ?2 = 0.690 × 10 N/m , α1 = 0.010 × 106 N/m2 , α2 =
6 2

?3 = ?0.012 × 106 N/m2 , α3 = which are taken from [33, p.498].

? ? ? ? ?2.0, ?

4.0, ?

? ? ? 1.3, ? ? ? ?

(61)

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

W (C) = W dev (Cdev ) =
A=1

w dev (λdev A )

(62)

together with the constraint ζ = λ1 λ2 λ3 = 1, where w dev is of the form w dev (λdev A ) = ?m dev αm [(λA ) ? 1]. m=1 αm
k

(63)

2 Here Cdev is the deviatoric part of C and (λdev A ) (A = 1, 2, 3) are the eigendev ? (C, Θ) values of C . For this model the modi?ed stored energy function W becomes

? (C, Θ) = W dev (Cdev ) + U (Θ) W where U is interpreted as a penalty function, which we take as U (Θ) = Kγ (Θ).

(64)

(65)

Here K > 0 is a penalty parameter and γ 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 men = 8 and use standard trilinear interpolation functions as the basis for the spaces Uh and Vh . For the quasi-incompressible case, we take mdis = 1 and use a constant interpolation function as the basis for the space Eh . Thus, the functions in Uh and Vh are piecewise linear and continuous, while those in Eh are piecewise constant and possibly discontinuous across element boundaries.

5.1 Tumbling Elastic Block Consider a cube of dimension l = 0.02m as shown in Figure 1. The cube is composed of a homogeneous, elastic material of Ogden type with parameters as given in (61) and density ρ0 = 1000kg/m3. We suppose the cube is initially at 20

rest and subject to tractions f1 and f2 that vanish after a short period of time. The traction f1 is a 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 uniform force distribution on the face X3 = l/2 de?ned by f2 = p(t)(0, f, f /2). Here f = 32N/cm2 and p(t) is given by p(t) =
? ? ? t,

0 ≤ t ≤ 0.005s t > 0.005s.

(66)

? ? 0,
F2

X3 X2 X1 F1

l

Fig. 1. Schematic of elastic cube showing applied loads.

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 ? (Θe )| < 5 × 10?10 for each time step. above and iterate until maxe |G n Figure 2 shows snapshots, at time increments of 0.005s, of the cube motion in the time interval [0, 0.1s]. The data in the ?gure corresponds to the compressible motion, however, the quasi-incompressible motion is visually identical. For clarity, we show the snapshots in the X2 -X3 plane only.

Fig. 2. Cube motion (left-to-right) in X2 -X3 plane at time increments of 0.005s.

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

experiences numerical blow-up for the time step ?t = 0.005s in the time interval of the computation.
0.02 0.0002 0.00015 0.015 0.0001 5e-05 0.01 dt=.005 0 -5e-05 0.005 dt=.0025 0 0 0.04 0.08 t 0.12 0.16 0.2 -0.0001 -0.00015 -0.0002 0 0.04 0.08 t 0.12 0.16 0.2 H J J2=J3=0 J1, dt=.005 J1, dt=.0025

Fig. 3. Energy and angular momentum of compressible cube for mid-point rule.
0.01 0.0002 0.00015 0.0075 0.0001 5e-05 0.005 0 -5e-05 0.0025 dt=.005 -0.00015 dt=.0025 0 0 0.04 0.08 t 0.12 0.16 0.2 -0.0002 0 0.04 0.08 t 0.12 0.16 0.2 -0.0001 H J J2=J3=0 J1, dt=.0025 J1, dt=.005

Fig. 4. Energy and angular momentum of compressible cube for conserving scheme.

Figure 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.0100s, 0.0125s]. For each scheme, we supplied the same initial conditions at t0 = 0.0100s and computed the solution up to a time t1 = 0.0125s. Denoting the solution at time t1 for a given scheme and time step by ?t1 ,?t , we de?ne the L2 -error in displacements as L2 Error = |?t1 ,?t ? ?t1 ,conv |2
1/2

?

(67)

where ?t1 ,conv is the solution at t1 computed using a time step of ?t = 2.5 × 10?7 s. As shown in Figure 5, the conserving scheme and the mid-point rule have nearly identical accuracy properties for this model problem. Figure 6 shows total energy and angular momentum time histories, computed using a time step ?t = 0.005s, for the conserving scheme applied to the quasiincompressible formulation. 22

0.001

0.0001

L2 Error

1e-05

1e-06

1e-07

Mid-Point Rule Conserving Scheme

1e-08 1e-06

1e-05

0.0001 Time Step

0.001

0.01

Fig. 5. L2 -error in displacements of compressible cube versus time step.

0.01 0.008 0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 0 0.02 0.04 t 0.06 0.08 0.1 H J

0 J1 -0.0001

-0.0002

-0.0003

-0.0004 J2=J3=0 -0.0005 0 0.02 0.04 t 0.06 0.08 0.1

Fig. 6. Energy and angular momentum of quasi-incompressible cube for conserving scheme with ?t = 0.005s.

5.2 Stretched Elastic Plate

Consider a thick rectangular plate of thickness w = 0.002m, length l = 0.01m and height h = 0.01m as shown in Figure 7. The plate is composed of a homogeneous, incompressible, elastic material of Ogden type with parameters as given in (61) and density ρ0 = 1000kg/m3. In this example, we suppose the plate is subject to a zero displacement boundary condition along the shaded region shown in the ?gure, is initially at rest and is subject to a traction f that vanishes after a short period of time. The traction f is applied to the face X3 = ?h/2, is uniform across the thickness, and varies in the X2 -direction according to
? ? ? p(t)(f /2, 0, ?f ), ? ? p(t)(0, 0, ?f ),
l ?2 ≤ X2 < 0

f (X2 , t) =

0 ≤ X2 ≤

l 2

(68)

23

where f = 5kN/cm2 and p(t) is given by p(t) =
? ? ? t,

0 ≤ t ≤ 0.004s t > 0.004s.

(69)

The zero displacement boundary condition is imposed along the surface de?ned by {X3 = h/2} ∩ {?l/8 ≤ X2 ≤ l/2}.

? ? 0,

l w

X3 h X2

X1

F

Fig. 7. Schematic of elastic plate showing applied load and zero displacement boundary condition.

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 8 shows snapshots of the plate motion in the interval [0, 0.0045s]. This motion was computed with the quasi-incompressible conserving scheme with a step size ?t = 0.0001s. We employed the Uzawa ? (Θe )| < 5 × 10?10 for each time iteration strategy and iterated until maxe |G n step.

Fig. 8. Motion of quasi-incompressible elastic plate computed with conserving scheme. The direction of time is left-to-right, top-to-bottom.

24

Figure 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.0040s, but the total angular momentum is not an integral due to the displacement boundary conditions. These features are clearly resolved by the conserving scheme.
0.01 0.008 0.006 0.004 0.002 0 -0.002 -0.004 -0.006 -0.008 -0.01 0 0.001 0.002 0.003 t 0.004 0.005 H J 1e-06 8e-07 6e-07 4e-07 2e-07 J1 0 -2e-07 -4e-07 -6e-07 -8e-07 -1e-06 0 0.001 0.002 0.003 t 0.004 0.005 J2 J3

Fig. 9. Total energy and angular momentum of quasi-incompressible plate for conserving scheme with ?t = 0.0001s.

5.3 Inversion of a Spherical Cap Consider a thick spherical cap with dimensions r = 0.1m, β = 0.008m and α = π/4 as shown in Figure 10. The cap is composed of a homogeneous, compressible, elastic material of Ogden type with parameters as given in (61) and density ρ0 = 100kg/m3 . We suppose the cap is initially at rest and subject to a system of equilibrated force distributions f1 and f2 that vanish after a short period of time. The force distribution f1 is applied to the convex face of the cap, acts in the negative X1 direction, is uniform in a disk centered about the X1 axis and has a resultant given by p(t)(?f, 0, 0) where f = 160kN and p(t) = ?
? ? ? t,

0 ≤ t ≤ 0.001s t > 0.001s.

(70)

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 this example we compare two time integration schemes: the implicit midpoint rule and the energy-momentum preserving scheme (10). In the spatial discretization of the cap we use 4 elements through the thickness and a quasiuniform discretization of each constant-radius surface into 128 elements, for a total of 512 elements. Figure 11 show snapshots of the cap motion at time increments of 0.0005s computed with the conserving scheme. 25

? 0,

β

F2 X3 α

F1

r

X1

F2

Fig. 10. Schematic of elastic cap showing geometry and applied loads.

Fig. 11. Snapshots of cap motion computed with conserving scheme. The direction of time is left-to-right, top-to-bottom.

Figures 12 and 13 show total energy and angular momentum time histories computed using a time step ?t = 0.0005s. 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.

6

Closing Remarks

In this article we have considered energy and momentum conserving time discretization schemes for general initial boundary-value problems in ?nitedeformation elastodynamics. Conserving schemes were developed in a weak form for both compressible and incompressible hyperelastic material models, implemented using ?nite element discretizations in space and applied to three example problems. Compared to the implicit mid-point rule, the conserving 26

7 6 5 4

1e-11

5e-12

J1=J2=J3=0 0 J

H

3 2 1 0 0 0.005 t 0.01 0.015 -1e-11 0 0.005 t 0.01 0.015

-5e-12

Fig. 12. Energy and angular momentum of cap for conserving scheme with ?t = 0.0005s.
7 6 5 4 0 H 3 2 1 0 0 0.005 t 0.01 0.015 -1e-11 0 0.005 t 0.01 0.015 J 5e-12 1e-11

J1, J2, J3

-5e-12

Fig. 13. Energy and angular momentum of cap for mid-point rule with ?t = 0.0005s.

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 e?cient 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 27

acknowledges the ?nancial support of the US National Science Foundation.

References
[1] U.M. Ascher, H. Chin & S. Reich (1994) “Stabilization of DAEs and Invariant Manifolds,” Numerische Mathematik, 67, 131–149. [2] D.P. Bertsekas (1982) Constraint Optimization and Multiplier Methods, Academic Press. [3] K.E. Brenan, S.L. Campbell & L.R. Petzold (1989) Numerical Solution of Initial-Value Problems in Di?erential-Algebraic Equations, Elsevier Science. [4] F. Brezzi & M. Fortin (1991) Mixed and Hybrid Finite Element Methods, Springer series in computational mathematics, v. 15, Springer-Verlag. [5] C. Chen & W. von Wahl (1982) “Das Rand-Anfangswertproblem f¨ ur Quasilineare Wellengleichungen in Sobolevr¨ aumen Niedriger Ordnung,” Journal f¨ ur die reine und angewandte Mathematik, 337, 77–112. [6] P.G. Ciarlet (1988) Mathematical Elasticity, Volume I: Three-Dimensional Elasticity, Elsevier Science. [7] C.M. Dafermos & W.J. Hrusa (1985) “Energy Methods for Quasilinear Hyperbolic Initial-Boundary Value Problems. Applications to Elastodynamics,” Archive for Rational Mechanics and Analysis, 87, 267–292. [8] D.J. Dichmann & J.H. Maddocks (1996) “An Impetus-Striction Simulation of the Dynamics of an Elastica,” Journal of Nonlinear Science, 6, 271–292. [9] D.G. Ebin & R.A. Saxton (1986) “The Initial-Value Problem for Elastodynamics of Incompressible Bodies,” Archive for Rational Mechanics and Analysis, 94, 15–38. [10] D.G. Ebin & S.R. Simanca (1990) “Small Deformations of Incompressible Bodies with Free Boundaries,” Communications on Partial Di?erential Equations, 15, 1589–1616. [11] D.G. Ebin & S.R. Simanca (1992) “Deformations of Incompressible Bodies with Free Boundaries,” Archive for Rational Mechanics and Analysis, 120, 61–97. [12] D.G. Ebin (1993) “Global Solutions of the Equations of Elastodynamics of Incompressible Neo-Hookean Materials,” Proceedings of the National Academy of Sciences, USA, 90, 3802–3805. [13] G. Fichera (1972) Existence Theorems in Elasticity, Handbuch der Physik, Volume VIa/2, Springer-Verlag. [14] M. Fortin & A. Fortin (1985) “A generalization of Uzawa’s Algorithm for the Solution of the Navier-Stokes Equations,” Communications in Applied Numerical Methods, 1, 205–208.

28

[15] C. F¨ uher & B. Leimkuhler (1991) “Numerical Solution of Di?erential-Algebraic Equations for Constrained Mechanical Motion,” Numerische Mathematik, 59, 55–69. [16] C.W. Gear, B. Leimkuhler & G.K. Gupta (1985) “Automatic Integration of Euler-Lagrange Equations with Constraints,” Journal of Computational and Applied Mathematics, 12 & 13, 77–90. [17] R. Glowinski & P. Le Tallec (1981) “Finite Elements in Nonlinear Incompressible Elasticity,” Finite Elements, Volume V: Special Problems in Solid Mechanics, edited by J.T. Oden & G.F. Carey, Prentice Hall. [18] R. Glowinski & P. Le Tallec (1989) Augmented Lagrangian and OperatorSplitting Methods in Nonlinear Mechanics, SIAM Studies in Applied Mathematics. [19] O. Gonzalez (1996) “Time Integration and Discrete Hamiltonian Systems,” Journal of Nonlinear Science, 6, 449–467. [20] O. Gonzalez & J.C. Simo (1996) “On the Stability of Symplectic and EnergyMomentum Algorithms for Nonlinear Hamiltonian Systems with Symmetry,” Computer Methods in Applied Mechanics and Engineering, 134, 197–222. [21] O. Gonzalez “Mechanical Systems Subject to Holonomic Constraints: Di?erential-Algebraic Formulations and Conservative Integration,” Physica D, accepted. [22] O. Gonzalez, J.H. Maddocks & R.L. Pego, “Multi-Multiplier Ambient-Space Formulations of Constrained Dynamical Systems: The Case of Linearized Incompressible Elastodynamics,” Archive for Rational Mechanics and Analysis, submitted. [23] M.E. Gurtin (1981) An Introduction to Continuum Mechanics, Academic Press. [24] T.J.R. Hughes, T. Kato & J.E. Marsden (1977) “Well-posed Quasilinear Second-order Hyperbolic Systems with Applications to Nonlinear Elastodynamics and General Relativity,” Archive for Rational Mechanics and Analysis, 63, 273–294. [25] T.J.R. Hughes (1987) The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice Hall. [26] W.J. Hrusa & M. Renardy (1988) “An Existence Theorem for the Dirichlet Problem in the Elastodynamics of Incompressible Materials,” Archive for Rational Mechanics and Analysis, 102, 95–117. [27] T. Kato (1976) “Quasi-linear Equations of Evolution with Applications to Partial Di?erential Equations,” in Hyperbolicity, edited by G. DaPrato & G. Geymonat, Centro Internazionale Matematico Estivo, II ciclo, Cortona, 125– 191. [28] B.J. Leimkuhler & S. Reich (1994) “Symplectic Integration of Constrained Hamiltonian Systems,” Mathematics of Computation, 63, 589–605.

29

[29] D.G. Luenberger (1984) Linear and Nonlinear Programming, Addison-Wesley. [30] J.H. Maddocks & R.L. Pego (1995) “An Unconstrained Hamiltonian Formulation for Incompressible Fluid Flow,” Communications in Mathematical Physics, 170, 207–217. [31] J.E. Marsden & T.J.R. Hughes (1983) Mathematical Foundations of Elasticity, Dover. [32] R.W. Ogden (1982) “Elastic Deformation of Rubberlike Solids,” Mechanics of Solids, 499–537. [33] R.W. Ogden (1984) Non-Linear Elastic Deformations, Ellis Horwood. [34] W.C. Rheinboldt (1984) “Di?erential-Algebraic Systems as Di?erential Equations on Manifolds,” Mathematics of Computation, 43, 473–482. [35] R.D. Richtmyer & K.W. Morton (1967) Di?erence Methods for Initial-Value Problems, Interscience. [36] M. Sabl? e-Tougeron (1988) “Existence pour un probl` eme de l’? elasto-dynamique Neumann non lin? eaire en dimension 2,” Archive for Rational Mechanics and Analysis, 101, 261–292. [37] S. Schochet (1985) “The Incompressible Limit in Nonlinear Elasticity,” Communications in Mathematical Physics, 102, 207–215. [38] J.C. Simo, R.L. Taylor & K.S. Pister (1985) “Variational and Projection Methods for the Volume Constraint in Finite Deformation Elasto-Plasticity,” Computer Methods in Applied Mechanics and Engineering, 51, 177–208. [39] J.C. Simo, J.E. Marsden & P.S. Krishnaprasad (1988) “The Hamiltonian Structure of Nonlinear Elasticity: The Material and Convective Representation of Solids, Rods, and Plates,” Archives for Rational Mechanics and Analysis, 104, 125–183. [40] J.C. Simo & R.L. Taylor (1991) “Quasi-Incompressible Finite Elasticity in Principal Stretches. Continuum Basis and Numerical Algorithms,” Computer Methods in Applied Mechanics and Engineering, 85, 273–310. [41] J.C. Simo & K.K. Wong (1991) “Unconditionally Stable Algorithms for Rigid Body Dynamics that Exactly Preserve Energy and Angular Momentum,” International Journal for Numerical Methods in Engineering, 31, 19–52. [42] J.C. Simo, N. Tarnow & K.K. Wong (1992) “Exact Energy-Momentum Conserving Algorithms and Symplectic Schemes for Nonlinear Dynamics,” Computer Methods in Applied Mechanics and Engineering, 1, 63–116. [43] J.C. Simo & N. Tarnow (1992) “The Discrete Energy-Momentum Method. Conserving Algorithms for Nonlinear Elastodynamics,” ZAMP, 43, 757–793. [44] J.C. Simo & N. Tarnow (1994) “A New Energy and Momentum Conserving Algorithm for the Non-Linear Dynamics of Shells,” International Journal for Numerical Methods in Engineering, 37, 2525–2550.

30

[45] J.C. Simo, N. Tarnow & M. Doblar? e (1995) “Non-linear Dynamics of Threedimensional Rods: Exact Energy and Momentum Conserving Algorithms,” International Journal for Numerical Methods in Engineering, 38, 1431–1474. [46] C.A. Truesdell & W.M. Noll (1965) The Nonlinear Field Theories of Mechanics, Handbuch der Physik, Volume III/3, Springer-Verlag. [47] K.C. Valanis & R.F. Landel (1967) “The Strain-Energy Function of a Hyperelastic Material in Terms of the Extension Ratios,” Journal of Applied Physics, 38, 2997–3002.

31

### 5.4 Conservation of Momentum in Two Dimensions

5.4 Conservation of Momentum in Two Dimensions_法律资料_人文社科_专业资料。SPH 4U 5.4 ? Conservation of Momentum in Two Dimensions Page 1 of 2 Recall: ...

### Conservation of Momentum

Conservation of Momentum Introduction Momentum is one of the central quantities in physics, and less intuitive to many people than the conservation of energy...