δv is given as,

ˆ

DΠL (φ, p)[δv] = DΠ(φ)[δv] + pDJ[δv] dV = 0 (6.31)

V

Substituting Equation (6.29) and using Equation (3.76) for DJ[δv] and

(5.47) for the derivative of the second Piola“Kirchho¬ stress gives, after

12 LINEARIZED EQUILIBRIUM EQUATIONS

some algebra, the principle of virtual work as,

ˆ

‚Ψ

DΠL (φ, p)[δv] = : DC[δv] dV + Jp · δv dV ’ δWext (φ, δv)

‚C

V V

(6.32a)

= S : DE[δv] dV + p · δv dv ’ δWext (φ, δv)

V v

(6.32b)

= σ : δd dv + pI : δd dv ’ δWext (φ, δv) (6.32c)

v v

= σ : δd dv ’ δWext (φ, δv) = 0 (6.32d)

v

where Equation (6.11b) for the relevant push forward operation and Equa-

tions (4.49) and (4.50b) have been invoked. Observe that Equation (6.32c)

clearly identi¬es the Lagrange multiplier p as the physical internal (hydro-

static) pressure.

In the context of a Newton“Raphson solution process, the governing

Equations (6.30) and (6.32) need to be linearized with respect to both vari-

ables p and φ in the direction of respective increments u and ∆p. Starting

with Equation (6.30) and using the linearization of J as given in Equation

(3.76) we have,

D2 ΠL (φ, p)[δp, ∆p] = 0 (6.33a)

D2 ΠL (φ, p)[δp, u] = δp · u dv (6.33b)

v

The linearization of Equation (6.32d) with respect to increments in p is

easily found using (6.32b) to give,

D2 ΠL (φ, p)[δv, ∆p] = ∆p · δv dv (6.34)

v

The obvious symmetry between Equations (6.34) and (6.33b) leads to a

symmetric tangent matrix upon discretization.

Finally, for the purpose of obtaining the derivative of (6.32d) in the

direction of an increment u in the motion, it is convenient to revert to a

material description and rewrite Equation (6.32d) as,

S = S + pJC ’1

DΠL (φ, p)[δv] = S : DE[δv] dV ’ δWext (φ, δv);

(6.35)

V

The linearization of this expression in the direction of u is obtained in the

13

6.6 VARIATIONAL METHODS AND INCOMPRESSIBILITY

same manner as in Section 6.3 to give,

D2 ΠL (φ, p)[δv, u] = DE[δv] : C : DE[u] dV

V

T

+ S : [( 0 δv] dV (6.36)

0 u)

V

where the tangent modulus is given in Section 5.5.2 as,

‚ ˆ

(S + pJC ’1 ) = C + C p

C=2 (6.37a)

‚C

‚S

ˆ

C=2 (6.37b)

‚C

’1 ’1

C p = pJ[C — C ’ 2I ] (6.37c)

Similarly, in the current con¬guration Equation (6.36) becomes,

σ : [( u)T

D2 ΠL (φ, p)[δv, u] = δd : c : µ dv + δv] dv (6.38)

v v

where the spatial tangent modulus is,

ˆ

c = J ’1 φ— [C];

ˆ ˆ

c = c +c p ; c p = p[I — I ’ 2i ] (6.39)

In conclusion the linearization of the governing Equations (6.30), (6.34), and

(6.38) is the combination of Equations (6.33), (6.34), and (6.38).

The above Lagrangian multiplier approach can be the basis for a success-

ful ¬nite element implementation provided that the interpolations of φ and

p are carefully chosen so as to avoid volumetric locking. It does, however,

su¬er from the limitation of the presence of additional pressure variables.

These two problem will be considered in the next two sections, and the

Lagrange multiplier approach will not be pursued further.

6.6.3 PENALTY METHODS FOR INCOMPRESSIBILITY

Penalty methods are an alternative to the Lagrangian multiplier approach

to incompressibility. There are two ways in which penalty methods can be

introduced in the formulation. A popular physically based approach, which

conveniently eliminates the pressure as an independent variable, is to con-

sider the material as being nearly incompressible, whereby a large value of

the bulk modulus e¬ectively prevents signi¬cant volumetric changes. A sec-

ond, less intuitive route is to perturb the Lagrangian functional given in

Equation (6.28) by the addition of a further “penalty” term that enables

the pressure to eventually be arti¬cially associated with the deformation,

thereby again eliminating the pressure variable. It will transpire that these

14 LINEARIZED EQUILIBRIUM EQUATIONS

two methods lead to identical equations. The perturbed Lagrangian func-

tional is,

12

ΠP (φ, p) = ΠL (φ, p) ’ p dV (6.40)

2κ

V

where κ is the penalty parameter and clearly ΠP ’ ΠL as κ ’ ∞.

