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

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

(113a)

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 "=

{-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)

(2.4b)

- 8X(~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

--

= S<x)g, g e Rm}.

(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 (*) +

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) :----"'.~(k1~
, n + l - - ( ' ] n~(k)
+l '
[~(k)]-l.

:=

--"

(3.7)

fo est[O2x + T02o~
-

"'

-- 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,,+,)- ~

(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 <k). =O

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,sx) 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,

gElat

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

~.

glact

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

R(')"

gelat

(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

= [ 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[32 ] ]

(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

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

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

10

17

(ti -- 0 . 0 3 )

(ti -- 0.15)

(ti - 2.65)

( ~ = 6.15)

0 . 8 9 7 e + 00

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 <

< 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

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

4
5

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:

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

2~

'

'

'

'

2O
$ fSSS

.S . . . . .

" ....

" " " " "

/"
./

/
/
i
]

AssumedStress FormulationI l l
DisplacementFormulation. . . . .

i
i
0o

i
i

;0

15

20

.,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<p.+,)
".+,~"i" : =

f,

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,

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

_t t

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

(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

(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.