Complementary mixed finite element formulations for elastoplasticity

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 74 (1989) 177-206 NORTH-HOLLAND

COMPLEMENTARY MIXED FINITE ELEMENT FORMULATIONS FOR

ELASTOPLASTICITY J.C. SIMO and J.G. KENNEDY Applied Mechanics Division, Stanford University, CA 94305, U.S.A.

R.L. TAYLOR Department of Cit~l Engineering, University of California, Berkeley, CA 94720, U.S.A. Received 17 December 1987 Revised manuscript recieved 11 October 1988 A global formulation of the principle of maximum plastic dissipation is systematically exploited to construct complementary mixed finite element formulations for elastoplasticity. Completely general return ~apping integration algorithms are obtained as Euler-Lagrange equations of a temporally discret~ed Lagrangian functional. The flow rule is no longer enforced point-wise, but rather at the element level in a weak fashion that couples all Gauss points within an element. The resulting algorithms can be linearized exactly in closed-form for an arbitrary functional form of the yield condition and hardening law. The present approach enables one to extend successful superconvergent quadrilaterals to elastoplastic analysis. Numerical simulations involvint; plane-stress elastoplasticity are presented.

Introduction

Over the past ten years, a general methodology for the formulation of integration algorithms for elastoplastic constitutive equations has been developed, which finds its point of departure in the radial return algorithm of Wilkins [1] for plane strain J2-flow theory. From a computational perspective, considerable work has been devoted to the accuracy analysis of this basic algorithm, see [2, 3], and its extensions to more general plasticity models, see [4-10]. Currently, these algorithms are viewed as a class of product formula algorithms emanating from an elastic-plastic operator split. Within this framework, consistency is enforced by means of a closest-point-projection of an elastic predictor state onto the elastic domain. In [11, 12], this point of view is exploited to develop general purpose algorithms capable of accommodating arbitrary functional forms of the flow rule, hardening law and yield condition. From a mathematical perspective, an independent development in a rather general setting starts with the work of Moreau [13], who coined the expression catching up algorithms. Unfortunately, it is fair to say that, with the exception of the work of Nguyen [14], this and subsequent formal developments, see [15] or [16], have passed largely unnoted in the computational literature. Implementations of existing return mapping algorithms rely crucially on the view of the discrete elastoplastic problem of evolution as a strain driven problem. As a result, within the 0045-7825/89/$3.50 © 1989, Elsevier Science Publishers B.V. (North-Holland)

178

J.C. Simo et al., Complementary mixed finite element formulations

framework of displacement-like finite element formulations, plastic loading is tested independently at each quadrature point of the element, and an independent return mapping algorithm is performed at each of the quadrature points for given incremental displacements. Effectively, such an algorithmic treatment corresponds to a strain space formulation of elastoplasticity, where total and plastic strains are the independent variables, while the stress is regarded as a dependent variable which is computed from the ~lastic strains by means of the stress-strain relations. In contrast to this view, a sizable body of the mathematical literature on plasticity has been concerned with the formulation of the elastoplastic problem with the stress field as the independent variable, the so-called stress problem, see [15-18, 21]. Mathematically, this problem is considerably simpler than the strain-problem which requires careful definition of the smoothness of the strain field, the so-called space of bounded deformations, and elaborate tools of convex analysis, see [23, 46] for a comprehensive review. In this paper, we explore an alternative algorithmic treatment in which stress and displacements are independent fields and are interpolated independently through a mixed variational formulation of the Hellinger-Reissner type. The central issue addressed here is not simply the derivation of a mixed formulation of the elastoplastic problem, indeed many such formulations exist, but rather the precise development and implicit implementation of the algorithmic structure associated with the discrete mixed problem. Our development hinges on a variational formulation of elastoplasticity based on the principle of maximum plastic dissipation, often credited to Von Mises (see [24]), and weU-known to be equivalent to normality. This principle plays a central role in the mathematical treatment of the elastoplastic problem by methods of convex analysis, as in [15, 16, 25-27, 45]. Following standard procedures, we obtain a complementary variational formulation of the continuum elastoplastic problem by a Legendre transformation that introduces the stress tensor as a basic independent variable. The novelty in the present approach lies in our subsequent algorithmic treatment of this functional, which is transformed into a time-discrete functional by means of an implicit backward-Euler difference scheme. We then show that the first variation of this discrete functional yields the weak form of the appropriate retui'n mapping algorithm formulated in stress space, along with the discrete algorithmic versions of the hardening law and the consistency condition. By introducing finite element interpolations for the stress and the displacement fields, we obtain a discrete optimization problem which can be solved by convex mathematical programming techniques. In particular, we consider the solution procedure by the classical Newton method in detail. This leads to the extension to the present context of the notion of consistent elastoplastic tangent moduli proposed in [6]. From a computational standpoint, the present formulation results in a finite element architecture that is substantially different from current algorithmic implementations. The fundamental difference is that the plastic return mapping algorithm can no longer be formulated independently at each Gauss point, in contrast with displacement-like methods. Here, the closest-point-projection iteration that restores consistency is performed at the global element level and involves all the Gauss points within the element. The present approach enables one to extend successful mixed finite element interpolations formulated within the framework of the Hellinger-Reissner principle to the elastoplastic regime. As an application, we consider an interpolation recently proposed by Plan and Sumihara [28], which is highly insensitive to mesh distortion, free from locking in plane-strain quasi-incompressible elasticity, and appears to yield superconvergent results in bending dominated problems.

J.C. Simo et al., Complementary mixed finite element formulations

179

1. Variational f o r m u l a t i o n o f plasticity

The variational formulation of plasticity considered in this paper relies crucially on the principle of maximum plastic dissipation, often credited to Von Mises (see [24, p. 60]), and subsequently considered by several authors, e.g. [29, 44]. For recent discussions within the context of finite strains see [30-33].

1.1. The principle of maximum plastic dissipation; local form We summarize the basic results needed in our variational formulation of plasticity. Further details in a somewhat more restricted context are given in [34]. According to ~ e principle of maximum plastic dissipation, the actual state is the one among all possible admissible states (i.e., states that satisfy the yield condition) that maximizes plastic dissipation. A precise formulation goes as follows. Assume the existence of a Helmholtz free energy function, ~(e t, e p, at), of the form ,,

ja,

(1.1)

Here, e t := V'u t is the strain tensor, e p is the plastic strain tensor, and at is a vector of suitable internal plastic variables, all evaluated at time t. Here, VSut denotes the symmetric part of the displacement gradient. Next, via the Legendre transformation1, introduce the internal variable vector qt, O(q,) := - ~(,~,) - ~,. q,,

(~.2a)

q~ = - O , , ~ ( a t ) ,

(1.2b)

so that a t - -aq@(¢).

Observe that qt = - O , ~ ( e t, e p, a t) may be interpreted as the thermodynamic force (affinity) conjugate to at. By definition, the local plastic dissipation is given by

a,) . ~_' Dp[e,, e,", ~,; ,~,P, ,~,] := - aci,(e,,0eEP, p

a~Ce,,0aeP, ,~,) • 4 , .

(1.3a)

Standard arguments exploiting the Clausius-Duhem inequality (along with the assumption that unloading is elastic) show that the stress tensor is given by o"t = O~/aEt-Vqt(et- eP). It then follows from (1.1-3a) that (1.3b) Let the yield function be defined in stress space by 0(o', q) = 0. The admissible set of stress states, then, in the convex set defined as 1 Introduction of this parameterization is essential to obtain a symmetric discrete return mapping algorithm and algorithmic elastoplastic tangent moduli in Section 3.

J.C. Simo et al., Complementary mixed finite element formulations

180