The stationary condition of the Functional (6.40) with respect to p now

becomes,

p

DΠP (φ, p)[δp] = δp (J ’ 1) ’ dV = 0 (6.41)

κ

V

Consequently, enforcing this equation for all δp functions gives an equation

arti¬cially relating p and J as,

p = κ(J ’ 1) (6.42)

Referring to Section 5.5.3, Equation (5.60), it is clear that this equation

represents a nearly incompressible material with the penalty number κ as

the bulk modulus. The use of this equation in a ¬nite element context, ei-

ther directly as a nearly incompressible material or indirectly as a perturbed

Lagrangian method, will lead to a formulation involving only kinematic un-

known variables.

The stationary condition of Functional (6.40) with respect to the motion

is identical to Equation (6.32) and gives the principle of virtual work, where

now the internal pressure in the Cauchy stresses can be evaluated using

Equation (6.42) or, in general for a nearly incompressible material, using

Equation (5.59). In conclusion, the linearized equilibrium equation is given

directly by Equation (6.15), where, as shown in Section 5.5.3, the tangent

modulus is now,

ˆ

c = c +c p +c κ (6.43a)

2ˆ

ˆ = 4 ‚ Ψ = 2 ‚S

ˆ

’1

ˆ

c =J φ— [C]; (6.43b)

C

‚C‚C ‚C

c p = p[I — I ’ 2i ] (6.43c)

dp

cκ = J I — I = κJI — I (6.43d)

dJ

6.6.4 HU-WASHIZU VARIATIONAL PRINCIPLE FOR

INCOMPRESSIBILITY

Neither the Lagrange multiplier method nor the penalty method as presented

above will result in a simple and e¬cient ¬nite element formulation. This is

15

6.6 VARIATIONAL METHODS AND INCOMPRESSIBILITY

because the kinematically restricted motion implicit in the ¬nite element dis-

cretization is, in general, unable to distort while simultaneously meeting the

incompressibility requirement at each point of the body. This phenomenon

manifests itself as a catastrophic arti¬cial sti¬ening of the system known as

locking. One common practical solution to this problem involves the use of

di¬erent discretizations for p and φ in the Lagrange multiplier approach.

Another ad hoc solution, associated with the enforcement of incompressibil-

ity using nearly incompressible constitutive equations, is to ¬rst separate

the volumetric and distortional strain energy as,

ˆ

Π(φ) = Ψ (C) dV + U (J) dV ’ Πext (φ) (6.44)

V V

and then underintegrate the volumetric term containing U (J).

A formal framework for the solution of locking problems for nearly in-

compressible materials that has potential for further developments is pro-

vided by a functional that permits the use of independent kinematic descrip-

tions for the volumetric and distortional deformations. This can be achieved

by introducing a three-¬eld Hu-Washizu type of variational principle as,

¯ ˆ ¯ ¯

ΠHW (φ, J, p) = Ψ (C) dV + U (J) dV + p(J ’ J) dV ’ Πext (φ)

V V V

(6.45)

¯

where J is a new kinematic variable representing the dilatation or volume

change independently of the motion and p is a Lagrange multiplier enforc-

¯

ing the condition that J = J. The opportunity of using independent dis-

¯

cretizations for φ and J that this functional a¬ords will introduce su¬cient

¬‚exibility to prevent locking.

¯

The stationary conditions of Functional (6.45) with respect to φ, J, and

p will yield the equilibrium equation and the constitutive and kinematic

relationships associated with the volumetric behavior. To this end, we ¬nd

the directional derivative of ΠHW in the direction δv as,

ˆ

‚Ψ

¯

DΠHW (φ, J, p)[δv] = : DC[δv]dV + pDJ[δv] dV

‚C

V V

’ DΠext (φ)[δv] (6.46)

Repeating the algebra used in Equation (6.32), the stationary condition with

respect to φ gives the principle of virtual work as,

¯

DΠHW (φ, J, p)[δv] = σ : δd dv ’ δWext (φ, δv) = 0 (6.47)

v

The stationary condition of the Functional (6.45) with respect to changes

16 LINEARIZED EQUILIBRIUM EQUATIONS

¯

in J gives a constitutive equation for p as,

dU

¯ ¯ ¯

DΠHW (φ, J, p)[δ J] = ¯ ’ p δ J dV = 0 (6.48)

dJ

V

Similarly, the stationary condition of (6.45) with respect to p gives a kine-

¯

matic equation for J as,

¯ ¯

DΠHW (φ, J, p)[δp] = (J ’ J)δp dV = 0 (6.49)

V

¯

Obviously if δp and δ J are arbitrary functions, Equation (6.49) gives

¯

J = J and Equation (6.48) gives p = dU/dJ. However, as explained, a

¬nite element formulation based on this procedure where the discretization

¯

of J is based on the same interpolation as φ confers no advantage, and