- Email: [email protected]
Penalty function formulations for incompressible nonlinear elastostatics
COMPUTER METHODS IN APPLIED MECHANICS NORTH-HOLLAND PUBLISHING COMPANY
PENALTY
AND ENGINEERING
35 (1982) 107-l 18
FUNCTION FORMULATIONS FOR INCOMPRESSIBLE NONLINEAR ELASTOSTATICS J.C. SIMO and R.L. TAYLOR
Department of Civil Engineering, University of California, Berkeley, CA 94720, U.S.A. Received
1 December
1981
The analysis of incompressible nonlinear elastic solids is considered by a penalty function approach. In particular the necessary compatibility conditions for developing penalty forms for isotropic nonlinear elasticity are addressed, and the commonly used form of strain energy functional leading to a penalty formulation of the incompressibility constraint extended. The Mooney-Rivlin model is used to show how previous developments can lead to physically meaningless situations. For a certain class of incompressible materials, which includes the important Mooney-Rivlin model, a new and particularly simple formulation is proposed.
1. Introduction The penalty function method, originally proposed by Courant [l], provides a simple and effective procedure of reducing a constrained minimization problem, to one without constraints. The method has a simple geometrical interpretation when placed in the general context of optimization theory in Banach spaces [2]. Furthermore, in contrast with most of duality methods, it can be justified without invoking any convexity assumptions on the functional to be minimized [2,3]. This remarkable fact, makes the method particularly attractive in problems in nonlinear elastostatics, where the total potential energy functional is rarely convex.’ In the context of nonlinear elastostatics, the penalty function method provides a procedure of enforcing the incompressibility constraint, without restraining the configuration space to isochoric deformations. For this class of problems, the application of the method hinges on a suitable extension to the compressible range, of the constitutive model for the given incompressible material. The most widely used extension [3-61, assumes as strain energy in the compressible range, the sum of the strain energy potential of the incompressible material plus a penalty term enforcing the incompressibility constraint. The purpose of this paper is to discuss the form of this penalty term and, moreover, show that this form of strain energy is not the only possible one leading to a penalty formulation of the incompressibility constraint. For a certain class of incompressible materials, to which the Mooney-Rivlin model belongs, ‘As a matter of fact, convexity references therein.
00457825/82/0000-0000/$02.75
implies
0
uniqueness
and therefore
1982 North-Holland
precludes
buckling
[12]. See also [8] and
108
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
an alternative form of the strain energy potential in the compressible range is proposed, which leads to a simpler form of the elasticity tensor. This form is, therefore, particularly useful in the context of a finite element formulation. 2. Isotropic nonlinear elastostatics Consider a hyperelastic, isotropic body B, assumed to be identified with its reference configuration B C W3, a bounded open set with smooth boundary aB. Let @ : B - R3 be any (finite) deformation, and a& C lJB the part of the boundary where ip is pcescribed. Similarly B, C dB is that part of the boundary where the Kirchhoff stress vector PN = t is prescribed. The con~guration space C is then the set of diffeomo~hisms defined by C = {@ : B -+ R’ 1det(V@) > 0 and @jasd= g}
(1)
where a& n &B, = 0 and a& U B3, = dB. Denote by F = V@ the deformation gradient and let C = F’F be the right Cauchy-Green tensor, with principal invariant Ii (i = 1,2,3). If we let J = det(F) then p. and p = pa/J are the densities in the undefo~ed and deformed configurations respectively. Since the body is assumed to be hyperelastic and isotropic the strain energy functional W: C-+R exists and is given by @-+
W(G) = W(F) =
f%(C) = W(I1,L, 13)
and the first Piola-K~rchho~ and Cauchy stress tensors, denoted by P and CTrespectively computed from the strain energy W according to
are
P = c3bvKw
(31
CFfPF.
(4)
and
Since I, = tr(C), 12= ;{I: - tr(C’)), and 1, = det(C) = J2, application of the chain rule and the Cayley-Hamilton theorem gives for the Cauchy stress tensor the constitutive equation [9]
where B = RF’ is the left Cauchy-Green given by
tensor, and the functions fii(ll+ 12, 4) (i = 0, 1.2) are
However, the response functions Bj (i = 0, 1,2) are not independent. key compatibility conditions [ 101
They are related by the
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
109
(7)
which can be easily checked using (6). Denote by A the fourth order elasticity tensor, functional by
defined in terms of the strain energy
cw
a’W
A=zim=dF with components AiGJ, with respect to the standard basis {E,) and {e,) in the configurations B and 0(B), respectively. Let V be the space of kinematically admissible variations, that is, the tangent space to the con~guration space C. Thus V is defined by
V=(v:B-+R3~vJa&=0}.
(9)
With this notation the condition Vu - (A VW) = VW - (A Vu) for any v, w E V is the necessary and sufficient condition for the potential W: C-tR to exit [8-91. This condition is equivalent to the symmet~ of efasticity tensor; i.e. l
l
It can be shown that the symmetry condition (10) holds, provided the compatibility conditions given by (7) hold. Therefore (7) are just the integrability conditions for P = JaF-” to be derivabie from the potential W. 2.1. Incompressibility-
The Mooney-Rivlin
model
In the case of an incompressible
hyperefastic material, the configurations @: B + R3 are restricted by the constraint J = det(F) = 1. The constrained configuration space Ci”, C C is then defined by (11) and the constrained
tangent space of admissibfe variations by (12)
The strain energy functional, denoted now by I@ : C&c+ R, depends only on the first and second invariant I1 and &. Then instead of (3) and (5) the appropriate expressions for P and o”
110
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
are [9,10] P = -$F+ + aW/dF
(13)
@ = --$I + fi*B + &IF
(14)
and
where 0 : B--f R is the hydrostatic pressure, a bounded function to be determined by the boundary conditions. The response functions & (CX= 1,2) depend only on the invariants II and &. They are given now by
and the integrability conditions (7) for the response functions @a (a = 1,2) reduce to a&iarp= -a&all.
(16)
Mooney-Rivlin model corresponds to the assumption fil = a, and & = -az, where u1 and a2 are given constants independent of I, and I,, al > 0 and a2 > 0 ]lO], By redefining the hydrostatic pressure to be p = j3 + @I + &., (14) takes the form The
a=-(p+(u,-a&l+ci
where the extra stress 5 is given by
The strain energy potential for this model is I@ = &l(Il-
3) -I-a,(& - 3)) .
3. Penalty formulation for nonlinear elastostatics In elastostatics the variational formulation of the incompressibility constraint by a penalty function procedure involves the extension of the constitutive model for the incompressible material to the compressible range. In the linearized theory this extension is immediate. The constitutive equation ET= -PI+ 2,~e, where e is the deviatoric part of the infinitesima1 strain tensor E, is extended to the compressibIe range by considering instead u = K div(u) + 2p.e. This is the same as adopting for the strain energy cfi in the compressible range a form w = l&‘+ $K{div(u)l’, where k = p tr(e * e) is the distortional part. A penalty formulation makes use of this compressible model, with the bulk modulus K as a penalty parameter [ll]. In the nonlinear theory, however, the situation is quite different. A first formulation [3-61
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
111
considers a strain energy potential of the form W(&, &,4) = ti(&, 12)f U(&)
(19)
in the compressible range. @(I,, IJ is the strain energy of the given incompressible material, and the term U(&) has the structure of a penalty function enforcing the incompressibility constraint. A suitable form for this term will be discussed in Section 3.1. The type of strain energy functional given by (19) is not the only possible one leading to a penalty formulation. Alternative expressions to (19) more convenient for computational purposes, will be considered in Section 3.2. In the sequel the following convention will be adopted: variables with a superimposed b-3 will always be associated with the incompressible model, while those with a superimposed ‘-’ will be associated with the constitutive model extended to the compressible range. 3.1. GeneraE fury ulati~~ Let us assume a relationship between the %train energies @(l,, 1J and c?/(I1, I,, &) of the form (19). The associated response functions &(I1, IZ) and Ba(l,, 1,, 13) are related by /% = f&/J,
b = k&J,
m
which follow from the definitions (6) and (15) and from (19). Since (2/J)~~(~/~~~)= #aY, it is convenient for simplicity to regard U in (19) as a function of J rather than of 1,. A suitable form for the term U(&) in (19) is provided by U(I,) = $A{G(J)}’ + cH(J)
(21)
where G : (0, a)+lR is a penalty function with the property that G(J) = 0 if and only if J = 1 and A > 0 is the penalty parameter. The constant c and the function W(J) with the property H(J) = 0 and H’(J) = 1 iff J = 1, guarantee that the reference configuration @ = Identity is stress free. The usual way in which the condition is enforced [6,7], amounts to taking H(J)
= J.
From (21) and the definition of ,&(fl, &, J3) given by (6) it follows that
(22) The response functions /!&(i = 0, 1,2) defined by (20) and (22) satisfy the compatibility equations (7), and completely characterize the constitutive equation (5) for the Cauchy stress tensor. The constant c is given by (23) for the reference configuration to be stress free. The corresponding expression for the first Piola-Kirchhoff tensor follows from (4) and the elasticity tensor A can then be computed using (8).
112
J.C. Simo and R.L.
Taylor, Penalty function formulations
for incompressible
nonlinear elastostatics
In order to discuss appropriate forms for the functions G(J) and H(J) in (21) it is convenient to consider separately the effect of the penalty term JA{G(J)}~.Thus, we rewrite (19) in the form
I@(&,Iz, 13)= $h{G(J)}2
(244 WW
-I- W" ,
w* = cH(J)+ Iv.
3.1.1. The choice of the penalty term In the context of any numerical procedure, like the finite element method, the condition A 3 ~0can only be enforced in a relative sense. It is therefore appropriate to view (24) in physical terms as the strain energy functional of a compressible material, which exhibits incompressible behavior as the parameter A + ~0. Accordingly, (24) should remain physically meaningful for the widest possible range of values of A. In particular, situations in which the stresses take a finite value as the volume tends to zero, or become zero as the volume tends to infinite, not only are physically meaningless, but for highly confined nearly incompressible bodies may lead to numerical instabilities of an iterative solution strategy. The behavior of the strain energy given by (24) can be easily examined as A varies in (0, cc) by considering the simple problem of the extension of a cube with unit length sides. We shall take the stretching ratio of its side, as the control variable. Thus, any possible deformation @ E C has a deformation gradient F = yI
(25)
with determinant J = y3. The associated first Piola-Kirchhoff perform the extension) is a hydrostatic pressure2
tensor (the force required
to
(26)
p = PW.
For a given value of A E (0, co), the following conditions are expected to hold: (c.1) The strain energy VSI+ ~0, as the stretching ratio y + 0 or y + a. (c.2) p(y) is a monotone function in (0, ~0) with p(1) = 0, p--f ---co as y -+ 0, and p + -1-00 as y++m. We shall consider in the sequel, the special case of the Mooney-Rivlin model. For this type of material the strain energy and first Piola-Kirchhoff stresses for the extension problem are w =
ih(f17(~“))2+ w* ,
W” = -(2az + a#YI(y”) + t(al(y2 - 1) + a,(?” -
1)), (27)
P = {2a2(y3 - y’H’(y”)) + a,(7 - y2H’(y3))}1. (a) For A = 0 we have the limiting case of no penalty term and I@ = W*. The choice H(J) = In(J) ‘Taking the stresses as control nonlinear elasticity. Seven solutions
(28) variables leads to the classical example which shows for a simple neo-Hookean material are possible [O].
lack of uniqueness
in
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear
ehStOStatiCS
113
produces the strain energy @ and first Piola-Kirchhoff stresses P plotted in Fig. 1 and Fig. 2, respectively (curves A = 0). Both I&’and P satisfying conditions (c.1) and (c.2). (b) For A > 0 the response is affected by the penalty term. Two common forms of penalty function are G(J) = (J - 1) [3-61 and G(J) = In(J) 171. The first one, together with H(J) = J 161, produces the family of curves plotted in Fig. 3. Notice that an instability appears around J = 0.5, and that P vanish when y = 0. A similar instability phenomenon occurs for G(J) = In(J), although around the very large value of J = lo3 [7]. However, the combined function {G(J)}* = (In J)” + (J - l)*
(29)
is a suitable form of penalty function and produces the family of curves plotted in Fig. 1 and Fig. 2 for different values of A. For each A E (0, a), @ and P satisfy conditions (cl) and (c.2). The functions H(J) and G(J) given by (28) and (29) do not introduce any special computational effort. In fact, the corresponding response function B0 is b0=${A(J(J-1)+lnJ)+c-12B2}.
(30)
3.2. Alternative formulation Instead
of considering
a relationship
x=1000 X=lOO A=
of the type given by (19) between
strain energy
x=0
IO
5C
T
4c 3c 0k
20
2 2
IO
u z 2
0
E L
-10
E c
IO 0 -10
tL 0
: 1, 0.5
-20 -30 -40
, EXTENSION
1.5
2.0
2.5
3.0
RATIO
Fig. 1. Strain energy in isotropic extension. Proposed model.
-5o_ -1 ”
us
1.5
1.0
EXTENSION
Fig. 2. Piola-Kirchhoff model.
2.0
2.5
3.0
RATIO
hydrostatic pressure. Proposed
114
J.C. Simo and R.L. Taylor, Penalty function formulations
0
05
IO
I 5
EXTENSION
Fig. 3. Piola-Kirchhoff
for incompressible
2.0
25
nonlinear elastostatics
30
RATIO
hydrostatic pressure. Other models.
potentials, attention will now be focussed in the constitutive equations (5) and (14) for the Cauchy stresses in the compressible and incompressible cases, respectively. We shall consider first the case of the Mooney-Rivlin model. Extensions to more general incompressible models will be examined later. Let g : (0, =J)-+R be a smooth real valued function, such that g(x) = 0 if and only if x = 1. Consider first the replacement of the term --PI in (t7a) by hg(J), where A > 0 is a real constant. Now let A -+ ~0; i.e., consider a sequence {A,,) such that lim.,, A,, = 00. Since p : B --$ R is a bounded function (say in the uniform norm) and g(J) = p/A, it follows that if we let A --) ~0,then g(J)+ 0. By the assumptions on g(J), J + I as A --, CQ. The described replacement of hydrostatic pressure -PI, shows that the incompressible case can be viewed as the limit of a sequence of compressible cases. However, in contrast with the linearized theory, this fo~ulation requires the modification of the term 5 in (17a). In fact, if ci remains unchanged one is led to the constitutive model
and substitution
into the general integrability
ial = i3&/aI,= 0 ,
conditions (7) yields
$zz= -ap,/ax, =0 .
Thus ri = 0, and we arrive to the remarkable
conclusion that the only possible hyperelastic
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
1JS
compressible model, compatible with the assumption of constants response functions & (a = 1,2) is the hydrostatic pressure a = g(J)I. Since the response functions PI and p2 can not be a pair of constants in the compressible range, although they must reduce to the Mooney-Rivlin constants aI and a2 whenever .I = I, we are led to the assumption
Substitution into the integrability conditions equations which have the simple solution
therefore, given by
a consistent extrapolation
LT =
{hg(J)
of the Mooney-Rivlin
- (al - a2)h(J)}I
+ y
and
h(J)
differential
model to the compressible
range is
(34)
B - T B-l
h : (0, co)--) R
where for more generality the function been introduced. Ilf we take g(J) = 4 $ (G(J))’
(7) yields a pair of ordinary
with the property
h(x)
= 1
iff x = 1 has
= -& H(J),
where the functions G(J) and N(J) have the same meaning as in the previous section, (34) corresponds to a strain energy potential of the form VV(Il, 12, Q
= U(I3) + i(al(Il
- 3) + a2(I2I;1
- 3))
where U(&) is again given by (21) with the constant c = a, - ct2. This strain energy potential has a structure different from that considered in (19) and discussed earlier. In contrast with the previous formulation, the resulting first Pio~a-Kir~hho~ tensor does not contain the invariant 12, and 13 only appears in /%. Thus, it leads to a simpler expression of the elasticity tensor A, more convenient for computational purposes. The form of strain energy given by (35) also leads to a penalty formulation of the incompressibility constraint. This point is examined next. 3.2.1. Pevtalty ~~~~~~ai~u~ for gh’lzeM~u~~~-~~~~~~ rnuciei Let us rewrite (35) in a form analogous to (24a), with W*
= (a2 - a&Y(J)
The total potential
-
$aj(il - 3)
energy functional
+- $at(Zf;’
corresponding
W*
given now by
- 3). to the strain energy
W*
is, for any
116
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
configuration
@ E C, given by
(37) where x is the vector from the origin to x = Q(X), and 6 are the body forces per unit of mass in the reference configuration. With this notation the variational formulation of the boundary value problem for an elastic material with strain energy given by (36) is then find @ E C such that
I, (G(JN2 dV} .
(38)
For isochoric deformations @ E C&c, I3 = I and (36) reduces to the strain energy for the Mooney-Rivhn model. By the assumptions on G(d), the condition c;PE Cl,,, is equivalent to G(~(~)) = 0. Therefore, the variational formulation of the boundary value problem for the Mooney-Rivlin model is given by find @ E C’i,, such that
II*(@) = mix {n*(A)} s.t. G(J(A)) = 1 .
If we let A + m, (38) is just the penalty formulation [l-3]. 3.2.2. ~xte~s~o~ to more general co~stit~tive models Let us examine for what class of incompressible response functions of the form
of the constrained
(39)
problem given by (39)
materials a simple relationship
between
analogous to that found for the Mooney-Rivlin model is in general possible. Substitution of (40) and (41) into the integrabihty conditions (7) yields, after some manipulation, the system of equations
(42)
The solution of (42) together with (16) is given by
h(J) =Ql/J,
c;bz(J)= 1/P-1
)
(43)
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
117
together with to a strain energy potential of the form
where A and 23 are constants. For this class of constitutive models exactly the same formulation as that used for the Mooney-Rivlin model is applicable. Eq. (44) includes some forms of strain energy commonly used in rubber elasticity, like ti = 1”<&- 3) + (I* - 3) .
(45)
4. Final remarks The penafty fo~ulation of the incompressibiIity constrain in finite elasticity amounts to considering an extend constitutive model which, although no longer incompressible, exhibits incompressible behavior in the limit as the penalty parameter tends to infinity. It has been shown in this paper that consideration of the physical signi~cance of this extended model leads to a proper choice of the penalty term. In particular, the proposed form has proper behavior over all the range of volume changes. Consequently, numerical instability should not appear in an iterative solution strategy based upon Newton’s method. Most importantly, the proposed form of the penahy function does not involve any additional computational effort. Finally, for a certain class of materials a simpler extension to the compressible range was proposed as an alternative to the one commonly used [3-71. This class of materials includes some frequently used in rubber elasticity, in particular the important Mooney-Rivlin model. For this constitutive model the proposed penalty formulation results in substantial savings over previous proposals in computing a tangent stiffness matrix.
Acknowledgment
Support for this research was provided by the Malaysian Rubber Producers Research Association of Brichendonbury, England, as part of a program of research on the use of rubber bearings for the seismic protection of buildings. This support is gratefully acknowledged.
References [l] R. Courant, Calculus of Variations and Supplementary Notes on Exercises (New York University, New York, 1962). [2] G.D. Luenberger, Optimization by Vector Space Methods (Wiley, New York, 1969). [3] J.T. Oden, A theory of penalty methods for finite element approximations of highly nonlinear problems in continuum mechanics, Comput. and Structures 8 (1978) 445449.
118
J.C. Simo and R.L. Taylor, Penalty function formulations
for incompressible
nonlinear elastostatics
[4] J.T. Oden, Finite Element for Nonlinear Continua (McGraw-Hill, New York, 1972). [5] D. Malkus, Finite element with penalties in nonlinear elasticity, Internat. J. Numer. Meths. Engrg. 16 (1980) 121-136. [6] D. Malkus and E.R. Fuller Jr., An isoparametric finite element model for large strain elastostatics, J. Res. Nat. Bureau of Standards, 86(l) (1981). [7] T.J.R. Hughes and E. Carnoy, Nonlinear finite element shell formulation for large membrane strains, Civil Engineering Lab., Port Hueneme, CA, Rept. N68305-0357-4185 (1981). [8] J. Marsden and T.J.R. Hughes, Topics in the mathematical foundations of elasticity, in: R.J. Knops, ed., Nonlinear Analysis and Mechanics: Heriot-Watt Symposium, Vol. II, Research Notes in Mathematics, Vol. 27 (Pitman, Belfast, 1978). [9] M. Gurtin, Topics in Finite Elasticity, SIAM Vol. 35 (SIAM, New York, 1981). [lo] C. Truesdell and W. Nell, The non linear field theories of mechanics, in: S. Flugge, ed., Handbuch der Physik, Vol. III/33 (Springer, New York, 1965). [ll] 0. Zienkiewicz, The Finite Element Method (McGraw-Hill, New York, 3rd ed., 1977). [12] R. Hill, On uniqueness and stability in the theory of finite elastic strain, J. Mech. Phys. Solids 5 (1957) 229-241.
PENALTY
AND ENGINEERING
35 (1982) 107-l 18
FUNCTION FORMULATIONS FOR INCOMPRESSIBLE NONLINEAR ELASTOSTATICS J.C. SIMO and R.L. TAYLOR
Department of Civil Engineering, University of California, Berkeley, CA 94720, U.S.A. Received
1 December
1981
The analysis of incompressible nonlinear elastic solids is considered by a penalty function approach. In particular the necessary compatibility conditions for developing penalty forms for isotropic nonlinear elasticity are addressed, and the commonly used form of strain energy functional leading to a penalty formulation of the incompressibility constraint extended. The Mooney-Rivlin model is used to show how previous developments can lead to physically meaningless situations. For a certain class of incompressible materials, which includes the important Mooney-Rivlin model, a new and particularly simple formulation is proposed.
1. Introduction The penalty function method, originally proposed by Courant [l], provides a simple and effective procedure of reducing a constrained minimization problem, to one without constraints. The method has a simple geometrical interpretation when placed in the general context of optimization theory in Banach spaces [2]. Furthermore, in contrast with most of duality methods, it can be justified without invoking any convexity assumptions on the functional to be minimized [2,3]. This remarkable fact, makes the method particularly attractive in problems in nonlinear elastostatics, where the total potential energy functional is rarely convex.’ In the context of nonlinear elastostatics, the penalty function method provides a procedure of enforcing the incompressibility constraint, without restraining the configuration space to isochoric deformations. For this class of problems, the application of the method hinges on a suitable extension to the compressible range, of the constitutive model for the given incompressible material. The most widely used extension [3-61, assumes as strain energy in the compressible range, the sum of the strain energy potential of the incompressible material plus a penalty term enforcing the incompressibility constraint. The purpose of this paper is to discuss the form of this penalty term and, moreover, show that this form of strain energy is not the only possible one leading to a penalty formulation of the incompressibility constraint. For a certain class of incompressible materials, to which the Mooney-Rivlin model belongs, ‘As a matter of fact, convexity references therein.
00457825/82/0000-0000/$02.75
implies
0
uniqueness
and therefore
1982 North-Holland
precludes
buckling
[12]. See also [8] and
108
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
an alternative form of the strain energy potential in the compressible range is proposed, which leads to a simpler form of the elasticity tensor. This form is, therefore, particularly useful in the context of a finite element formulation. 2. Isotropic nonlinear elastostatics Consider a hyperelastic, isotropic body B, assumed to be identified with its reference configuration B C W3, a bounded open set with smooth boundary aB. Let @ : B - R3 be any (finite) deformation, and a& C lJB the part of the boundary where ip is pcescribed. Similarly B, C dB is that part of the boundary where the Kirchhoff stress vector PN = t is prescribed. The con~guration space C is then the set of diffeomo~hisms defined by C = {@ : B -+ R’ 1det(V@) > 0 and @jasd= g}
(1)
where a& n &B, = 0 and a& U B3, = dB. Denote by F = V@ the deformation gradient and let C = F’F be the right Cauchy-Green tensor, with principal invariant Ii (i = 1,2,3). If we let J = det(F) then p. and p = pa/J are the densities in the undefo~ed and deformed configurations respectively. Since the body is assumed to be hyperelastic and isotropic the strain energy functional W: C-+R exists and is given by @-+
W(G) = W(F) =
f%(C) = W(I1,L, 13)
and the first Piola-K~rchho~ and Cauchy stress tensors, denoted by P and CTrespectively computed from the strain energy W according to
are
P = c3bvKw
(31
CFfPF.
(4)
and
Since I, = tr(C), 12= ;{I: - tr(C’)), and 1, = det(C) = J2, application of the chain rule and the Cayley-Hamilton theorem gives for the Cauchy stress tensor the constitutive equation [9]
where B = RF’ is the left Cauchy-Green given by
tensor, and the functions fii(ll+ 12, 4) (i = 0, 1.2) are
However, the response functions Bj (i = 0, 1,2) are not independent. key compatibility conditions [ 101
They are related by the
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
109
(7)
which can be easily checked using (6). Denote by A the fourth order elasticity tensor, functional by
defined in terms of the strain energy
cw
a’W
A=zim=dF with components AiGJ, with respect to the standard basis {E,) and {e,) in the configurations B and 0(B), respectively. Let V be the space of kinematically admissible variations, that is, the tangent space to the con~guration space C. Thus V is defined by
V=(v:B-+R3~vJa&=0}.
(9)
With this notation the condition Vu - (A VW) = VW - (A Vu) for any v, w E V is the necessary and sufficient condition for the potential W: C-tR to exit [8-91. This condition is equivalent to the symmet~ of efasticity tensor; i.e. l
l
It can be shown that the symmetry condition (10) holds, provided the compatibility conditions given by (7) hold. Therefore (7) are just the integrability conditions for P = JaF-” to be derivabie from the potential W. 2.1. Incompressibility-
The Mooney-Rivlin
model
In the case of an incompressible
hyperefastic material, the configurations @: B + R3 are restricted by the constraint J = det(F) = 1. The constrained configuration space Ci”, C C is then defined by (11) and the constrained
tangent space of admissibfe variations by (12)
The strain energy functional, denoted now by I@ : C&c+ R, depends only on the first and second invariant I1 and &. Then instead of (3) and (5) the appropriate expressions for P and o”
110
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
are [9,10] P = -$F+ + aW/dF
(13)
@ = --$I + fi*B + &IF
(14)
and
where 0 : B--f R is the hydrostatic pressure, a bounded function to be determined by the boundary conditions. The response functions & (CX= 1,2) depend only on the invariants II and &. They are given now by
and the integrability conditions (7) for the response functions @a (a = 1,2) reduce to a&iarp= -a&all.
(16)
Mooney-Rivlin model corresponds to the assumption fil = a, and & = -az, where u1 and a2 are given constants independent of I, and I,, al > 0 and a2 > 0 ]lO], By redefining the hydrostatic pressure to be p = j3 + @I + &., (14) takes the form The
a=-(p+(u,-a&l+ci
where the extra stress 5 is given by
The strain energy potential for this model is I@ = &l(Il-
3) -I-a,(& - 3)) .
3. Penalty formulation for nonlinear elastostatics In elastostatics the variational formulation of the incompressibility constraint by a penalty function procedure involves the extension of the constitutive model for the incompressible material to the compressible range. In the linearized theory this extension is immediate. The constitutive equation ET= -PI+ 2,~e, where e is the deviatoric part of the infinitesima1 strain tensor E, is extended to the compressibIe range by considering instead u = K div(u) + 2p.e. This is the same as adopting for the strain energy cfi in the compressible range a form w = l&‘+ $K{div(u)l’, where k = p tr(e * e) is the distortional part. A penalty formulation makes use of this compressible model, with the bulk modulus K as a penalty parameter [ll]. In the nonlinear theory, however, the situation is quite different. A first formulation [3-61
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
111
considers a strain energy potential of the form W(&, &,4) = ti(&, 12)f U(&)
(19)
in the compressible range. @(I,, IJ is the strain energy of the given incompressible material, and the term U(&) has the structure of a penalty function enforcing the incompressibility constraint. A suitable form for this term will be discussed in Section 3.1. The type of strain energy functional given by (19) is not the only possible one leading to a penalty formulation. Alternative expressions to (19) more convenient for computational purposes, will be considered in Section 3.2. In the sequel the following convention will be adopted: variables with a superimposed b-3 will always be associated with the incompressible model, while those with a superimposed ‘-’ will be associated with the constitutive model extended to the compressible range. 3.1. GeneraE fury ulati~~ Let us assume a relationship between the %train energies @(l,, 1J and c?/(I1, I,, &) of the form (19). The associated response functions &(I1, IZ) and Ba(l,, 1,, 13) are related by /% = f&/J,
b = k&J,
m
which follow from the definitions (6) and (15) and from (19). Since (2/J)~~(~/~~~)= #aY, it is convenient for simplicity to regard U in (19) as a function of J rather than of 1,. A suitable form for the term U(&) in (19) is provided by U(I,) = $A{G(J)}’ + cH(J)
(21)
where G : (0, a)+lR is a penalty function with the property that G(J) = 0 if and only if J = 1 and A > 0 is the penalty parameter. The constant c and the function W(J) with the property H(J) = 0 and H’(J) = 1 iff J = 1, guarantee that the reference configuration @ = Identity is stress free. The usual way in which the condition is enforced [6,7], amounts to taking H(J)
= J.
From (21) and the definition of ,&(fl, &, J3) given by (6) it follows that
(22) The response functions /!&(i = 0, 1,2) defined by (20) and (22) satisfy the compatibility equations (7), and completely characterize the constitutive equation (5) for the Cauchy stress tensor. The constant c is given by (23) for the reference configuration to be stress free. The corresponding expression for the first Piola-Kirchhoff tensor follows from (4) and the elasticity tensor A can then be computed using (8).
112
J.C. Simo and R.L.
Taylor, Penalty function formulations
for incompressible
nonlinear elastostatics
In order to discuss appropriate forms for the functions G(J) and H(J) in (21) it is convenient to consider separately the effect of the penalty term JA{G(J)}~.Thus, we rewrite (19) in the form
I@(&,Iz, 13)= $h{G(J)}2
(244 WW
-I- W" ,
w* = cH(J)+ Iv.
3.1.1. The choice of the penalty term In the context of any numerical procedure, like the finite element method, the condition A 3 ~0can only be enforced in a relative sense. It is therefore appropriate to view (24) in physical terms as the strain energy functional of a compressible material, which exhibits incompressible behavior as the parameter A + ~0. Accordingly, (24) should remain physically meaningful for the widest possible range of values of A. In particular, situations in which the stresses take a finite value as the volume tends to zero, or become zero as the volume tends to infinite, not only are physically meaningless, but for highly confined nearly incompressible bodies may lead to numerical instabilities of an iterative solution strategy. The behavior of the strain energy given by (24) can be easily examined as A varies in (0, cc) by considering the simple problem of the extension of a cube with unit length sides. We shall take the stretching ratio of its side, as the control variable. Thus, any possible deformation @ E C has a deformation gradient F = yI
(25)
with determinant J = y3. The associated first Piola-Kirchhoff perform the extension) is a hydrostatic pressure2
tensor (the force required
to
(26)
p = PW.
For a given value of A E (0, co), the following conditions are expected to hold: (c.1) The strain energy VSI+ ~0, as the stretching ratio y + 0 or y + a. (c.2) p(y) is a monotone function in (0, ~0) with p(1) = 0, p--f ---co as y -+ 0, and p + -1-00 as y++m. We shall consider in the sequel, the special case of the Mooney-Rivlin model. For this type of material the strain energy and first Piola-Kirchhoff stresses for the extension problem are w =
ih(f17(~“))2+ w* ,
W” = -(2az + a#YI(y”) + t(al(y2 - 1) + a,(?” -
1)), (27)
P = {2a2(y3 - y’H’(y”)) + a,(7 - y2H’(y3))}1. (a) For A = 0 we have the limiting case of no penalty term and I@ = W*. The choice H(J) = In(J) ‘Taking the stresses as control nonlinear elasticity. Seven solutions
(28) variables leads to the classical example which shows for a simple neo-Hookean material are possible [O].
lack of uniqueness
in
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear
ehStOStatiCS
113
produces the strain energy @ and first Piola-Kirchhoff stresses P plotted in Fig. 1 and Fig. 2, respectively (curves A = 0). Both I&’and P satisfying conditions (c.1) and (c.2). (b) For A > 0 the response is affected by the penalty term. Two common forms of penalty function are G(J) = (J - 1) [3-61 and G(J) = In(J) 171. The first one, together with H(J) = J 161, produces the family of curves plotted in Fig. 3. Notice that an instability appears around J = 0.5, and that P vanish when y = 0. A similar instability phenomenon occurs for G(J) = In(J), although around the very large value of J = lo3 [7]. However, the combined function {G(J)}* = (In J)” + (J - l)*
(29)
is a suitable form of penalty function and produces the family of curves plotted in Fig. 1 and Fig. 2 for different values of A. For each A E (0, a), @ and P satisfy conditions (cl) and (c.2). The functions H(J) and G(J) given by (28) and (29) do not introduce any special computational effort. In fact, the corresponding response function B0 is b0=${A(J(J-1)+lnJ)+c-12B2}.
(30)
3.2. Alternative formulation Instead
of considering
a relationship
x=1000 X=lOO A=
of the type given by (19) between
strain energy
x=0
IO
5C
T
4c 3c 0k
20
2 2
IO
u z 2
0
E L
-10
E c
IO 0 -10
tL 0
: 1, 0.5
-20 -30 -40
, EXTENSION
1.5
2.0
2.5
3.0
RATIO
Fig. 1. Strain energy in isotropic extension. Proposed model.
-5o_ -1 ”
us
1.5
1.0
EXTENSION
Fig. 2. Piola-Kirchhoff model.
2.0
2.5
3.0
RATIO
hydrostatic pressure. Proposed
114
J.C. Simo and R.L. Taylor, Penalty function formulations
0
05
IO
I 5
EXTENSION
Fig. 3. Piola-Kirchhoff
for incompressible
2.0
25
nonlinear elastostatics
30
RATIO
hydrostatic pressure. Other models.
potentials, attention will now be focussed in the constitutive equations (5) and (14) for the Cauchy stresses in the compressible and incompressible cases, respectively. We shall consider first the case of the Mooney-Rivlin model. Extensions to more general incompressible models will be examined later. Let g : (0, =J)-+R be a smooth real valued function, such that g(x) = 0 if and only if x = 1. Consider first the replacement of the term --PI in (t7a) by hg(J), where A > 0 is a real constant. Now let A -+ ~0; i.e., consider a sequence {A,,) such that lim.,, A,, = 00. Since p : B --$ R is a bounded function (say in the uniform norm) and g(J) = p/A, it follows that if we let A --) ~0,then g(J)+ 0. By the assumptions on g(J), J + I as A --, CQ. The described replacement of hydrostatic pressure -PI, shows that the incompressible case can be viewed as the limit of a sequence of compressible cases. However, in contrast with the linearized theory, this fo~ulation requires the modification of the term 5 in (17a). In fact, if ci remains unchanged one is led to the constitutive model
and substitution
into the general integrability
ial = i3&/aI,= 0 ,
conditions (7) yields
$zz= -ap,/ax, =0 .
Thus ri = 0, and we arrive to the remarkable
conclusion that the only possible hyperelastic
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
1JS
compressible model, compatible with the assumption of constants response functions & (a = 1,2) is the hydrostatic pressure a = g(J)I. Since the response functions PI and p2 can not be a pair of constants in the compressible range, although they must reduce to the Mooney-Rivlin constants aI and a2 whenever .I = I, we are led to the assumption
Substitution into the integrability conditions equations which have the simple solution
therefore, given by
a consistent extrapolation
LT =
{hg(J)
of the Mooney-Rivlin
- (al - a2)h(J)}I
+ y
and
h(J)
differential
model to the compressible
range is
(34)
B - T B-l
h : (0, co)--) R
where for more generality the function been introduced. Ilf we take g(J) = 4 $ (G(J))’
(7) yields a pair of ordinary
with the property
h(x)
= 1
iff x = 1 has
= -& H(J),
where the functions G(J) and N(J) have the same meaning as in the previous section, (34) corresponds to a strain energy potential of the form VV(Il, 12, Q
= U(I3) + i(al(Il
- 3) + a2(I2I;1
- 3))
where U(&) is again given by (21) with the constant c = a, - ct2. This strain energy potential has a structure different from that considered in (19) and discussed earlier. In contrast with the previous formulation, the resulting first Pio~a-Kir~hho~ tensor does not contain the invariant 12, and 13 only appears in /%. Thus, it leads to a simpler expression of the elasticity tensor A, more convenient for computational purposes. The form of strain energy given by (35) also leads to a penalty formulation of the incompressibility constraint. This point is examined next. 3.2.1. Pevtalty ~~~~~~ai~u~ for gh’lzeM~u~~~-~~~~~~ rnuciei Let us rewrite (35) in a form analogous to (24a), with W*
= (a2 - a&Y(J)
The total potential
-
$aj(il - 3)
energy functional
+- $at(Zf;’
corresponding
W*
given now by
- 3). to the strain energy
W*
is, for any
116
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
configuration
@ E C, given by
(37) where x is the vector from the origin to x = Q(X), and 6 are the body forces per unit of mass in the reference configuration. With this notation the variational formulation of the boundary value problem for an elastic material with strain energy given by (36) is then find @ E C such that
I, (G(JN2 dV} .
(38)
For isochoric deformations @ E C&c, I3 = I and (36) reduces to the strain energy for the Mooney-Rivhn model. By the assumptions on G(d), the condition c;PE Cl,,, is equivalent to G(~(~)) = 0. Therefore, the variational formulation of the boundary value problem for the Mooney-Rivlin model is given by find @ E C’i,, such that
II*(@) = mix {n*(A)} s.t. G(J(A)) = 1 .
If we let A + m, (38) is just the penalty formulation [l-3]. 3.2.2. ~xte~s~o~ to more general co~stit~tive models Let us examine for what class of incompressible response functions of the form
of the constrained
(39)
problem given by (39)
materials a simple relationship
between
analogous to that found for the Mooney-Rivlin model is in general possible. Substitution of (40) and (41) into the integrabihty conditions (7) yields, after some manipulation, the system of equations
(42)
The solution of (42) together with (16) is given by
h(J) =Ql/J,
c;bz(J)= 1/P-1
)
(43)
J.C. Simo and R.L. Taylor, Penalty function formulations for incompressible nonlinear elastostatics
117
together with to a strain energy potential of the form
where A and 23 are constants. For this class of constitutive models exactly the same formulation as that used for the Mooney-Rivlin model is applicable. Eq. (44) includes some forms of strain energy commonly used in rubber elasticity, like ti = 1”<&- 3) + (I* - 3) .
(45)
4. Final remarks The penafty fo~ulation of the incompressibiIity constrain in finite elasticity amounts to considering an extend constitutive model which, although no longer incompressible, exhibits incompressible behavior in the limit as the penalty parameter tends to infinity. It has been shown in this paper that consideration of the physical signi~cance of this extended model leads to a proper choice of the penalty term. In particular, the proposed form has proper behavior over all the range of volume changes. Consequently, numerical instability should not appear in an iterative solution strategy based upon Newton’s method. Most importantly, the proposed form of the penahy function does not involve any additional computational effort. Finally, for a certain class of materials a simpler extension to the compressible range was proposed as an alternative to the one commonly used [3-71. This class of materials includes some frequently used in rubber elasticity, in particular the important Mooney-Rivlin model. For this constitutive model the proposed penalty formulation results in substantial savings over previous proposals in computing a tangent stiffness matrix.
Acknowledgment
Support for this research was provided by the Malaysian Rubber Producers Research Association of Brichendonbury, England, as part of a program of research on the use of rubber bearings for the seismic protection of buildings. This support is gratefully acknowledged.
References [l] R. Courant, Calculus of Variations and Supplementary Notes on Exercises (New York University, New York, 1962). [2] G.D. Luenberger, Optimization by Vector Space Methods (Wiley, New York, 1969). [3] J.T. Oden, A theory of penalty methods for finite element approximations of highly nonlinear problems in continuum mechanics, Comput. and Structures 8 (1978) 445449.
118
J.C. Simo and R.L. Taylor, Penalty function formulations
for incompressible
nonlinear elastostatics
[4] J.T. Oden, Finite Element for Nonlinear Continua (McGraw-Hill, New York, 1972). [5] D. Malkus, Finite element with penalties in nonlinear elasticity, Internat. J. Numer. Meths. Engrg. 16 (1980) 121-136. [6] D. Malkus and E.R. Fuller Jr., An isoparametric finite element model for large strain elastostatics, J. Res. Nat. Bureau of Standards, 86(l) (1981). [7] T.J.R. Hughes and E. Carnoy, Nonlinear finite element shell formulation for large membrane strains, Civil Engineering Lab., Port Hueneme, CA, Rept. N68305-0357-4185 (1981). [8] J. Marsden and T.J.R. Hughes, Topics in the mathematical foundations of elasticity, in: R.J. Knops, ed., Nonlinear Analysis and Mechanics: Heriot-Watt Symposium, Vol. II, Research Notes in Mathematics, Vol. 27 (Pitman, Belfast, 1978). [9] M. Gurtin, Topics in Finite Elasticity, SIAM Vol. 35 (SIAM, New York, 1981). [lo] C. Truesdell and W. Nell, The non linear field theories of mechanics, in: S. Flugge, ed., Handbuch der Physik, Vol. III/33 (Springer, New York, 1965). [ll] 0. Zienkiewicz, The Finite Element Method (McGraw-Hill, New York, 3rd ed., 1977). [12] R. Hill, On uniqueness and stability in the theory of finite elastic strain, J. Mech. Phys. Solids 5 (1957) 229-241.