~,,q := {(d', ~) ER6 x Rq I 4~(&, #)~<0}.

(1.4)

Denoting by o', and q, the actual states of stress and internal variables, the principle of maximum plastic dissipation is equivalent to the statement that, for given ( ~ , &t),

Dp(~r,, q,; ~p, ,~,)I> Dp(,~, #; ~,~, ,~,) v(~, #)E ~ , ,

(1.5)

or, equivalently,

(~rt, q,) = arg[ max {DP(&, ~; ~P, &t)}]. (o',#)EKuq

(1.6)

The fundamental significance of this principle lies in the following proposition.

PROPOSITION 1.1. Maximum plastic dissipation implies normality of the flow rule in stress space, a potential form of the hardening law, and loading/unloading conditions in KuhnTucker form. PROOF. First one recasts (1.6) into a minimum problem merely by reversing the sign, that is

(~rt, q,)= arg[ man. {-DP(~, ~; ~P tit)}]. (&,#)EK~,q

(1.7)

Next, one reformulates (1.7) as an unconstrained problem by means of the Lagrange multiplier technique. Accordingly, consider the cone K p "= {By E

L'(a) 18~~>0}.

(t.8)

Introduce the Lagrangian functional associated with (1.7) defined as

lip "= -DP( ~, #; ~P, ~,) + Y6(#, #),

(1.9)

where # E K p. The standard Kuhn-Tucker optimality conditions for (1.7) are then formulated in terms of the Lagrangian (1.9), and are given by (see, e.g., [35, p. 314] or [36, p. 724]) a[p

0~r

~" ÷ #,0.6(,T,, q,)=o

(1.1Oa)

a~_p Oq -

&' + ~,,Oq~(O',,q,) = 0 ,

(1.10b)

1,t~>0,

0(o't,q,)~0

and ~,tff(¢t, q t ) - 6 .

(1.10c)

These conditions yield the statement of normality of the flow rule, the potential form of the hardening law, and the loading/unloading conditions. (The consistency condition ~ = 0 also needs to be added). I-1

J.C. Simo et al., Complementary mixed finite element formulations

181

R E M A R K 1.1. It follows from (1.2b) that dtt = -O:qO(qt)#t. Consequently, the local hardening law (1.10b) can be rephrased in the equivalent form

4,=-.i,,[02,o(q,)l-'oj,(o.,,

q,),

(1.11)

which serves to define the generalized hardening moduli h(~r, q). If one postulates a quadratic form for O(qt) , i.e., O(q,) = ½q,.D-lq,, then [~2q0]-1-D is constant.

1.2. Complementary variational formulations of elastoplasticity The complementary variational formulation of elastoplasticity is obtained through t w o steps: (i) A global formulation of the principle of maximum plastic dissipation, followed by (ii) a Legendre transformation that eliminates the strain tensor e, in favor of the stress tensor ~rt as an independent variable.

1.2.1. Notation Let us start by introducing the space of kinematically admissible variations (virtual displacements) defined for simplicity a s 2 v : = {,~.

n->R' l H E [ H * ( n ) ] ~ ;

(1.12)

~/la.n = 0 } ,

where N ~<3 is the spatial dimension, HI(/~) denotes the space of functions with derivatives bounded in energy, and 0 . n C On is the part of On, the boundary of the body/] C R N, where the displacement field if specified as u,[o n = ~. In addition, we let 0~,/~ C On to be part of the boundary where the stress vector ~s specified as ¢,n]o~a = t. Next, we introduce the total energy functional, at time t, •

ill

•

D

II(u,, et, e~, at):= fn ~(e,, e~, at) d/] + IIext(u,),

(1•13a)

where no,t(u,) := -

pb . u, d n - f, ul] l. u, d r

(1.13b)

is the potential energy of the external loading and pb is the body force•

1.2.2• Complementary energy functional A mixed formulation in which the strain tensor e, is replaced by the stress tensor ut as an independent variable is obtained by introducing the Legendre transformation of q ' ( e , - eP), i.e.,

• (t,-

-x(¢,) +

[e,-

(1.14)

where x(crt) is the complementary stored energy function• We then have the standard 2 This assumption is not realistic. A more appropriate functional framework is the space of bounded deformations (see Temam [23]). J

J.C. Simo et al., Complementary mixed finite element formulations

182

relations V ~ ( E , - e,p) = ~r,,

Vx(~rt) = e , - ..~P,.

(1.15)

By substitution of (1.2a) and (1.14) ~nto (1.1) and inserting the resu?~ into II(ut, ~.t, ~.P, at) defined by (1.13a) we obtain the complementary functional l l ( u , o ' , e P q)t "=

L

{-X{,~:)-O(q)-o~'q+o':[VSu-e,P]},dn+Ilext(U,),

(1.16)

where u t -- I] • V, (o't, qt) EO(~gq and r~text(~t) is defined by (1.13b).

1.2.3. Variational characterization of plastic response As shown below, a variational characterization of the evolution of plastic flow in stress space is accomplished by introducing the Lagrangian functional associated with the plastic dissipation over the entire body ~2 C R 3. From the local expression (1.3b) for the plastic dissipation, by appending the constraint through a Lagrange parameter as in Proposition 1.1, we obtain

"= fo

,,, + q, .

-

(1.17)

( ,,, , q , ) ] t i n ,

where ~,, E K p is the Lagrange parameter, in physical terms the plastic consistency parameter, and K p is the proper positive cone defined by (1.8).

2. Discrete variational problem To derive discrete governing equations we first construct the discrete complementary functional by integrating the dissipation functional Lp over the time interval [0, t,,+l] C R+ with a backward-Euler difference scheme. Then, discrete variational equations are obtained as Euler-Lagrange equations of the discretized functional. Explicitly, the discrete complementary functional for elastoplasticity is defined as

ftn +i

l=l(u, o', e p, q)[. "= l](u, o', e p, q)l.+t +

jr,,

Lp d~',

(2.1a)

which is a discrete statement of the first law of thermodynamics for isothermal quasi-static problems. That is, (2.1a) may be interpreted in physical terms as follows Potential energy[, = Potential energyl, +t + Dissipationl~ +t.

(2.1b)

One should also note that the second law of thermodynamics (or Kelvin's dissipation inequality) implies Dissipationl~ +t ~0. Making use of a backward-Euler difference scheme, we define

183

J.C. Simo et al., Complementary mixed finite element formulations

f nt'+t L P d T "----" f~ [On+ 1" ( ~"n P+ l- £ P ) - A T n + l t ~ ( O . n + l , q n + l ) + q n + l . ( O l ' + l - O l n ) ] d , ( "~ jt

(2.2) where ATn+I : = '~n+l -- q/n" By substituting (2.2) into (2.1a) the potential energy II" at time t" becomes a function of the state variables at t'+l, which we shall denote by the symbol 0-n+,Accordingly, substitution of (2.2) and (1.16) into (2.1a) yields the explicit expression

~''+1"----f~[--X(°''+I)- O(qn+l)- a " ' q n + l + fo [~r"+l:(VSu'+l- ep)] da

A~'+I~/)'+I] d a (2.3)

+ IIext(Un+l) •

Enforcement of the stationarity conditions for (2.3) through application of the directional derivative formula 3 leads to the following discrete Euler-Lagrange equations for the variables (o', q) E K,,q, AT E K p and u with u - ti E V:

t"

t"

Dff.'+l " 1~ --- JD [ O n + 1" VSl~ - pb " ~] d a - Ja

gn

DO_'+, • 8 ~ - f~ 8o': [(V'u - e P )

dr

=

(2.4a)

0

(2.4b)

- 8¢X(~r'+,) - A%+,8~4}'+,1 d a = O,

DB.,,+1" 8q = fn Bq.[-Oq@(q'+ 1) "~"{~qO(q~) DB.'+ I " By = fo

[" ~

(2.4c)

-- AT,,+ lOq~tn+l] d a = O,

(2.4d)

~(o',,+,,q,,+z)SydKt •O.

These four variational equations hold for any ~ E V , 8~rE[L2(n)] 6, 8 q E [ L 2 ( n ) ] q and By E K p.

R E M A R K S 2.1 P and the internal variables a'+l at the current (1) Note that the plastic strain tensor e'+x time t.+, are no longer present in variational equations (2.4). Equation (2.4b-d) may be viewed as the weak form of the following closest point projection algorithm formulated in stress space: e'P+ 1 -- f "p "1" A ~ . . ' + l O t r ~ ( O ' n + l ,

q'+l)

,

~lf'+ 1 ---1 Or" "b A T n + 1 0 ~ q 0 ( O ' ' + l , q ' + l )

,

(2.5)

in which (eP+1, an+,) are eliminated by means of the complementary relations E'+I = v s u . +1 - 0 . x ( ~ . + 1 ) '

a.+~ = - 0 , O ( q . + , )

3The directional derivative formula takes the form D f ' B m : = ~d I~--of.(... ¢z + a8o~,...).

(2.6)

J.C. Simo et al., Complementary mixed finite element formulations

184

A main implication discussed below is that within the present stress based framework one can no longer obtain perfectly plastic behavior, since in the absence of hardening the elastoplastic problem becomes ill-posed in stress space. (2) If equations (2.4b-d) are enforced point-wise one recovers the displacement model. Within the stress based variational framework presented in Section 3, the consistency condition and hardening law are enforced point-wise, but the flow rule is enforced in a weak sense. Alternatively, the consistency condition and hardening law may be enforced in a weak sense as discussed in the Appendix. A related approach within the restricted framework of J2-plasticity is recently discussed in [37]. Weak enforcement of the flow rule, on the other hand, leads to a closest-point-projection iteration which is no longer performed independently at each Gauss point as in typical displacement implementations, but rather is performed in a global element iteration that couples all the quadrature points within an element, as is discussed below. (3) By resorting to a Hu-Washizu type of variational formulation, it is possible to enforce point-wise the flow rule, hardening law and consistency condition while, at the same tim~, obtaining a weak formulation of the stress-strain relations. In fact, this is the proper variational setting of assumed stress methods widely used in current inelastic calculations. We refer to [38] for a detailed account.

3. Mixed finite element formulation

The mixed finite element formulation discussed below will be based on discontinuous stress interpolations over the elements that define the finite element discretization. Our motivation for this approach is quite explicit. Although the developments that follow are general, the actual implementation o~tlined in Section 4 for plane stress is based on a discontinuous stress interpolation proposed by Pian and Sumihara [28]. This interpolation does not exhibit "locking" behavior in the incompressible limit for plane strain and is highly insensitive to mesh distortion.

3.1. Approximation; discrete Kuhn- Tucker conditions We consider discontinuous stress interpolations defined by the finite dimensional approximating subspace s

--

I

= S
(3.1)

Here, S : n e - , R 6 × R m are (6 x m) prescribed local functions. An explicit construction of S h is discussed in Section 4. As discussed in Remark 2.1(2), within the four field variational framework provided by (2.3) and (2.4), q,+~, AT,+~ ~ K p are additional fields which may be defined point-wise over the element or interpolated through the weak forms (2.4c, d). Here, we consider the former approximation scheme and replace (2.4d) with the standard discrete Kuhn-Tucker conditions

AT.+

(3.2a) (3.2b)

J.C. Simo et al., Complementary mixed finite element formulations

185

Sim~'llarly, point-wise enforcement of (2.4c) yields the discrete hardening law

-o,o(e.+,) + o,O(q~)= a~+,o,#,(~r.+,, q.+,) .

(3.3)

3.2. Structure of the general return mapping Within the context of the discontinuous approximation (3.1) for o', a general solution strategy proceeds as follows. Equations (2.4b), (3.2a) and (3.3) are solved at the element level to obtain o'n+l, A3',+1 and qn+l for a given strain VSUn+l which, subsequently, are substituted into the momentum balance equation (2.4a). To formulate the discrete problem resulting from (2.4b), (3.2a) and (3.3) after substitution of the interpolation (3.1), we introduce the following notation r

Ee(~n+l) :--- JQest(x)O~rX(S(X)#~n+I)d a

(3.4a)

,

L(/3n+1, qn+1,3'n+i):= fa, St(x)a~¢Orn+1,qn+1)3'n+1(X)d a ,

(3.4b)

Etrial n+1:mfae st(x)[Vsun+I -

(3.4~)

ep] d~,

where, hereafter, 3',+1 is used in place of the more appropriate symbol A3'n+ 1. In the case of plastic loading within the element (i.e., 3'n+1> 0), equations (2.4b), (3.2a) and (3.3) then yield the discrete problem R ( ~ , q, 3')n+1

: --'--

-Ee(~n+l)- L ( ~ n + l , qn+l,3'n+1)+

o~O(q.+,) - a,O( q.) + v~+,o,¢(~+,, q~+,) = o, 3'.+xO(S/3n+'J,qn+1)= 0,

l~trial ----0 ~"n+l ,

(3.5a) (3.5b) (3.5c)

where 3', +1 ~ 0 and 0n +1 ~<0. The system (3.5) constitutes the appropriate return mapping for the elastoplastic problem. As a consequence of the discontinuous stress interpolation (3.1), (3.5) may be solved (locally) at the element level for/3n+1, 7n+1 and qn+l. A general algorithm for the solution of (3.5) for an arbitrary form of the yield function O(zr, q) is given below.

3.3. Numerical solution strategy; general yield condition The solution of (3.5) at the element level may be accomplished by a systematic application of Newton's method as follows. We regard (3.5a, c) as a nonlinear system of equations in/3n+ 1 and Y,,+I in which qn+l is treated as a function of both /3n+1 and 7n+1 through (3.5b). A separate Newton iteration for qn+l is then performed on (3.5b) once /3n+1 and 3'n+1 are determined. The linearized system corresponding to (3.5a, c) is obtained through a systematic application of the directional derivative formula and (at each time step t,+l) takes the form

o = R ¢*) + ~d I~--o R(P(*) + a A~(*~, V(*) + a A3'(*))

(3.6a)

J.C. Simo et al., Complementary mixed finite element formulations

186

or, equivalently 0 - R(,) _ [ ~ ( k ) ] - I

A~[](k) _ fo e StN(k)

(3.6b)

AT(*)(X) d ~ .

In addition we have

0= &(*)+ ~aadI,~=o~('B(k)+ a A~ (k), T (*) +

a

AT(k) )

(3.6c)

or, equivalently, 0-- 0 (k) + N(k)'s Aft (k) - [0qc~(k)]t [02qqO .4_ T(k)C~2q~(k)]-l[Oq~(k)] AT(k) '

(3.6d)

where the superscript (k) refers to the iteration index in the local Newton iteration and =(~,+,, A ( ' ) (k) :----"'¢.~(k÷1~ , n + l - - ( ' ] n~(k) +l ' [~(k)]-l.

s

:=

--"

(3.7)

fo est[O2¢x + T02o~ -

2

"'

-- T 2Oo.q~[C~ 2 2q q e "["

TO2qq~]-1 02qtro](k)s d~'~ ,

(k)

+

At each step (k) of the Newton iteration, the linear system (3.6b, d) is solved for Aa(~') "-'/"n + 1 and (k) A separate Newton iteration for ~,,÷: ,.(,) is then performed. Details of the step-by-step AT,+~. solution algorithm for this general problem can be found in [39]. Observe that the symmetry of -.+t='<*)in (3.7) leads to a symmetric return mapping algorithm on the element level as well as symmetric algorithmic elastoplastic tangent moduli (see Section 3.5) on the global level. 3.4. Numerical solution strategy; model problem

To illustrate the numerical solution of system (3.5), we consider a simple model problem.

3.4. I. Formulation of model problem We consider plane stress elastoplasticity with a Von-Mises yield criterion and linear isotropic hardening formulated in terms of the equivalent plastic strain q = #P, i.e., ~(o', q) := f(o') - ~ 3~ g ( q ) ,

(3.8a) (3.8b)

g( q) := orr + H ' q .

(3.8c)

,[2,

Here ¢rr > 0 and H' > 0 are material constants. Following [7], we define the function f(zr) as

f(o.):-~[o.teo.] 1/2 ,

-1 0

2 0

(3.9)

Observe that for this choice 0~,f= [o "t e ~r] -~/2 P o', and, by construction, ]]0~,fH = 1. Consequently, (3.8b) becomes linear in ~, and takes the form

1. C. Simo et el., Complementary mixed finite element formulations

187

Note that the Legendre transformation (1.2) is unnecessary in the ease of a single scalar internal variable since a symmetric algorithm emanates directly from a formulation in q. We complete our model problem by further assuming that the elastic response is linear, that is X(lflr) " = I o ' t c - l o

",

(3.10)

where C := [O2~X(o')] -1 is the elasticity tensor, assumed to be constant (and point-wise stable, i.e. g' t> llgll 2, for some constant a >0).

3.4.2. Return mapping algorithm Because of the linearity assumption on the elastic response, Ee(/3.+ 1)= Hell.+1 is linear in gl,,+1, with H e defined as

H" "= fn, s ' ( x ) c - l s ( x )

(3.11)

da .

In addition, the discrete yield function 4}(er.+1, q. +1) becomes linear in 3'. +1 (recall that ,/. +1 is used in place of the more appropriate symbol Ay.+ 1) due to the linearity of the hardening law. As a result, in the ease of plastic loading within the given element (i.e.T.+I >0), (3.5a-c) yield the following discrete return mapping problem, Rn+l

"ffi - H •~1.+1 - L(~.+I, T;,+I) + 17trial --,,+~ ~

~.+1 := f(er,,+,)- ~

0

(3.12a)

,

(3.12b)

K ( q . ) - 3ZH'T,,+; = 0 ,

where backward-Euler integration has been used on (3.8b) to yield the following linear discrete hardening law q . + l - qn = Vr] %+1.

(3.13) .(k)

Linearization of (3.12a) about the state (fl, T, q).+l yields n+t m

[

'sdal' ' Aft
H -1

.I.+1

(3.14a)

where ]-1 : = H e + f

J-r n+lJ

ja,

~t~2

[(k) ~. (k) l'-~

~, v t r ~ J n +

l..,Tn+

l lat j d~'~ .

(3.14b)

The system (3.12) may be solved effectively for/3,+1 • R m and T,+I E R by the algorithm summarized in Boxes 1 and 2. Observe that this algorithm can be phrased as an elastic predictor/plastic corrector solution scheme. However, in sharp contrast with displacement formulations, see for instance [6, 7], the return mapping algorithm cannot be formulated independently at each Gauss point. If yielding occurs at any Gauss point, the return mapping in Box 2 is conducted over the entire element. The following remarks concerning the procedure summarized in Boxes 1 and 2 should be noted.

J.C. Simo et a!., Complementary mixed finite element formulations

188

Box 1 Integration algorithm. Elastic predictor. Model problem (i) Initialization (k = 0):

fo s,(x)c,s¢x) d£~ ,

H e :=

~ t r i a l :__. S ~"n + 1 fn,

'

(x)[

VS

u,,+l - eP] d•

(ii) Calculate elastic predictor: #(0)

__ j g l e - l ~ t r i a l aa a.~n+ 1 •

(iii) Check for yielding at the element:

f(S(xs)#'°)) IF ~¢°)(Xs) ~<0 for all x 8 E ~2, THEN Elastic Process; EXIT ELSE Plastic Process; GO TO Box 2 ENDIF.

R E M A R K S 3.1

(1) In Boxes 1 and 2, the notation (.)(k), Xt8 and N8 refer to the following: (.)(k) . __/. ~(k) k In+l

,

xs := x E n e at Gauss point g, N s := number of Gauss points in/~,,

f~:= domain associated with element e. (2) In addition, the notation (@) stands for the Macauley bracket of @, i.e., (@)= @ if @>0, and (O}ffi0 if ~ < 0 . (3) Integration over ~, is evaluated via Gauss quadrature of the form

fn

x ( x ) d n := ~ [X(x)LsJsWg,

(3.15)

g--1

where W. is the Gauss weight and Js is the Jacobian of the isoparametric mapping at x s, i.e. J8 :ffi det I0qt(~)/0~]xg, where W(~) is defined in Section 4. (4) The crucial role played by the discrete Kuhn-Tucker condition in step (iv) of Box 2 should be noted. (5) The present complementary formulation is ill-posed for perfect plasticity (H' =0), as suggested by the factor 1/H' appearing in Box 2. 3.5. Displacement approximation; residual and tangent operator

We derive the appropriate expressions for the residual force vector and the tangent stiffness

J.C. Simo et al., Complementary mixed finite element formulations

189

Box 2 Integration algorithm. Plastic corrector. Model problem (iv) Enforce discrete Kuhn/Tucker conditions to compute 3,(k): 1

'Y~)=~H'--'-; ('~'k'(xe))'

g=l""'Ne'

I,~,:= { g e {1,2,. . . ,Ng} [g(k)(xs)>O) .

(v) Compute constitutiveresidualand check convergence:

L (k):- 2

[st(x)a~ffk)Ls~(sk)']sWg,

gEla©t

R(k) 2_ _He~(k) -- L(t) + j~trial IF

IIR"°II

>

,

T O L THEN

(vi) Compute consistent tangents and get Ap(k): [~(,)]-t

:= H • +

2 (k)s]..~. (k) :,w,, [s t a..f

~.

g¢lact

&.8(.) = (,..--(k))-, + IH--"; ~' [$t(O'f(k))(a"f(k))t$]'J~We

R(')"

gela©t

(vii) Update/3 (k) and compute ~(k+ ~):

#(k+,)

=

p(h) + Ap("),

~(k)(Xe):-- f(S(xs)p ~ ) - ~

K(q,,(X.)),

g = 1,.. ., N e ,

SET k = k + 1 and GO TO (iv) ELSE (viii) SET

q.+,(=.) = q.(~.) + ~

e,,+,(x.}

"v. ,,(~.),

?,,+,(xe)a.f,,+,( .)

ENDIF, EXIT.

matrix emanating from the weak form of momentum balance (2.4a). Assume a C°(n) displacement approximation defined by

v" :- { h e [C°(n)]N I

re.

}

= ~, NA(x)SdA, 8d A e. R N C I7. A=1

By introducing the [N × (N × N e n ) ] matrix of shape functions N := [Nl(x)l"" where 1 is the (N × N) unit matrix, and setting 8 d t " - - [ S d t l " • . 8 d Nt ~ , ] , we obtain

(3.16)

fcNen(X)l],

190

J.C. Simo et al., Complementary mixed finite element formulations

~hJo" = N(x)Sd,

V s~ h JOe... B(x)Bd,

(3.17)

where, following standard notation (see e.g. [40]), we have denoted by B "= V'N. 3.5.1. Residual force vector With the above displacement interpolation the weak form of momentum balance (2.4a) takes the form E

D[.+, . ~ h = ~ a h ' where a~ = 8dt[R~-f~].

(3.18)

ell

The vector f~ gives the element contribution to the external loading (emanating from H,,t(u)). On the other hand, the element internal force vector, R e, is given by

Re "= fl~ Bt(x)o'n+l

d~'~ -" G t p n + l

,

(3.19)

where we have set G "- fn. St(x)B(x) d n .

(3.20)

Note that the element internal force vector R, becomes completely determined from (3.19) once/]~.~ is obtained, from a given strain e,+~ ffiV'u, by the return mapping algorithm in Box 1 and Bvx 2. 3.5.2. Consistent tangent stiffness matrix The element tangent stiffness matrix is obtained by (exact) linearization of the weak forms (2.4) through a systematic application of the directional derivative formula. This, in turn, involvea the linearization of (3.19) and the linearization of the algorithm in Boxes I and 2 as follows. Here, we regard all the fields (t r, q, y) as functions of uhE V h. Then, for any incremental displacement field AuhE V h such that Au h := N(x)Adn+l, the directional derivative of the discrete weak form Ge(d,+~) is computed as d__

d~ ~-o

a.(d.+, +.

Ad.+,) = 8d'

ad.+1

' • Ad.+,

(3.21)

Thus, the element tangent stiffness matrix is given by K~'=Gt[OlJn+l/dn+l]. The matrix [~t3n+l/~d,+l] i,~ obtained by linearization of the discrete system (3.5). To this end, from (3.4c) and (3.20) we first observe that Etrial n+l ~--"G d ~ + l w f o e S(x)e p dl~

~

~ik7 trial a-*n+ |

ads+' = G .

(3.22)

Y.C. Simo et al., Complementary mixed fini:e element formulations

191

For simplicity, hereafter we restrict our attention to the model problem. Differentiation of (3.12a, b) then yields ,-,- 1 ~(~n+l ~t. Jne f S t Ou.(~tn+1 (~dn+l #7,,+~ d n = 6;, =,,+x #d.+ t [a,,.f]t$

a/],+l

8d.+t

(3.23a)

2 07,+1 = 0 , 3 H' 8d.+~

(3.23b)

where ~+1 is defined by (3.14b). Solving (3.23b) for [Oy~+llad.+l] and substituting in (3.23a) yields

a.s,,+, = [=._, + 3 f,,. S t(oo.(k,,+l)(a~.ck,,+,) tS d13 1-1 G.

(3.24)

Finally, by inserting (3.24) into (3.21) one obtains the expression of the element tangent stiffness in the form S'(a,,¢k,,+x)(O=,ck,,+x)tS dfJ]-I G.

(3.25)

R E M A R K 3.2. The tangent stiffness for (3.5) is found by differentiation of (3.5a, c), again considering q. +1 as a function of both/3,, +1 and % +1 via (3.5b), and proceeding along the same lines as above to obtain ([39])

• ~--"o[-,~ n + l -- fo

-1 ~n+ i s t~T'n+l~rn+ t 1s d D

G, (3.26)

¢.+, : - 1/[(0.¢)'(0',,e + where ~.+t and N.+t are defined in (3.7). Notethat Ke in (3.26) is symmetric for arbitrary

functional forms of the yield criterion (3.2) and hardening law (3.3). Global expressions for the residual force vector and the tangent stiffness matrix are obtained through the standard assembly procedure.

4. Application and numerical simulations We consider stress interpolations for a quadrilateral element with bi=linear isoparametric interpolation for the displacement field. The stress interpolation is a modification of that proposed by Pian and Sumihara [28]. For linear isotropic elasticity the resulting element appears to be optimal. Its performance under distortion and in the nearly incompressible limit is excellent. The proposed modification enables one to achieve all the computational advantages reported in [41]. The performance of this element, 11owever, appears to be superior to that obtained with the interpolation proposed by these authors.

192

J.C. Simo et al., Complementary mixed finite element formulations

4.1. Stress interpolation for a quadrilateral element

Let g~e C R2 be a typical quadrilateral element with nodes x A E R 2 (A = 1,2,3, 4), and denote by qt. g] --->g~e the isoparametric mapping, where g] : = [ - 1, 1] x [ - 1, 1]. We have t~ ~ b ---> ~ , ( ~ ) = ~, f c ~ ( ~ ) x ~ ,

(4.1)

A=I

where ~A(~) are the standard bilinear interpolation functions on the unit square/]. Let .~ "= qt(0) be the centroid of /~e and /~ the Jacobian of the isoparametric mapping at I!he centroid, i.e.,

o,t,(~) ]

[Pl P~]

(4.2)

Define ~rh(~) as the stress tensor tr convected by the constant mapping #, that is ~.h := l~-lO.l~-t. We consider a stress interpolation, discontinuous between elements, and defined in terms of five parameters as

,, ----'/'oh + A(~ 2 -

1 ~2)[0

(4.3)

where %h, A, B are constant and ~' are defined by the expression fa, ( ~ ' - ~ ' ) d ~ = O. The interpolation for the stress tensor ¢ is defined by mere transformation (push-forward) of (4.3) with F. Replacing ~.h by the alternative constant matrix tr~ : - l ~ . ~ t , the final result may be expressed in reduced matrix notation as (recall that x = ~ ( ~ ) )

0.11 0 .22 0.12

h

= [ I 1 ( 6 ~ - ~2)a,

I(¢'

- ~')a2]g-S(x)O,,

(4.4)

where I : = Diag[1, 1, 1] is the (3 x 3) unit matrix, tS~ is the (5 x 1) vector of stress parameters, and a t , a2 are constant (5 × 1) vectors defined in terms of the components of 1~ as

a, := [(~',y (~)" (~:,)l-'-'', (4.5) R E M A R K S 4.1 (1) The stress interpolation outlined above is obtained from that proposed by Plan and Sumihara [28], by replacing fi with (~:i- ~ ) , i = 1, 2. For linear elasticity where the tangent elasticities C are constant, this modification results in a block diagonal weighted compliance matrix H ~-1, i.e.,

"vol(/'Je)C -1 H~-J _

0[3×2 ] ]

(4.6)

J.C. Simo et aL, Complementary mixed finite element ~ormulations

193

where O-1 is a (2 x 2) matrix. Thus, He can be inverted in closed form. For the elastoplastic case this computational advantage is lost, even if a~,, f(o') is constant, as in the case of the Von Mises yield condition. (2) It can be easily seen that the interpolation outlined above passes the mixed patch test in the form recently discussed by Zienkiewicz et al. [42]. In particular, the one element test with minimum displacement constraints gives 2 × 4 - 3 = 5 displacement degrees of freedom, which is the number of free stress parameters/3.

"4.2. Numerical simulations, extension of a perforated plate Two numerical examples are considered to demonstrate the applicability of the new variational formulation of elastoplasticity as well as the reliable performance of the corresponding numerical implementation. All calculations are performed on a Convex C1 computer by inplementing the algorithm shown in Boxes 1 and 2 in an enhanced version of the nonlinear finite element computer program FEAP, developed by R.L. Taylor, and described in [40, Chapter 24].

4.2.1. Exter&~ion of a perforated plate A perforated plate characterized by the plane stress/,_-flow material model in Section 3.4 with isotropic hardening is considered first. Displacement controlled boundary conditions are imposed, as depicted in Fig. 1 and the material properties are: Young's modulus E - 70, Posson's ratio v - 0 . 2 , uniaxial tensile yield stress K0 -0.243, hardening modulus H ' - 0 . 2 . Figure 2 shows a comparison of the load-displacement curves (load calculated from nodal reactions), predicted from the mixed formulation and the standard displacement formulation, each using an identical mesh of 4-node quadrilaterals. The plane stress version of the classical radial return employed here is due to Simo and Taylor [6]. The mixed formulation provides more accurate load-displacement predictions than the displacement formulation as seen by comparison of the solutions for the 72 element mesh and the 722 element mesh. Aside from the added compliance of the mixed solution, however we observe that the load-displacement curves are in good agreement for the two formulations. .4 t t t

t t.~ a

lllll

18

I_ r

J I0

~1

Fig. 1. Geometry, mesh and loading configuration for a perforated plate.

J.C. Simo et al., Complementary mixed finite element formulations

194 2.2

S.

2 . 0 ""

1.8-

1.6-

o,

1.4-

1.2 "

1.0 "

(a) Oispl. method, 722 elmts • • • • • ,, (b) Mixed method, 722 elmts . . . . (¢) Displ. method, 72 elrn;s {d) Mixed ,method,72 elmts

0,8 "

0.6 0

u

u

t

2

4

6

8

Displacement

Fig. 2. Load-displacementcurves for perforated plate with: (a) 722 element mesh, displacementformulation; (b) 722 elementmesh, mixedformulation; (c) 72 elementmesh, displacementformulation;(d) 72 elementmesh, mixed formulation.

Figures 3-6 compare the evolution of the plastic zones predicted by the mixed and displacement formulations, again using an identical mesh of 4-node quadrilaterals for each. Specifically, the figures contain contour plots of ~ values, where •=

f(o')

Ko : = [K(q)]q_o.

(4.7)

X.V~'~ K0

The legend .shown in Fig. 3 applies to all of Figs. 3-6. In each figure, the plastic zone is characterized by ~k ~ 1.0. Again, the mixed formulation gives somewhat better (more flexible) results than the displacement formulation. With a 72 element mesh, the mixed formulation is better able to predict the semicircular elastic region neighboring (x, y) = (10, 0) as shown for = 0.15. Aside from this added accuracy, the yield zones predicted by the two formulations are in good agreement. As expected, the rate of convergence exhibited by the global Newton iteration is excellent, as shown in Tables 1-4. In discussing the rate of global convergence for this problem, it is important to recognize two facts. Firstly, Tables 1-4 suggest that the radius of covergence for this problem is such that quadratic convergence is attained only for energy norms below 10 -6. Secondly, as the mesh becomes more refined for either formulation, the rate of global convergence experiences some degradation, as shown in Tables 1 and 2 or 3 and 4. Consequently, since the mixed problem is more flexible than the displacement problem for a

J.C. Simo et al., Complementary mixed finite element formulations

l-'-1

195

~ < .90 .90 < ~)< 1.0

!Ii i Fig. 3.

Yield zone e v o l u t i o n for d i s p l a c e m e n t f o r m u l a t i o n with 722 e l e m e n t s . ti = 3.15. (d) ti = 4.65. (e) ti = 6.15.

1.o (a) ~-0.15.

(b) t i - 1 . 6 5 .

(c)

Table 1 G l o b a l N e w t o n i t e r a t i o n e n e r g y n o r m c o n v e r g e n c e r a t e for d i s p l a c e m e n t f o r m u l a t i o n , 722 e l e m e n t mesh L o a d step

Iteration 1 2 3 4 5 6 7 8

1 ( ~ -- 0.03) 0.905e + 0.214e 0.124e 0.831e O. 133e O.190e~ m

00 04 05 08 13 24

5 (~2 - 0.15) 0.929e + 00 0.120e - 04 0.449e - 06 0.A.A~A.e-- 08 0.680e - 12 0 . 2 2 8 e - 19 0 . 2 2 6 e - 30 n

10 (ti --- 2.65) 0.258e + 0.150e 0.827e 0.303e 0.730e 0.535e0.717e0.406e -

03 02 04 05 07 11 19 28

17 (ri -- 6.15) 0.257e + 0.996e 0.238e 0.149e 0.916e 0.211e0.177e0.135e -

03 03 04 06 07 09 14 24

J.C. Simo e: al., Complementary mixed finite element formulations

196

[~

~ < .90 .90 < ~ < 1.0

BI II

1.o

Fig. 4. Yield zone e v o l u t i o n for d i s p l a c e m e n t f o r m u l a t i o n with 72 e l e m e n t s . (a) ii - 0.15. (b) Ii -- 1.65. (c) ii = 3.15. (d) ti = 4 . 6 5 . (e) ti = 6.15.

Table 2 G l o b a l N e w t o n iteration e n e r g y n o r m c o n v e r g e n c e rate for d i s p l a c e m e n t f o r m u l a t i o n , 72 e l e m e n t m e s h Load step

Iteration 1 2 3 4 5 6 7

1

5

10

( ~ = 0.03)

( ~ = 0.15)

( ~ ffi 2.65)

0.270e 0.111e 0.119e 0.233e 0.161e 0,127e --

+ -

00 04 05 09 16 30

0.294e 0.104e 0.177e 0.305e 0.258e 0.217e --

+ -

00 04 07 12 21 31

0.817e 0.138e 0.705e 0.623e 0.124e 0.240e 0.662e

+ -

02 02 04 06 09 16 29

17 ( ~ ffi 6.15) 0.717e + 0.628e 0.4~ ~e 0.637e 0.269e 0.273e -

02 03 05 09 16 29

J.C. Simo et al., Complementary mixed finite element formulations

197

$ < .90 .90 < $ < 1.0

/

1.o

Fig. 5. Y i e l d z o n e e v o l u t i o n f o r m i x e d f o r m u l a t i o n w i t h 722 e l e m e n t s . (a) I i - - 0 . 1 5 . ( b ) t i - 1.65. (c) t i - 3.15. ( d ) ti - 4.65. (e) ri = 6.15. Table 3 G l o b a l N e w t o n i t e r a t i o n e n e r g y n o r m c o n v e r g e n c e r a t e f o r m i x e d f o r m u l a t i o n , 722 e l e m e n t m e s h Load step

Iteration

1

5

10

17

(ti -- 0 . 0 3 )

(ti -- 0.15)

(ti - 2.65)

( ~ = 6.15)

1

0 . 8 9 7 e + 00

2

0 . 2 2 0 e - 04

3 4 5 6 7 8 9 10

0.175e 0.449e 0.320e 0.559e m m ~ --

-

05 07 13 23

0.921e 0.217e 0.196e 0.883e 0.719e 0.125e 0.483e ~ ~ m

+ 00 - 04 - 05 - 07 - 10 - 15 --27

0.256e + 0.831e 0.683e 0.753e 0.871e 0,275e 0.176e 0.192e 0.578e0.255e -

03 03 02 04 05 06 07 11 19 26

0.245e 0.185e O. 156e 0.189e 0.951e 0.375e 0.210e 0.164e --

+ -

03 03 04 05 07 11 18 25

J.C. Simo et al., Complementary mixed finite element formulations

198

r-]

~ < .9o .9o <

B

< 1.o

¢> 1.0

Fig. 6. Yield zone e v o l u t i o n for mixed f o r m u l a t i o n with 72 e l e m e a t s . (a) ii = 0.15. (b) 6 = 1.65. (c) 6 = 3.15. (d) 6 - 4.65. (e) 6 - 6.15.

Table 4 G l o b a l Newton iteration e n e r g y n o r m c o n v e r g e n c e rate for mixed f o r m u l a t i o n , 72 e l e m e n t m e s h L o a d step 1

5

Iteration

( 6 = 0.03)

(6 - f). 15)

1 2

0.264e + 00 0.222e - 04 0.2b0e - ~ 0.332e - 08 0 . 2 2 8 e - 14 0.172e - 26

3

4 5

6

0.288e 0.164e 0.165e 0r 645e 0.956e 0.116e

+ -

00 04 06 1] 20 29

10 (6 - 2.65) 0.800e 0.970e 0.797e 0,755e 0.222e 0.320e

+ -

02 03 04 06 09 16

17 (6 - 6.15) 0.779e 0.196e 0.139e 0.278e 0.210e 0.564e

+ -

02 03 05 09 16 27

J.C. Simo et al., Complementary mixed finite element formulations

199

given mesh, the mixed problem converges slower than the displacement problem when the same mesh is used. One should note, however, that within a global Newton iteration, the mixed formulation typically requires less computational effort than the displacement formulation, since the mixed method performs the plastic return mapping only once within a plastic element, whereas in the displacement method the return mapping is performed four times (once at each Gauss point). 4.2.2. Bending o f a tapered beam In this example we consider a tapered beam with elastoplastic response characterized by plane stress J2-flow with isotropic hardening. The elastic version of this problem is often referred to in the literature as Cook's Membrane Problem". We impose load controlled" boundary conditions, as depicted in Fig. 7, and assume the following material properties:

1

Fig. 7. Geometry, mesh and loading configuration for a tapered beam. 2~

,

,

•

'

'

'

'

2O $ fSSS

.S . . . . .

/

" ....

" " " " "

/" ¢./

/ / i ]

s

AssumedStress FormulationI l l DisplacementFormulation. . . . .

]

i i 0o

i i

I

;0

I

,

15

,

I

20

,

i

.,5

~'0

35

Elements per Side - n

Fig. 8. Displacement of the top, right corner versus the number of elements per side in the finite element mesh for the (a) assumed stress formulation, and the (b) displacement formulation. Equal numbers of elements in the horizontal and vertical directions are ,~oa tt,.,.,,,,.h....~ . . . . . . .

,s,

a ~ww~n,ev~L.

200

J.C. Simo et al., Complementary mixed finite element formulations

Young's modulus E = 70, Poisson's ratio v =0.3333, uniaxial tensile yield stress K0 =0.243, hardening modulus H' = 0.2. The loading is increased in increments of AF = 0.1 until a final value of F = 1.8 is reached, which corresponds to a state in which nearly the entire specimen is plastic (except for two small elastic zones in the lower left and upper fight comers). Figure 8 shows the vertical displacement of the top fight node plotted versus number of elements per side, at a load level of F = 1.8, computed with the proposed mixed for~nulation and the standard 4-node bilinear isoparametric displacement formulation with the usual "strain driven" return mapping algorithm. As demonstrated in Fig. 8, for this problem, the displacement formulation exhibits a significant degradation in accuracy over the mixed formulation. In particular, the mixed formulation solution with 64 elements is more accurate . than the displacement formulation solution with 1024 elements.

5. Concluding remarks

We have presented a discrete mixed variational formulation of plasticity based on the principle of maximum plastic dissipation. The relevant variational principle is simply the discrete statement of the first law of thermodynamics over an arbitrary time-step [tn, tn+~]. Explicitly, we express the potential energy at time t n as the sum of the potential energy at time tn+~ plus the plastic dissipation during [t~, t~+~]. Integration of the dissipation over [t~, t,,+~] with a backward-Euler algorithm yields a discrete Lagrangian whose Euler-Lagrange equations are precisely the momentum balance, the flow rule, the hardening law and the consistency condition, where the latter three characterize a general closest-point-projection algorithm. This results in a discrete optimization problem, formulated in stress space, the solution of which has been considered in detail. Concerning this solution procedure, the following remarks are noteworthy. (i) Within the context of mixed finite element formulations with independent stress interpolations, the plastic return mapping algorithm can no longer be formulated independently at each Gauss point. The closest-point-projection that restores plastic consistency must be performed at the global element level, and involves all the Gauss points within the element. This is in sharp contrast with traditional displacement-like formulations for which the return mapping algorithm is performed independently at each quadrature t~oint. (ii) The notion of consistent elastoplastic tangent moduli, [6], also carries over to the present stress-space algorithmic framework. As in displacement-like formulations, these moduli are obtained by consistent linearization of the return mapping algorithm. We have shown that a closed-form expression can be obtained for arbitrary functional forms of the yield condition and hardening law. (iii) A main motivation for the present study is the extension to the elastoplastic regime of recently proposed mixed formulations. The present developments are applicable to any mixed finite element formulation of the Hellinger-Reissner type. In particular, as an application, we have considered an interpolation recently proposed by Pian and Sumihara [28], which is highly insensitive to mesh distorsion and appears to yield super'convergent results for bending dominated problems. (iv) The present framework carries over to the elastoplastic analysis of plates and shells. Recent work, [43], indicates that successful mixed interpolations for the membrane and

J.C. Simo et al., Complementary mixed finite element formulatio,,.s

201

bending fields can be developed within the context of a Hellinger-Reissner variational setting. As demonstrated in the numerical simulations, use of the traditional displacement formulation of elastoplasticity over the present mixed formulation may lead to significant degradation in the accuracy of the calculations for certain classes of problems, particularly those entailing mesh distortion.

Appendix A. Weak enforcement of the consistency condition

As noted previously, within the variational framework provided by (2.3) and (2.4), the consistency parameter A~, E K p can be spatially interpolated within the element. To this end, we introduce the approximating set K p' -- {A~ h E L 2 ( n ) ! A~/ h[ne = F t ( x ) Y e , eye ~ R s

and Ye ~ 0 } ,

(A.I)

where we have used the standard notation % ~>0 if and only if all the components are positive, i.e. ~,~ I>0 for A-~ 1 , 2 , . . . ,s. In addition, to ensure that KPhc K p we require that the prescribed set of functions FA : n,--, R+, (A - 1 , . . . , s), is positive, that is,

FA(X)>~O V x e n ,

and

Afl,..~.,s.

(A.2)

Since by assumption YeA~>0 for A = 1 , 2 , . . . , s, (A.2) implies that AT h ffi FA(X)7 ~ ~>0, so that A7 E K p are required. The stresses again are represented by the discontinuous interpolations defined by (3.1).

A. 1. Discrete Kuhn- Tucker conditions We recall that the consistency parameter A3,.+t E K ph, and the yield function 0. +t : - 0(or. +1, q. +1) (with zr. +1 E S h) must satisfy the Kuhn-Tucker conditions

(A.3) The discrete counterpart of these conditions, which result from the spatial approximation (2.3) and (2.4) is obtained as follows. First, introduce the notation f

6.+,(x):=

r(x)c,(s(x)IJ.+,, q.+,)d•.

(A.4)

Clearly, since ~b.÷~~<0 and FA(x)>~O,for A = 1,2,... ,m, it follows that ~,;÷i~<0. (Recall that this means ~.A+I(X) ~<0 for A = 1,... ,m). In summary we have the following discrete Kuhn-Tucker conditions o

~n+l~ 0 '

~n+l~O '

~tn+1~n+1 ~__0 ,

which govern approximation (A.1) and (A.2).

(A.5)

202

J.C. Simo et al., Complementary mixed finite element formulations

REMARK A.1. It appears that requirement (A.2) on the approximation may be relaxed to the weaker condition

fn FA(x)dfJ>~O

(A.6a)

for A = 1 , 2 , . . . , s .

This condition also leads to the same discrete Kuhn-Tucker conditions (A.5), as the following estimate shows:

OA"= fa. F+~(x), df~ ~<[max:~a~ok(x)] faeFA(X)df~ <-0.

(A.6b)

The last inequality holds because of (A.6a) and the fact that 0 ~<0 for all x E n,.

REMARK A.2. It is interesting to observe that the variational setting emanating from (2.4) also determines a discrete--interpolatedmyield function given by (~n+l(X) ,'~- rt(x)y-l~n+i

Y "~- fn e r(x)rt(x) d a ,

,

(A.7)

where 0,,.~ is defined by (A.4). To see this, simply note from (2.4d) that fa, Ay.+,0.+,(x)dn= 7:+,

fo.rr' dar-' fa"

'Io

= v.+,

q).+l dl'~

q.+,)dn e

- fn. A%+,ck(o',+~, q,+,) d/].

(A.8)

Thus, ~,.l(x) satisfies (2.4d).

A.2. Structure of return mapping Withi:a the context of the discontinuous approxima~tions (3.1), (A.1) and (A.2) for the stress tensor cr and the consistency parameter Ay, respectively, a general solution strategy proceeds as follows. Equations (2.4b, d) are solved at the element level to obtain stresses o',,+1 for a given strain VZu.,l which, subsequently, are substituted into the momentum balance equation (2.4a). To formulate the discrete problem resulting from (2.4b, d) after subsitution of the interpolations (3.1), (A.1) and (A.2), we introduce the following notation E~(/3"+1) "- fn, S'(x)O,, X( S(x) lJ,,. ~) d n ,

L
• =

f,

s

St(x)[V su.+ 1 - e.p] d£~ .

(A.9)

J.C. Simo et al., Complementary mixed finite element formulations

Equations (2.4b, e

203

d) then yield the discrete problem

i~?trial (p..,)-L(p.+,)v.., +=.., =o,

t

7.+1~.+1 =0,

(A.IO)

lp.+1 ~>0.

A general algorithm can be developed for the solution of problem (A.10) for an arbitrary form of the yield function ~(o', q) and generalized hardening moduli h(o', q). For concreteness, however, we shall again restrict ourselves to the particular form discussed in the model problem.

A.3. Numerical solution strategy To illustrate the numerical solution of system (A.10), we again consider the model problem discussed in Section 3.4, where, for linear elastic response, Ee(gl~.l) = Hell,,+1 is linear in /3,. t, with H e defined in (3.11) and Y :=

fQer(x)rt(x)

d/].

(A.11)

The discrete yield vector 0~+ 1 defined by (A.4) also becomes linear in Y~+l due to the linearity of the hardening law, i.e., ~,,+1 :ffi fo, l'(x)[f(S(x)[3"+t) -

X/]~K(eP)] d,O - ~H'Y1~,+~.

(A.12)

The discrete system (A.10) now reduces to a nonlinear equation in/3,.t E R m, which can be solved systematically using Newton's method. The step-by-step return mapping algorithm closely resembles that in Box 1 and Box 2 and may be found in [39]. The tangent stiffness for this case then becomes -1

Ke = Gt r ~._ 1

1

_t t

L'"+t = ~H; L"+IY L"+I

G

(A.13)

where G is defined in (3.20).

A.4. Interpolation of consistency parameter for a quadrilateral To complete the foregoing finite element interpolation, we consider the following choice for

r(x), ¢eh,

(A.14)

where ~A E S], A = 1,2,3, 4, are the Gauss points in natural coordinates, 8(.) denotes the delta function and ~ = qr,l(x). Clearly, (A.14) is the simplest choice that satisfies the condition rA(X)30, for x ffi qt(~), as required. In fact, Y becomes a diagonal matrix since

204

J.C. Simo et al., Complementary mixed finite element formulations

YAs=f~ ~ ( ~ - _Ea)8(g- ~s)S(~)d~ 1 d~2 =/~A8 J(ga),

(A.15)

where J(~A) is the Jacobian of the isopara_metric mapping evaluated at ~:A and a 2 × 2 quadrature is used (W8 = 1). Also note that ~ defined by (A.12) becomes ~A ~- ~A

-

2

t

(A.16)

Therefore, denoting by fA := f(S(~:,~)/3), the consistency parameter is determined simply as ~/A = (fA - V ~ K(e.P))/}H'. One should note that, in the case that a 2 x 2 quadrature is used over 1~e(Ws = 1), the delta function interpolation (defined by (A.14)) as well as the bilinear interpolation of A~,.+I between Gauss points results in a formulation identical to that in Box 1 and 2.

Acknowledgment We are indebted to M.S. Rifai for his help and comments on the subject of this paper, and his suggestion for the simulation in Fig. 8. S!lpport for this research was provided by a grant from the National Science Foundation. J.G. Kennedy was supported by a Fellowship from the Shell Development Company. This support is gratefully acknowledged.

References

[11 M.L. Wilkins, Calculation of elastic-plastic flow, in: B. Alder, S. Fernback, M. Rottenberg, eds., Methods of Computational Physics 3 (Academic Press, New York, 1964).

[2] R.D. Krieg and D.B. Krieg, Accuracies of numerical solution methods "forthe elastic-perfectlyplastic model, ASME J. Pressure Vessel Tech. 99, 1977.

[31 H.L. Schreyer, R.L. Kulak and J.M. Krame,r, Accurate numerical solutions for elastic-plastic models, ASME J. Pressure Vessel Tech. 101 (1979) 226-334.

[41 I.S. Sandier and D. Rubin, An algorithm and a modular subroutine for the cap model, Intemat. J. Numer. and Analyt. Methods Geomech. 3 (1979) 173-186.

[51 M. Ortiz, and E.E Popov, Accuracy and stability of integration algorithms for elastoplastic constitutive equations, lnternat. J. Numer. Methods Engrg. 21 (!085) 1561-1576.

[6] J.C. Simo and R.L. Taylor, Consistent tangent operators for rate independent elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 48 (1985) 101-118.

[7] J.C. Simo and R.L. Taylor, A return mapping algmithm for plane stress elastoplasticity, Internat. J. Numer. Engrg. 22 (3) (1986) 649-670.

[81 J.C. Simo and M. Ortiz, A unified approach to finite deformation elastoplasticity based on the use of hyperelastic constitutive equations, Comput. Methods Appl. Mech. Engrg. 49 (1985) 221-245.

[9] J.C. Simo, J.W. Ju and R.L. Taylor, Softening response, completeness condition, and numerical algorithms for the cap model, Internat, J. Numer. Methods Engrg., to appear. [I0] H. Matthies, A decomposition method for the integration of the elastic-plastic rate problem, lnternat. J. Numer. Methods l~ngrg., to appear.

[111 J.C. Simo and T.J.R. Hughes, General return mapping algorithms for rate independent plasticity, in: C.S. Desai, E. Krempl, P.D. Kiousis, T. Kundu, eds., Constitutive Equations for E.ngineering Materials (Elsevier, Amsterdam, 1987) 221-231.

J.C. Simo et al., Complementary mixed finite element formulations

205

[12] J.C. Simo, J.G. Kennedy and S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/ unloading conditions and numerical algorithms, Internat. J. Numer. Methods Engrg. 26 (1988) 2161-2195. [13] J.J. Moreau, Evolution problem associated with a moving convex set in a Hilbert space, J. Differential Equations 26 (1977) 347. [14] Q.S. Nguyen, On the elastic-plastic initial-boundary value problem and its numerical integration, Internat. J. Numer. Methods Engrg. 11 (1977) 817-832. [15] H. Matthies, Problems in plasticity and their finite element approximation, Ph.D. Thesis, Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA, 1978. [16] P. Suquet, Existence et regularit6 de solutions des equations de al plasticit6 parfaite, (Existence and regularity of solutions to perfect plasticity equations), These de 3e Cycle (Universit6 de Paris, 1978). [17] G. Duva,~Jt and J.L. Lions, Les Inequations en Mechanique et en Physique (Dunod, Paris, 1972). [181 C. Johnson, Existency theorems for plasticity problems, J. Math. lures Appl. 55 (1976) 431-444. [19] C. Johnson, On finite element methods for plasticity problems, Numer. Math. 26 (1976) 79-84. [2O] C. Johnson, A mixed finite element for plasticity, SIAM J. Numer. Anal. 14 (1977) 575-583. [21] C. J ofinson. On plasticity with hardening, J. Math. Anal. Appl. 62 (1978) 325-336. [22] H. Matthies and G. Strang, The solution of nonlinear finite element equations, Intemat. J. Numer. Methods Engrg. 14 (11) (1979) 1613-1626. [23] R. Temam, Mathematical problems in plasticity (Gauthier-Villars, Paris, 1985) (Translation of 1983 French original edition). [24] R. Hill, The Mathematical Theory of Plasticity (Oxford University Press, Oxford, 1950, last ed. 1983): [25] G. Maler, Quadratic programming and theory of elastic-perfectly plastic structures, Meccanica 3 (1968) 265-273. [26] O. Maier, A matrix structural theory of piecewise linear elastoplasticity with interacting yield planes, Meccanica 5 (1970) 54-56. [27] J.J. Moreau, Application of convex analysis to the treatment of elastoplastic systems, in: P. Germain and B. Nayroles, eds., Application~ of Methods of Functional Analysis to Problems in Mechanics (Springer, Berlin, 1976). [28] T. Pian and Sumihara, Rat,,~nal approach for assumed stress finite elements, Internat. J. Numer. Methods. Engrg. 20 (9) (1984) 1685--1695. [29] J. Mandel, Contribution theorique al'6tude de i'ecrouissage et des lois de I'~couleme~,t plastique, Proceedings llth International Congress Appl. Mech. (1964) 502-509. [3o1 J. Lubliner, A maximum-dissipation principle in generalized plasticity, Acta Mech. 52 (1984) 225-237. [31] J. Lubliner, Normality rules in large-deformation plasticity, Mech. of Mater. 5 (1986) 29-34. [32] J.C. Simo, A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decompositiu.~: Part I. Continuum formulation, Comput. Methods Appl. Mech. Engrg. 66 (2) (1986) 199-219. [33] J.C. Simo, A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part II. Computational aspects, Comput. Methods App|. Mech. Engrg. 68 (1) (1988) 1-31. P41 J.C. Simo and T. Honein, Variational Formulation, Discrete Conservation Laws and Path-Domain Independent Integrals for Elasto-Vis~plasticity, to appear. [351 D.G. Luenberger, Linear and Nonlinear Programming (Addison-Wesley, Reading, 1984). [36] G. Strang, Introduction to Applied Mathematics (Wellesley-Cambridge Press, Wellesley, 1986). [37] P. Pinsky, A finite element formulation for elastoplasticity based on a three field variational principle, Comput. Methods Appl. Mech. Engrg. 61 (1) (1987)41-60. [38] J.C. Simo and T.J.R. Hughes, On the variational foundations of assumed strain methods, J. Appl. Mech. 53 (1) (1986) 51-54. [39] J.C. Simo, J.G. Kennedy and R.L. Taylor, Complementary mixed finite element formulations for elastoplasticity, Report UCB/SEMM-87/09, University of California, Berkeley, CA (1987). [4O] O.C. Zienkiewicz, The Finite Element Method, Third Edition (McGraw-Hill, London, 1977). [41] T. Belytschko and W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 54 (1985) 279-301.

206

Y.C. Simo et al., Complementary mixed finite element formulatior~

[42] O.C. Zienkiewicz, R.L. Taylor and S. Nakazawa, The patch test for mixed methods, Internat. J. Numer. Methods Engrg. 23 (1986) 1873-1883. [43] J.C. Simo, D.D. Fox and M.S. Rifai, On a geometrically exact shell model. Part II. The linear theory; Computational aspects. Comput. Methods Appl. Mech. Engrg. 73 (1) (1989) 53-92. [44] J. Mandel, Thermodynamics and plasticity, in: J.J. Delgado et al., eds., Foundations of Continuum Thermodynamics (MacMillan, New York, 1974). [45] G. Strang, H. Matthies and R. Temam, Mathematical and computational methods in plasticity, in: S. Nemat-Nasser ed., Variational Methods in the Mechanics of Solids, (Pergamon, Oxford, 1980). [46] R. Temam and G. Strang, Functions of bounded deformation, Arch. Rational Mech. Anal. 75 (1980) 7-21.