Limit Analysis of Collapse States

Limit Analysis of
Collapse States
Edmund Christiansen
Department of Mathematics and Computer Science
Odense University
Campusvej 55
DK-5230 Odense
Denmark
HANDBOOK OF NUMERICAL ANALYSIS, VOL. IV

Finite Element Methods (Part 2)-Numerical Methods for Solids (Part 2)
Edited by P.G. Ciarlet and J.L. Lions
( 1996 Elsevier Science B.V. All rights reserved
193
Contents
PREFACE 197
CHAPTER I. Introduction 199

1. The physical problem 199
2. The mathematical model 202
3. Remarks on the development of limit analysis 206
CHAPTER II. The Duality Problem 209
4. Why so different from elasticity? 209

5. Existence of the collapse solution 213
5.1. The stress space X 216
5.2. The flow space U 217
5.3. The duality of limit analysis 219
5.4. The reduced problem 225
6. Generalizations and variations 227
6.1. Boundary flow constrained to a subspace 227
6.2. Presence of a pre-load 228
6.3. Bounded yield set 228
6.4. Yield conditions for granular materials 229
7. Green's formula in the collapse state 230
CHAPTER III. Discretization by Finite Elements 235

8. Formulation of the discrete problem 235
9. Elements for divergence-free flow 244
10. Bounds and convergence 247
CHAPTER IV. Solution of the Discrete Problem 255

11. General remarks 255
12. Linear programming methods 259
12.1. Solution with the simplex method 262
12.2. Solution with the primal affine scaling algorithm 265
12.3. Solution with the dual affine scaling algorithm 267
13. Convex programming methods 272
14. Divergence-free elements 276
14.1. Solution of the primal problem 277
14.2. Solution of the dual problem 278
15. Concluding remarks 282
15.1. On the solution to the test problem 284
195
196 E. Christiansen
CHAPTER V. Limit Analysis for Plates 285

16. The continuous problem 285
17. Duality of limit analysis 288
18. Discretization of the plate problem 291
19. Solution of the discrete problem 295
REFERENCES 303
LIST OF SYMBOLS 307
SUBJECT INDEX 311

Preface
Limit analysis is the analysis of the simplest possible model for the collapse of
a material subject to a static load distribution. The model involves a duality
problem in infinite dimensional mathematical programming, formally similar to
the min-max problem of game theory. It involves:
(a) Mathematical analysis within the framework of Sobolev spaces.
(b) Numerical analysis based on the finite element method.
(c) Nonlinear optimization methods to solve the discrete problem formu-
lated in (b).
The object of limit analysis is to improve the understanding of plastic materi-
als, or rather, our models for plastic material behaviour, specifically equilibrium
equation and yield condition. From a numerical point of view the ambition is to
make limit analysis useful for applications as a supplement to the elasticity anal-
ysis, which is now a standard tool in engineering. The model of linear elasticity
results in linear elliptic boundary value problems, which are so well behaved,
both mathematically and numerically, that realistic problems can be handled in
great detail. However, even small forces may cause unbounded elastic stresses,
clearly in conflict with the physical reality, and the very question of collapse is
beyond the model of elasticity.
As is the case in elasticity theory, the numerical analysis is closely related
to the mathematical analysis of the model: Both employ the theory for partial
differential equations using Sobolev spaces. However, since one of the funda-
mental concepts for plastic materials, the yield condition, is given as a pointwise
condition, i.e., in terms of L , the spaces involved are different from those in
elasticity theory. The minimum energy principle in elasticity is replaced by a
saddle point problem for a bilinear form in stresses and plastic flow. Unique-
ness and regularity, known from the theory of elliptic boundary value problems,
do not hold. For the numerical analysis this means that mixed finite elements
must be used, and there is a delicate balance between the discrete spaces for
stresses and plastic flow. Also the convergence results for the numerical solution
are weaker than in the elliptic case.
The solution methods for large nonlinear optimization problems are based on
a sequence of steps, each of which involves the sort of sparse numerical linear
algebra that is used in elasticity problems. Also in this respect limit analysis
stands on the shoulders of the elastic model.
The conclusion of this presentation is that limit analysis is ready for applica-
tions. Realistic problems can be solved. However, we need more computational
197
198 E. Christiansen
experience on such problems to help us choose between the possible combina-

tions of finite element spaces and several possible optimization methods for the
discrete problem.
CHAPTER I
Introduction
1. The physical problem
The plastic analysis of solids must be seen as a supplement to the elastic analysis.
The linear elastic model results in elliptic boundary value problems, which are
essentially well behaved from a mathematical and numerical viewpoint, whereas
the models of plasticity are intrinsically harder. The plastic analysis will proba-
bly always lag behind the elastic analysis. However, in view of the quality which
practical structural analysis has reached, it is apparent that the linear elastic
model is insufficient: In standard applications even small forces may cause un-
bounded elastic stresses and thus require special treatment and, by definition,
the problem of collapse is beyond the model of elasticity. It is a safe, and not
at all new, prediction that nonlinear analysis of solids and structures will move
from the research table to practical engineering applications, but since the sur-
vey by CoHN [1979] most progress has still been at the research level. One goal
still to be achieved is that the finite element programs, used every day by the
working engineer, will be able to perform a plastic collapse analysis in addition
to the linear elastic analysis, which is now standard.
Another immediate goal for computations in limit analysis is to supply a feed-
back to the model. The development of good yield criteria for different materials
has been impeded by the lack of consequential computations to compare with
experiment. It is impossible to refine a model without the ability to compute the
behaviour, which the model predicts. With the recent progress in mathematical
optimization and computing power such computations have become possible.
Limit analysis is the simplest way to perform plastic analysis of a solid. The
material is assumed to be rigid, perfectly plastic, which means that small forces
will be neutralized by stresses in the material without any deformation taking
place. These stresses are not uniquely determined by the equilibrium equation,
which in the rigid plastic case is non-elliptic. If the forces are too large, so that
the material cannot sustain stresses to neutralize them, a plastic, i.e., permanent,
deformation will occur, and the material will flow for as long as the forces remain
and geometric changes can be ignored. Hence the deformation of the material
is described by flow or displacement rates, and not by displacements. We shall
not try to justify the rigid, perfectly plastic model here, but refer to, e.g., SAVE
[1979].
199
200 E. Christiansen CHAVrER
F_
ox
FIG. 1.1. Test problem in plane strain.
Within the model of rigid, perfectly plastic materials the problem of limit
analysis can be formulated as follows: Consider a fixed force distribution on
the continuum, consisting of volume and surface forces. What is the maximum
multiple of this force distribution, which the solid can sustain without collapsing,
i.e., without plastic flow occurring? And, at the moment of collapse, what are
the fields of stresses and flow in the material? In particular, what is the plastified
region, i.e., the region of the continuum where the stresses are at the limit, and
where local deformation may occur?
Clearly the model of limit analysis is based on a very idealized situation
corresponding to a rigid, perfectly plastic material and static or slowly increasing
loads. It is the simplest model, which makes it possible to analyse quantitatively
the plastic collapse of a material. For a justification of the model of limit analysis
we refer to, e.g., KALISZKY [1975] or SAVE [1979].
An alternative approach is the incremental analysis in an elastic, perfectly
plastic model, but as pointed out by SAVE [1979, Section 3.1], the validity of
both models requires that deformations for loads smaller than the limit load do
not alter the initial geometry of the material, and in that case the incremental
method offers no improvement. The model of limit analysis describes only the
fields for stresses and flow at the collapse moment, i.e., when plastic deformation
first occurs. It does not give information about these fields, when deformation
has changed the initial geometry. It is, in a way, a flash photo of the collapse
moment.
In order to visualize the difference between the elastic and the plastic analysis
of a solid we shall consider a specific example. It will be used as a reference
test problem later.
EXAMPLE 1.1 (A test problem in plane strain). Consider a homogeneous and iso-
tropic block of material, e.g., a metal, as shown in Fig. 1.1. The block is assumed
very wide in the z-direction, and quadratic in the x -y plane. All forces and
deformations take place in the x - y plane, reducing the problem to 2 space di-
mensions. This is the plane strain model. The external load consists of a uniform
tensile force in the x-direction applied at the ends x = +1. To make it more in-
SECTION 1 Introduction 201
FIG. 1.2. Elastic and plastic solutions for the test problem in Fig. 1.1.
teresting we make symmetric thin cuts at x = 0 half way to the mid-plane y = 0

(see Fig. 1.1). Due to symmetry the solution need only be computed in the first
quadrant, and only this part is shown in Fig. 1.2.
Figure 1.2 shows on the left side the elastic deformation with the tensile force
suitably scaled and, on the right side, the plastic limit solution, both computed
using the finite element method. The plastic solution needs explanation: The
collapse flow field u has been computed at the grid points and multiplied
by a time constant. The resulting displacement field is then plotted. Since the
collapse velocity is determined only at the instant of collapse, this displacement
may never be seen, but there is general agreement that this is the best way
to visualize the flow field. It gives a better idea of the local deformations at
collapse than the velocity field itself does. We see that the elastic deformation
is smooth and distributed over the solid, while the plastic flow field is of a
different nature. We can identify two rigid regions separated by a relatively
narrow interface, which seems to contain the entire collapse mechanism. The
stresses are not shown in the figure, although they have been computed as part
of the solution. The elastic stress field may be computed from the deformation.
The exact stress field is known to have a singularity at the tip of the cuts and
to be smooth everywhere else. The plastic collapse stress tensor o'(x) is much
harder to find. It is not uniquely determined outside the plastified region, i.e.,
the part of the material, where there is deformation taking place. In the chapter
on computation we shall see several ways to visualize the collapse stresses, also
for the problem in this example.
The qualitative difference between the elastic and the plastic behaviour of
a solid is reflected in the mathematics and the numerics of the problem. The
plastic model does not have the nice mathematical properties of the elastic
model, and numerically the solution is much harder to compute. All useful
202 E. Christiansen CHAP'TER I
numerical solution methods in limit analysis for a continuous medium known

to this author are based on discretizations of a variational principle using the
finite element method. As is the case for elliptic problems the numerical analysis
of the finite element method applied to limit analysis is closely related to the
mathematical analysis of the continuous problem, so we have to start there.
After formulating the continuous problem we shall attempt to point out some
of the steps, mainly analysis of corresponding formulations for discrete problems
and structures, that have led the way.
2. The mathematical model
We first formulate the problem in its most basic form. Subsequently it will be
indicated how the analysis applies also to more general boundary conditions and
forces. We concentrate on the mathematical formulation to be analysed later.
For a closer link to the mechanics we refer to, e.g., KACHANOV [1971], MARTIN
[1975, Chapter 9] or ZAVELANI-ROSSI [1979].
Let (2 denote the bounded domain in space occupied by a rigid, perfectly
plastic continuum. The boundary A02 consists of two parts, S and T: 02 = S u T.
On S the boundary is fixed, i.e., u = 0 on S, where u = u (x) denotes the plastic
displacement rate at each point x in 12, the closure of 2. On T there are given
surface forces g(x), x C T, (possibly zero), and inside the material there is given
a field of volume forces f (x), x c n (e.g., zero or gravity forces). This pair of
force fields (f, g) is kept fixed. The problem of limit analysis is to find the limit
multiplier, i.e., the upper limit A*of the set of real multipliers A,such that the
body of material can "carry" the force distribution (Af, Ag).
Associated with u is the strain rate tensor = (u) defined by
ij(u) =21 xj duxi+)+ (2.1)

This is the notation normally used for the strains themselves, and not their time
derivatives. Since displacements and strains do not occur in limit analysis (only
the corresponding rates), we ask the reader to accept this notation. It can also
be justified by the fact that the time scale does not occur in the model.
Given a force distribution (f, g) a virtual displacement rate u in the contin-
uum will be associated with the external work rate
F(u) = If u + f g-u. (2.2)
fl T
The stresses inside the material are given by the symmetric stress tensor or =
or (x). The classical form of the equation of equilibrium between or and the force
distribution (f,g) defined above is
-V.o=
=f in, (23)
(2.3)
v.r =g on T,
SETION 2 Introduction 203
where V denotes the divergence operator, here applied to the columns of or,
and v is the outward normal to 12. The coordinate form of Eq. (2.3) is
3 otij
-E ix =fi inn fori=1,2,3,
j=l
3
vjoij = gi on T for i = 1,2,3.
j=3
In rigid plasticity there is a priori no relation given between the strains and
the stresses. In the collapse state, however, such a relation will be deduced.
We shall always interpret Eq. (2.3) as being equivalent to the more precise
variational form Eq. (2.6), which we now derive: To a given pair of stresses -
and displacement rates u we associate the internal work rate or energy dissipa-
tion rate
a(r, ) = f(oa, E(U))
f=ul~rijCj(u) (2.4)
n 1,J
=- (V.r).u + f (v.r)-u. (2.5)

2 T
For the moment we ignore technical details and assume that or and u are
sufficiently smooth. The second equality then follows from Green's formula.
The equation of virtual work (rate) states that the external work rate defined
by Eq. (2.2) and the internal work rate defined by Eq. (2.4) and Eq. (2.5) must
be equal for all admissible u:
a(or, u) = F(u) for all u satisfying the boundary condition on S. (2.6)
Equating Eq. (2.2) and Eq. (2.5) we see that the classical form of the equilibrium
equation Eq. (2.3) is equivalent to the generalized form Eq. (2.6), at least for
sufficiently smooth fields ff and u. Note that the boundary condition on the
stresses (the second part of Eq. (2.3)) is only implicit. In the terminology of
elliptic boundary value problems the boundary condition on the displacement
field is an essential boundary condition, while the condition on the stress field
is a naturalboundary condition as in the Neumann problem.
There is a limit to the strength of a plastic material. In this model it means
that the stress tensor must satisfy a yield condition in addition to the equilibrium
equation Eq. (2.6). The yield condition is a pointwise constraint of the form
o(x) E B(x) for all x 2. (2.7)
If the material is assumed to be homogeneous B(x) is independent of x and
will be denoted B. For the mathematical theory it is no complication to allow
204 E. Christiansen CHAPTER I
the yield condition to depend on x. About the set of admissible stresses B(x) at
each point we make the following assumptions:
(a) There exists e > 0, independent of x, such that ii lIiji A E implies oa c
B (x).
(b) B (x) is a convex subset of the space of stress tensors.
(c) B(x) is a closed subset of W9.
These three conditions are crucial for the mathematical and numerical analysis
of the model. The theory of convex optimization can be applied to prove that the
concept of a collapse state is well defined. (a) states that the material can sustain
any stress field, if only the stresses are sufficiently small at every point, i.e., there
are no points or directions of "zero-strength" for the solid. (b) corresponds
to adequate physical assumptions, while condition (c) is purely mathematical
without physical significance. It states that a limit of admissible stresses is itself
admissible.
It is important to note that B(x) is not assumed to be bounded. In three-
dimensional problems and in the two-dimensional plane strain model B is un-
bounded for many materials. It is bounded in plane stress and for plate bending
problems. An unbounded yield set means that certain types of stresses may be
arbitrarily large. In order to understand this, let us consider a specific classical
yield condition:
The Von Mises yield condition in three space dimensions for a homogeneous
material can be written as follows
(11t - 22) 2 + (22 - 033) 2 + (033 - 11)2 + 6(o122 + '223 + o'321) < 2oO2, (2.8)
co being the yield stress in simple one-dimensional tension. Define the trace of
a tensor t as
tr(r) = E ii
and let I denote the identity tensor (ij). Then we define the deviator of o as
the tensor
a = - tr(O) I, (2.9)
where N is the spatial dimension (2 or 3) of the problem. Clearly o satisfies the

yield condition Eq. (2.8), if and only if orD does, but since 0oD has trace zero,
Eq. (2.8) puts a bound on its single components. Typically the yield condition
will be of the form
rD(x) E B(x) for all x E , (2.10)
where B(x) is a bounded subset of 9. Such yield conditions, which can be
formulated exclusively in terms of the deviator rD of the stresses, are insensitive
to the addition of tensors of the form to = PI, where is any scalar function
defined on 12. Hence the material can sustain any such tensor field, no matter
how large .4 is. These tensors correspond to hydrostatic pressure in the material,
so using a yield condition of this form is equivalent to assuming that the material
SECTION 2 Introduction 205
does not yield under hydrostatic pressure. We shall see later that there is a
corresponding dual property of the plastic flow u: It is divergence free, or
incompressible: V.u = 0. We do not have to require this property. It follows
automatically, but it is important to be aware of this dual relationship in the
mathematical and numerical analysis of the problem.
We can now formulate the static principle of limit analysis: Find the limit
multiplier A* of a given force distribution (f,g) that the solid can "carry". In
other words, we want to maximize A subject to the constraints that there exists
a stress field ao-, which is in equilibrium with (Af, Ag) and which satisfies the
yield condition Eq. (2.7). With the above notation this may be written
A* = max{A 3a Eo-
K: a(o-,u) = AF(u) for all u with u = O on S},
(2.11)
where K denotes the set of stress fields that satisfy the pointwise yield condition
at every point. Details will be made precise later.
Since GREENBERG and PRAGER [1952] it is customary to term a stress field
or statically admissible with a given force distribution (f,g), if it satisfies the
equilibrium equation Eq. (2.3) or equivalently Eq. (2.6). If, in addition, the
stress field satisfies the yield condition, it is called safe. Similarly a velocity field
u is termed kinematically admissible, if it satisfies the boundary condition u = 0
on S.
The kinematic principle of limit analysis is formulated in terms of the plastic
flow: We define the energy dissipationrate associated with u as follows
D(u) = max a(o-, u). (2.12)
o-EK
The kinematic principle states that the limit multiplier A* is the minimum of
D(u) taken over all u with external work rate normalized to unity:
A*= min D(u). (2.13)
F(u)=1
The physical reasoning behind the kinematic principle is as follows: If A > A*,
then there is a u with F(u) = 1, such that A > D(u). This implies that for all
stress tensors satisfying the yield condition we have that A = AF(u) > a(o-, u).
Hence, for this particular flow field a, no admissible stress tensor can neutral-
ize the work rate of the external forces (Af, Ag), and the material will yield.
Conversely, if A < A*, no such u exists.
The static and kinematic principles are formulated independently of each
other. However, we shall see that they are dual mathematical programming
problems. In particular, Eq. (2.11) and Eq. (2.13) give the same value for A*, and
the extremal values are attained for some stress field -*C K for the problem
Eq. (2.11) and some admissible velocity field au* for the problem Eq. (2.13).
The pair (o*, u*) is a saddle point for the bilinear form a(o-, u) and can be
interpreted as the fields for stress and flow in the collapse state. More precisely,
206 E. Christiansen CHAPTER
the following holds for all or satisfying the yield condition and all u satisfying
the boundary conditions and normalized to F(u) = 1:
a(o, u*) < A* = a(&r, u*) < a(o*, u). (2.14)
This variational principle, formulated for discrete structures or continuous
solids, is as fundamental to limit analysis as the minimum energy principle is to
the elliptic problems of elasticity.
In order to simplify the notation we have made a couple of assumptions
in the above formulation: The boundary condition on u either prescribes u
completely, or leaves it free. A more general boundary condition may allow the
surface to move only in certain directions, e.g., tangential to the surface, and
be subject to surface forces in these directions. Also, there may be a fixed force
distribution such as gravity, sometimes referred to as a pre-load or a dead load,
on the continuum in addition to the load (f, g), for which we want to find the
maximal multiplier. The analysis of the next chapter is readily modified to cover
these cases as well.
3. Remarks on the development of limit analysis
Like so many other numerical methods limit analysis has its origin in clever
hand calculations. Suppose that for some value A we are able to guess or con-
struct a stress field, which satisfies the yield condition and is in equilibrium with
the forces (Af, Ag). Then, by the static principle, A is a lower bound for the limit
multiplier A*. Similarly, by the kinematic principle, for any flow field u, which
satisfies the boundary condition and has the work rate normalized to F(u) = 1,
D(u) will provide an upper bound of A*. Thus the two principles can be used
to establish bounds for the collapse multiplier, simply by guessing stress or flow
fields. This approach is frequently referred to as lower bound methods, respec-
tively upper bound methods. In the lower bound method it may be difficult to
verify the equilibrium equation for the stresses, while the calculation of D(u)
is the hard part when using an upper bound method. If the yield condition is
of the form Eq. (2.10), then D(u) is equal to +oo unless u is divergence free.
Only in rare cases is it possible to determine the exact limit multiplier this way
by getting the same value for the upper and lower bound.
The static principle was formulated in a special case already by KAZINCZY
[1914] and KisT [1917]. Precise physical arguments for both the static and kine-
matic principles (often referred to as "theorems" in the literature of mechanics)
based on the convexity of the yield set were given in GVOZDEV [1938]. More
details about the early history of limit analysis are given in NEAL [1956, Chap-
ter 3]. It must be emphasized that all early formulations and analysis are given
for simple structures and not for a general continuous medium.
In 1951-52 the theory expanded in two directions: In DRUCKER, PRAGER and
GREENBERG [1952] the static and kinematic principles were established for con-
tinuous media, and in CHARNES and GREENBERG [1951] the connection between
SECnON 3 Introduction 207
the static and kinematic principles for trusses on one side and linear program-
ming duality on the other side was pointed out. It was later proved for frames
in CHARNES, LEMKE and ZIENKIEWICZ [1959]. The realization of the connection
to duality theory seemed to surface simultaneously in several places and im-
mediately made a major impact. The application of linear programming and
the simplex method to limit analysis made it possible in general to find the ex-
act limit multiplier for finitely determined structures. This started an explosive
development in the plastic analysis of structures as reflected in a multitude of
papers. Many problems for simple structures can adequately be formulated with
piecewise linear yield conditions, and nonlinear optimization algorithms was not
ready for use in this context yet.
Limit analysis for continuous media needed more time to develop. Compu-
tations had to await the success of the finite element method. It appears that
the first descriptions of the finite element method applied to limit analysis was
ARGYRIS [1967, Chapter IV] (bound methods for discrete structures), HAYES
and MARCAL [1967] (upper bound method in plane stress), and HODGE and BE-
LYTSCHKO [1968] (upper and lower bound methods for plates). Finite element
function spaces were used to represent stresses or flow field, and the associated
computed bound was optimized over the discrete space. Thus the exact limit
multiplier for the discretized structure was determined, and this value was then
a bound (lower or upper) for the limit multiplier for the continuous problem.
In both cases nonlinear optimization methods were used, but of course only
for rather coarse elements. Convergence, as the element size tends to zero,
was not discussed. At this time the development again exploded with papers
using the finite element method to discretize continuum problems combined
with convex and linear programming. At first only bound methods were used
with finite element representations of stresses or flow, but not both. In CAPURSO
[1971] and ANDERHEGGEN and KNOPFEL [1972] mixed finite elements were intro-
duced to discretize stresses and flow simultaneously, in both cases using linear
programming, where the duality between the static and the kinematic principle
was established. Without a duality theorem the mathematical foundation of the
mixed method is missing, a point emphasized by CASCIARO and DICARLO [1974,
p. 172].
Mathematical respectability, in form of a proof within a rigorous mathemati-
cal frame of the duality between the static and kinematic principles, took even
longer. NAYROLES [1970] gave a formulation of the problem as a duality prob-
lem in infinite-dimensional mathematical programming. The duality was proved,
and also the existence of a collapse mechanism in the duality between o with
components in L and with components in the dual space (LO)*. In order to
recover u from in this formulation it was necessary to work within the quo-
tient space modulo rigid displacements. The spaces and boundary conditions do
not quite match the physical problem. In the fundamental work DUVAUT and
LIONS [1972] the rigid perfectly plastic problem is introduced as a limit case of
the elasto-plastic model.
MOREAU [1976] and CHRISTIANSEN [1976] deduced from analysis of the duality
208 E. Christiansen CHAPTER I
that the strain rates e (u) must be measures. The space of velocities u with
precisely this property was introduced in MATHIES, STRANG and CHRISTIANSEN
[1979]. It is called the space of bounded deformation and is denoted BD(12).
The properties of BD were analysed in SUQUET [1978, 1979], TEMAM and STRANG
[1980a] and TEMAM [1981b]. These results were used to analyse the mathematical
problem of plasticity in TEMAM [1981a, 1983], TEMAM and STRANG [1980b] and
SUQUET [1981] among others.
In the attempt to prove the existence of collapse flow CHRISTIANSEN [1980a]
and SUQUET [1981] independently concluded that even BD(f2) is not quite large
enough. The collapse flow may have discontinuities along the boundary of 12,
and the associated energy dissipation is not measured if the strain rate tensor is
defined only by distributional derivatives, as in BD(12). A proof of duality and
existence of collapse fields incorporating this is given in CHRISTIANSEN [1986].
CHAPTER II
The Duality Problem
4. Why so different from elasticity?
In this section we try to point out some of the mathematical characteristics of

limit analysis without referring to Sobolev space theory. The purpose is to make
it possible to understand the nature of the collapse solution without knowing
the technical mathematical details.
Figure 1.2 in the introduction is an attempt to illustrate the difference between
elastic and plastic behaviour of a solid. This section is an attempt to point out the
profound theoretical differences between the underlying mathematical models.
In elasticity, both linear and nonlinear, the solution is uniquely determined
as being the displacement field minimizing an energy functional. The collapse
solution in limit analysis is a triple (A*, o,*, u *), where Ao*provides the maximum
in the static principle Eq. (2.11), u* gives the minimum in the kinematic principle
Eq. (2.13), both values being equal to A*. This is the saddle point problem
Eq. (2.14) for the bilinear form a(o', u) in Eq. (2.4). We shall see that existence
holds in a generalized sense only. Uniqueness of the collapse mechanism simply
does not hold in general, either for stresses or for the flow. And at the moment
of collapse there may be large rigid regions of the material with zero strain
rates, a phenomenon not seen in elasticity due to the nature of solutions to
elliptic boundary value problems.
Concerning the existence problem it turns out that u * does not necessarily
exist in the ordinary Sobolev space sense, i.e., as a locally integrable function
in 12, which satisfies the boundary conditions (see, e.g., SHOEMAKER [1979]).
This is illustrated by the following example, which is similar to Example 2.3 in
CHRISTIANSEN [1982].
EXAMPLE 4.1 (An example in antiplane shear). Consider an infinite layer of con-
stant thickness in the y-z-plane between x = -1 and x = 1 as indicated in
Fig. 4.1. Only forces in the y-direction are allowed, and they may only depend
on x. The problem is then completely described by the following scalar functions
on the interval -1 < x 1:
rt = xy =y (X),
u= = Uy(X),
209
210 E. Christiansen CHAPTER
y1.
I
u(-1)=0 I
I.(X) I Tg(l)=a
I
__L_ Ix
a 1
11- 4,4,
I
FIG. 4.1. A problem in antiplane shear.
f = fy = fy(x),
g = gy = gy(x).
The layer is kept fixed at x = -1 by imposing the boundary condition u(-1) = 0.
We now solve the limit analysis problem for this example. The rate of internal
work, Eq. (2.4), is
I I
a(o, u) = oru' dx = - u'u dx + o(1)u(1). (4.1)

-1 1
As yield condition we use o-l < 1. We consider the following external forces:
f(x) -1 for [xI <a,
lo for Ix > a,
g(l) = a,
for some a, 0 < a < 1. Thus one side of the layer at x = -1 is fixed. There is a
uniform gravity-like force pulling down a central section of the layer, -a < x <
a, and we are pulling up the other side of the layer at x = 1. The equilibrium
equation, Eq. (2.6), for the pair (Af, Ag) is
a(a, u) = AF(u) = A (- u dx + au(1)) (4.2)
for all u with u(-1) = 0. Combining Eq. (4.1) and Eq. (4.2) we get the classical
form
(x) = fA for xI < a,
0 for xj > a, (4.3)
r(1)= Aa.
The static principle takes the form
maximize A
subject to Eq. (4.3) and or1 1
SECTION 4 The duality problem 211
The solution is easily seen to be
A*=
a
with the uniquely determined collapse stresses
-1
for -1 x < -a,
o*(x) = x/a for-a x < a,
+1 for a < x 1.
The dual solution u* must satisfy Eq. (2.13) and Eq. (2.14):
A* = a(o-*, u*) = D(u*)
or
J
1
*u* dx = sup fiou*
1
dx= IuIM. (4.4)
I|U*/IjM denotes the norm of u*' as a measure on the closed interval [-1, 1], i.e.,
the dual norm to the continuous functions on [-1,1]. It is equal to the total
variation of u* on the interval [-1, 1]. Note that u*' is equal to the distributional
derivative of u* in the open interval ] - 1,1[ plus a contribution from the end
points of the interval. (In Eq. (4.1) and Eq. (4.4) we use the integral notation
also to denote the formal duality between a measure and a continuous function
on a closed interval.) Hence o* must be equal to the sign of u*', whenever u*' is
nonzero (considered as a measure on [-1, 1]). This implies that u* must satisfy
the following:
u* is decreasing for -1 < x < 0,
au* is constant for -a < x < a, (4.5)
u* is increasing for 0 < x 1.
(The intervals in the condition overlap, because u(-a) must be greater than
or equal to the constant value in the open interval -a < x < a and similar for
u(a).) We see that the loaded region of the layer [-a, a] moves as a solid block.
In addition u* is normalized to satisfy
a
1 = F(u*) = u* dx +au*(l) = a (-2uo +u*(1))

-a
or
u*(1) = 2u0 + -, (4.6)

a
where u0 denotes the constant value of u* on [-a, a]. This leaves room for any
collapse flow, which satisfies Eq. (4.5) and Eq. (4.6) in addition to the boundary
condition u(-1) = 0. Thus there are solutions for each u0 < 0.
212 E. Christiansen CHAPTER II
For a < 1 there are both smooth and non-smooth collapse flows. For a = 1
there is only a one-parameter family of dual solutions, and every collapse flow
has a discontinuity at x = -1, if uo yA 0, or at x = 1, if uo # -a - l , or both. The
discontinuity at x = -1 appears to violate the boundary condition u(-1) = 0,
but the corresponding work rate is included in the internal work rate a(ar, u).
On the other hand, the discontinuity at x = 1 causes a technical problem: The
corresponding work rate is ignored in a(ao, u), if u' is the derivative of u in the
distributional sense. We must define u in such a way that a discontinuity along a
loaded part of the boundary is significant in the energy dissipation rate. Hence
the collapse flow must be defined as a pair (u0, uT), where u is a function
defined on the loaded part of the boundary, and such that Eq. (4.1) holds for all
stresses. More generally Un and uT are related through the equation of virtual
work (rate) Eq. (2.4) and Green's formula Eq. (2.5).
An analogous example in plane stress is given in CHRISTIANSEN [1982, Sec-
tion 3]. These examples show several important and unpleasant aspects of limit
analysis, in addition to being some of the rare cases, which can be solved by
hand. The strains in the collapse state are measures in the closed interval [-1, 1],
as emphasized by MOREAU [1976]. This means that the flow field is a function
of bounded variation. In several dimensions this generalizes to the space BD
of bounded deformation introduced in MATrHIES, STRANG and CHRISTIANSEN
[1979]. For a vector field u to be of bounded deformation only the strain rates
sij(u) are required to be measures, not all partial derivatives of u. We have
also established by simple calculations that the collapse flow does not always
exist as an integrable function, but may have values on the loaded part of the
boundary, which are neither given by boundary conditions, nor are they the
limit of values from the inner of the domain occupied by the solid. If u has
a discontinuity along the part of the boundary, where the boundary condition
u = us is given, then the associated work is included in the internal work rate
Eq. (2.5) by adding a term of the form
J (u - Us).(0.-). (4.7)
This is the relaxation method used in TEMAM and STRANG [1980b]. If, however,
the discontinuity occurs at the loaded part of the boundary, the flow at the
boundary, denoted by uT in the example, must be treated as a part of the un-
known solution. It is very satisfactory that the proof of existence of a collapse
solution (in the next section) produces the pair (u 0 ,uT) as a purely mathe-
matical object. The existence of the collapse flow as such a pair is suggested
in CHRISTIANSEN [1976, 1980a] and with a slightly different notation in SUQUET
[1980, 1981]. The existence problem is pinpointed and dealt with in a different
way in TEMAM and STRANG [1980b] and TEMAM [1981a].
Another difficulty, not illustrated by the example, is that the stress tensor
need not be continuous in 2. Since the strain rates are measures, this means
that the expression in Eq. (2.4) for the internal work rate is not necessarily well
SETION 5 The duality problem 213
defined. As we shall see in the next section, the expression Eq. (2.5) is always
well defined.
The example also illustrates non-uniqueness of the collapse flow. An example
showing non-uniqueness of the collapse stresses is given in CHRISTIANSEN [1982,
Example 3.1]. We see that the solution contains at least one rigid block of
material without deformation, and there may be discontinuities. In some cases
arbitrarily small perturbations (from a = 1 to a < 1) may change a problem with
only discontinuous solutions to a problem which also has infinitely many smooth
solutions.
From a mathematical point of view the above example is one-dimensional.
True two-dimensional antiplane shear problems are studied in STRANG [1979a,b],
and in STRANG and TEMAM [1982]. The solutions are similar to those above. Of
particular interest is the following example suggested by Robert V Kohn.
EXAMPLE 4.2 (Two dimensional antiplaneshear). Consider an infinite pipe of
constant square cross section in the x-y-plane: 2 = {(x,y) I lxI 1, Iyl 1}.
The yield condition is 11
oll < 1, where ot = (z, ryz). Apply zero volume forces
f = f(x, y) = 0 in n and the following surface forces:
1 y2, X = 1,
2
g= gz = -i+y , x=-l,
0, y = 1.
The limit multiplier is A*= 1. The collapse flow is given by u = = inside
12, while the flow on the boundary uT consists of a -function at (1,0) and a
8-function of opposite sign at (-1, 0). Hence we must allow 8-distributions in
the surface flow.
The two-dimensional problem in antiplane shear has been analysed in more
general terms in DEMENGEL [1989, Section 2.1]. In particular, conditions are
given, under which the collapse flow will exist in the space BD(Q2).
Mathematically it is a complication that solutions do not exist in standard
Sobolev spaces. From a numerical point of view it is bound to cause conver-
gence problems that neither uniqueness nor regularity of the collapse fields can
be expected. If the collapse flow involves discontinuities along the boundary,
finite element spaces with the same property will give more accurate solutions.
The fact that the solution may contain rigid blocks suggests that adaptive finite
elements should be used. On the other hand, the nature of the solution may ex-
plain the early success of bound methods: If the rigid regions can be estimated
or guessed, then good bounds may be obtained, even with hand calculations.
5. Existence of the collapse solution
In this section we give a precise mathematical frame for the duality problem of
limit analysis. The existence theorem is the one given in CHRISTIANSEN [1986],
but the presentation is somewhat different. The formulation is also closely re-
lated to the formulation in TEMAM and STRANG [1980b]. The static formulations
214 E. Christiansen CHAPTER 11
are equivalent. The two formulations give the same limit multiplier and the ex-
istence of a collapse stress tensor. The main difference is our generalization of
the velocity field to a pair (u n, u T), which is significant for the existence of a
collapse velocity field. In SUQUET [1980, 1981] UT is introduced as the limit of a
flow defined outside D.
12 is assumed to be a bounded domain in N, N = 3 or N = 2, with Lipschitz
continuous boundary 0 as defined in NEAS [1967]. Then the outward normal
v to 12 exists almost everywhere in an. Let S and T be open disjoint subsets
of 012, such that S U T = Afl. Assume that their common boundary is Lipschitz
continuous in local coordinates on 012. (For N = 3, S and T meet in a Lipschitz
continuous curve. For N = 2 they meet in a finite number of points.) On S the
flow u is prescribed, while T is free and subject to loading (see below).
The first step is to specify the spaces for stresses and plastic flow. The min-
imum problem of elasticity results in Hilbert space analysis within the frame
of Sobolev spaces with p = 2. In limit analysis the duality between or and u
is non-symmetric, and there is nothing canonical about the choice p = 2. The
starting point is the following facts:
* The yield condition on the stresses is given in terms of L .
* The strain rates are bounded measures in 12.
* The internal work rate is given by Eq. (2.4) and Eq. (2.5).
* The equilibrium equation, Eq. (2.3) or Eq. (2.6), must be satisfied withf and
g in appropriate spaces.
As we shall see subsequently, these facts, combined with the results in TEMAM
and STRANG [1980a] and TEMAM [1981b], settle the question of which Sobolev
spaces to work with in limit analysis. Within these spaces the duality between the
static and kinematic principles holds, and the saddle point in Eq. (2.14) exists.
In the previous section it was demonstrated that the distributional derivatives
of the flow u are insufficient to describe the total energy dissipation. However,
it is still preferable to reserve the notation (u) for the strain rates in the usual
distributional sense. Hence formula (2.1) is specified as follows: For a locally
integrable vector field u in 12 the associated strain rate tensor (u) is defined
to have the following distributional derivatives as its components
ij I4uix + aixi) (5.1)
Let M(12) denote the space of bounded measures in 12, i.e., the distributions,
which are continuous with respect to the maximum norm on the test space
CO(12). The space of vector fields of bounded deformation BD(12) introduced
in MATrHIES, STRANG and CHRISTIANSEN [1979] is defined as
BD(2) = {u [L1(12)]N ] sij(u) M(Q), 1 < i,j < N}. (5.2)

BD(12) is independent of the choice of coordinate system. For a coordinate-free
definition, see TEMAM and STRANG [1980a].
BD(Q) is a Banach space with the norm

N N
IIUIIBD() = E-
l11uiILg() + E IIij(u)iM(). (5.3)
i=l i,j=l
The following essential properties about BD(12) are proved in TEMAM and
STRANG [1980a]. There the results are formulated with the assumption that 0/2
is of class C1 , but as pointed out in CHRISTIANSEN [1986] the proofs given are also
valid when a2 is assumed to be Lipschitz continuous only. The distinction is
of importance, since objects occupying geometries like, e.g., a cube are obvious
applications.
THEOREM 5.1 (The trace theorem for BD). Assume that f1 C RN is bounded
with Lipschitz continuous boundary. There exists a unique continuous linearop-
erator y from BD(12) onto [L' (a2)]N, such that
y(u) = ulI a1 Vu C BD(12) n [C ( f2 )]N . (5.4)
For each i, j, and for every p e C' (n) the following generalized Green'sformula
holds
Ja +Ui +2 (5.5)
where v is the unit outward normal on 02, and yi(u) is the ith component of
y(u).
THEOREM 5.2 (Imbedding theorem for BD). Assume that 2 C RN is bounded

with Lipschitz continuous boundary. Let q = N/(N - 1). Then BD(12) is con-
tinuously imbedded in [Lq(2)]N. In other words:
BD(n2) C [L(12)]N (5.6)
and
N
>E IIUillq(n) < CIIUIIBD(n) Vu C BD(2). (5.7)
i=1
THEOREM 5.3. Assume that 12 C RN is bounded with Lipschitz continuous

boundary. If the N components of u are distributionsin 12, and if eij(u) M (12)
for 1 < i, j < N, then u belongs to BD(2).
For the proofs of the above three theorems see TEMAM and STRANG [1980a].
We now define the stress space for the main case with unbounded yield set.
Modifications for the (simpler) case with bounded yield set will be described
later.
216 E. Christiansen CHAPrE II
5.1. The stress space X
Let p be fixed in the interval N < p < oo. The space I of stress tensors consists
of all symmetric N x N tensors or with the following properties
or [L (l)] N2, (5.8)
or E [L()] N2,
(5.9)
V.or C [LP(f2)]N, (5.10)
N .
Pa E [WI-l/PP(T)] (5.11)
Recall that orD is the deviator of or defined by Eq. (2.9). Z is equipped with
the norm
Ilorll = i[IDIIL(n) + IIVoIILP(n) + IIV'oW-llpPP(T)- (5.12)
The yield condition is assumed to be of the form
aD(x) B(x) x E f2,

2
B (x) c RN closed, bounded uniformly in x and convex for all x C 1.
(5.13)
B (x) must have non-empty interior in the uniform sense that there exists a a > 0
independent of x, such that the following implication holds:
Iljl < Vi,j or E B(x). (5.14)
We also assume that B(x) can be defined by inequalities which are piecewise
continuous in x.
Yield conditions of this form cover standard cases, such as the Von Mises and
Tresca conditions in three-dimensional problems and in plane strain, but not
plane stress problems or conditions such as the Coulomb yield condition for
granular materials. The modifications for these cases are discussed in the next
section.
REMARK 5.1. The choice of Ll in Eq. (5.8) is arbitrary. Any Lq with q in the
interval 1 < q < oo would be equivalent. This follows from Theorem 5.4 below.
REMARK 5.2. The condition in Eq. (5.9) is consistent with the yield condition in
Eq. (2.10).
REMARK 5.3. Condition (5.10) is motivated by the first term in Eq. (2.5) and
Theorem 5.2. p = N is conjugate to q = N/(N - 1), but since 2 is bounded we
are free to choose p > N.
REMARK 5.4. It is proved in TEMAM [1977] that Eq. (5.8) and Eq. (5.10) imply
the existence of v -o on 12Q
in a weaker sense than Eq. (5.11). Condition (5.11)
is motivated by the second term on An2 in Eq. (2.5) and the observation that
the flow u on the loaded part of the boundary T may be a measure on the
boundary (Example 4.2). The formal integral of the second term in Eq. (2.5)
is well defined, if or is continuous on T, which follows from Eq. (5.11)
for p > N. (By Sobolev space theory functions in Wl-l 1/PP(Q2) are restrictions
to 12)of functions in WP (12), and these are continuous for p > N. See, e.g.,
NEAS [1967].)
The requirements in the definition of I restrict more than meets the eye:
THEOREM 5.4. Let 12 C RN be bounded with Lipschitz continuous boundary. Let

p > N and assume that
or E [L1(12)]N2 , D [L (2)]N 2, Vor e [LP(f2)]N.
Then
2
E [Lq(Q)]N Vq < oo
and
I tr(,o)llLq/R < Cq([ll0 IIL + IIVorlILP) -

The norm on the left-hand side is the quotient norm modulo constants in L q ().
PROOF. See Lemma 2.1 in CHRISTIANSEN [1986]. [
5.2. The flow space U

We formulate the problem for the homogeneous boundary condition u=0 on
S. The case of more general boundary conditions will be dealt with in the next
section. The space U of displacement rates consists of all pairs (u Q,UT) of the
form
(u12,UT) E BD(Q) x [(W1-1 /PP(T))*]N (5.15)
which satisfy the following condition: The linear form
r-- - (V.or).u + I(v.o).uT (5.16)

1 T
defined for symmetric tensors oar [WIP(f2)]N2 must be continuous with respect
to the maximum norm on or. Thus (u a , uT) in Eq. (5.15) belongs to U if and
only if there is a symmetric tensor of measures on , L EG[(C())*]N, such
that for all symmetric o- C [Cl(I)] N 2 the following equality holds:
(ij, 'ij)C(T) xC() = f(Vo)U2 + J(J.).UT. (5.17)

i,j=1 12 T
218 E. Christiansen CHAPrER II
U is an incomplete space equipped with the norm

Il I = Ilu IBD(n) + IluTII(W'-1/p,p(T))*. (5.18)
In the second integral of Eq. (5.16) we have kept integral notation for the
duality W-l1/PP(T) x (W-l1/PP(T))*. In Eq. (5.17) we have tried to emphasize
the generalized nature of the strain rate tensor Ix(u) by using the more correct
duality notation.
In view of Theorem 5.2 and the trace theorem for the Sobolev space WP(f2),
Eq. (5.16) defines a continuous linear form om [WlP(f2)]N2 with its "own"
norm. The tensor p = Ip(u) in Eq. (5.17) is the generalized strain rate tensor
associated with the pair u = (u Q,u T). Clearly, p (u) and the distributional strain
rate tensor e(u I) defined by Eq. (5.1) are equal as distributions in 2. The
difference between Iu(u) and (u) is characterized by the following theorem.
THEOREM 5.5. Assume that n2 C RN is bounded with Lipschitz continuous

boundary. Let u E U with = pI(u) and e = E(u) defined by Eq. (5.17) and
Eq. (5.1) respectively. Then the following equality holds for all symmetric tensors
o, E [C1 ()]N 2 :
N
E(ij - eij, i)C(xC(n) = f(Vo)'(u - (u)) -(V ) y(U)
ij=l T S
(5.19)
where y denotes the trace operator in Theorem 5.1. In particular, Iu(u) and
(u) are equal as measures in 2 if and only if
y(ufl)=uT on T and y(u a )=OonS. (5.20)
PROOF. Equation (5.19) follows immediately from Eqs. (5.5) and (5.17). o
REMARK 5.5. Example 4.1 and Example 4.2 and further examples in CHRIS-
TIANSEN [1982] show that it is necessary to consider the flow on the loaded
part of the boundary separately in order to prove existence of a collapse flow
field. In many applications it will hold that y(u 1) = u T, in which case the col-
lapse flow falls within the frame of TEMAM and STRANG [1980b] and KOHN and
TEMAM [1983]. Intuitively, this will happen, when the surface shear forces are
not too strong (according to Eq. (5.49) the normal components of y(u ) and
uT are always identical for the collapse flow). For the special case of antiplane
shear, this is proved in DEMENGEL [1989, Corollary 2.2].
REMARK 5.6. If we define the internal energy dissipation rate by Eq. (5.17),
then the boundary condition u = 0 on S is included in the definition of U in the
same "relaxed" way as in TEMAM and STRANG [1980b]. The second term on the
right-hand side of Eq. (5.19) is identical to Eq. (4.7), when us = 0.
BD(f2) is a subset of U through the imbedding u - (, Y(u)IT) for all u E

BD(2). In this case the generalized strain rate tensor is given by
N N
(ijij)(c())*xC(Q) = (ij, ij)M()xC) J(-O')y

|- (u)
ii,ij=1
2
for all symmetric or [l()]N . If u satisfies the boundary condition in the
sense that y(u) = 0 on S, then the internal energy dissipation rate is given by
the classical expression.
It is also instructive to identify a vector of bounded measures m on T with an
element of U through the imbedding m - (0, m) c U. The generalized strain
rate tensor t is now given by the expression
N
E (AXij, 'i)c ()xC() = (m, v .)M(T)xC(T7)

i,j=l
for all symmetric or E [C(1)]N2 .
5.3. The duality of limit analysis

On the product space Z x U the internal work rate is defined by the right-hand
side of Eq. (5.17):
a(or, u) =- (V.o').ua + f(vor).u T (o, u) x U. (5.21)

n T
By Theorem 5.2 this is a continuous bilinear form on the product space X x U.

The volume and surface forces are defined as a pair (f ,g), where
f [LP(f2)]N and g [W-/1/PP(T)]N (5.22)
for any p > N. By Theorem 5.2 the work rate for the external forces is well
defined as the following continuous linear form on U:
F(u) = /f.u+ fgu. (5.23)

? T
In the special case, where T = 0Q2, i.e., there is no boundary condition on u,

(f, g) must in addition satisfy the usual compatibility condition (the total load
must be zero):
If + g= (5.24)
REMARK 5.7. Since p > N it follows from Eq. (5.22) that on smooth parts of T
only continuous surface forces are allowed in this formulation (e.g., Remark 5.4).
This restriction appears to be of the same nature as the exclusion of point loads
220 E. Christiansen CHAPrER II
in the corresponding plate model. If, for example, a piecewise constant load g
is the limit of a sequence of continuous loads by use of an approximating unit,
then the corresponding limit multipliers converge, and the limit may be inter-
preted as the limit multiplier for g. This argument also justifies discretization
and numerical solution of such problems. Physically, point loads and step loads
are idealizations. In the elastic model they are admitted, but in the plastic model
they are not.
We shall need the following theorem.
THEOREM 5.6. Let 1] C IRN be bounded with Lipschitz continuous boundary. Let
S and T be open disjoint subsets of 012 with S u T = a01 and S n T Lipschitz
continuous in local coordinateson 012. Letf and g be given as in Eq. (5.22). Then
there exists at least one symmetric or E [W1'P(f2)]N2, such that the equilibrium
equation on classicalform is satisfied:
-Vo-r =f in 12 and v. =g on T. (5.25)
PROOF. Using the regularity of T and Theorem 2.3.9 in NEIAS [1967] it can
be seen that g may be extended to an element of W-l/PP(0Q), such that
Eq. (5.24) is satisfied. In NECAS [1966] it is proved that the system of elasticity is
W 1, (Q2)-coercive. Hence the Neumann problem corresponding to the equations
of equilibrium withf and g has a unique elastic solution u [W 2'P(f2)]N. The
corresponding elastic stress tensor can be used as the solution. [
Finally it is convenient to introduce a short notation for the set of admissible

stress tensors:
K = {o C I oD(x) E B(x) a.e. in f2} (5.26)
or e - is admissible, if it satisfies the yield condition everywhere in 2, except
possibly for a subset of measure zero.
The static principle of limit analysis, Eq. (2.11), now takes the form
A*= sup{A I 3]o C K: a(or, u) = AF(u) Vu U} (5.27)
=sup inf a(a, u). (5.28)
o-EK F(u)=1
The equality between Eq. (5.27) and Eq. (5.28) follows from a standard algebraic
argument: The inner infimum in Eq. (5.28) is, for fixed or, the infimum of a
linear functional a(or, ) over an affine hyperplane. This value is different from
-oo only if the linear functional is constant on the hyperplane. In this case,
a(or, ) is proportional to the linear functional F, which defines the hyperplane,
i.e., a(or, ) = F(u) Vu c U for some real A.
The kinematic principle is
A*= inf D(u)= inf supa(,u), (5.29)

F(u)=l F(u)=1 ocK
SECION 5 The duality problem 221
where
D(u) = sup a(r,u) Vu E U. (5.30)
acK
D(u) is called the total energy dissipation rate associated with u.

The duality theorem of limit analysis states that the values of Eq. (5.28)
and Eq. (5.29) are equal, and that the extremal values are attained as maxima,
respectively minima.
THEOREM 5.7 (Duality theorem of limit analysis). Assume that 12 C RWN is bound-
ed with Lipschitz continuous boundary. Let -Ybe defined by Eqs. (5.8-5.12), U
by Eqs. (5.15-5.16), and let K satisfy Eq. (5.26) and Eqs. (5.13-5.14). The bilinear
form a(o, u) on X x U is given by Eq. (5.21), and F by Eqs. (5.22-5.23). Then
max min a(, u)= min maxa(r, u). (5.31)
a-EK F(u)=l F(u)l1 oA-K
If oar maximizes the left-hand side, and u * minimizes the right-hand side, then
(o*, u *) satisfies the inequality
a(or, u*) < a(or*, u*) = a(o*, u) = A* (5.32)
for all orE K and all u U with F(u) = 1. Here A* is the value of Eq. (5.31).
In particular, (o*, u*) is a saddlepoint for a(or, u).
PROOF. We first prove the duality part of Eq. (5.31) by proving the equality
sup inf a(or, u)= min sup a(or,u). (5.33)
aoCK F(u)=1 F(u)=1 aEK
We apply Theorem 2.1 in CHRISTIANSEN [1980a] to the bilinear form -a(o, u)

on x U and the convex sets K C X and C C U given by
C = {u U I F(u) = 11.
C is closed in U, and due to Eq. (5.14) or = 0 is an interior point of K. Equa-
tion (5.33) will follow by establishing the following two implications (a) and
(b):
(a) If l is a continuous linear functional on U, which satisfies infu,,c (u) >
-oc, then there exists o E X, such that
(u))=a(or,u) Vu CEU. (5.34)
(b) If 0 is a continuous linear functional on X,which satisfies supK 0(I(or) <
so and the condition
a(o, u) = O Vu E U r '(o) = 0,
then there exists u U such that
(o) =a(, u) Va . (5.35)
If I in (a) is bounded from below on the affine hyperplane given by F(u) = 1,
then by the argument following Eq. (5.28), we have
t(u) = tF(u) u U
for some real number t. Let oa c [WlP(/2)]N2 be a symmetric stress field satis-
fying Eq. (5.25) as stated in Theorem 5.6. Then
P(u) = tF(u)=a(to, u) Vu E U.
This proves that the implication (a) holds.
Now let tP satisfy the conditions in (b). By definition of a(or, u) the following
implication holds for 0r:
(V-or = Oin and v-or = on T) (or) = 0. (5.36)
Consider the linear map
2
D: [W1,p(2)]N sym - [LP(f) W-l/p,P(T)]N
defined on the symmetric tensors in [W 1 P( 2 )]N2 by

D(r) = (-V ao, v a).
By Theorem 5.6, D is onto and hence an open map. Condition (5.36) implies
that the null space of D is contained in the null space of P. Hence, there exists
a continuous linear functional -q on [LP(/) x W 1- 1/PP(T)]N, such that
(or) = 7(D(o)) Vo E [W1P(Q)]Nsym (5.37)

71 may be identified with a pair (u , U T) of the form
T) N
(UO,U [Lq(2) x (WI-IPP(T))*] , (5.38)
where lip + 1/q = 1. Equation (5.37) now reads
(O) fJ(V.O). u + f (VO)-UT (5.39)

fQ T
for all symmetric or E [Wl'P(/2)]N2. We claim that Eq. (5.39) holds for all E
: Let or E X. From Theorem 5.6 applied tof = V or E [LP(Q2)]N and vPo- E
[W-ll/PP(T)]N we conclude that there exists a symmetric Oro E [W1p(2)]N2 ,
such that V (or- o0) = 0 in 2 and v.(o - or0) = 0 on T. By Eq. (5.36) P(o-) =
P(oro), and thus Eq. (5.39) holds for or, as well as for o.
It remains to prove that u a belongs to BD(2). Part of the assumption (b) on
P is supIEK 'P(or) < oo. By Eq. (5.14) K contains a ball with respect to the L -
norm on oa, so is continuous on [Wlp(f2)]N2 sYm equipped with the maximum
norm. Hence, there exists bounded measures /-ij = 1/ji E C*(/), such that
N
(S(O) = E ( [ ij
, 'ij)C*(Q)xC(,)
(5.40)
i,j=l
for all symmetric o E [W1"p(f 2 )]N2 . Considering first ar with compact support in
12 we conclude from Eqs. (5.40) and (5.39) that
[ij = ij (u ) as distributions in (2
(5.41)
and hence that u E BD(d2). Using again Eqs. (5.40) and (5.39) we see that
(un,uT) C U, and that
(a) = a(a, u) Var E.
This proves condition (b) and hence also Eq. (5.33).
It remains to verify that the maximum on the left-hand side of Eq. (5.31) is
attained. Let A* be the value of Eq. (5.31), or equivalently of Eq. (5.27), and
let ark E K be a sequence, which satisfies
a(ok,U) = AkF(u) Vu U, Ak A*. (5.42)
The sequence of deviators rID~ is bounded in L
and hence there is a subse-
quence, also denoted Oak, such that a-D is weak* convergent in L- as the dual
of L 1 . Let or E [L(1)]N2 be the limit. In particular, 0rD converges to o in the
distributional sense. If u has components in C (dI) and V u = 0, then
a(oa ,u) = a(ak,u) = AkF(u).

From Eq. (5.42) it follows that ar satisfies the following condition:
If u [C0 (n)]N with V.u = 0,

then (Vo-a + A*f, u)~,gx = 0.
' x! denotes the duality between distributions and testfunctions. Equa-
tion (5.43) is a familiar property in the analysis of the Navier-Stokes equation.
From TEMAM [1977, Proposition 1.1, p. 14] we conclude that there exists a dis-
tribution 0 E @'(2), such that (in the theory of Navier-Stokes equation q is
the pressure)
Vo- + Abf = V q = V.(4 pI) in @'(2). (5.44)
Using NECAS [1966, Theorem 1, p. 108] we see that (p C L 2 (d2) (in fact by The-
orem 5.4 S e L q(d2) for all q < oo). Let
ar* = a - I.
The components of or* are in L 2 (f2). or* satisfies the yield condition, since its
deviator is or. Also V-a = -A*f has components in LP(2). Hence, ar* is a
maximizing stress tensor for Eq. (5.27), if we can prove that v yr*= A*g on T.
By standard Sobolev space theory, , or* exists in W- 1 /2 2 ( 99), and the following
Green's formula holds for any u c [C l ( Yi)]N with u = 0 on S:
fL Oisij(u)=- (J(V.o*).u
+/(v'ar*)-U
n i, n T
=A* f .u + f (v.-r*).u. (5.45)

2 T
224 E. Christiansen CHAPMER II
If in addition u satisfies V-u = 0, the following holds:
J EC 0ijeij() = fE UijEij(U) = lim AkF(u)

Q iJ n i
=A* /f.u +A* fg.u. (5.46)

2 T
From Eqs. (5.45) and (5.46) we conclude that the following implication holds
for o-*:
(u E [WI'2 ()]N, u = O on S, V-u = ) X f(r.o'* - A*g).u = 0.

T
In FOIAS and TEMAM [1978] it is proved that
{y(U) I C [W'l2 (n)]N, VUV = 0} = {V [Wl/ 2 '2 (an)]N if V = 0},
where y denotes the trace operator from [W1,2(Q)]N onto ]Wi/22(012)]N. Thus
the null space of the functional
v- (vor* - A*g).v
T
defined on the space
{v G [W1/ 22 (dQ9)]N Iv = 0 on S}
contains the null space of the functional
V -I J vV.
T
We conclude that
v.a* - A*g= tv = tvPI on T
for some real constant t. But since S in Eq. (5.44) is only determined up to
an additive constant, we may replace so by qp + t. The resulting o* maximizes
Eq. (5.27), concluding the proof. El
From Eq. (5.32) it follows that

a(r*, u*) = D(u*) = maxa(o, u*). (5.47)
crCK
Now suppose for the moment that a(uo, u) is given by the expression (2.4). We
recall that Eq. (2.4) is not well defined on x U, but according to the results
in Section 7 Eq. (2.4) may be interpreted as a(o-, u), if u is divergence free and
thus, in particular, for u = u*. Then Eq. (5.47) implies that at each point x in
SECION 5 The duality problem 225
2, where the collapse strain tensor * = (u*) is nonzero, the collapse stress
tensor o* must be at the yield surface at a point, such that e*(x) is normal to
the yield surface at o*(x). Conversely, if at some point x C 12 the collapse stress
tensor or* is not at the yield surface, then the collapse strain rates are zero,
and there is no local deformation in the collapse state. This is the principle of
complementary slackness in limit analysis. The subset of 12, where or* is at the
yield surface and where local deformations may occur, is usually denoted the
plastified region or with some ambiguity the plastic region. This region need not
be uniquely determined for a given problem.
5.4. The reduced problem
We have earlier noted that D(u) defined by Eq. (5.30) equals +oo, if u is
not divergence free. If, on the other hand, u is divergence free, then a(o', u) =
a(aoD, u) for all stress fields. This suggests that the saddle point problem of limit
analysis could be formulated equivalently within a product of smaller spaces
than Z and U. This observation turns out to be useful both mathematically and
numerically. Two complications must be handled, though: From a mathemati-
cal viewpoint, a(oD, u) is not necessarily well defined, and from a numerical
viewpoint, the condition V-u = 0 increases the number of constraints, rather
than reducing the size of the problem. The latter problem can be handled by
using the finite element technique developed for the Navier-Stokes equation.
The former is dealt with below.
We must specify what it means that u e U is incompressible. From Eq. (5.17)
we see that D(u) = +oo, unless u satisfies the following condition
(tr(/,), q)c()*xC(Q) = 0 V E Cl(n)
or equivalently
N
tr(/) = ii = 0 in C(n)*, (5.48)
i=1
where = /,(u) is the generalized strain rate tensor associated with u =
(u , UT). From Theorem 5.5 we see that Eq. (5.48) is equivalent to the fol-
lowing conditions
V u =0 in 1 (as distribution),
v.(uT - y(u')) = 0 on T, (5.49)
v y(ua)=0 onS.
This is the generalization of incompressible flow to the pair (u Q, u T). In addition
to Ua being divergence free in the classical sense the normal component of the
flow must be continuous along the boundary.
We conclude that the value of the kinematic principle, Eq. (5.29), remains
unchanged, if U is replaced by the subspace
U0 = {u E U I tr(,/) = 0 in C(n)*} (5.50)
For u E U0 and o- E [cl(Q)]N2 we get from Eq. (5.17)

a(o, u) =a(orD, ).
Hence U0 is in duality with the quotient space
o = X/{(pII p c L1 (c)}. (5.51)

The corresponding yield set in - is
Ko = K/{pI I p E L(d)}. (5.52)
The components of X0 will be denoted &. Note that a(&, u) is only well defined
on o0x U0 through an appropriate choice of representative of &.
We may now formulate the duality problem of limit analysis on the product
space L0 x U0. The static formulation of the "reduced" problem is
sup inf a(, ,u) (5.53)
&-EKo EUo
F(u)=1
or equivalently
sup{A 3 E Ko: a(&,u) = AF(u) Vu U0}. (5.54)
The kinematic form is
inf Do(u) (5.55)
uEUo
F(u)=1
with
Do(u) = sup a(&,u) = D(u) u E U. (5.56)
a-EKo
The relation between the static formulations Eq. (5.27) and Eq. (5.54) is
expressed by the following theorem, which is of independent interest. It is es-
sentially proved, although not explicitly stated, in TEMAM and STRANG [1980b].
THEOREM 5.8. Assumptions as in Theorem 5.7. If & E o satisfies
a(&, u) = F(u) Vu Uo
then there is a representative of &, such that or is in equilibrium with (f,g):
a(o, u) = F(u) Vu E U
or equivalently
-V.r =f in f2 and V.o = g on T.
PROOF. This was essentially proved in the existence part of the previous proof.
Equation (5.43) with A* = 1 holds for any representative of O. We now construct
or exactly as we constructed or*. []
THEOREM 5.9 (Reduced duality of limit analysis). We make the same assump-
tions as in Theorem 5.7. Let Uo be given by Eq. (5.50) and o by (5.51). Then
u E U minimizes Eq. (5.29) if and only if it minimizes Eq. (5.55), and E K
maximizes Eq. (5.27) if and only if it belongs to a class & E Ko, which is optimal
for Eq. (5.54). In particular,the duality between Eqs. (5.54) and (5.55) holds and
the value is identical to A*, the value of Eq. (5.31).
PROOF. By definition Do(u) = D(u) for u E Uo. If u *E U minimizes Eq. (5.30),

then A* = D(u*) < o, which by Eq. (5.17) implies that tu(u*) satisfies
Eq. (5.48). In that case u* belongs to U0 and hence minimizes Eq. (5.55). It
follows that Eqs. (5.29) and (5.55) have the same optimal solutions.
If & is feasible for Eq. (5.54) with some value of A, then it follows immediately
from Theorem 5.8 that & has a representative or, which is feasible for Eq. (5.27)
with the same value of A. Since the opposite direction is obvious the proof is
complete. E]
REMARK 5.8. In computations with divergence-free flow we may choose the dis-
crete stresses in the obvious space of representatives
2 = or I tr(or) = }.
With a suitable choice of finite element spaces a(or, u) will be well defined and
computable on the discrete subspace of X x U0.
6. Generalizations and variations
In the previous section we made some simplifying assumptions for the sake of
clarity. We shall now point out the changes for the more general case. Also the
duality theorem was proved for yield conditions formulated in the deviator rD
alone. We shall discuss other types of yield conditions.
6.1. Boundary flow constrainedto a subspace
Until now, we have assumed the boundary to be either completely fixed, or
free and subject to loading. In applications we may face more general boundary
conditions, such as v.u = 0 for a part of the boundary, which can move only
in directions tangent to the surface. denotes the outward normal to 2. Other
examples are v.u = c, where c is an unknown part of the solution. Also it may
be required that the tangential component u - (u -v) vanishes. The special
case, where the conditions are of the form ui = 0 on part of the boundary with
one of the coordinate axes as normal is treated in CHRISTIANSEN [1980a]. This
case occurs, when a problem is reduced by symmetry.
We shall not try to formulate the duality theorem in such generality that the
above mentioned examples of boundary conditions are covered. In the appli-
cations, which have come up, it has been obvious how to formulate and prove
the theorem. It is no additional difficulty to handle such boundary conditions
numerically.
228 E. Christiansen CHAPTER I1
6.2. Presence of a pre-load

There may be other loads present than the load for which the limit multiplier is
to be determined. A typical example is the force of gravity. In the literature this
is called a pre-load or a "dead load". The duality theorem of limit analysis is
easily modified to cover also this case: Let (f 0,go) be the pre-load. The equation
of virtual work is now
a(o, u) = AF(u) +Fo(u) u E U, (6.1)
where Fo(u) is defined by
Fo(u) = Jfo-u + JgoT U (6.2)

a T
We define the following affine-linear form

al(o, u) = a(o-, u) - Fo(u) V(o, u) E X x U. (6.3)
The duality between the static and kinematic principles is now of the same form
as in the previous section with al(or, u) instead of a(o, u). It is again a matter
of simple verification to see that the duality theorem still holds.
6.3. Bounded yield set

In the previous section the duality theorem was proved for yield conditions of
the form Eq. (5.13), which are appropriate for example for isotropic materials.
In particular, the Von Mises condition, Eq. (2.8), which applies to metal-like
materials, is of this form. Also in the plane strain model with such materials the
yield condition is of the form Eq. (5.13).
We now consider the case of a bounded yield set. We assume that the yield
condition is of the following form:
o(x)E B(x) Vx E Q,
B(x) C RN closed, bounded uniformly in x and convex for all x . (6.4)
Also Eq. (5.14) is still assumed to hold. In the two-dimensional plane stress
model the yield condition is of this form for all materials. For example the Von
Mises condition in plane stress is
2
aX + y2y_ a ++3a 2 < 2 (6.5)
This condition clearly implies that each component of belongs to LO(12).
Hence in the definition of the stress space conditions (5.8)-(5.9) may be
replaced by the condition
N 2.
OrE [L (2)] (6.6)
The norm Eq. (5.12) on is replaced by the norm
Io'l1 = III'1L (f) + IIVo(TIILP(L ) + I1V'o'llW1-V/P,P(T) (6.7)
The duality theorem of limit analysis now holds without further due.
THEOREM 6.1. Assume that 12 c RN is bounded with Lipschitz continuous

boundary. Let _ be defined by Eq. (6.6) and Eqs. (5.10)-(5.11), and U by
Eqs. (5.15)-(5.16) (as before). Let
K = C I o(x) B(x) a.e. in f2},
where B(x) satisfies Eqs. (6.4) and (5.14). Then the conclusions of Theorem 5.7
hold.
PROOF. All steps in the proof of Theorem 5.7 can be repeated without mod-
ification, except for the existence of a maximizing stress tensor o*, which is
now much simpler: There is a sequence {fk} in K satisfying Eq. (5.42). This
sequence (and not only its deviators) has a weak* convergent subsequence, and
the limit is easily seen to be a maximizing stress tensor. c
In limit analysis with bounded yield set the collapse flow field is not divergence
free. The total energy dissipation is finite for any flow field. Mathematically this
is a much simpler analysis, and according to experience it is also easier to handle
numerically. This could be the reason why the case of bounded yield set has
been over-emphasized in the literature on limit analysis.
6.4. Yield conditionsfor granularmaterials

We now turn to a type of yield condition, which is more difficult to handle than
both of the above mentioned, but which is too important for applications to be
ignored. It applies to granular materials like soil and concrete (not reinforced).
Characteristic for these materials is that they can resist compression, but not
tension, or at least only to a lesser extend. This means that stress tensors of the
form o = I are admissible only if 4 (x) < o for all x in 2, where o0o > 0.
The classical yield condition for granular materials goes back to COULOMB
[1773]. A more recent reference is KALISZKY [1975]. On any section with normal
n the shear stress rn and the normal stress or, must satisfy the inequality
/OLn + |In[ < c, (6.8)

where , and c are non-negative constants. /Ais the coefficient of internalfric-
tion, and c is a measure of cohesion in the material. If c = 0 the material is
without cohesion. In a frictionless material A = 0, and in this case the Coulomb
condition becomes equivalent to the Tresca yield condition, which bounds the
maximum shear stress at every point of the material. The Tresca yield condition
can be written in the form Eq. (5.13), (KALISZKY [1975, p. 66]), so the case of a
frictionless material has been covered in Section 5.
The Coulomb condition may also be formulated in terms of the principal
stresses o-I, o-I and orIm (the eigenvalues of the stress tensor). Let p be given by
, = tan SD. Then the Coulomb yield condition can be written as three inequali-
ties:
[ol - aIiI + ( + all) sin q0 - 2ccos (p < 0,
1OIrI - 'II I + (II + o-III) sin - 2c cos < 0, (6.9)
A- I I + (Il + mII)
sin SO- 2ccos p < 0.
In Eq. (6.9) an inequality can only be active, if it involves the maximum and
minimum principal stresses.
From a computational viewpoint it is a complication that the principal stresses
must be determined at each node of the discretization, before the yield condition
can be applied. In the case of plane strain the extremal principal stresses are
known, and the Coulomb condition may be expressed in rectangular coordinates
of our choice:
/(x - rjyy) 2 + 4o2y + (x x + oyy) sin (q - 2c cos p< 0. (6.10)
For other yield conditions of this type we refer to KALISZKY [1975, Section 8.3].
Common for them is the property that the cone
{ I (X) < aOrx C 2}
is contained in the set of admissible stress tensors. Hence a necessary condition
for the total energy dissipation rate to satisfy D(u) < o is
V-u 0 in 12. (6.11)
In particular, the collapse flow u* must satisfy Eq. (6.11). In analogy with
Eq. (5.48) the precise form of this semi-incompressibility condition for u E U
is
N
tr(/J) = E/ii >0 in C(Q)*, (6.12)
i=1
where p is the generalized strain rate tensor for u = (u , u T).

Theorem 5.7 does not hold as stated for such yield conditions. If the condition
(5.14) holds, the proof of Eq. (5.33) is valid without change. Hence the duality
between the static and the kinematic principles holds, and a minimizing flow field
exists for the kinematic principle. However, the existence proof of a collapse
stress field does not hold. To our knowledge this analysis has not been made.
If the condition (5.14) does not hold, arbitrarily small forces may cause yield.
This seems to be an extreme case without much practical interest. For the
Coulomb yield condition it occurs, if c = 0 in Eq. (6.8), and in this case the
limit multiplier is either 0 or +oo for any load distribution.
7. Green's formula in the collapse state
In this section we shall try to make sense of the expression in Eq. (2.4) for
the internal energy dissipation rate, or equivalently, analyse the validity of
Green's formula for the collapse fields (*,u*). On X x U, a(o, u) is defined

by Eq. (5.21), which is identical to Eq. (2.5). By Eq. (5.17) a generalized Green's
formula holds for u U and or with components in C1 (12). However, the left-
hand side of Eq. (5.17), which is the generalized form of Eq. (2.4), is not defined
for (, u) E x U; 1tij is a measure, but oij need not be continuous.
Recall that (u) denotes the classical (distributional) strain rate tensor de-
fined by Eq. (5.1). We first analyse e(u) under the extra assumption that the
divergence V-u belongs to Lq(12) for some q in the interval 1 < q < oo. Let u
satisfy
u BD(n), Vu E Lq(2) (7.1)
and let or be symmetric and satisfy
oa E [L2(12)]N2 , TD E [L-(n)]N2, V.o E [LN(n)]N. (7.2)
Define the following distribution, which we shall denote by (D, ED(u)):
|(OD D()) _N/ tr(v)(V U)S - o

(V A)-u(P
-/ Aor, u V ) V EC(n) (7.3)
We have used the notation
(a, u V) . = ijUi ax'

ij
In Eq. (7.3) each term on the right-hand side is well defined by Theorem 5.4
and Theorem 5.2. It must be emphasized that in general "(aD, ED(u))" is not
an inner product as suggested by the notation. It is a convenient, but inaccurate
notation, which is justified only when or and u are sufficiently smooth. In that
case Eq. (7.3) is a consequence of Green's formula. The following theorem is
proved in KoHN and TEMAM [1983, Theorem 3.2]. It is the closest to a regularity
theorem for the collapse state we get.
THEOREM 7.1. Assume that 2 C RN is bounded with Lipschitz continuous

boundary. Let 1 < q < oo and let u and a satisfy Eqs. (7.1) and (7.2) respec-
tively. Then the following statements hold:
(a) The distribution (rD, eD(u)) defined by Eq. (7.3) is a bounded measure
on 12 and it satisfies
|(0",
CD(U))(PO~
<,p1EDIO| D(1)0"I SOC (S2) (7.4)
n n
(b) If the boundary of 12 can be defined locally by C 2 -functions, then there is
a uniquely defined distribution "(o. v, u)" on 02, which satisfies
/an (, u), J(,'D

l2
CD(U)) + |J(V.u) tr(o,)
n
+ J(vr,).u ~s+ (o, u v) VD E C((2). (7.5)

2 n2
Furthermore (. v, u) belongs to the dual of C' (32), and Eq. (7.5) holds for all
cP
E C'().
(c) If {Uk} is a sequence in BD(Q) with V-uk C Lq(d2), which converges to
u in the following sense
Uk -- U in [LNIN-I(2)]N,
V.uk -V.u in Lq(2), (7.6)
fn IheD(Uk)li -- , f ED(u){n
then
I (C 0 D(Uk))(P "f(D, D(U))P Vc E C(). (7.7)

12 12
(d) If { oak } is a sequence satisfying Eq. (7.2), and which converges to ar in

the following sense
ak f in [L2(2)]N 2 ,
Vaok - V in [LN(f2)]N, (7.8)
IlakIllL < C, c independent of k
then
(ea ED(U))qO (SD, ED(U))(P V, E c 0 (h). (7.9)

Q 12
PROOF. See KOHN and TEMAM [1983, Theorem 3.2]. o
In Eq. (7.5) (. v,u) is not in general an inner product as suggested by

the notation. However, if a E Z, and u E U with y(u) = 0 on S and V-u
Lq(2), q > 1, then by Eq. (5.11) and Theorem 5.1 we may conclude that
(erD, eD(U)) satisfies the identity
J(orD,e(U))S = ( or, U Vqo) - Jtr(o)(Vu)p

f
- (V Or).U qp + J( v).y(u)p Vp E c 1(Q). (7.10)

1 T
If (o-, u)= (*, u*) is the saddle point in Theorem 5.7, then V-u = 0 and
E(u) = D(u). If furthermore y(u) = 0 on S (the boundary condition without
relaxation), then we may apply Eq. (7.10) with q0 = 1 and get
/(o, e(u)) = f(o, (u))

n Q
=- J(V (r)u + /(o v)(U). (7.11)

1 T
This is Green's formula in the case /(u) = e(u), or equivalently (u12 ) = 0

on S and y(u 2 ) = UT on T, where gu(u) is the generalized strain rate ten-
sor. Anticipating the more general case, we now state the above result with
different assumptions, moving the regularity conditions from the domain Q to
the pair (, u). The technical details are not essential, and our assumptions
may be stronger than necessary. We shall assume that one of the following two
conditions holds:
or E W 1'P(T3) for some > 0, (7.12)
where Ta = {x E f2 I dist(x, T) < S} is a neighbourhood of T in 2.
y(u 12 ) W (l/r)r'(T) for some r > 1. (7.13)
THEOREM 7.2. We make the same assumption as in Theorem 5.7. Let (or, u) =
(a*, u*) E Z x U be the saddle point in Eq. (5.32). Assume that y(u a) 0 on =
S and that Eq. (7.12) or Eq. (7.13) holds. Then (or, e(u)) = (D, eD(u)) is a
bounded measure on 2, which satisfies Eq. (7.4).
PROOF. See CHRISTIANSEN [1986, Corollary 3.4]. E
In the general case we need a similar result for the generalized strain rate
tensor /L(u) for u = (u2, UT). We replace condition (7.13) by the analogous
condition on u T:
uT C w-(l/r),r(T) for some r > 1. (7.14)
THEOREM 7.3. Make the same assumptions as in Theorem 5.7 and let (or, u) =
(a*, u *) be the saddle point in Eq. (5.32). We define (o-, /) = (D, /D) weakly
by the equation
/(rD ~D)(p = - (V. O)U fo -/ UJ t X V + ( V) u Tq
(7.15)
for all p E C1 ( 12). If condition (7.12) or condition (7.14) holds, then (, ) is a
bounded measure on 2, and it satisfies
(a,,a)f|
< 11-D~ |JL I(P/
|LI V C0 (12). (7.16)
PROOF. See CHRISTIANSEN [1986, Theorem 3.1]. o
Applying Eq. (7.15) for or smooth and p = 1 shows that (or, is) is an extension
of the left-hand side of Eq. (5.17) and that Green's formula holds in the collapse
state.
In Theorem 7.3 (*, u*) can be replaced by any other pair (or, u) as long as
u is divergence free in the sense that Eq. (5.48) holds.
If the conclusion of Theorem 7.3 holds, then it follows from Eq. (5.47) that
fJ(o,/ (u*))= sup

orEC(1)nK
f ij
oEijbij(u*).
This is the generalized form of the principle of complementary slackness. If Ao*

and u* can be interpreted pointwise, then or* (x) must at each point x maximize
the expression
E ijij
i ,j
subject to the yield condition at the point x.

CHAPTER III
Discretization by Finite Elements
8. Formulation of the discrete problem
The duality problem of limit analysis is well suited for discretization by the finite
element method. In its simplest form Yh C Z and Uh C U are finite-dimensional
subspaces, and K C Z is replaced by Kh C Yh. In general Kh need not be con-
tained in K. Then the discrete version of the problem is:
A =max{A | 3' h E Kh: a(0rh, Uh) =AF(Uh) VU h E Uh}
= max min a(Orh,Uh) (8.1)
ohEK h F(Uh)=l
= min max a(o'h, Uh)

F(uh)=1 rhEKh
= min Dh (Uh)
F(uh)=l
The above duality is proved in Theorem 8.1.

It is natural to choose mixed finite elements to discretize stresses and flow
simultaneously. If only o- is discretized in Eq. (2.11) or Eq. (5.27) the equilibrium
equation must be satisfied exactly for some O'h C h, which can be the case only
for special loads. If only u is discretized in Eq. (2.13) or Eq. (5.29), we will get
D(u) = +oo, unless u is divergence free.
Let 3-h be a triangulation of 2 (CIARLET [1978, 1991]) and let Xh and Uh be
corresponding finite element spaces representing X and U. The discrete fields
for or and u will be denoted h and uh respectively. Frequently we will use
different finite element functions to define the single coordinates of o-h and h,
but for practical reasons they should be based on the same triangulation. We
ignore the complications which arise when 12 is a domain with curved boundary.
These are not specific to limit analysis, so we may refer to CIARLET [1978, Section
4.4] or CIARLET [1991, Chapter VI].
It is convenient to think in terms of conforming elements, i.e., Zh c Z and
Uh C U. However, non-conforming elements (see, e.g., STRANG and Fix [1973,
Chapter 4] or CIARLET [1991, Chapter 5]) may be used, if only the energy func-
tions a(O-h, Uh) and F(Uh) are well defined. Computational experience indicates
that piecewise constant elements for O'h paired with piecewise bilinear elements
for Uh is a useful combination.
235
236 E. Christiansen CHAPTER III
Due to the essential boundary condition on u, and in order to compute

a('h, Uh) by Eq. (2.4), it is convenient to choose globally continuous element
functions for Uh. Then the partial derivatives of h are integrable functions
(CIARLET [1991, Theorem 5.1]). Since the collapse fields o* and u* are typically
discontinuous there may not be obtained significantly better accuracy by using
more smooth element functions.
For simplicity we will assume that the energy functions a(ah, Uh) and F(Uh)
can be computed exactly on 1h and Uh by the expressions (2.4) and (2.2). The
effect of using numerical integration is not specific to limit analysis, so we may
refer to CIARLET [1991, Chapter IV].
The yield condition o-(x) E B(x) will be imposed only through the nodal
values of Orh(x): If {x>} are the nodal points for the value of 0 'h(x), then the
discrete set of admissible stresses Kh is defined as
Kh = {'h E h I h(Xv) E B(x) Vx,}. (8.2)
Assume for the moment that B(x) = B is independent of x. Then Eq. (8.2) im-
plies that Oh(x) E B Vx C 2 for piecewise linear and multilinear elements (lin-
ear n-simplices and n-rectangles of degree 2 in the notation of CIARLET [1991,
Chapter II]). However, for other element types, such as quadratic or bi-quadratic
elements, the yield condition may be violated between nodes, even if it is sat-
isfied at all nodes. We shall return to this problem in Section 9. Examples of
element combinations in plane strain are:
* Linear elements for both stresses and velocity. The nodes are the vertices of
the triangles.
* Bilinear elements (rectangles of degree 2) for both stresses and velocity. The
nodes are the vertices of the rectangles.
* Bilinear elements (rectangles of degree 2) for velocity and piecewise constant
elements for stresses. The nodal values for the stresses are the values in each
rectangle.
* Linear elements for velocity and piecewise constant elements for stresses. The
nodal values for the stresses are the values in each triangle. This combination
has too many degrees of freedom for the stresses compared to the velocity.
In our experience such lack of balance causes ill-conditioning. This element
combination cannot be recommended.
* Cl-bi-cubic rectangles (Bogner-Fox-Schmit elements, see CIARLET [1991, p.
92]) for the stream function and bi-quadratic rectangles for the stresses. This
combination is used to implement divergence-free flow. The nodal values for
the flow are the values of the stream function, its gradient (i.e., the flow) and
the mixed second-order derivative at the vertices. The nodal values for the
stresses are values of the stress components at the vertices, at the midpoints
of the edges and at the center of the rectangles.
Since the notation for the discrete problem is somewhat complicated, we start
by giving a specific example.
SECTION 8 Discretizationby finite elements 237
EXAMPLE 8.1 (Linear elements for both h and uh in plane strain). Let 3-h be a
triangulation of 2 C Ri2 into triangles with node set NX. Associated with each
node v E X is a piecewise linear scalar function , with value equal to 1 at
node v and equal to 0 at all other nodes. Each of the three components in the
stress tensor, 11ll,c022 and 12, is given by a scalar function, which is a linear
combination of node functions of type 4,. Hence associated with each node
v E X there are the following three basis functions for the discrete stress space
;h:
The discrete stress tensors h E Vh may now be written
ah = (evlx v + 2r2 + 612 12

vex
Similarly there are two basis functions for the discrete velocity space Uh asso-
ciated with each node Ee X:
= ], p2 = ]
Each uh E Uh may be written as a linear combination
IIT ~ +
Uh = E ( 7 A11)'
We shall now describe the construction in the example in more general terms.
Let h and Uh be the finite element representations of X and U respectively. Let
, and Xu denote the corresponding set of nodes for the single components of
a h E Xh and uh E Uh respectively. If several degrees of freedom are associated
with the same node (as is the case for example for cubic Hermite elements
defined in CIARLET [1991, p. 84]), then the node is counted once for each degree
of freedom. Let , be the scalar basis function associated with the node v E
A;,. Based on , we define one basis function for h corresponding to each
independent component of the stress tensor. In full three-dimensional problems
there are six basis functions based on 0,:
i ,1o 0o O 1 = 0 O etc. (8.3)
In plane problems there are only three basis stress functions for each node,
corresponding to the three independent components of a symmetric 2 x 2 tensor.
The complete basis for h is the set of all basis tensor element functions:
Ii i V E N }.
238 E. Christiansen CHAPrER III
Every discrete stress tensor orh CE h may be written as a linear combination of

these basis functions:
'h = SE E iivji (8.4)
vctV i<j
The coordinate e? is the nodal value of the component -h,ij corresponding to

the node v.
Similarly, for every node /L E fu we define a basis function for Uh corre-
sponding to each velocity component at that node. If C denotes the scalar
basis function associated with the node g we get the following functions for Uh
in three-dimensional problems:
," O, ,3' =i O (8.5)
Disregarding for the moment boundary conditions on u the complete basis for
Uh is
{k Ilk=1,2,3, CtEV}.
Any discrete velocity uh E Uh may now be written
Uh E E I p ,k (8.6)
AeJVu k
where ,k is the nodal value of the velocity component Uh,k corresponding to

the node /z.
After insertion of the expressions (8.4) and (8.6) the energy functions may
be expressed as
F(Uh) = E E 715F(4,) (8.7)
AC(-\N k
and
a(h,,u,h)= E E E E '71 a(/v', ,) (8.8)

vaffo ij MENc. k
Thus F(Uh) is a linear form and a(o'h, Uh) a bilinear form in the "long vectors"
(~J) and (o), which means that F corresponds to a vector and a to a matrix. It is
a computational necessity (at least for the classical algorithms in numerical linear
algebra) that the coordinates of these vectors are given by one-dimensional
numberings, i.e., by one-to-one mappings
n = n(v, i, j) E {1, 2,..., N} (8.9)
and
m = m(A", k) 1, 2,..., Ml. (8.10)
With such numberings we may write the nodal values as

xn = f/J, where n = n(v, i, j) (8.11)
and
Ym = 7, where m = m(/p, k). (8.12)
The vectors x = (x.)n=l E RN and y = (y)mI E RM are in unique correspon-
dence with the discrete fields O'h cE h and uh E Uh respectively. By insertion in
Eqs. (8.7) and (8.8) we get
M
F(uh) = Eymbm = bTy, (8.13)
m=l
where
bm = F(IkI), m = m(p, k) (8.14)
and
M N
a(oh, Uh) = E ymxnamn = yTAx = xT(ATy). (8.15)
m=l n=l
A is the M x N matrix with entries

k )
amn = a(', f (8.16)
with m = m(A, k) and n = n(v, i, j) given by Eqs. (8.9) and (8.10). These num-
berings will typically be chosen in order to obtain a certain matrix structure.
As usual, when working with the finite element method, the components of
Eqs. (8.14) and (8.16) are computed through integration element by element
and then assembled.
The matrix product notation bTy is used for the ordinary inner product in
RN. Also the components of "long" vectors and matrices are denoted by indices
m, n instead of i, j to avoid confusion with the use of i and j as indices to space
coordinates, such as in -ij and 'j.
With the notation Eqs. (8.13) and (8.15) the discrete version of the duality
problem can be formulated as a finite-dimensional mathematical programming
problem in the variables x E RN and y C RM. Inserting Eqs. (8.13) and (8.15) in
Eq. (8.1) gives
A = max{A I 3x E Kd: Ax = Ab} (8.17)
= max min yTAx
XEKd bTy=l
= min max xT(A T y)

bTy=l XGKd
= min Dd(y). (8.18)

bTy=l
We have introduced the notation
Dd(y) = max xT(ATy) (8.19)
xEKd
and
Kd = {x E RN I h E Kh}, (8.20)
where x and h are related through Eqs. (8.4) and (8.11), and Kh is defined
by Eq. (8.2). (The subscript "d" is mnemonic for "discrete".) We shall always
refer to the discrete static formulation (8.17) as the primalproblem and to the
discrete kinematic formulation (8.18) as the dual problem.
The nonlinear constraints x E Kd in the primal problem is of a particularly
simple structure: Recall that the vector (x,) is just a numbering of the vari-
ables 'J through Eq. (8.11). There is one constraint for each node v, and this
constraint only involves the "short" vector (i) with six components in three-
dimensional problems. This sparse and perfectly banded structure of the nonlin-
ear constraints is well suited for nonlinear optimization algorithms, such as the
projected Lagrangian algorithm described in MURTAGH and SAUNDERS [1982].
EXAMPLE 8.2 (Linear elements in plane strain). Consider a problem in plane
strain with the variables
=
[i o12l
2 and u=
7
12 fT22 U2
The Von Mises yield condition in plane strain is

(fi 1 1 - 22)2 + 4f122 < 4 2 (8.21)
Let v be any numbering of the nodes, i.e., the vertices in the triangulation of
X, and let the numbering in Eq. (8.9) be defined by
n(v, 1,1)= 3(v - 1) + 1,
n(v, 2, 2) = 3(v - 1) + 2,
n(v, 1, 2) = 3(v - 1) + 3.
This corresponds to ordering the nodal values of h by taking the triplets
(fil, 22,
22, 2 ) sequentially vertex by vertex. With this ordering the yield con-
dition on h is implemented as follows:
2 2
(X3(vl)+l - 3 ((v-1l)+
2) + 4(X 3 (-1l)+ 3 ) 32 V.
Kd in Eq. (8.20) is the set of x E RN, which satisfy this condition. In particular,
N equals three times the number of vertices in the triangulation, and the yield
condition is implemented as N/3 nonlinear inequalities, each involving three
variables.
We now prove the duality between the discrete formulations Eqs. (8.17) and
(8.18).
THEOREM 8.1 (Discrete duality theorem). Let A be an M x N matrix and b C

1EM. Let Ko C RN be convex and compact, and assume that x = 0 is an interior
point of Ko. Let V C RN be a subspace and define Kd = Ko + V. Then the prob-

lems in Eqs. (8.17) and (8.18) have the same value, and there exists a saddle point
(x*, y*) for the bilinearform yTAx:
(y*)TAx < (y*)TAx* = Ah = yTAx* (8.22)
for all x E Kd and y E RM with bTy = 1.

Furthermore, if b V RA = {Ax I x CE N}, then the common value of Eqs. (8.17)
and (8.18) is zero.
PROOF. In the case b E RA there exists x E RaN such that b = Ax. Then the
proof goes as the proof of Theorem 5.7, although much simpler in this finite-
dimensional case.
If b V RA, the value of Eq. (8.17) is clearly equal to 0. Now NAT = (RA)-
where NAT and (RA)- denote the null-space of the transpose of A, and the or-
thogonal to RA, respectively. Hence there is y NAT, such that bTy = 1, proving
that also the value of Eq. (8.18) equals 0 in this case. []
The duality in Eqs. (8.17)-(8.18) is just a reformulation of Eq. (8.1). If we

insert Eqs. (8.13) and (8.15), Eq. (8.22) may be expressed in the discrete stresses
and flow corresponding to x and y:
a(O'h, Uh) (h, Uh) = A = a(a', Uh) (8.23)
for all a-h E Kh and Uh E Uh with F(uh) = 1.

The exceptional case b RA in Theorem 8.1 has an important implication:
We must choose the finite element spaces 25h and Uh in such a way that b E RA.
Otherwise the discrete limit multiplier A will equal zero, independent of the
element size.
With one exception, to be pointed out below, the condition b C RA should be
ensured by satisfying the following consistency condition:
The matrix A defined by Eq. (8.16) must have full row rank M.
This consistency condition is independent of the orderings of Eqs. (8.9) and

(8.10). It may equivalently be formulated as a condition on the discrete spaces
,h and Uh:
a(ah, Uh) = 0 Vah E h must imply that Uh = 0 (8.24)

Loosely speaking there must be sufficiently many degrees of freedom for the
stresses relative to the velocity.
In the very special case with no boundary conditions at all on u, rigid body
motions are admissible, and Eq. (8.24) cannot be satisfied. In this case the con-

dition b E RA = (NAT) is simply the discrete analogue of the condition (5.24).
Unfortunately the consistency condition (8.24) is not sufficient to guarantee
good solutions for the collapse fields for stress and flow. If there are too many
degrees of freedom for the stresses O'h relative to the velocity uh, the collapse
stresses will be poorly determined, and vice versa. In our experience the balance
between primal variables x and dual variables y is quite delicate. In order to
illustrate this point we return to the example in antiplane shear.
EXAMPLE 8.3 (Antiplane shear). This is a modification of Example 4.1. We use
the same geometry and notation. Now both sides of the layer are kept fixed
by imposing the boundary condition u(-l) = u(1) = 0. Hence there are only
volume forces, and they are given by
f( -1 for - <x<a,
fx)= 1 l for a < x < 1
for some constant a in the interval 0 < a < 1. The static principle is
A*= max{A I r'(x) = -Af(x), JIo(x) 11
with the solution
* 2
+a
and
o*(x) = 1 - A*x - a I Vx E [-1,1].
The collapse flow is
*(x_:C-1/(1 + a) for - 1 < x < a,
U(* = o for a < x < 1.
Both o-* and u* are uniquely determined. The region -1 < x < a is moving
"down" as a solid block. The rest does not move. The entire deformation takes
place at x = -1 and x = a.
Now discretize the problem as described earlier in this section. Use linear
elements for both h and Uh. To fix ideas let a = and h = . In this case the
exact collapse stress o-* actually belongs to the discrete space h, SO we are
tempted to expect that it coincides with the discrete collapse stresses ah*. Not
so! The discrete equilibrium equation in Eq. (8.1) is
a(oh, Uh) = AF(Uh) Uh E Uh
or equivalently (recall that u(-1) = u(1) = 0)

1 1
-f ChW=A f Vk,
-1 1
where is any of the "hat-functions" that constitute the basis of the linear
element space Uh. If ah alternates from node to node between two values, then
o- will simply change sign from interval to interval. For such a oah we will have
1
~ 0 v4,
for all basis functions qk and thus for all Uh E Uh. Due to the boundary condi-
tion u(-1) = u(1) = 0 there are no basis functions for Uh associated with the
nodes x = -1 and x = 1. This means that any such alternating Oh is in discrete
equilibrium with the zero load, in addition to the constant stresses.
Now let (A*, or*, u*) be the exact solution found previously. o-* belongs to
2
h for a= 1 and h= 1 and is in equilibrium with the load A*f (both ex-
act and discrete). Let h E 2h be given by h(-1)=Oh(O) = oh(l)= 1 and
Och(-) = oh(1) = -1, so that it is in discrete equilibrium with zero. Then for
small positive values of t the stress function a* + th will still be in discrete
equilibrium with A*f, and it will satisfy Io* + th I 1 - t < 1. It follows that
A*/( - t) is an admissible load multiplier for the discrete problem and that
h > A*. The discrete collapse stress function o-h* will satisfy Ioh*(x)l = 1 in the
intervals adjacent to the points x = -1 and x = , and in these intervals defor-
mation may take place. The discrete collapse flow is no longer uniquely deter-
mined, and the plastified region may expand from two points for the continuous
problem to intervals of positive length in the discrete case.
The point of the above example is that the discrete problem may have more
unwanted features than the corresponding continuous problem. It may display
a higher degree of non-uniqueness, both for stresses and flow, and it may spread
out an otherwise narrow and well defined plastified region. In our experience
these features play a significant role in two-dimensional problems.
The negative features in the example do not occur, if the piecewise linear el-
ements for oih are replaced by piecewise constant elements. The set of stresses,
which are in discrete equilibrium with the zero load, is a one-dimensional sub-
space consisting only of the constant functions as in the continuous case. We
shall return to this element combination in Example 10.1 in a different context.
It is an important property of the exact collapse flow field u * that it is diver-
gence free: V.u* = 0. We shall now examine in which sense this condition is
satisfied by the discrete collapse flow. u satisfies the condition
Dh(Uh) = max a(aTh, Uh) = A < oo, (8.25)
oahEKh
where Kh is given by Eq. (8.2). Inserting Oh = KOh in Eq. (8.25) yields
J PhVuh = 0 V/h, (8.26)
where Oh denotes any finite element function of the type used to represent the
single components of the stress tensors ah.
The exact interpretation of Eq. (8.26) depends on the finite element spaces
used. With piecewise linear elements for uh and piecewise constant elements
for O'h, and hence for Oh, Eq. (8.26) means that V-uh = 0 exactly, since it is
constant on each triangle. (As mentioned earlier this element combination has
too many degrees of freedom for stresses relative to velocity.) If other elements
for u h are combined with piecewise constant stresses, we can conclude only that
V.u has zero average over each element. When piecewise linear elements are
244 E. Christiansen CHAPER III
used for both ah and Uh, Eq. (8.26) implies that V u has zero average over the
union of the elements adjacent to each vertex. In general we may conclude that
V.u * is zero in some average sense only. In some of the applications we shall see
that this zero mean may result from undesirable and non-physical fluctuations
in the numerical solution.
In the next section we shall discuss the possibility of imposing discrete incom-
pressibility as constraints.
9. Elements for divergence-free flow
In this section we shall discretize the reduced duality problem in Eqs. (5.54)-
(5.55). It only applies to yield conditions of the form in Eq. (5.13). The expected
gain is a smaller discrete problem considering only divergence-free flow and
stresses with trace zero, o-rii = 0, cf. Remark 5.8. We also expect an easier
discrete problem, since the set of admissible stress tensors is bounded in the
reduced problem.
Finite element spaces for divergence-free flow are known from the numerical
solution of the Navier-Stokes equation and are described in TEMAM [1977, Sec-
tion 4.4-4.5]. They can be obtained in two ways: For certain elements discrete
(but not exact) incompressibility can be obtained by removing degrees of free-
dom. Unfortunately, the resulting space does not have a standard finite element
basis. Application of these elements corresponds to imposing Eq. (8.26) as a
set of linear constraints on Uh. TEMAi [1977, Section 4] gives specific examples
of velocity spaces in the plane: Piecewise quadratic elements (over triangles)
satisfying Eq. (8.26) for h piecewise constant and a slightly larger space sat-
isfying Eq. (8.26) for obh piecewise linear. However, the obvious combinations
with piecewise constant and piecewise linear elements respectively, do not sat-
isfy condition (8.24), since -h has smaller dimension than Uh (see Table 15.1).
To our knowledge computations in limit analysis with this approach have not
been reported, but it certainly deserves further investigation.
The other approach is based on the observation that the divergence of the
curl of any vector field is zero. This approach is briefly described in CmRIs-
TIANSEN [1991]. Instead of choosing a finite element space for the flow u itself,
we discretize the vector field gI to Y'h, i.e., the components of tPh belong to
some finite element space. The discrete flow uh is then given by the curl of 'Ph:
Uh = V x rh . (9.1)
One advantage of this method is that V-Uh = 0 is satisfied exactly. Care must
be taken that h satisfies the boundary conditions. An example is given below.
Now we pay the price for the expected gain: In order to get Uh continuous,
the components of fh must be of class C1, i.e., have continuous derivatives.
In the plane this implies the use of elements like Argyris triangle or Bell's
triangle (CIARLET [1991, Section 9]). If the geometry permits a triangulation
into rectangles (in two or three dimensions) the tensor product of the standard
cubic Cl-elements (Bogner-Fox-Schmit rectangle, CIARLET [1991, Section 9]) is

a very convenient choice.
With Cl-elements for 'Ph the dimension of the discrete flow space Uh is so
large that piecewise constant or linear elements cannot be used for the stresses.
Condition (8.24) cannot possibly be satisfied. More degrees of freedom are
required for the stresses. This can be obtained by using quadratic or bi-quadratic
elements (CIARLET [1991, Section 6-7]). With these elements the yield condition
may be violated between nodes, even though it is imposed at the nodes. As
earlier mentioned, this may be considered part of the discretization error, and we
have no reservations against such elements. Having specified the finite element
spaces the rest is a matter of numbering variables and computing integrals as
described in Section 8. A specific example will be given below.
In plane problems this approach is particularly attractive: 'Ph has only one
non-zero component, which we shall denote Pth, and Eq. (9.1) reduces to
Uh( 0 th , _ th ) (9.2)
The stress tensor of trace zero satisfies 0'2 2 = -o'l and may be identified with
the vector (o-, o-2):
I(-> I (9.3)
With this notation the Von Mises yield condition in plane strain, Eq. (8.21),
becomes
12 + 22 < a
2
(9.4)
The left-hand side is the Euclidean norm of (01, -2 ). This is the simplest possible
yield condition.
EXAMPLE 9.1 (Example 1.1 revisited). Consider the plane strain problem in Ex-
ample 1.1. h is given by bi-cubic Cl-elements over rectangles (actually
squares). There are four basis functions associated with each vertex /u, cor-
responding to the following nodal values of h at the node tz (see CIARLET
[1991, p. 92]):
', at at 02 tF
'aX1' ax2' O
XiX' 2
We shall denote the corresponding four basis functions o0, q 10, qiO and C1.
If V, denotes an arbitrary vertex, then the bi-cubic Cl-basis function </' 3 is
completely characterized by the following identities:
((v) = ., '9 4( = '0, ' (Vo 0,
C1
t(VI) = 0 ad 10(V , a=
O(Vvv = O,
0 O)
9x, a 1 aX2, .,
246 E. Christiansen CHATER III
a 0a 0
(VV) =O, a (v)=0, ax-(v)=08, a
(V ) = O (V) = ' (V = ', aal 0
'tla \ X-2 v/l aX2
Recall that = (xI, x2 ) is just the product of two cubic C1 -element functions
in xl and x 2, respectively.
With a notation analogous to Eq. (8.6) Th may be written
1 1
=
/>h E E (9.5)
,tEX a=0 ,8=0
The boundary condition on u comes from the symmetry:

u2 =0 forx 2 =0,
ul=O forx =0 and 0O<x 2 <a,
where x2 = a marks the tip of the cut. By Eq. (9.2) this boundary condition may
be implemented by the requirement
qth(Xl,X 2) =O for (x 1,x 2) E {x2 = O}U {x =0, O<x2 <a}
or equivalently
700 =0 for all nodes C/-
E{x2 = 0}U {xl =0, O x 2 a}.
The values of Uh are easily extracted from Eq. (9.5) using Eq. (9.2):
Uh(Vt-) = (01 _10)
It turns out that the space of piecewise bi-quadratic elements for the stress
components rl = a 1l = -22 and o2 = 012 has a dimension which is compatible
with the space for Uh. The nodes are the vertices, the midpoints of the sides of
the rectangles and the midpoint of the rectangles (see CIARLET [1991, p. 77]).
Associated with each node is a basis function characterized completely by being
equal to 1 at "its own" node and 0 at all other nodes. Let S, denote this set of
nodes, and let 4, be the associated bi-quadratic basis functions. Corresponding
to each node we have the following two basis functions for the reduced space
of stresses in Eq. (9.3):
1 I[V O 1 .2 V
The discrete stress tensor may now be written
ah = E (=Vb + P2+V) (9.6)
The yield condition (9.4) is imposed as follows:
(el)2 ()2 0 VV E V,. (9.7)

We now introduce numberings of the variables and "long" vectors x and

y as in Eqs. (8.9)-(8.12) and compute the components of A and b defined as
in Eqs. (8.13)-(8.16). Note that the integrals in Eq. (8.16) involve second-order
derivatives of 'h, but since Th is made of C 1 -elements, the second-order deriva-
tives are perfectly integrable. Also note that all integrals can be computed as
products of one-dimensional integrals. Each integral is simple, but there is quite
a manifold of them. They may with advantage be computed using a program
for symbolic calculation.
The solution of the discrete problem and the results for this example will be
described in Section 14.
10. Bounds and convergence
In this section we examine the relationship between the solutions of the ex-
act continuous problem in Eq. (5.31) and the approximate discrete problem in
Eq. (8.1). Let it be said immediately that the situation is far from as satisfactory
as for elliptic problems. There is no regularity theorem for the solution; in fact,
the collapse fields may not even be continuous. There is no energy norm, which
is minimized by the solution, and there is no equivalent of C6a's lemma (CIAR-
LET [1991, Theorem 13.1]). Still it must be possible to obtain better convergence
results than those reported here, but that remains to be done.
According to our experience the convergence of A*will be linear:
IA - A* < ch (10.1)
for some constant c. We shall briefly say that A* - A*is O(h). This convergence is
not improved by choosing elements of higher order, since the collapse fields are
not smooth. Fortunately the linear convergence also seems to be independent
of our ability to prove it.
In numerical analysis error estimates like Eq. (10.1) are frequently combined
with the possibility of using Richardson extrapolation to improve the results
or to estimate the actual size of the difference A - A*. This is based on the
following condition, which is stronger than Eq. (10.1):
A* - A*
lim A = c. (10.2)
h-o h
It is well known that conditions like Eq. (10.2) can be checked from computed
values of A*, i.e., A* need not be known. In limit analysis Eq. (10.2) typically
does not not hold. Consider again the example in antiplane shear:
EXAMPLE 10.1 (Example 8.3 revisited). We now solve the problem in Example
8.3 using piecewise constant elements for h and piecewise linear elements for
Uh. The equation for discrete equilibrium with the load Af is
1 1
h U
h = A uh VUh C Uh
-1 -1
or equivalently
I 1
iO-h =A/ffQ, i= 1,...,N-, (10.3)

~~1 1
where i is the linear "hat-function" at node xi = -1 + ih, h = 2/N, and N is

the number of intervals. Let si denote the (constant) value of ah in the interval
xi- 1 < x < xi for i = 1,.. .,N. Then Eq. (10.3) can be written
Xi+l
S,-Si+l=f, i= l,...,N-1.
Note that this determines -r h uniquely up to an additive constant in contrast

to Example 8.3. Let the index k be determined by the location of a, such that
Xk < a < Xk+1, and let d = a - Xk. (Recall that f(x) = -1 for -1 < x < a and
f(x) = +1 for a < x < 1, and that a > 0.) Then the discrete equilibrium equation
takes the form
Si+l - Si = Ah for i < k,
- 1
Sk+l - Sk = Ah (2hd - d2),
Sk+2 - k+l = -h 1 (2h(h - d) - (h - d) 2 ),

Si+l - Si = -Ah for i >k+1.
The maximum value for A is achieved with sl = -1 and Sk+l = 1:
A*
= 2
h l +a+d-h-h ld2'
The continuous problem was solved in Example 8.3:
* 2
1+a'
The error on A* is
A*-A' =2h 1- (d/h) + (d/h)2

A~2h
-= (1 + a)(1 + a + d - h - d2/h) '
We see that
3
lim inf h A* - lim sup h A
h-O h h 0 h
Hence Eq. (10.1) is satisfied, but Eq. (10.2) is not. The error primarily depends
on the location of singularities relative to the finite element grid.
It is our experience that the behaviour illustrated in the example is typical in
limit analysis. Equation (10.1) usually holds, but Eq. (10.2) does not.
Error estimates on the discrete limit multiplier A are based on the saddle
point properties of the exact and the discrete solutions. In the following
SEcInoN 10 Discretization by finite elements 249
(A*, v*, u*) denotes a solution to the continuous problem, i.e., a triple satisfy-
ing Eq. (5.32), and (Ah, o-(, u*) denotes a solution to the discrete problem in
Eq. (8.1), i.e., a triple satisfying Eq. (8.23).
THEOREM 10.1. Let (A*, or*, u*) and (Ah, a,, u) denote an exact and a discrete
solution, respectively. If or* satisfies the exact yield condition, i.e., o* E K, then
a(oh - oa, uh) < Ah - A* < a(ao, Uh - u*) (10.4)
for all oh C Kh and all Uh e Uh with F(Uh) = 1.
PROOF. In Eq. (5.32) insert or = oa and u = u and subtract from Eq. (8.23).
The condition o- e K in the above theorem will automatically be satisfied, if

Kh C K. Since the yield condition is only imposed at the nodes, this condition
may fail to hold for two reasons: The yield condition may vary between nodes,
and certain element functions may achieve their maximal value between nodes.
Problems caused by changes in the yield condition with x E f2 should be treated
on an ad hoc basis, so consider the case of a homogeneous material. Then
Kh C K will hold for piecewise constant, linear and multilinear elements: For
these elements the convex yield condition will be satisfied everywhere, if it is
satisfied at the nodes. For quadratic and cubic elements this is not the case. The
effect of this on the discrete limit multiplier A depends on the specific elements
for both O'h and uh and on the collapse fields. Under quite weak conditions this
effect can be seen to be of order O(h) (bounded by a multiple of h), where
h is a linear measure of the element size. Hence the effect may be considered
equivalent to a discretization error.
In order to prove convergence of A from Eq. (10.4) we need stability condi-
tions, which bound o* and uh in some norms and regularity conditions, which
guarantee that or* and u* may be approximated by oh h and Uh E Uh, re-
spectively. Examples of such conditions are
IIUIIBD() < C Vh (10.5)
paired with the regularity condition
or* w1 '00(f2) (10.6)
or
|E(u*)llc(n) <c Vh (10.7)
paired with a weaker regularity condition on r*. We may then write the left-
hand side of Eq. (10.4) as
h
-A > fX (hrij - Oij) -(i ). (10.8)
We shall always assume that A-h is a regular finite element triangulation as

defined in CIARLET [1991, Chapter III]. Also, we restrict ourselves to the case
u E BD(2), i.e., Iu(u) = e(u).
THEOREM 10.2. Assume continuous elements for both h and Uh. If Eqs. (10.5)
and (10.6) hold, then
A; - A* > -ch (10.9)
for some positive constant c independent of h.
PROOF. We use Eq. (10.6) and CARLET [1991, Theorem 16.1] combined with
Sobolev's imbedding theorem to conclude that
min Ioh - -*I c(f) < ch
OhE&h
for some constant c. Since oh in Eq. (10.8) is arbitrary, Eq. (10.9) follows from
Eq. (10.5) by definition of the norm on BD(n2) (see Eq. (5.3)). [1
THEOREM 10.3. Assume continuous elements for Uh, but allow alsopiecewise con-
stant elements for o-h. Assume that Eq. (10.7) holds, and that o-r* is "piecewise
W 1,' " in the following sense: n2 can be divided into a finite number of regions
12
k, such that
(a) the components of o-* belong to L(fl) and to Wl'(Q2k) for each 2k;
(b) the total volume of elements overlapping with more than one J2 k tends to
zero as O(h), as h - 0O.
Then there is a positive constant c such that
Ah - A* > -ch
independent of h.
PROOF. We estimate the right-hand side of Eq. (10.8). The integral over elements
contained in one f2 k is O(h) by the argument in the previous proof. The integral
over the remaining elements is bounded by a constant times their total volume,
which by assumption is O(h). E
In order to obtain the opposite bound we need a stability condition on o:

|II-h|XL < C1 and IIVO;hILc < C2 Vh. (10.10)
THEOREM 10.4. Assume continuous elements for oh and uh. If u* is "piecewise

W' " as defined in the previous theorem and Eq. (10.10) holds, then
Ah - A* ch. (10.11)
PROOF. We write the right-hand side inequality in Eq. (10.4) as
Ah * - (V'h)(Uh - U*) + (Oh)(Uh - U).

12Q~ T
SECTION 10 Discretization by finite elements 251
The proof now goes much as the previous proof. [
Condition (10.5) is rather weak, because Iu 1IBD is dominated by the total

energy dissipation rate D(u) defined by Eq. (5.30). (uhIBD need not be dom-
inated by the discrete form Dh(Uh) defined in Eq. (8.1), so Eq. (10.5) is not
guaranteed to hold.) On the other hand condition (10.6) implies that or* is con-
tinuous, clearly a restrictive condition, which can only be checked a posteriori,
and even that may not be trivial. Condition (10.7) is very restrictive. The prob-
lem with these conditions is that they are formulated in terms of either or or
u. Good conditions should reflect the fact that singularity in one is paired with
regularity in the other.
These error estimates are special cases of the following theorem, which fol-
lows directly from Eq. (10.4).
THEOREM 10.5. Assume that Il111 and 11-112 are norms on Z and U, respectively,
such that
la(o, u)l < I1111 llul112 V(o,u) E x U
and let (Ah, ah, u ) be a sequence of discrete solutions with oh E K.
(a) If II hll < C 1, then
A -A* C min ]Uh-U*112.
hF(uh)=l
(b) If IIUh112 < C2, then
A - A* -C 2 min ]orh-T-orll
OhCKh
The above error estimates are one-sided. Clearly the best one-sided estimates
are of the form A > A*or A*< A*, which brings us back to the so-called bound
methods mentioned in Section 3. If only (or u) is discretized, there results
a one-sided approximation of A*. In some cases this may also be obtained with
the mixed finite element methods advocated here:
A discretization h X Uh is a lower bound method if Kh C K (Kh is defined
by Eq. (8.2)), and if the following implication holds for oh C Xh:
If a(h, Uh) = F(Uh) VUh Uh, then a(oh, U) = F(u) Vu C U.
In fact it suffices that the implication holds for or, which will then be an ad-
missible stress tensor in exact equilibrium with the load multiplier A, implying
A* < A*. This property depends not only on the discretization, but also on the
actual forces (f,g).
A discretization Xh x Uh is an upper bound method, if the discrete energy
dissipation rate Dh(Uh) is exact on Uh, i.e., if the following equality holds:
max a(orh, Uh) = maxa(, Uh) Vuh G Uh.
0hGKh oAK
In this case it immediately follows from Eqs. (5.29) and (8.1) that A* > A*. It is
easy to see that if the yield condition is independent of x E 12, then a discretiza-
tion with piecewise linear elements for Uh and piecewise constant elements for
252 E. Christiansen CHAFER III
arh is an upper bound method (see CHRISTIANSEN [1981, Theorem 6]). This was
explicitly seen in Example 10.1.
We do not know conditions, which guarantee convergence of the collapse
fields (or*, uh) to (or*, u ). The fields will in general not be uniquely determined,
and the degree of uniqueness need not be the same for the continuous and the
discrete problems. The second best is to guarantee that if a sequence of discrete
solutions (or;, u ) converges in some sense, then the limit is a solution to the
continuous problem.
THEOREM 10.6. Consider a sequence 0a of h-values satisfying Xhl C Zh2 C and

Uhl C Uh2 C U for h2 < h. Assume thatKh = K n hh E and let C = {u E
U I F(u) = 1}. Let (Ah, o-, u ) be a sequence of solutions to the discreteproblem
(8.23) for h E X. Assume that A - hAand that (o;,u ) converges to (o, u o ) E
K x C in the following weak sense:
lima(o, u) = a(r , u) u E U, (10.12)
lima(o, u;) = a(o, u ) or cE . (10.13)

h-O
Assume furthermore that (aO,u) satisfies the following weak approximation

property:
Vu E C Ve >0 3h E C 3Uh E C n Uh: la(, uh - u) < , (10.14)
Vo-r B Ve > O h c 3oh K Zh: a( oh - , U0 ) < . (10.15)
Then (A , ro, u 0) is a solution to the continuous problem (5.32).
PROOF. For fixed hi E e and any ohl c K n hl, it follows from Eqs. (10.13) and
(8.23) that
a(orh,, U) = limAa(rh,uh)
h-0
S limAh
h-0O
= AO
From Eq. (10.15) we conclude
a(o, u ) < A0 Vo E K.
Similarly it follows from Eqs. (10.12) and (8.23) that for h i C C and any u h, E
SECTION 10 Discretization by finite elements 253
C n Uh,
a(o ,Uhl) = im a(o,* Uhl) = imh = AO,
h-*0
which combined with Eq. (10.14) implies
a(o , u) = A u C.

Hence the triple (A , or, u) satisfies Eq. (5.32). [1
CHAPTER IV
Solution of the Discrete Problem
11. General remarks
Using a mixed finite element method the problem of limit analysis was dis-
cretized to the dual forms Eqs. (8.17)-(8.18):
A = max{A I 3x E Kd: Ax = Ab} (11.1)
= min{Dd(y) bTy = 1}, (11.2)
where A, b, Kd and Dd are defined by the discretization Eqs. (8.2)-(8.20). Recall
that in general Kd is convex in RN and contains the sum of a ball and a subspace
as illustrated in Example 8.2. Then Dd(y) is finite, only if ATy is orthogonal to
this subspace.
We shall rewrite the discrete problem in standard mathematical programming
notation. The discrete static formulation Eq. (11.1) will be referred to as the
primal problem:
max A,
(P): Ax = bA, (11.3)
x E Kd.
The discrete kinematic formulation Eq. (11.2) is the corresponding dualproblem
(D): minDd.
bTy 1. (11.4)
The problem in Eqs. (11.3)-(11.4) is an interesting mathematical optimization

problem in its own right. Even for moderate size grids it is non-trivial, and
for realistic problems and grids it is a real challenge. Not only is it large and
sparse; in our experience it is also ill conditioned. We do not claim to have a
"perfect" method yet. In fact it is not obvious which of our options, if any, will
survive as "the best method". However, it does appear that we are able to solve
realistic problems in limit analysis with reasonable accuracy. At least the limit
multiplier and one of the collapse fields can be determined, but we cannot always
determine collapse fields for both stresses and plastic flow simultaneously.
We shall illustrate the methods on a classical test problem in plane strain.
Three-dimensional problems can also be solved, but the plane strain problem
255
256 E. Christiansen CHIAPTER IV
X2 I
U
1
a =
U1=0 I
.'
I A
U2 =U L
FIG. 11.1. The test problem reduced by symmetry.
with its two-dimensional convenience offers the typical difficulties of the gen-
eral case: Unbounded yield set, ill conditioning, non-uniqueness, etc. Our test
problem has been used in laboratory tests and reported in text books on ma-
terial science. It was presented to this author by professor McLintock, MIT, in
1975 as a challenge to the numerical techniques in limit analysis. In this sense
the problem serves two purposes: Numerical limit analysis has come of age, and
realistic problems can be solved. Furthermore, the problem is well suited to
reveal the virtues and shortcomings of potential solution methods.
EXAMPLE 11.1 (The reference test problem). The problem is described in Exam-
ple 1.1. In the plane strain model a rectangular block of material is given ex-
ternal symmetric thin cuts of variable depth. The volume forces are zero, and
the surface forces consist of a uniform pull at the ends. For symmetry reasons
the problem is reduced to the first quadrant (see Fig. 11.1). Actually, due to
non-uniqueness non-symmetric collapse fields exist, but by reflection and super-
position it is easy to see that there will always be a symmetric collapse solution.
The geometry of the reduced problem is seen in Fig. 11.1. The boundary con-
ditions are
u2 =0 forx 2 0,
=
ul = 0 for xl = 0 and 0 < x2 < a.
The Von Mises yield condition for a homogeneous material in plane strain is
normalized to
2 (11.6)
1(O-n -T22) + 22 < 1.
The problem is analysed in MARTIN [1975, Section 12.3] and KACHANOV [1971,
p. 195]. The following bounds are given for the case a = , i.e., the cuts extend
two thirds to the center of the bar:
4
2
<A* < 3
These bounds are based on exact upper and lower bound solutions. We find the
true value for this case to be near A* = 0.92 (see Table 12.1 and Table 14.1).
SECTION 11 Solution of the discrete problem 257
Only computations with uniform finite element meshes are performed. It will
be obvious from the results that the use of adaptive mesh refinement is even
more useful in limit analysis than for elastic problems, due to the rigid regions.
However, since there is ill-conditioning and non-uniqueness in limit analysis,
caution must be taken not to influence the collapse solution by choosing a mesh
with built-in preferences.
It should be emphasized that the results reported here, with one exception
which will be pointed out, have been obtained with general purpose linear
programming or nonlinear optimization software.
The following example shows a complete discretization of the test problem
using piecewise constant element functions for the stresses and bilinear elements
over rectangles for the velocity. This element combination has provided the best
results so far. Unfortunately it can only be applied in simple geometries. The
more generally applicable linear-linear element combination over triangles also
works, but in our experience it does not determine the stress field in the collapse
state quite as well.
EXAMPLE 11.2 (The test problem with constant-bilinearelements). Consider the

test problem of Example 11.1. We describe the discretization with piecewise
constant element function for the stresses and bilinear elements over rectangles
for the flow.
Let -h be a uniform division of the domain
= {(Xl,X 2) 21 O Xl
x < L,0<x2 < 1}
into squares with side h = 1/N. Each component of the discrete stress tensor is
a linear combination of the functions 4,, where do, is equal to 1 in the element
(square) with index v and equal to 0 elsewhere. In the notation of Section 8 the
"nodes" for the discrete stress space -h are the squares of -h The basis for Yh
consists of the tensors
1l R[V 0] 22 [ 0 12 [0 v
V, , = ' V 01
The discrete tensors tOh C h are of the form
= (4l 11 +2 222
++ 212). , (11.7)
where 6j' is the constant value of h,ij in the square with index v. The yield
condition (11.6) takes the form of the following explicit set of nonlinear in-
equalities:
1( 52
_ 22)21+ (2)2 < VV e (11.8)
The node set Nu for the flow is the set of vertices of 3-h' Let C denote the
258 E. Christiansen CHAPIER IV
bilinear element function with value 1 at vertex , E JYVand 0 at all other nodes.
The basis for the discrete flow space consists of the vector functions
The discrete flow uh E Uh can be written
uh= E3 (4+' 1iqt),

(q1 2 )(1 (11.9)
where k is the value of Uh,k at the vertex with index Ak.Due to the bound-
ary condition in Eq. (11.5) the following components n k are fixed to zero, or
equivalently the corresponding basis functions fk are absent in Eq. (11.9):
2 = 0 for vertices t with x2 = 0,

(11.10)
l = 0 for vertices / with xl = 0 and 0 < x 2 < a.
In this example there are no volume forces,f = 0, while the surface forces are
given by g = (1, O0)at the right-hand boundary with xl = L and g = 0 elsewhere.
Hence the terms in Eq. (8.7) are:
F(fI<l) = h for vertices with xl = L and 0 < x2 < 1,
F(1)= h for the vertex with xl = L and x2 = 0,
F(p&) = h for the vertex with xl = L and x2 = 1,
F(p ) = 0 for all other vertices,
F(2,2) = 0 for all vertices.
These are the components of b in Eqs. (11.3)-(11.4).

The integrals in Eq. (8.8) are easily computed using that the bilinear basis
functions 4i (x,y) are of the form '(x, y) = rl(x)k(y), where ifj is a standard
linear element function in one variable. Recall that
2
a(of, = E J a.
i,j=1n
Since each basis function O 'j is only different from zero in element v, there is
nothing to assemble in the stiffness matrix Eq. (8.8), so we may write down the
matrix elements directly. (Special for piecewise constant elements.) Let v be
any of the squares in J-h and let be the vertex of the lower left-hand corner.
Then we find:
a(-b1, aqi) =-h, a(,l', p2) = 0,
1 2
a(2, m ) =0, a(v2, ) =-2
a( V I~L ) =-1h,
2, q'1 2 a(02, f,2 ) = -1h.
V I~f~t 2
SECToN 12 Solution of the discrete problem 259
The values for the other three vertices are the same except for permutations
and sign variation. As usual this part of the discretization is easier to do than
to describe.
It remains to number the variables ((i) and (o/k), or equivalently to choose
the order of columns and rows in the matrix A in Eq. (8.15) or Eq. (11.3). Since
the length satisfies L > 1 the elements and vertices are both ordered "bottom-
up, left-to-right" to obtain smallest band width of A. For the same reason the
stress components 'ij are numbered element by element, i.e., the mapping (8.9)
is given by
n(v,1, 1) = 3(v - 1) + 1,
n(v, 2, 2) = 3(v - 1) + 2,
n(v,1, 2) = 3(v - 1) + 3,
where v is the element number. The column corresponding to A, i.e., the b
vector found above, is added as the last column, i.e., the variable A gets the
number n = 3N 2 L + 1. Similarly the velocity components are numbered by the
mapping (8.12), which we define as follows:
2
m(gt, 1) = (At - 1) + 1,
2
m(At, 2) = (gt - 1) + 2.
Modern mathematical optimization software permits incorporating the bound-
ary conditions on u by fixing the variables in Eq. (11.10). Otherwise the number-
ing m(A, k) must be modified to leave out the corresponding rows (see CHRIS-
TIANSEN and KORTANEK [1990]).
This completes the setup of the convex programming problem Eqs. (11.3)-
(11.4). The following sections will describe several ways to solve it.
Other element combinations, such as bilinear-bilinear or linear-linear, for
which we shall report results, are implemented with small modifications to the
above example. Preferably one should try to program the logic and leave the
details to the program, as is the standard now for elasticity problems.
12. Linear programming methods
Traditionally the discrete problem has been solved using linear programming.
The main reasons are:
* Many structural problems can adequately be formulated with linear yield
conditions.
* Good software, based on the simplex method, was available for linear pro-
gramming, long before similar software for convex programming was com-
petitive.
* The duality between the static and kinematic formulations was first estab-
lished within the linear programming frame.
260 E. Christiansen CHAPTER IV
The price to pay for these virtues is, at first glance, modest. A new error is
introduced, when the convex yield condition x c Kd is replaced by a piecewise
linear condition Bx < c. This error is easily controlled and within the typical
tolerance of applications.
After linearization of the yield condition the problems (11.3)-(11.4) turn into
a pair of dual LP-problems:
max A,
(Pin) Ax= bA, (12.1)
Bx < c.
The equality constraints correspond to the equilibrium equation, while the in-
equalities are the linearization of the yield condition. The dual of the problem
in Eq. (12.1) is
min c T Y2,
ATy = BTy 2,
(Diin): bTy = 1, (12.2)
Y2 ) 0.
Of course Eq. (12.2) can also be obtained directly from Eq. (11.4) by inserting
the linear yield condition Bx < c in the definition of Dd(y).
Either of these forms may be fed to your favourite LP-solver. With various
problems and LP-codes both forms have turned out slightly better conditioned
than the other, so we have no recommendation as to which form should be
used. Adding linear constraints in the yield condition corresponds to adding
rows in Eq. (12.1) and columns in Eq. (12.2). This should favour the dual form
Eq. (12.2), but our experience is not unanimous.
Before reporting on computational results we shall try to describe in detail
how the linearization of the Von Mises yield condition in plane strain may be
performed.
EXAMPLE 12.1 (Linearized Von Mises condition in plane strain). The normalized
Von Mises yield condition in plane strain is
2
El: 4( 11 - 2 2) + 2 -2 1. (12.3)
El is an ellipse in the (al - c022, r12) plane. It is shown in Fig. 12.1, which also
shows the approximation of the ellipse by the intersection of 16 half-planes. Let
this approximating convex polygon be denoted Elin- Clearly E1 C Elin. Let E3
denote the inflated ellipse (for /3 > 1):
2 2
Ep: 4(01Io - 2 2) + 022 . (12.4)
For a given number of lines we want to choose these lines in such a way that
Elin C Ed for /3 as small as possible. Then the following relation between the
R '7 1.
___ TI _6
_
II2< 4
aE - aYY
FIG. 12.1. Linearization of the Von Mises yield condition in plane strain.
exact limit multiplier A*, corresponding to El, and the approximate value Ali,
computed from Elin, will hold:
,Ahi A<* A-<

. (12.5)
The optimal choice of K lines can be constructed with a simple trick (CHRIs-
TIANSEN and KORTANEK [1991, p. 52]): Inscribe the unit circle in a regular K-
polygon. Then rescale by a factor 2 in the direction of the 0o' - 022 axis. The
circle maps into the ellipse E1 , and the regular polygon maps into the optimal
circumscribed polygon. It is convenient to let the number of lines K be divisible
by 4 and include the four lines parallel to the axes. The error is estimated by
Eq. (12.5) with the following values of /3:
/3 = for K = 4,
/3 = 1.0824 for K = 8,
= 1.0196 forK=16.
The linear constraints are explicitly given by
2~rk 2'rrk
(0J11 - 22) COS K + 212 sin K < 2, k =0,..., K - 1. (12.6)
For K = 16, Eq. (12.5) gives an error interval of about 2%. This exceeds the
needs for most applications, but is adequate for testing. In three space dimen-
sions it will clearly be necessary to settle for less.
At any given point of the material at most two of the constraints in Eq. (12.6)
can be active. Based on physical considerations we expect that most of the
constraints can be omitted in large regions of the material. Unfortunately this
seems to hold only for certain solution methods.
EXAMPLE 12.2 (LP-formulationfor the test problem). Consider again the test
problem in Example 11.1 discretized with constant-bilinear elements as in Ex-
ample 11.2. In the LP-formulations (12.1)-(12.2) A and b are the same as found
in Example 11.2. The yield condition (11.8) is replaced by the linear inequalities
(12.6):
os 2 + 2 sin K 2 (12.7)
for k = 0, ... K - 1 and all squares v in 3-h. These are the constraints Bx c in
Eq. (12.1). The relation between x and (>?), i.e., the numbering (8.9) of variables
is specified in Example 11.2. This completes the setup of the LP-problem (12.1).
We shall add some efficiency considerations.
Traditional LP-codes may handle the problem (12.1) more efficiently if we
help to minimize fill-in during the solution process. Each row in Eq. (12.7) only
involves the variables associated with one element. Hence fill-in will be reduced,
if this row is numbered adjacent to the velocity components associated with
the vertices of that element. Similarly, if the problem must be brought to LP-
standard form, the same consideration should be given to the slack variables.
However, modern LP-codes for sparse problems may do this automatically.
In the remainder of this section we shall present results for the test problem
obtained with three linear programming algorithms: The simplex method, the
primal affine scaling algorithm and the dual affine scaling algorithm. The order
of presentation is the chronological order of appearance and of application in
limit analysis. It also happens to be the reverse order of performance in our
evaluation.
12.1. Solution with the simplex method

Solutions for the test problem in Example 11.1 using the simplex method are
reported in CmSTIANSEN [1981]. Figure 12.2 visualizes the collapse flow field by
the method described in Example 1.1: The velocity is multiplied by a suitable
small time, and the resulting deformation is plotted. The result is shown for
a = and a = 2, both with h = .
The solution in Fig. 12.2 is obtained with piecewise bilinear elements for both
stresses and flow and with 8 lines in the linearized Von Mises yield condi-
tion (K = 8 in Eq. (12.6)). The flow field and the resulting deformations are
physically acceptable, although there are some non-physical fluctuations in the
xl-coordinate in the upper left corner, which must be attributes of the solution
method. In our experience such fluctuations in those parts of the material, which
are expected to be rigid, tend to appear in the numerical solutions. Neverthe-
less, the solution clearly indicates rigid non-deformed regions in the upper part
and at the lower left corner separated by a relatively narrow zone, where all
the deformation occurs. This is the plastified (or "plastic") region.
The collapse stresses are not indicated on Fig. 12.2. They are not physically
acceptable for this method. The value of shows strong variation from node to
node. The nodes where o' h is at the yield surface are distributed over the whole
SEcnON 12 Solution of the discrete problem 263
FIG. 12.2. Collapse deformation obtained with bilinear elements for stresses and flow. Linearized
yield condition, simplex method.
domain in unsystematic fashion. We know from the principle of complementary

slackness (see discussion following Eq. (5.47)) that ao must be at the yield
surface in the plastified region, and that is indeed the case for this discrete
solution. However, unsystematically scattered "plastified nodes" in the rest of
the domain can only be attributed to the numerical solution method. There are
two reasons for this deficiency:
(a) With this discretization, i.e., the choice of finite element spaces for
stresses and flow, the discrete stresses are undetermined by the equilibrium
equation Ax = Ab and the variational principle. Loosely speaking, there are too
many primal variables (degrees of freedom for h) relative to the number of
dual variables (degrees of freedom for Uh). This can be changed by using alter-
native finite element spaces, but then of course the opposite situation may (and
to some extend will) occur.
(b) The simplex method, being an extreme point method, has a strong pref-
erence for solutions with variables at their bounds (non-basic variables) and
inequalities at the limit of their range. Even if a whole face of the feasible
set is optimal, the simplex method will choose an extreme point. This "arbi-
trary" selection has no mechanism to prevent non-physical fluctuations from
node to node. In a manner of speaking, nature handles non-uniqueness far bet-
ter than the simplex method does. Based on our experience with interior-point
LP-methods and with convex programming methods (see later) we maintain
that this phenomenon is mainly caused by the simplex method, and only to a
lesser extent by the fact that we have replaced convex constraints by piecewise
linear constraints. The extreme point nature of the simplex method, responsible
for several of its virtues, is in this context a serious drawback.
Figure 12.3 shows the collapse velocity fields for a = and h- , obtained
with bilinear-bilinear elements (left) and constant-bilinear elements (right).
-==
Ir i V -r
~N_
r_
\
!N,
N
,\
7Y_ : I
-N1
SIZ_
K I FI 1 \ _181
. I,
K_
IN
'N x I II N . ;s - \
-V-i:~~~~~~~N
\N
'I- 'K 7K
V XN X
v 'K 7K
ij
t
N
7~\ --I 1 I I
_ V V\X L - - -
II _- N-'
~iK7VI
-
I - \- 4- I -
f'
--S
FIG. 12.3. Collapse velocity obtained with bilinear-bilinear elements (left) and constant-bilinear
elements (right). Linearized yield condition, simplex method.
(Hence Fig. 12.3, (left) corresponds to Fig. 12.2, right, only with a different value
of h.) The constant-bilinear element combination has fewer primal variables
(3N 2 compared to 3(N + 1)2 before) with the same number of dual variables
(2(N + 1)2 minus ((N + 1) + N + 1) lost to the boundary condition). However,
the collapse stresses are still not sufficiently well determined to reveal physical
information. On the other hand, the velocity field now also shows fluctuations:
The velocity alternates from node to node between zero and values, which seem
to reveal physical information. The human eye will easily "recognize" three rigid
regions, but that is an interpretation. This is even more obvious in Fig. 12.4 with
the solution for L = 2, a = and h = 1.
With bilinear-bilinear elements the discrete collapse flow u appears to be
unique (different codes give identical solutions), while ol is not. With constant-
bilinear elements neither u, nor ow- is uniquely determined.
In spite of the above, the constant-bilinear element combination cannot be
rejected. With other optimization algorithms it will provide us with the best
collapse fields, and for all algorithms (also the simplex method) it results in an
optimization problem, which is more sparse and better conditioned than any
other element combination. This means that computations with finer grids are
possible.
Due to the poorly determined stresses it was not possible to eliminate con-
straints from the linearized yield condition Eq. (12.6) and thereby reduce the
size of the LP-problem. The largest problems solved with the simplex method
in CHRISTIANSEN [1981] are with a 24 x 24 grid and length L = 1. The limiting
factor with sparse LP-codes is neither computing time nor memory size. For
fine grids the problem becomes so ill-conditioned that the LP-program cannot
improve the objective value or verify optimum.
Based on the limited experience reported here we reluctantly conclude that
-r -L
4X
L'
'- NZ
\4 7-
-e I -e i-e
--
'x- N
\\
NN 1
N
--
N C;-
'x
--e --
\\
1 --e --t
N4 N'
-r -t
V_
'x
N
N
- --e
-e -e
'x
'v
x~ i-e
-- -t
-- l - -- B
l1
|
- - -
-
L_
I --
-e -L1
-e --L
FIG. 12.4. Collapse velocity for L = 2 obtained with constant-bilinear elements. Linearized yield
condition, simplex method.
after playing a key role in the development of limit analysis the simplex method
is not well suited for problems in limit analysis for continuous media.
12.2. Solution with the primal affine scaling algorithm

We claimed above that the extreme point nature of the simplex method is re-
sponsible for its inability to determine the collapse fields for the test problem.
In the wake of KARMARKAR [1984] several classes of interior point methods
have surfaced. The simplest of these, the primal affine scaling algorithm, can be
thought of as an affine version of Karmarkar's polynomial time algorithm, but
was suggested as early as DIKIN [1967]. The algorithm takes the LP-problem on
standard from, maintaining the equality constraints and keeping all variables
strictly positive. In the process of convergence the feasible point is forced to-
wards an optimal point, which may be an extreme point or, if a face is optimal,
a point on an optimal face. This is the difference to the simplex method; the
solution will only be an extreme point, if it is necessary in order to be optimal.
In KARMARKAR [1984] focus was on the efficiency of the interior point methods
compared with the simplex method. In the context of limit analysis we are more
interested in the qualitative properties of the method as consequence of the
interior point versus extreme point principle.
For a description and analysis of the primal affine scaling algorithm we re-
fer to BARNES [1986], VANDERBEI, MEKETON and FREEDMAN [1986], KORTANEK
and MIAOGEN [1987], TSUCHIYA [1991]. The application to our test problem is
reported in CHRISTIANSEN and KORTANEK [1990, 1991].
This algorithm makes it possible to determine collapse fields for both stresses
and flow. With the piecewise constant-bilinear element combination and 16 lines
i
0 U)~~~~
~ii~'ll~'llllii
I
I i
I
iliT7
[
-1
ri' '' II "I

-i
I !, ''
I ~MIIIiII
I 0
1
i
]
2 T
I CFL all~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
I
iI
D1
- I
-
rK
II
.
: l
4_
T
E222 -L4
-44-
R
,
I
iI I I I I II I III
I I I I
I_
I. 2.4 .222i
- -
co
EI
I- I
~ 0
Oioo
a a
. . . . .
I22A:
-,
I I I I I I I I I I I I I I I
FIG. 12.5. Collapse fields obtained for a= and a= , = and constant-bilinear elements.
Linearized yield condition and the primal affine scaling algorithm.
(K = 16 in Eq. (12.6)) only the linear constraints with numbers 0-6 and 15 in
Fig. 12.1 are active. The other half of the constraints may be removed (they
are checked a posteriori), a considerable reduction in problem size. Figure 12.5
shows the collapse fields for K = 16, L =1, h = 21, a = and a = 2. The stress
field is visualized as follows: If an element is considered plastified, the number
(referring to Fig. 12.1) of the active line or lines (at most 2) is written inside the
element. For example Tll - 022 = 2 in squares with the number 0, and a12 = 1 in
squares with the number 4. In elements without a number the stresses are not at
the linearized yield surface. Since the algorithm generates a sequence of interior
points, the yield surface is never reached exactly. It is therefore necessary to
interpret a constraint to be active if the computed stress tensor is close enough.
More precisely, we say that the linear constraint Eq. (12.7) is -active, if
(11- 22 2irk
cos -
12 . 2rrk
+ 22 sin > 2 - 6. (12.8)
K K
Figure 12.5 shows the 8-active lines for = 10- 3 . The picture is almost the same
for 8 = 10 2 and 8 = 10 - 4 .
The collapse fields indicated in Fig. 12.5 are physically acceptable. There are
no plastified elements scattered outside a well-defined region, and there are only
very small fluctuations in the velocity. The values for the collapse multiplier Ah
were computed for N < 30, h = 1/N. They coincide with the results for the
method in the next subsection and will be given later.
For the linear-linear element combination over triangles the stresses are not
so well determined. Plastified nodes are scattered outside the plastified region,
and no linear constraints can be removed. On the other hand the collapse ve-
locity is without any fluctuations.
As is the case for the simplex method the limiting factor for the problem
size, which can be solved, is the deterioration of accuracy near the solution.
The primal affine scaling algorithm is particularly sensitive to loss of accuracy.
In each iteration the scaled gradient must be projected onto the null space of the
linear constraints in order to satisfy the equality constraints in the LP-standard
form. Loss of accuracy in this projection accumulates and is not recovered by
the algorithm. This is a weak point of the primal affine scaling algorithm.
It is to be expected that a linearization of convex constraints will cause ill
conditioning near the optimal solution, and a natural conclusion is to abandon
the linear approach altogether. However, we shall see that the ill-conditioning
is also present with convex constraints. And in the next subsection we present
another linear method, which seems to handle the problem extremely well.
12.3. Solution with the dual affine scaling algorithm

We find the following results surprising both from a mathematical programming
and from a limit analysis point of view. The reason is that an algorithm, which
should not work at all on this problem, provides the best results so far.
The dual affine scaling algorithm is essentially the same as the primal algo-
rithm just described, applied to the dual of the standard LP-form. Consider the
general LP-problem on standard form
minc Tx,
Ax = b, (12.9)
x 0.
The dual problem is
max bTy,
(12.10)
ATy c
or equivalently
max bTy,
ATy +s = C, (12.11)
S >0,
where s is the vector of slack variables. Equation (12.11) is of the same form
as Eq. (12.9), except that only some of the variables, namely s, are restricted
to be nonnegative. The components of y are free. The primal affine scaling
algorithm can be modified to handle this form: Only the slack variables are
subject to scaling, and the algorithm generates a sequence of points (y, s) with
the s-components strictly positive.
An obvious advantage of the form (12.11) over (12.9) is that it suffices to
compute a new y at each iteration, which satisfies ATy < c. Then we may define
s = c - ATy > 0, and the linear constraints will be satisfied with no accumulated
loss of accuracy. Hence the modified primal algorithm applied to Eq. (12.11) is
less sensitive to round-off error than the primal algorithm applied to Eq. (12.9).
Now, if some variables in the original problem are free (i.e., unrestricted
in sign; by duplication of variables such problems may still be written in the
form (12.9)), then a corresponding number of inequalities in Eq. (12.10) will be
inequality pairs of the form
aTy < Cj, aify > Cj.
The corresponding slack variables will necessarily be zero. There are no feasible
points for Eq. (12.11) with s strictly positive. We say that the feasible set for
Eq. (12.11) has empty interior, and in this case the primal algorithm simply does
not apply.
In limit analysis most of the variables are free, so how can the dual affine scal-
ing algorithm be applied? Like the simplex method, the interior point methods
need a "phase 1", where an interior feasible starting point is found. This is done
by solving an associated problem with an objective function, which will force the
variables towards feasibility. In ADLER, KARMARKAR, RESENDE and VEIGA [1989]
the objective function for the original problem is included with some weight in
phase 1 in order to produce a good starting point, such that fewer iterations
will be necessary in phase 2. This idea has been improved in ANDERSEN [1993]
to the extent that the optimal solution (to the original problem) will be found
in phase 1, if the feasible set has empty interior. So phase 2 never starts!
The dual algorithm just described is very efficient in general (ANDERSEN
[1993]). It is still amazing that, for the linearized problem of limit analysis,
phase 1 alone works so well (ANDERSEN and CHRISTIANSEN [1995]). It is more
accurate, more efficient in use of cpu-time and memory than other methods,
including smooth convex optimization algorithms.
Results for the test problem obtained with (phase 1 of) the dual affine scaling
algorithm are seen in Figs. 12.6-12.8. The piecewise constant-bilinear element
combination is used, and all 16 lines of Eq. (12.7) are included. With this algo-
rithm it turns out that there is no significant saving by removing lines. Figure 12.6
shows the case a = with h = ii, i.e., a 198 x 198 grid. (All other algorithms,
which we have tried, including convex methods, have had serious problems with
ill-conditioning no later than h = .) The top part shows the collapse deforma-
tion, and on the bottom part the stress field is visualized. The plastified elements
are indicated by a "dot" instead of line numbers, for obvious reasons. The ac-
tive lines are as one would expect from Fig. 12.5. The field for both stresses and
flow are informative, although the rigid regions show the same fluctuations in
the flow as observed earlier. With the linear-linear element combination these
fluctuations are not seen, but the stresses are not so well determined, and the
problem is more ill-conditioned allowing only more coarse grids.
On Figs. 12.7 and 12.8 corresponding solutions for the cases a = 2, h = 15
and L a-2 ,h= h are shown. The "extra" band of plastified elements
and the "extra" deformation layers for L = 2 are unexpected, but it must be
recalled that the exact collapse field is probably not unique.
SECrnON 12 Solution of the discrete problem 269
FIG. 12.6. Collapse fields for L = 1, a = and h = 18 obtained with the linearized yield condition
and the dual affine scaling method.
270 E. Christiansen CHAFPTER IV
FIG. 12.7. Collapse fields for L = 1, a = and h = 10 obtained with the linearized yield condition
FIG. 12.8. Collapse fields for L = 2, a = and h = obtained with the linearized yield condition
Clearly much effort used to compute the rigid regions could be saved by using
adaptive elements. That would almost certainly influence the computed collapse
fields, and numerical experiments are needed.
The computed values for the limit multiplier are shown in Fig. 13.2 (open
triangles) together with values from other methods. The values with the same
yield condition, but the linear-linear element combination, are also shown (black
triangles), but only for coarse grids. For finer grids they coincide with the values
from the constant-linear element combination. They should converge to the
same limit, which also appears to be the case. On the other hand, the values
272 E. Christiansenl CHAPrER IV
TABLE 12.1. Results for L = 1, a = , and 2 obtained with the dual affine scaling algorithm. The
computed convergence order and the values from extrapolation to order 1 are also shown.
- 2 3
Ah eeorderextr. A order extr. h order h xtr.
12 0.9826 1.2117 1.4826
24 0.9584 0.9342 1.1776 1.1434 1.4400 1.3973
36 0.9493 0.79 0.9311 1.1656 0.90 1.1415 1.4250 0.90 1.3949
48 0.9448 1.01 0.9312 1.1596 1.03 1.1417 1.4177 1.11 1.3960
60 0.9421 0.97 0.9311 1.1561 1.07 1.1421 1.4134 1.03 1.3962
72 0.9403 1.19 0.9315 1.1537 0.82 1.1416 1.4104 0.73 1.3952
84 0.9391 1.08 0.9316 1.1519 0.97 1.1415 1.4082 1.08 1.3954
96 0.9381 1.08 0.9317 1.1506 0.96 1.1414 1.4066 1.00 1.3954
150 0.9359 1.10 0.9319 1.1473 1.02 1.1415 1.4025 0.90 1.3951
198 0.9349 1.06 0.9320
computed with the linearized yield condition with 16 lines are the upper bound
of an error interval of length 2%, independent of h, relative to the values for
the exact convex Von Mises condition, also shown in the figure (open circles).
Tables 12.1 and 12.2 list some computed values for Ah, experimental con-
vergence orders based on these values and the extrapolated values to order 1.
(For definition and method see CHRISTIANSEN and PETERSEN [1989].) Both the
computed orders and the extrapolated values indicate that Eq. (10.1) holds, but
that Eq. (10.2) does not. The values for order and extrapolation in the table are
computed from the Ah-values in the table (although we have computed many
more values). If we had used all Ah-values, it would have been more apparent
that Eq. (10.2) does not hold. Larger intervals between the h-values has the
effect of emphasizing the dominating linear term in the error.
13. Convex programming methods
We now return to the convex programming problem in Eqs. (11.3)-(11.4). In this

section we concentrate on the form (11.3), which is a smooth convex program-
ming problem: The objective function is linear, and the constraints consist of
sparse linear equalities (equilibrium) and very sparse smooth nonlinear inequal-
ities (the yield condition). The algorithm of RoBInSON [1972] and, in particular,
the implementation of MURTAGH and SAUNDERS [1982], commercially available
in the program package MINOS, is designed to solve such problems efficiently.
The discretization of the test problem with the piecewise constant-bilinear
TABLE 12.2. Same as in Table 12.1 for L = 2.

-1
h a= 3 a=33
A
h* order extr. A* order extr.
12 1.2012 1.6667
24 1.1713 1.1414 1.6662
36 1.1607 0.89 1.1394 1.6614 1.6517
48 1.1552 0.93 1.1389 1.6588 0.80 1.6510
60 1.1519 0.93 1.1386 1.6570 0.42 1.6498
69 1.1502 0.99 1.1386 1.6557 0.56 1.6491
99 1.1465 0.87 1.1382 1.6537 0.61 1.6484
Ah order extr.
10 1.4722
20 1.4481 1.4240
30 1.4374 0.50 1.4161
40 1.4316 0.73 1.4140
50 1.4280 0.92 1.4136
60 1.4256 1.06 1.4138
70 1.4239 1.01 1.4138
100 1.4209 0.93 1.4137
element combination is described in Example 11.2. MINOS solves the test prob-
lem more efficiently and with at least the same accuracy as any simplex code,
which we have tried on the linearized problem. (One of these was the LP-mode
of MINOS.) Were it not for the results obtained with the dual affine scaling
algorithm described in the previous section the LP-approach would hardly be
justified. Now the situation is not clear.
MINOS was applied both to the constant-bilinear and the linear-linear el-
ement combinations. In contrast to our experience with the linearized yield
condition, the discretization with linear elements does not give rise to a more
ill conditioned problem for MINOS than the constant-bilingear combination.
We could compute with the same element size (30 x 30 grid). The values for
the limit multiplier are equally good, the collapse stress field is almost as good,
and the collapse deformation field has no fluctuations in the rigid regions.
The collapse fields for the test problem with piecewise linear elements for both
274 E. Christiansen CHArER IV
FIG. 13.1. Collapse fields for the test problem with a = , h = obtained with MINOS for the
exact Von Mises condition and linear-linear elements.
stresses and flow, obtained with MINOS, are shown in Fig. 13.1. The stresses are
indicated by drawing the vector ((o-xx - yy),) rXy) at each node. The defor-
mation field is displayed as described earlier. The stress field is not as well
determined as with the interior point methods for the linearized problem. The
stress tensor fluctuates in a non-physical manner from node to node in the
(expected) rigid regions. Again these fluctuations are caused by the extreme
point nature of the algorithm. MINOS solves a sequence of linearly constrained
problems operating with a basis as in the simplex method. The nonlinear ob-
jective function (in our case the nonlinear terms come from the linearization
of the constraints) eventually may force variables to be superbasic, i.e., neither
basic nor equal to zero. Since the optimization problem in limit analysis is ill
conditioned and has a very "flat" optimum, the algorithm cannot force enough
variables to be superbasic. It would be very interesting to try a projected La-
grangian algorithm as in MINOS, with the principle of basic variables replaced
by the idea from interior (or exterior) point methods.
The values for the limit multiplier obtained with the exact Von Mises condi-
tion are seen in Fig. 13.2 (open circles).
EXAMPLE 13.1. The linear-linear element combination was used to discretize
the plane strain problem indicated in Fig. 13.3. The triangulation is obtained by
drawing a diagonal in each quadrangle. The structure is subject to a uniform
SEcnoN 13 Solution of the discrete problem 275
.*
h.
1.4
r~~~~~~~~~ 8
WV 8
1.2
VB S
1.0 -
0.8 -
I' I I I h
40 20 10 6 4
FIG. 13.2. Computed values for the limit multiplier for the test problem with L = 1, a = , a =
and a = : V linearized Von Mises condition (12.7) with K = 16, and constant-bilinear elements. 2
linearized Von Mises condition (12.7) with K = 16, and linear-linear elements. O exact Von Mises
condition and constant-bilinear elements. * exact Von Mises condition and the elements described
in Section 14.
vertical load on the top surface. The weight of the material itself is ignored, but
might just as easily have been included as described in Section 6. The boundary
conditions are u = 0 at the lower edge of the frame, and u = 0 at the left and
right sides of the frame. Thus the structure is assumed to be part of a larger
structure, which is obtained by reflection in the left and the right side of the
frame. The resulting problem was solved by MINOS for two different yield
conditions: The Von Mises condition (11.6) and the Coulomb condition (6.10).
The discrete collapse stress field for the Von Mises condition is indicated in
Fig. 13.3 by the same method as in Fig. 13.1. We can identify regions with com-
pression, tension and shear, but there appears to be non-physical fluctuations in
the stress field, causing the stress tensor to be at the yield surface at too many
nodes. The corresponding collapse flow field is visualized in the usual way in
Fig. 13.4.
The discrete problem with the Coulomb yield condition turns out to be more
ill conditioned than with the Von Mises condition, and more computational ex-
perience is clearly necessary. The collapse solution for the Coulomb condition
is seen in Fig. 13.5. In order to avoid the non-physical fluctuations, the stresses
have been computed at the center of each diagonal. The stress tensor at this
point will be at the yield surface, only if it is at the yield surface, and has almost
I
FIG. 13.3. Original geometry and collapse stresses obtained with MINOS for the exact Von Mises
condition and linear-linear elements.
identical values, at the adjacent nodes. This has the effect of "purifying" the
plastified region and must be thought of as an interpretation, not as a compu-
tational result.
14. Divergence-free elements
In Section 9 we described how the reduced form with divergence-free flow and
trace-zero stress tensor could be discretized for the test problem in Example
11.1. The discrete problem is still of the form Eqs. (11.3)-(11.4), but now Kd is
bounded. The dual property is that the objective function Dd(y) in Eq. (11.4)
is now finite for all y RM. However, Dd(y) is not differentiable (we shall give
the precise form shortly), so standard algorithms do not apply efficiently to the
problem in Eq. (11.4).
The finite element spaces for o'h and the stream function 'Th and their bases
are described in Example 9.1. The matrix elements in Eq. (8.16) are computed
element by element as usual. In order to illustrate the simple, but rather te-
SECION 14 Solution of the discrete problem 277
FIG. 13.4. Collapse flow for the problem in Fig. 13.3.
dious procedure, let us indicate how the integral over one element (square) is
computed in the innermost of ten nested loops:
Loop 1-2: Lower left vertex of the square is specified.
Loop 3-6: Node for oh (i.e., a pair of columns) is specified.
Loop 7-8: Node for 'h is specified.
Loop 9-10: Type of basis function (i.e., a row) is specified. Now the contri-
bution to two matrix elements is computed.
The loops 3-10 only consist of two steps each. In the following we assume
that A and b have been computed. The condition x c Kd is identical to Eq. (9.7).
14.1. Solution of the primal problem

The primal problem Eq. (11.3) can be solved by MINOS (MURTAGH and SAUN-
DERS [1982]) as in the previous section. The problem is formally the same. The
constraint matrix A and b are different, and Kd defined by Eq. (9.7) is bounded.
The collapse fields for the test problem with a = obtained with this approach
are visualized in Fig. 14.1. Since the velocity is not linear on the grid lines for
this discretization the deformed grid appears curved between nodes. The stresses
are indicated by drawing the vector (L, o-2 ) (see Eq. (9.3)) at each node, where
FIG. 13.5. Collapse fields for the Coulomb yield condition.
it is at the yield surface: 2o-+ o-2 = 1. We see that the plastified region is not
well localized. When comparing the results with other element combinations,
it must be taken into consideration that this discretization has four times as
many degrees of freedom per element for both stresses and velocity. Hence
these results for an N x N grid should be compared with earlier results for a
2N x 2N grid.
The values for the collapse multiplier Ah are seen in Table 14.1. With these
more smooth finite element functions the error is still O(h), but in addition
Eq. (10.2) appears to hold. This may justify the use of smooth elements, even if
it does not change the order of the error. For general geometry this will require
use of Argyris' triangle or Bell's triangle.
14.2. Solution of the dual problem

Our main motivation for trying finite elements for divergence-free flow was to
solve the reduced dual problem, formally identical to Eq. (11.4):
minDd(Y), (14.1)
bTy = 1.
FIG. 14.1. Collapse fields for divergence-free elements, a = and h = 1, obtained with MINOS.
In Eq. (14.1) Dd(y) is finite for all y E RM after the discretization in Example
9.1. As mentioned earlier Dd(y) is not differentiable, but a closer look reveals
a structure for which there may still be an efficient algorithm: According to
Eq. (8.19) we have
Dd(y) = max XT(ATy),
xcKd
where
N
xT(ATy)= EXn(ATy)n.
n=1
Recall that the vector (x") is just a linear ordering of the variables (), which
define the discrete stress tensor through Eq. (9.6). If we keep the indices (v, i)
for the components of x and ATy we have
x T (AT y) Y)
By insertion of the normalized yield condition
(41)+ ( <)2
1 for all nodesv
TABLE14.1. Results for L = 1, a = , and a = obtained with

MINOS for the element combination described in Section 14.
h- 1 a a= 2
A* order extr. A order extr.
3 1.0555 1.5309
6 0.9898 0.9242 1.4775 1.4242
9 0.9670 0.93 1.9214 1.4493 0.18 1.3929
12 0.9559 1.06 0.9224 1.4336 0.68 1.3862
15 0.9493 1.05 0.9229 1.4238 0.86 1.3846
h- 1 a=
Ah order extr.
4 1.2598
6 1.2193 1.1383
8 1.1973 0.77 1.1314
10 1.1839 0.94 1.1303
12 1.1750 1.02 1.1305
14 1.1687 1.03 1.1307
we get
(( (Ay)1) + ( ATY) ) 2
Dd(y)= E
VV V/
vety
Dd(y) is a sum of norms of two-dimensional vectors. Each term only involves

two rows of AT, i.e., two columns of A. These are precisely the two columns
associated with the node v and corresponding to the value of the two compo-
nents of oh at that node. Now let Av denote the M x 2 matrix, which consists of
these two columns, i.e., the two columns corresponding to the primal variables
0 and 2. Then we get the following expression for the objective function in
the dual problem (14.1):
Dd(y) =E (aTy) + (ATY) 2

= E IiATyI. (14.2)
The norm is the ordinary Euclidian norm in R2.

SEcoN 14 Solution of the discrete problem 281
FIG. 14.2. Collapse fields for a = and h = obtained by solving the minimization problem (14.3).
The problem (14.1) may now be written

min
bTy
=1
VC
IIAryl. (14.3)
vX'
This belongs to a class of problems for which a second-order method was devel-
oped in CALAMJAI and CONN [1980, 1982] and OVERTON [1983]. The sum of norms
in Eq. (14.3) is not differentiable, when one of the terms 11*ATyll vanishes. Typi-
cally, in particular in limit analysis, many of these terms are zero at the optimal
vector y, a severe obstruction for the efficiency of standard algorithms.
The problem is solved as follows: When a term IIATyII becomes small, it
is tentatively set to zero. If this does not increase the value of the objective
function, the term is kept equal to zero by adding the linear constraint ATy = 0,
and the term is deleted from the objective function. In the resulting subspace for
y the objective function is now smooth, until either another term gets small and
the above step is repeated, or a minimum on the present subspace for y is found.
Thus a finite sequence of smooth convex linearly constrained problems is solved
282 E. Christiansen CHArPTER IV
by a second-order method. When a minimum subject to a number of constraints

of the form A y = 0 is found, there are two possibilities: Either we have found
the minimum for Eq. (14.3), or one (or more) of the constraints ATy 0 must be
removed again. (Unfortunately a term put equal to zero during the algorithm
need not be zero in the optimal solution.) An optimality check determines
which of the two situations has occurred. The principle of the method has a
neat interpretation in limit analysis: The nodes v for which ATy = 0 constitute
the rigid region for the discrete problem. So the algorithm is trying to decrease
Dh(Uh) = Dd(y) by expanding the rigid region with the option to reduce it later,
if need be.
There is a complication in limit analysis, which is not seen in the problems
that originally motivated this method. The set of columns from the matrices Av
in the constraints is linearly dependent. This causes problems for the optimality
check and the elimination of linear constraints.
Computing the limit load by minimizing a sum of norms has been applied
to two problems in limit analysis with bounded yield condition in OVERTON
[1984] and GAUDRAT [1991]. The above mentioned rank deficiency problem was
handled on an ad hoc basis. The general method described in this Section was
applied to our reference test problem (to be published) using the algorithm in
OVERTON [1983], but results were obtained only for coarse grids (12 x 12), due
to the rank deficiency problem. Very recently ANDERSEN [1995] has designed
an efficient large scale algorithm for minimizing a sum of norms. With this
algorithm we obtained results for much finer grids, confirming the first-order
convergence observed in Table 14.1.
The collapse fields for a 90 x 90 grid are shown in Fig. 14.2. Recall that this
should be compared with earlier results for a 180 x 180 grid. With this method
it was also possible to solve the classical punch problem from HILL [1950, p. 254]
of punching a die into a medium occupying a half space. To our knowledge this
problem has not previously been solved using numerical methods.
15. Concluding remarks
If the test problem is indicative of the behaviour of solution methods, as we think

it is, limit analysis is ready to be used along with the standard elastic analysis.
Plane problems can certainly be solved, and three-dimensional problems can be
approached with confidence.
Surprisingly fine grids can be handled by the dual affine scaling algorithm
with the linearized yield condition, both with the constant-bilinear element
combination over rectangles and with linear-linear elements over triangles. Also
the convex problem can be solved with MINOS, using linear-linear elements, for
more coarse, but still adequate, grids. These methods are ready to be applied.
The approach with elements for divergence-free flow looks good, but must be
extended to more general element shapes.
For the tested algorithms, except the dual affine scaling algorithm and the
SECrION 15 Solution of the discrete problem 283
TABLE 15.1. Dimensions of discrete spaces for some natural element combinations.
Elements or
h - uh dim(Xh) dim(Uh) ratio Expe-
rience
1 constant-bilinear [ 3N 2 2(N + 1)2 3 good
2 bilinear-bilinear 5 3(N + 1)2 2(N + 1)2 3 good
3 linear-linear A 3(N + 1)2 2(N + 1)2 3 good
2
4 constant-linear Z 6N 2(N + 1)2 3 bad
5 constant-quadratic A 6N 2 2(2N + 1)2 3
6 biquadratic-bicubic C1 [ 2(2N + 1)2 4(N + 1)2 (2) good
7 bilinear-bicubic C1 [ 2(N + 1)2 4(N + 1)2 (1) bad
8 const.-quadr. "V = 0"2(2N +1)2 2(2N + 1)2 2
2
9 linear-quadr. V = 0" A 2(N +1)2 + 2(2N 1)2 - 4N 1
2 2
10 cubic-Argyris A 18N + 12N + 2 9N + 14N + 6 (2)
11 type 3'-Argyris A 14N 2 + 12N + 2 9N 2 + 14N + 6 (14)
12 quadratic-Bell's A 2(2N + 1)2 6N 2 + 12N+6 (4)
algorithms for Eq. (14.3), the factor that limits the problem size, and hence the
grid size, is the deterioration of accuracy near the optimal solution. That is our
excuse for not reporting use of time and memory for the algorithms. We refer
to the papers cited. It appears that the development of speed and memory size
of computers has made it possible to solve problems so large that floating point
accuracy becomes a bottleneck, at least for some applications.
Table 15.1 lists some of the more natural element combinations and the di-
mensions of the discrete spaces for the test problem for L = 1 and h = 1/N.
The boundary conditions on u have been ignored, so the true dimension of Uh
is O(N) smaller than indicated.
For the divergence-free Cl-elements in rows 6, 7, 10, 11, 12 the column
dim(Uh) gives the dimension for the discrete stream function T"h (see Example
9.1). Row 8 indicates constant elements for the stresses combined with quadratic
elements for the flow satisfying (see TEMAM [1977, Section 4.4])
V'U h = VT C 3-h.
Then Eq. (8.26) is satisfied for any uh E Uh. Row 9 indicates linear elements
for stresses and quadratic elements for velocity satisfying
J Vh =O VT h,
T
where & denotes any linear function on the triangle T. In particular Eq. (8.26)
is satisfied for any uh C Uh.
The element combinations 5, 7, 8 and 9 in the table cannot satisfy condi-
tion (8.24), since dim(Vh) < dim(Uh). Hence they should not work. This has
been confirmed for row 7. To our knowledge experience with the combinations
in rows 5, 8-11 has not been reported. Based on the dimensions the element
combinations 10 and 11 should work.
15.1. On the solution to the test problem

After spending so much effort finding approximate solutions to the test problem
it seems reasonable to offer some comments on it.
We expected the collapse flow to consist of three rigid regions sliding along
narrow zones, possibly slip-planes. Only near the boundary x2 = 0, where the
three regions meet, was a more complicated flow expected. The fine-grid solu-
tions with the dual affine scaling method suggest a more interesting flow with
deformation zones of finite width. The collapse flow may even be continuous.
In order to test the various solution methods with respect to uniqueness,
localization of plastified region etc., they were also applied to the case with
no cut at all. A simple solution is easily found by hand: Constant stresses and
deformation (not flow) and A*= 2, but the fields are highly non-unique. All
methods gave Ah = 2 with various collapse fields, most of them similar or equal
to the simple exact solution.
The problem with these solutions, both exact and discrete, is that they are
not observed in real life. There we do not see uniform collapse fields, but a
phenomenon called necking: The deformation is localized to a region where the
material shrinks in the x 2 -direction and stretches in the xl-direction. In order
to find such a solution numerically the following boundary condition was added
to the test problem with L = 1:
u 2 =0 forxl=l1. (15.1)
Restricting the flow u cannot possibly decrease the limit multiplier A*. In fact
it was unchanged, A*= 2, as expected. The new collapse flow field displayed
necking (it was forced to), but since A* was not changed the new flow field was
also a solution to the original problem without the extra boundary condition.
We conclude that the numerical method can find the "right" solution, but in the
absence of uniqueness the numerical solution method will usually not have the
same preferences as nature. There are more surprises in limit analysis than in
elasticity.
CHAPTER V
Limit Analysis for Plates
16. The continuous problem
Also in the plate model the problem of limit analysis can be formulated as
an infinite-dimensional mathematical programming problem. Formally the main
difference is that there is one more derivative in the bilinear form for the internal
work rate.
Consider a plate occupying the area Q2 in the xl - x2 plane. We only consider
transversal loads (0, 0, f), f = f(xl, x 2 ). Let u = U3 be the transversal displace-
ment rate. At a distance X3, measured with sign, from the mid-plane of the plate
we get (HODGE [1959]):
au
i =-x 3x,ai = 1,2,
and hence
eij = -X3 i, j = 1, 2
In the plate model the other components of the strain rate tensor are neglected.
This implies, in particular, that transversal shear is a violation of the model.
The internal work rate may be written
l 3=-H/2 1i,J
=- | mij
fOx
dxd2,
a(mio, = x3ij , ij = dx 3 1, 2
where
H/2
H/2
are the bending moments. H is the thickness of the plate and may depend on
285
286 E. Christiansen CHAPER V
(X1 ,X2 ). Hence the internal work rate is expressed in the bending moments
m = (mij)2il and the transversal displacement rate u by the bilinear form
2 '2 u
b(m, u) = E mii 0 dx 1 dx 2 . (16.1)
The yield condition on the stress tensor can immediately be translated to

the tensor m of bending moments. For example, the Von Mises condition for a
homogeneous plate is
ml21 - m 11m22 +22 + 3m12 M (16.2)
The Tresca yield condition is
max (Imi, Imil, mi - miIl) < Mo, (16.3)
where mi and mI, are the principal moments.
The work rate of the transversal load f(xl, x 2) may be written as a linear form
in u:
F(u) = fu dx1 dx 2. (16.4)
The equilibrium equation for the bending moments is the equation of virtual
work rate. The expressions in Eqs. (16.1) and (16.4) must be equal for all u
satisfying the boundary conditions:
b(m,u) = F(u) for all u. (16.5)
We shall consider only homogeneous boundary conditions, since the non-
homogeneous case is not particular to limit analysis. The boundary 012 consists
of finitely many sections, each of which is either supported (u = 0) or non-
supported, and either clamped (u/ov = 0) or non-clamped. v denotes the out-
ward normal to 12, and au/dv the derivative of u in that direction. This gives
four possible combinations for each section of the boundary.
Let So C n2 denote the supported part of the boundary:
u =O on So. (16.6)
Condition (16.6) is an essential boundary condition, to be imposed explicitly.
The corresponding naturalboundary condition on To = 02\So is, in the homo-
geneous case with no load or twisting moments applied at the boundary:
v-(V-m) = , (v.m.t) =0 on To = 012\So. (16.7)

X is the unit vector with angle +iT/2 to the normal . Being a natural bound-
ary condition, Eq. (16.7) is not imposed explicitly. Handling of inhomogeneous
natural boundary conditions is mentioned below.
Similarly, let S1 C 0a denote the clamped part of the boundary (supported
or not):
au on (168)
= 0 on S1. (16.8)
SECION 16 Limit analysis for plates 287
Equation (16.8) is a natural boundary condition and is not imposed explicitly.

The corresponding essential boundary condition is in the homogeneous case
m, = Pv(m-v) = 0 on T1 = l2\S . (16.9)
If m and u are sufficiently smooth, Green's formula can be used repeatedly
to rewrite the internal energy dissipation rate Eq. (16.1) as follows:
b(m,u) =-Emij dx 1 dx 2
=/(V.m).Vu dxld 2 - /(v.m).Vu ds

n an
= f(V.m).VU d d 2 -
n
J
ana
(vm) ( +T ) ds
=f
(V.m).Vu
=
V
dxl dx2 +
X2
J
f(m).
+ Tl
A-(vmT) u ds (16.10)
a Oa
= J(V.m).Vu dx1 dx2 (16.11)
=- JV(V.m) udx dx2 + f .(V.m) uds

Q an
-| V.(V.m) u dxl dx2. (16.12)
In addition to Green's formula we have applied Eqs. (16.8)-(16.9) in

Eq. (16.10) and Eqs. (16.6)-(16.7) in Eq. (16.11) and in Eq. (16.12).
Inserting Eq. (16.12) into Eq. (16.5) gives the classical equilibrium equation
for plates
02mll 0+ 2 02m22
V.(V.m) = +2 + a = -f(16.13)
ax1 0x10X2 ax2
Conversely, if we impose the essential boundary condition Eq. (16.9) on m,
and require that Eq. (16.5) holds for all u satisfying Eq. (16.6), then for suffi-
ciently smooth m and u we conclude that Eqs. (16.7) and (16.8) hold in addition
to Eq. (16.13). We also see that inhomogeneous natural boundary conditions
on the moments, i.e., nonzero boundary loads, are handled by including these
loads in the expression (16.4). Inhomogeneous boundary conditions on u are
irrelevant for a rigid plastic plate.
From a mathematical and numerical point of view it is preferable to use the
most symmetric form Eq. (16.11) with one derivative on both m and u. Then
288 E. Christiansen CHAFFER V
the essential boundary conditions are meaningful, and standard finite elements
can be used.
The problem of limit analysis for plastic plates may now be formulated in
complete analogy with Section 5. Let K denote the set of bending moment
tensors m, which satisfy the yield condition and the boundary condition (16.9),
and let U be the space of transversal displacement rates u, which satisfy the
boundary condition (16.6). Then the static principle of limit analysis takes the
form
A*= max{A I 3m e K: b(m, u) = AF(u) Vu E U}
=max min b(m,u), (16.14)
mEK F(u)=l
where F(u) is defined by Eq. (16.4), and b(m, u) by Eq. (16.11).

The kinematic principle can be written
A*= min max b(m,u)
F(u)=l mcK
= min D(u), (16.15)

F(u)=l
where
D(u) = sup b(m, u). (16.16)
mEK
In the next section we prove that, under certain conditions, the duality be-
tween the static and kinematic formulations holds, and that collapse solutions
for m and u exist.
17. Duality of limit analysis
We shall formulate a mathematical frame for the duality problem for plastic
plates. The result is not as general as one could wish, since it requires the entire
boundary of the plate to be supported. We believe this restriction to be mainly
of proof technical nature.
Assume that the plate is supported at the entire boundary (and only there)
u = O ona
A, (17.1)
i.e., So = d01 in Eq. (16.6). It may be clamped along part of or all of its boundary,
as indicated by Eqs. (16.8) and (16.9). Where it is not clamped, i.e., it is simply
supported, we impose the essential boundary condition (16.9).
Let p > 2, such that W 1'P (2) is continuously imbedded in L (2) (in fact in
C(2)), and let lip + 1/q = 1. The space X of bending moments consists of all
symmetric 2 x 2 tensors m = (mij) with the following properties
mij E W 1 'P(n), (17.2)
m, = vP(m.v) = O on a2\S1, (17.3)
where S1 is the clamped part of the boundary. X is equipped with the WlP-norm.
SECTION 17 Limit analysis for plates 289
The space Y of transversal displacement rates consists of those u E W 1'q(2)

for which there exists a symmetric 2 x 2 tensor of measures = (ij), [zij E
C*(2), such that
(ijmij)x = (V(Vm), u)pxW1,q Vm E X (17.4)

i,j=l
Y is equipped with the Wl'q-norm. The measures [ij are the generalized second-
order derivatives of u. This generalization is necessary in order to incorporate
the energy dissipation rate from a hinge (i.e., a discontinuity in the derivatives
of u) along the clamped part of the boundary. The natural boundary condition
(16.8) is incorporated through Eq. (17.4).
On X x Y the bilinear form b(m, u) is defined by the expression in Eq. (17.4).
The transversal load f must satisfy
f G W- 1'P(2) = (W1(2))
*
(17.5)
and the external work rate is given by the linear form
F(u) = (f,U)w,,XW,q Vu E Y. (17.6)
THEOREM 17.1. Assume that f2 is bounded with Lipschitz continuous boundary.

For any f E W- P(12) there exists m E X, such that
V.(V.m)=f in 2.
PROOF. We may find a solution 0 E W1'P(12) to the Dirichlet problem

A =f in2,
0 =0 on 02.
Let ml = m2 2 = 0 m 1 2 = 0. Then
V.(V.m) = A =f in 2
and
v.(m.v) = Orv =0 on 01. []
The yield condition is assumed to be of the form

m(x) E B(x) x E 2,
B(x) closed, convex and bounded uniformly in x E 2, (17.7)
38>O: mijl <8 i,j mEB(x) xE 12.
Let K be the set of admissible moment tensors:
K={m E X m(x)
I c B(x) Vx E }. (17.8)
THEOREM 17.2 (Duality theorem for plates). Assume that 2 is bounded with
290 E. Christiansen CHAPrER V
Lipschitz continuous boundary. Let X, Y, b(m, u), F(u) and K be defined as

above. Then the duality between Eq. (16.14) and Eq. (16.15) holds. More pre-
cisely:
sup min b(m, u) = min sup b(m, u). (17.9)
mEK F(u)=1 F(u)=1 mCK
The supremum on the left-hand side of Eq. (17.9) is attainedin the following weak
sense: There exist m* C [L- (2)]3, which satisfies m*(x) E B(x) almost everywhere
in 12 and
V.(V.m*) E W-1'P(l),
b(m*,u) = A*F(u) Vu E Y,
where A* is the value of Eq. (17.9). If u* minimizes the right-hand side o
Eq. (17.9), we have the following saddle point inequality for (m*, u*):
b(m, u*) < A* = b(m*, u*) = b(m*, u) (17.10)
for all m E K and all u E Y with F(u) = 1.
PROOF. The proof is similar to the proof of Theorem 5.7, only considerably
simpler because of the boundary condition u = 0 on M02 and the bounded yield
set. We apply Theorem 2.1 in CHRISTIANSEN [1980a] to the bilinear form -b(m, u)
on X x Y and the convex sets K C X and C C Y given by
C = {u E Y IF(u)= 1).
C is closed in Y, and due to the third condition in Eq. (17.7), and the fact that
p > 2, m = 0 is an interior point of K. Equation (17.9) will follow by establishing
the following two implications (a) and (b):
(a) If W is a continuous linear functional on Y, which satisfies infauc 't(u) >
-0o, then there exists m E X, such that
q(u) = b(m, u) Vu E Y.
(b) If <I is a continuous linear functional on X, which satisfies supmeK (m) <
oo and the condition
b(m, u) = 0 u E Y => (m) = 0,
then there exists u E Y, such that
0q(m)=b(m,u) Vm X.
If t in (a) is bounded from below on the affine hyperplane C, then is
proportional to F, i.e., there exists a real number t, such that
q(u) = tF(u) Vu E Y.
From Theorem 17.1 we conclude that there is m E X, such that
F(u) = b(m,u) Vu E Y.
Hence (a) holds.
Now let 0 satisfy the conditions in (b). Due to Eq. (17.7) K contains a ball for
the maximum-norm, and since (P by assumption is bounded on this ball there
exist measures tii E C*0(), such that /12 = 21 and
0(m)= Z(ij,mij)c*xc Vm E X. (17.11)

i,j
The second condition on implies that the null space of the map
T: m - V.(V.m)
is contained in the null space of 4. Since by Theorem 17.1 T maps X onto
W- 1 ,P(2)and thus is an open map, we conclude that may be factorized over
the null space of T:
(m) = (, T(m))wq xw ,, Vm X
for some u E (W-1 P(12))* = Wq(Q2). Comparing with Eq. (17.11) we see that
u E Y. Thus (b) holds, and Eq. (17.9) is proved.
Let (mk) be a maximizing sequence for the left-hand side of Eq. (17.9), i.e., a
sequence in K, which satisfies
b(mk, u) = kF(u) Vu E Y, Ak -- A*,
where A*is the value of Eq. (17.9). By weak* compactness of K in Lo there
exists m*, which due to the convexity of B(x) satisfies m*(x) E B(x) almost
everywhere in 2. Since m* is also a distributional limit we have
V'(V.m*) = *f E W-11P(f2).
This completes the proof. [
18. Discretization of the plate problem
The discretization of the plate problem is formally the same as in Section 8.

However, the bilinear form b(m, u) and the essential boundary conditions on
both m and u strongly suggest that continuous finite element functions are used.
An obvious choice is piecewise linear elements for both m and u, or in simple
geometries piecewise bilinear elements. The bilinear form b(mh, Uh) is computed
by the expression (16.11) with first order derivatives on both m and u.
b(m, u) = (V.m).Vu dx1 dx2
/l ftmll Au aml
2 Au anl2 au &m
22 0u \
+ dX
J ax o1 Xl x- a 1
O X2x Ox2 d 1 .2
(18.1)
Let 3h be a triangulation of 2 in the sense of CIARLET [1991, Chapter II].
Based on this triangulation we choose finite element spaces (e.g., piecewise
292 E. Christiansen CHAPTER V
linear element functions) for u and each component of m. Let Xh and Yh de-
note the resulting spaces for the discrete bending moments mh and transver-
sal displacement rate Uh, respectively. With continuous finite element functions
Xh C X and Yh C Y will always be satisfied (CIARLET [1991, Theorem 5.1]),
where X and Y are the spaces defined in Section 17. Problems concerning
curved boundary will not be discussed, since they are not specific to limit anal-
ysis. We refer to CIARLET [1991, Chapter VI].
The yield condition (17.7) on mh will be imposed only at the nodes: If {x }I
are the nodal points for the value mh(x), then the discrete set of admissible
moment tensors Kh is defined as
Kh = {mh X h I mh(X) B(xv) Vxv}. (18.2)
As pointed out in Section 8, Kh C K may not be satisfied for higher-order ele-
ments.
As a matter of convenience we assume that the work rate for the transversal
load f given by Eq. (16.4), or formally by Eq. (17.6), can be computed exactly
for Uh Yh. The effect of using numerical integration is not specific to limit
analysis, and we refer to CIARLET [1991, Chapter IV].
The discrete problem for plates may now be written
A = max{A 13mh c Kh: b(mh, Uh) = AF(Uh) VUh G Yh} (18.3)
= max min b(mh, Uh)
mhEKh F(uh)=l
= min max b(mh, h)

F(uh)=l mhcKh
= min Dh(Uh), (18.4)

F(uh)=l
where
Dh(Uh) max b(mh, h). (18.5)
mh Kh
As in Section 8 we shall indicate how the discrete problem may be written in

a form suitable for mathematical programming. Let rm and rXu be the node sets
for the finite element representation of mh and Uh, respectively. With first-order
derivatives on both m and u in the bilinear form b(m, u) these node sets will
typically be the same, except for the boundary conditions Eqs. (16.6) and (16.9).
Let (o be the scalar basis function associated with the node v cE rm. Then there
are the following three basis functions for Xh associated with the node v:
[ 01] 2 [ O2 i [ V] (18.6)
Every discrete moment tensor mh G Xh may be written as a linear combina-

tion of the basis functions:
mh = E (1k-1 + +2 2) (18.7)
ve.
SECION 18 Limit analysis for plates 293
The coefficient vij is the value of mh,ij corresponding to the nodal value v.
For each node c authere is just one basis function C for Yh, since Uh is
a scalar function. Hence each Uh c Yh is of the form
Uh = E .1 ,L (18.8)
Equations (18.7) and (18.8) must be modified according to the boundary con-
ditions.
The yield condition (18.2) can now be written
[2 e22] B(x,) v E m, (18.9)
while the energy functions become

F(uh) = E 1rF(q), (18.10)
u.
b(mh, Uh)= E E ( bl/b((qll /) + 2.b(1 2

, 'q/p)
VCJ'm p.CA;
+ ,22Lb(v, ')). (18.11)
After insertion of Eq. (18.1) the single terms in Eq. (18.11) become
b( 1 , )= / t dx dx 2 , (18.12)
Q
b(ov A)t ,/ ) dX (18.14)
12
These are the elements of the matrix corresponding to b(mh, Uh). As usual they
should be assembled from the contributions from each element separately.
As in Section 8 we shall assume that the coefficients 6f~ and q,. have been
ordered in a linear numbering:
n = n(v, i, j) e {1, 2,..., N}, (18.15)
m = m(u) E 1, 2,..., M}. (18.16)
With such an ordering the coefficients (/) and (h.) may be identified with
vectors x = (xn) E RN and y = (Ym) RM:
xn = f where n = n(v, i, j), (18.17)
Ym = M, where m = m(.u). (18.18)
The work rate functions Eqs. (18.10) and (18.11) may now be written
M
F(Uh)= EYmbm = bTy, (18.19)
m=l
where
bm = F(q), m = m(F)
and
M N
b(mh, Uh) = Ey m Xn amn yTAx = xT (ATy) (18.20)

m=l n=l
A is the M x N matrix
amn = b(vj, l,'), m(A.),
m n = n(v, i,j).
In the variables x and y the discrete problem (18.3)-(18.4) for plate bending
is formally identical to the general discrete problem (8.17)-(8.18). For later
reference we shall rewrite it here. Let
Kd = {x E RN I mh Kh}, (18.21)
where x and mh are related through Eqs. (18.7) and (18.17), and Kh is defined
in Eq. (18.2). Then
Ah =max{A I 3x E Kd:Ax = Ab} (18.22)
= max min yTAx
xEKd bTy=l
= min min xT(ATy)

bTy=l xcKd
= min Dd(y), (18.23)

bTy=l
where
Dd(y) = max xT(ATy). (18.24)
xCKd
The discrete duality theorem, Theorem 8.1, applies, so we may conclude that
the duality between Eqs. (18.22) and (18.23), or equivalently between Eqs. (18.3)
and (18.4) holds, and that there exists a saddle point (m*, u*):
b(mh, uh ) < Ah = b(m, uh ) = b(m, Uh) (18.25)

for all mh E Kh and Uh E Yh with F(Uh) = 1.
In complete analogy with Theorem 10.1 we have:
THEOREM 18.1. Let (*, m*, u*) denote an exact solution, i.e., a triple satisfying
Eq. (17.10) and let (h, mh, uh) be a discrete solution satisfying Eq. (18.25). If mh
satisfies the exact yield condition m GCK, then
SECTION 19 Limit analysisfor plates 295
b(mh - m*, uh) < A - A* < b(m', Uh - u*) (18.26)

for all mh E Kh and all Uh E Yh with F(uh) = 1.
PROOF. In Eq. (17.10) insert m = m and u = u* and subtract from Eq. (18.25).
Equation (18.26) is formally identical to Eq. (10.4). However, the bilinear

form b(m, u) is of higher order than a(o-, u) in Eq. (10.4), so a greater variety
of norms may be used, and a convergence order greater than 1 is quite possible.
Following CHRISTIANSEN [1980b, Section 4] we shall give the following general
estimates, which are direct consequences of Eq. (18.26).
THEOREM 18.2. Assume that II11 and 11'112 are norms on X and Y respectively,
such that
jb(m, u)l < Ilmli 1lull2 V(m, u) CX x Y,
and let (A, m;, uh) be a sequence of discrete solutions with m* E K.
(a) If |Imn1l
h1 Cl, then
A - A* <C1 F(uh)=1
min IUh-u*112.
(b) If UhII1
2 < C2 , then
A - * > -C 2 min Ilmh - m*ll
mh Kh
19. Solution of the discrete problem
We shall refer to the discrete static formulation (18.22) as the primal problem
maxA,
(P): Ax= Ab, (19.1)
XEKd
and to the discrete kinematic formulation (18.23) as the dual problem
(D): minDd(Y) (19.2)

by = 1.
These problems are formally identical to Eqs. (11.3) and (11.4). However, the
differences in the linear constraint matrix A and the yield set Kd, which is now
bounded, are significant. In our experience the plate bending problem is consid-
erably easier to handle numerically. Already HODGE and BELYTSCHKO [1968] and
RANAWEERA and LECKIE [1970] used convex programming methods to compute
limit loads for plates (bound methods), and in CHRISTIANSEN and LARSEN [1983]
the conclusion is that the linear programming approach is not competitive with
nonlinear methods. Since then the power of the software available for nonlin-
ear optimization has increased significantly. For example MINOS (MURTAGH
and SAUNDERS [1982]) can be applied immediately to the form Eq. (19.1) for
general plate bending problems. This approach is ready for application.
Alternatively, consider the dual problem (19.2). The objective function Dd(y)
is given by Eq. (18.24):
Dd(y) = max xT(ATy). (19.3)
xcKd
Recall that x = (x,) is a numbering of the physically more meaningful variables

('J) defined by Eqs. (18.7) and (18.17). In these variables we may write
xT (AT y) E (( (A Y)V+ + 2 TAy) +2 (ATY)v2) (19.4)
In this expression
(ATy)j= E b(obi, )71K, (19.5)
where the relation between y = (Ym) and (, ) is given by Eq. (18.8).

In order to be explicit, consider the normalized Von Mises yield condition:
m2 - m1 1m 2 2 + m22 + 3m 2 X 1. (19.6)
Expressed in the variables (Y) the condition is
(1)2 _ /1g22 + (22)2 + 3(,12)1 < 1 for all nodes v. (19.7)
Thus Dd(y) in Eq. (19.3) may be computed by maximizing the sum in Eq. (19.4)
subject to the constraints (19.7). Since the constraints do not couple different
nodes v, each term in the sum (19.4) is maximized. For convenience let z=
ATy E I N , or in the notation of Eq. (19.5) z ii = (ATy)ij. Then each term in
Eq. (19.4) has the maximal value
max (nl z ll + 22Z2 + 12zv2),
where e = (E11, 522, 512), and Q is the symmetric positive definite matrix
1 -124o.
QJI0'
2 - 1 22 4
0 3 0 0 1
If we ignore boundary conditions on the moments, Dd(y) is explicitly given by
the expression
Dd (y)-- E ? TQ
p
I-Z
1 /4 ((Z1)2 + ZZ22 + (22)2) + (Z12)2 (19.8)
=XE IIlCzvlI (19.9)
where C is the Cholesky factor of 3Q -1:
C = X O] . (19.10)
0 1
Now let A, be the M x 3-matrix consisting of the three columns from A

associated with the node v, i.e., (a,)ij = a,i,j for all nodes v and (i,j) =
(1, 1), (2, 2), (1, 2). Then z, = ATy for all nodes v, and Eq. (19.9) may be written
Dd( ) = l CAVyl. (19.11)
Hence Dd(y) is the sum of Euclidian norms of the vectors (CA T)y. The dual
problem (19.2) is to minimize this sum subject to the single linear constraint
bTy = 1.
With the dual problem written in the form
min 1 (CA ) Y (19.12)

y
bTy=l 3 A
it can be solved with the algorithm in OVERTON [1983].
CHRISTIANSEN and LARSEN [1983] solved the dual problem (19.2) with the ob-
jective function given by the expression (19.8) using the algorithm GOLDFARB
[1969]. This is an adaptation for linear constraints of the variable metric method
and assumes differentiability of the object function, a requirement not met by
the problem (19.2). For the test problem reported here, the differentiability
problem was taken care of on an ad hoc basis, and the method can not be
expected to work in general. The two more recent methods discussed above
should be used.
In the following we report some of the results from CHRISTIANSEN and LARSEN
[1983]. For a full discussion we refer to that paper. The results have also been
obtained using the two more general methods described above: MINOS applied
to the form (19.1) and Overton's algorithm applied to the form (19.12). For the
case of a uniform load it is easy to obtain good accuracy, and the convergence
order of Ah as h tends to zero appears to be greater than 1. Figure 19.1 shows
the collapse velocity for a simply supported square plate with a uniform load.
The transversal displacement rate u* has been multiplied by a suitable time
scale, in order to visualize the deformation. The tensor of bending moments is
at the yield surface at every node, i.e., the whole plate is plastified.
FIG. 19.1. Collapse deformation for the simply supported square plate with uniform load.
TABLE 19.1. A for square (L = 1) and rectangular (L = 2) plates with uniform

load.
h- 1 L=1 L=2
Simply supported Clamped Simply supported Clamped
8 24.6574 42.0719 14.8363 26.1332
12 24.8634 43.0202 14.8924 26.4403
16 24.9280 43.4108 14.9116 26.5689
20 24.9590 43.6139 14.9204 26.6367
24 24.9765 43.7358 14.9253 26.6775
28 24.9873 43.8157 14.9283 26.7042
32 24.9945 43.8716
36 24.9995 43.9124
40 25.0031 43.9433
The concept of a point load is in violation of the plate model, since shear
in the z-direction cannot be neglected. Technically speaking, a point load does
not belong to the space W-l P(12) of admissible loads for plastic plates. From a
physical point of view a point load is the limit of loads, which are concentrated
on a small finite area, which shrinks into a point, while keeping the total load
fixed. The limit is independent of the actual sequence of loads. The sequence of
SECTION 19 Limit analysisfor plates 299
FIG. 19.2. Collapse deformation for the simply supported square plate with a point load at the center.
9 b
0
C
a 0
d
8
0 a
U-
7 as"@
@OO
a.
e
A* A ,** *A_ A I
A
6 f
*a
05 1 15 h
FIG. 19.3. A for the simply supported plate with various discretizations of a point load: (a) Discrete
point load. (b) Uniform load concentrated at a central square of side 2h, i.e., over the four elements
adjacent to the point. (c) Uniform load concentrated at a circle of radius h. (d) Uniform load
concentrated at a square of side hv tilted 45 degrees. (e) Uniform load concentrated at a square
of side h. (f) Same as (a), but for the rectangular plate of length L = 2.
13
[0 0 I(c)
F3
12
0 a (d)
11
no
Do
O
10 0 0 ____
0U
7 } * *
* *** * *
(a)
l l l l b
.05 .1 .15 h
Fig. 19.4. Same as Fig. 19.3, (a)-(d), for the clamped plate. Values for the rectangular case
coincides with (a).
collapse fields need not converge, and Theorem 17.2 does not apply. Even with
these limitations the concept of a point load is considered useful for applications.
There are several ways to discretize a point load. One of them is a discrete
point load at one node, but keeping in mind that a point load is a limit case, any
sequence of discrete loads, which shrink with the grid size, will give the same
limit of collapse multipliers. This can be used to obtain better approximations
than with a discrete point load.
Figure 19.3 shows the discrete collapse multipliers for five discretizations of a
central point load applied to the square and rectangular simply supported plate.
Clearly the choice of discretization is important. The limit as h -- 0 appears to
be the same, slightly less than seven, for the square and the rectangular plate.
SEcnON 19 Limit analysisfor plates 301
The marks L and U on the A-axis indicate the lower and upper bound for Ah
given in CASCIARO and DICARLO [1974].
The corresponding values for the clamped cases are shown in Fig. 19.4. The
values for the rectangular plate (L = 2) with a discrete point load are identical
to those for the square plate ((a) in Fig. 19.4). The results strongly indicate that
the exact limit multiplier A*is the same in the four cases considered: Clamped
and simply supported, square and rectangular plate.
Acknowledgement
I wish to give credit to the following students at Odense University, who have
contributed to the computational experience reported here: Knud Dalgaard
Andersen, Michael Bj0rnskov, Allan Porse Kristiansen, Christian Nissen, Erling
Kristian Purkaer.
References
ADLER, I., N. KARMARKAR, M.G.C. RESENDE and G. VEIGA (1989), An implementation of Kar-
markar's algorithm for linear programming, Math. Programming44, 297-335.
ANDERHEGGEN, E. and H. KN6PFEL (1972), Finite element limit analysis using linear programming,
Internat. J. Solids and Structures 8, 1413-1431.
ANDERSEN, K.D. (1993), An infeasible dual affine scaling method for linear programming, COAL
Bulletin 22, 20-28.
ANDERSEN, K.D. (1995), An efficient Newton barrier method for minimizing a sum of Euclidean
norms, SIAM J. Optim., to appear.
ANDERSEN, K.D. and E. CHRISTIANSEN (1995), Limit analysis with the dual affine scaling algorithm,
J. Comput. Appl. Math., to appear.
ARGYRIS, J.H. (1967), Continua and discontinua, in: Proceedings Conference on Matrix Methods in
Structural Mechanics, Wright-Patterson Air Force Base, OH 1965.
BARNES, E.R. (1986), A variation on Karmarkar's algorithm for solving linear programming prob-
lems, Math. Programming36, 174-182.
CALAMAI, P.H. and A.R. CONN (1980), A stable algorithm for solving the multifacility location
problem involving Euclidean distances, SIAM J. Sci. and Statist. Comput. 1, 512-526.
CALAMA, P.H. and A.R. CONN (1982), A second order method for solving the continuous multi-
facility location problem, in: G.A. WATSON, ed., Numerical Analysis: Proceedings Ninth Biennial
Conference, Dundee, Scotland, Lecture Notes in Mathematics 912 (Springer, Berlin), 1-25.
CAPURSO, M. (1971), Limit analysis of continuous media with piecewise linear yield condition, Mec-
canica 6, 53-58.
CASCIARo, R. and A. DICARLO (1974), Mixed FE. models in limit analysis, in: J.T. ODEN, ed., Com-
putational Methods in Nonlinear Mechanics (Texas Inst. for Computational Mechanics, Austin,
TX), 171-181.
CHARNES, A. and H.J. GREENBERG (1951), Plastic collapse and linear programming, Bull. Amer.
Math. Soc. 57, 480.
CHARNES, A., C.E. LEMKE and O.C. ZIENKIEWICZ (1959), Virtual work, linear programming and
plastic limit analysis, Proc. Roy. Soc. London Ser. A 251, 110-116.
CHRISTIANSEN, E. (1976), Limit analysis in plasticity; theory and approximation by finite elements,
PhD. Thesis, MIT, Cambridge, MA.
CHRISTIANSEN, E. (1980a), Limit analysis in plasticity as a mathematical programming problem,
Calcolo 17, 41-65.
CHRISTIANSEN, E. (1980b), Limit analysis for plastic plates, SIAM J. Math. Anal. 11, 514-522.
CHRISTIANSEN, E. (1981), Computation of limit loads, Internat J. Numer. Meth. Engrg. 17, 1547-
1570.
CHRISTIANSEN, E. (1982), Examples of collapse solutions in limit analysis, Utilitas Math. 22, 77-102.
CHRISTIANSEN, E. (1986), On the collapse solution in limit analysis, Arch. Rational Mech. Anal. 91,
119-135.
CHRISTIANSEN, E. (1991), Limit analysis with unbounded convex yield condition, in: R. VICHNEVETSKY
and J.J.H. MILLER, eds. ProceedingsIMACS 13th World Congress on Computational andApplied
Mathematics (Criterion Press, Dublin), 129-130.
CHRISTIANSEN, E. and K.O. KORTANEK (1990), Computing material collapse displacement fields on
a CRAY X-MP/48 by the Primal Affine Scaling Algorithm, Ann. Oper. Res. 22, 355-376.
303
304 E. Christiansen
CHRISTIANSEN, E. and K.O. KORTANEK (1991), Computation of the collapse state in limit analysis
using the LP Primal Affine Scaling Algorithm, J. Comput. Appl. Math. 34, 47-63.
CHRISTIANSEN, E. and S. LARSEN (1983), Computations in limit analysis for plastic plates, Internat.
J. Numer. Meth. Engrg. 19, 169-184.
CHRISTIANSEN, E. and H.G. PETERSEN (1989), Estimation of convergence orders in repeated Richard-
son extrapolation, BIT 29, 48-59.
CIARLET, PG. (1978), The Finite Element Method for Elliptic Problems (North-Holland, Amster-
dam).
CIARLET, P.G. (1991), Basic error estimates for elliptic problems, in: P.G. CIARLET and J.L. LIONS,
eds., Handbook of Numerical Analysis, Vol. II (North-Holland, Amsterdam), 17-351.
COHN, M.Z. (1979), Introduction to engineering plasticity by mathematical programming, in: M.Z.
COHN and G. MAIER, eds., Proceedings 1977 NATO Advanced Study Institute, University of Wa-
terloo (Pergamon Press, New York).
COULOMB, C.A. (1773), Essai sur une application des rgles de maximis et minimis quelques
problbmes de statique relatifs a l'architecture, Mim. Math. Phys. Acad. Roy. Sci. 7, 343-382.
DEMENGEL, F (1989), Compactness theorems for spaces of functions with bounded derivatives and
applications to limit analysis problems in plasticity, Arch. Rational Mech. Anal. 105, 123-161.
DIKIN, I.I. (1967), Iterative solution of problems of linear and quadratic programming, Soviet Math.
Dokl. 8, 674-675.
DRUCKER, D.C, W. PRAGER and H.J. GREENBERG (1952), Extended limit design theorems for con-
tinuous media, Quart. Appl. Math. 9, 381-389.
DAUTrr, G. and J.L. LIONS (1972), Les Inequations en M6canique et en Physique (Dunod, Paris).
[English translation: Inequalities in Mechanics and Physics (Springer, Berlin, 1976)].
FOIAS, C. and R. TEMAM (1978), Remarques sur les equations de Navier-Stokes stationnaires et les
ph6nomenes successifs de bifurcations, Ann. Scuola Norm. Sup. Pisa C. Sci. (4) 5, 29-63.
GAUDRAT, V.E (1991), A Newton type algorithm for plastic limit analysis, Comput. Methods. Appl.
Mech. Engrng. 88, 207-224.
GOLDFARB, D. (1969), Extension of Davidon's variable metric method to maximization under linear
inequality and equality constraints, SIAM J. Appl. Math. 17, 739-764.
GREENBERG, H.J. and W. PRAGER (1952), On limit design of beams and frames, Trans. Amer. Soc.
Civil Engrg. 117, 447.
GvozDEv, A.A. (1938), The determination of the value of the collapse load for statically inde-
terminate systems undergoing plastic deformation, in: Proceedings of the Conference on Plastic
Deformations Dec. 1936, Akademiia Nauk SSSR, Moscow-Leningrad. [Internat J. Mech. Sci. 1,
322 (1960)].
HAYES, D.G. and P.V. MARCAL (1964), Determination of upper bounds for problems in plane stress
using finite element techniques, Internat. J. Mech. Sci. 9, 245-251.
HILL, R. (1950), The Mathematical Theory of Plasticity (Oxford University Press, Oxford).
HODGE, P.G. (1959), Plastic Analysis of Structures (McGraw-Hill, New York).
HODGE, P.G. and T. BELYTSCHKO (1968), Numerical methods for the limit analysis of plates, J. Appl.
Mech. 35, 796-802.
KACHANOV, L.M. (1971), Foundations of the Theory of Plasticity (American Elsevier, New York).
KALIszK, S. (1975), Keplekenysegtan. Elmelet s mernoki alkalmazdsok (Akad6miai Kiad6, Bu-
dapest). [English translation: Plasticity. Theory and Engineering Applications (Elsevier, Amster-
dam, 1989)].
KARMARKAR, N.K. (1984), A new polynomial-time algorithm for linear programming, Combinatorica
4, 373-395.
KAzrNczY, G. (1914), Bemessung von statisch unbestimmten Konstruktionen unter Bericksichtigung
der bleibenden Formanderungen, Betonszemle 4, 5, 6.
KIST, N.C. (1917), Leidt een sterkteberekening, die uitgaat van de evenredigheid van kracht en
vormverandering, tot een goede constructie van ijzeren bruggen en gebouwen? Inaugural Disser-
tation, Polytechnic Institute, Delft, Holland.
KoHN, R. and R. TEMAM (1983), Dual spaces of stresses and strains, with applications to Hencky
plasticity, Appl. Math. Optim. 10, 1-35.
References 305
KORTANEK, K.O. and S. MIAOGEN (1987), Convergence results and numerical experiments on a linear
programming hybrid algorithm, European J. Oper. Res. 32, 47-61.
MARTIN, J.B. (1975), Plasticity (MIT Press, Cambridge, MA).
MATTHIES, H., G. STRANG and E. CHRISTIANSEN (1979), The saddle point of a differential program, in:
R. GLOWINSKI, E. RODIN and O.C. ZIENKIEWICZ, eds., Energy Methods in Finite Element Analysis
(Wiley, New York).
MOREAU, J.J. (1976), Application of convex analysis to the treatment of elasto-plastic systems, in:
P. GERMAIN and B. NAYROLES, eds., Applications of Methods of FunctionalAnalysis to Problems
of Mechanics, Lecture Notes in Mathematics 503 (Springer, Berlin), 56-89.
MURTAGH, B.A. and M.A. SAUNDERS (1982), A projected Lagrangian algorithm and its implemen-
tation for sparse nonlinear constraints, in: A.G. BUCKLEY and J.-L. GOFIN, eds., Algorithms for
Constrained Minimization of Smooth Nonlinear Functions, Math. ProgrammingStudy 16 (North-
Holland, Amsterdam), 84-117.
NAYROLES, B. (1970), Essai de th6orie fonctionnelle des structures rigides plastiques parfaites, J.
Mdcan. 9, 491-506.
NEAL, B.G. (1956), The Plastic Methods of Structural Analysis (Chapman and Hall, London).
NE6AS, J. (1966), Sur les normes quivalentes dans Wpk(1) et sur la coercivit6 des formes formelle-
ment positives, in: Equations aux Derivtes Partielles (Presses de 'Universite de Montreal,
Montr6al).
NECAS, J. (1967), Les Methodes Directes en Theorie des Equations Elliptiques (Masson, Paris).
NECAS, J. and I. HLAVACEK (1981), Mathematical Theory of Elastic and Elasto-PlasticBodies: An
Introduction (Elsevier, Amsterdam).
OVERTON, M.L. (1983), A quadratically convergent method for minimizing a sum of Euclidean
norms, Math. Programming27, 34-63.
OVERTON, M.L. (1984), Numerical solution of a model problem from collapse load analysis, in: R.
GLOWINSKI and J.L. LIONS eds., INRIA Conf Comput. Methods Engrg. Appl. Sci. (North-Holland,
Amsterdam), 421-437.
RANAWEERA, M.P. and EA. LECKIE (1970), Bound methods in limit analysis, in: H. TOTTENHAM and
C. BREBBIA, eds., Finite Element Techniques in Structural Mechanics (Stress Analysis Publishers,
Southampton).
ROBINSON, S.M. (1972), A quadratically convergent algorithm for general nonlinear programming
problems, Math. Programming3, 145-156.
SAVE, M. (1979), Fundamentals of rigid-plastic analysis and design, in: M.Z. COHN and G. MAIER,
eds., Proceedings 1977 NATO Advanced Study Institute, University of Waterloo (Pergamon Press,
New York).
SHOEMAKER, E.M. (1979), On nonexistence of collapse solutions in rigid-perfect plasticity, Utilitas
Math. 16, 3-13.
STRANG, G. (1979a), A family of model problems in plasticity, in: R. GLOWINSKI and J.L. LIONS, eds.,
Computing Methods in Applied Sciences and Engineering, Lecture Notes in Computer Sciences
704 (Springer, Berlin), 292-305.
STRANG, G. (1979b), A minimax problem in plasticity theory, in: M.Z. NASHED ed., Functional
Analysis Methods in Numerical Analysis, Lecture Notes in Mathematics 701 (Springer, Berlin),
319-333.
STRANG, G. and G.J. Fix (1973), An Analysis of the Finite Element Method, (Prentice Hall, Engle-
wood Cliffs, NJ).
STRANG, G. and R. TEMAM (1982), A problem in capillarity and plasticity, Math. ProgrammingStud.
17, 91-102.
SUOUET, P.M. (1978), Sur un nouveau cadre fonctionnel pour les quations de la plasticity, C.R.
Acad. Sci. Paris Ser. A 286, 1129-1132.
SUQUET, P.M. (1979), Un espace fonctionnel pour les quations de la plasticity, Ann. Fac. Sci.
Toulouse 1, 77-87.
SUQUET, P.M. (1980), Existence and regularity of solutions for plasticity problems, in: S. NEMAT
NASSER, ed., Variational Methods in the Mechanics of Solids (Bergman Press, New York).
306 E. Christiansen
SUQUET, P.M. (1981), Sur les equations de la plasticity: Existence et rigularit6 des solutions, J. Mec.
20, 3-39.
TEMAM, R. (1977), Navier-Stokes Equations (North-Holland, Amsterdam).
TEMAM, R. (1981a), Existence theorems for the variational problems of plasticity, in: M. ATTEIA,
D. BANCEL and I. GUMOWSKI, eds., Nonlinear Problems of Analysis in Geometry and Mechanics
(Pitman, London).
TEMAM, R. (1981b), On the continuity of the trace of vector functions with bounded deformation,
Appl. Anal. 11, 291-302.
TEMAM, R. (1983), Problemes Mathematiques en Plasticit6 (Gauthier-Villars, Paris).
TEMAM, R. and G. STRANG (1980a), Functions of bounded deformation, Arch. Rational Mech. Anal.
75, 7-21.
TEMAM, R. and G. STRANG (1980b), Duality and relaxation in the variational problems of plasticity,
J. Mic. 19, 493-527.
TsUCHIYA, T. (1991), Global convergence property of the affine scaling methods for degenerate
linear programming problems, Math. Programming52, 377-404.
VANDERBEI, R.J., M.S. MEKETON and B.A. FREEDMAN (1986), A modification of Karmarkar's algo-
rithm, Algorithmica 1, 395-407.
ZAVELANI-RossI, A. (1979), Collapse load analysis of discretized continua, in: M.Z. COHN and G.
MAIER, eds., Proceedings 1977NATO Advanced Study Institute, University of Waterloo (Pergamon
Press, New York).
List of Symbols
Continuum mechanics notation
12 domain occupied by material.

closure of 12.
the boundary of 2.
N spatial dimension of problem (and of 12).
S fixed part of the boundary n12.
T loaded part of the boundary 0d2.
outward normal to 12.
f =f(x) field of volume forces in 12.
g = g(x) field of surface forces on T.
u = (x) field of plastic flow (displacement rates).
U space of plastic flow fields u.
Uo subspace of U with divergence zero.
Un flow in the interior of D2.
T
u flow on the loaded boundary T C 12.
E = E(U) strain rate tensor with components Eqij
A- = (x) field of stresses with components o-ij.
O', (WII, o'inI the principal stresses.
tr(or) trace of oa, i.e. the sum of diagonal elements.
rD deviator of the tensor ar.
I the identity tensor with components (ij).
space of stress tensor fields or.
0 quotient space /{qJI o E Ll(12)}.
elements of the quotient space o0.
B(x) set of admissible stress tensors at point x c 12.
r(x) B(x) yield condition at point x c 12.
K set of stress fields satisfying the yield condition at every point.
K0 quotient set K/{qIo I E L(12)} C 0.
F(u) work rate for external forces.
a(o, u) internal work rate in the material.
D(u) total energy dissipation rate associated with u.
A load multiplier.
A*, or*, * A, or, u as above in the collapse state.
Vu divergence of vector field (V is also applied to tensors).
U V tensor with components uivj.
307
308 E. Christiansen
m bending moments in plate problems.

u transversal displacement rate for plates.
b(m, u) internal work rate for plates.
mni, mnI principal bending moments.
X the space of plate bending moments.
Y the space of transversal displacement rates for plates.
f the transversal load distribution for plates.
Discrete variables
In general subscript 'h' indicates discrete variables.
Uh space of discrete displacement rate fields u h

space of discrete stress tensors o h.
Kh discrete stress tensors 'h, satisfying the yield condition at all
nodes.
Dh(uh) discrete analog of D (u), usually different from D (u h).
A*
h the discrete limit multiplier.
finite element triangulation of 2.
set of nodes for the single components of Orh.
set of nodes for the single components of uh.
basis functions for -h corresponding to the node v.
coefficients of 0Th w.r.t. the basis {4qj}.
basis functions for Uh corresponding to the node .
coefficients of Uh w.r.t. the basis k.
N
bTXy vector, whose coordinates are { j} in some linear ordering.
= mX= vector, whose coordinates are {/k} in some linear ordering.
vector with coordinates b, = F(ak).
bTy the discrete expression for F(u h).
A M x N 'stiffness' matrix with entries amn = a(4O'j, fk).
yTAx the discrete expression for a(O-h, Uh).
Dd(Y) the discrete expression for Dh(uh).
Kd {X i 1 N I Oh Kh}
the discrete stream function.
A, the matrix consisting of the columns of A, which are associated
with the node v.
Function spaces and norms
Cm(n) m times continuously differentiable functions in 2.

Cm() subspace of Cm(2) with continuous derivatives in 12.
CO (2) subset of C 0 (12) with compact support in 12.
List of Symbols 309
L P (Q) the measurable functions f in a2 with f IflP < cx.

Ilflp (Jfn Ifp)l/P (the norm in Lp(2)).
W m P(n) functions in Lp(f2) with all derivatives of order < m in LP().
II[llwmP norm on WnP(f2).
WoP () closure of Co(Q2) in W m 'P(Q).
M (n) the space of bounded measures in 2.
C(Q)* = (Co(F) : the space of bounded measures in Q2.
BD(f2) the space of bounded deformation in 2.
y(u) trace ("restriction") of u E BD(n2) to [Ll(0i2)]N.
VN vector functions with N components from a space V.
Il llv norm on the vector space V.
V* the dual of the normed space (V, Illlv)
(, )vXv* duality between V and V*.
Subject Index
O(h)-notation, 247 External work rate, 202, 219, 289
Antiplane shear, 209, 213, 218, 242, 247 Finite elements, 235, 283
- Argyris' triangle, 244
Bending moments, 285, 286 - Bell's triangle, 244
- discrete, 292 - biquadratic, 236, 245, 246
- space for, 288 - bicubic, 236, 245
Bound method, 251 - bilinear, 235, 236, 257, 262-265, 291
- lower, 206, 251 - Bogner-Fox-Schmit, 236, 245
- upper, 206, 251 - conforming, 235
Boundary condition, 202,217,227, 236, 256,284, - constant, 235,236,243,247,252,257, 262,264,
286 265
- essential, 203, 236, 286-288 - linear, 236, 240, 242-244, 247, 252, 274, 276,
- natural, 203, 286, 287, 289 291
Bounded deformation - quadratic, 236, 244, 245
- imbedding theorem, 215 Friction, 229
- space of, 208, 212, 214
- trace theorem, 215 Green's formula, 203, 212, 215, 223, 230, 231,
233, 234, 287
Clamped boundary, 286, 288, 297, 298, 301
Hydrostatic pressure, 204
Cohesion, 229
Collapse flow, 201
Incompressible flow, 205
- space for, 217
Internal friction, 229
Collapse solution, 209
Internal work rate, 203, 214, 219, 285
Collapse state, 200, 201, 203, 205, 212, 225, 231
Complementary slackness, 225, 234, 263
Kinematic principle, 205, 206, 214, 220, 225
Consistency condition, 241
- discrete, 255, 295
Coulomb yield condition, 216, 229, 275
- plates, 288
Kinematically admissible, 205
Dead load, 228
Deviator, 204, 216 Limit multiplier, 202, 205, 206, 220, 248, 261
Displacement rate, 199, 201, 202, 297 - computed values, 275
- discrete, 238, 292 Linear programming, 207, 259
- space for, 217, 289 Lipschitz continuous boundary, 214
- transversal, 289 Lower bound method, 206, 251
Divergence-free flow, 205,224,225,243,244,276
Dual problem, 240, 255, 267, 278, 295 MINOS, 272, 274, 276, 279, 296
Duality, 207, 214, 221, 235, 240, 288, 294
- reduced, 225, 226, 244 Necking, 284
Energy dissipation rate, 205, 221 Perfectly plastic, 199

Equilibrium equation, 202, 203, 210, 214, 220, Plane strain, 200, 204, 216, 230, 236, 240, 245,
286, 287 255, 256, 260
311
312 E. Christiansen
Plane stress, 216, 228 Superbasic variables, 274

Plastic region, 225 Supported boundary, 286, 288, 298
Plastified region, 200, 201,225, 243, 262, 266, 268 Surface forces, 202, 219
Pre-load, 228
Primal problem, 240, 255, 277, 295 Test problem, 200, 256, 257, 261
Principal moments, 286 Trace of tensor, 204
Principal stresses, 229 Tresca yield condition, 216, 229, 286
Regularity conditions, 249 Upper bound method, 206, 251

Richardson extrapolation, 247
Rigid, 199 Volume forces, 202, 219
Von Mises yield condition, 204, 216, 228, 240,
Saddle point, 205, 209, 214, 221, 241, 290, 294 245,256,260, 274, 276,286, 296
Safe stress field, 205
Simplex method, 207, 259, 262 Work rate
Simply supported, 288, 298, 299 - external, 202, 219, 289
Slack variables, 267 - internal, 203, 214, 219, 285
Stability condition, 249
Static principle, 205, 214, 220, 242 Yield condition, 203, 204, 214, 216, 228, 236
- discrete, 255, 295 - antiplane shear, 210
- plates, 288 - Coulomb, 216, 229, 275
Statically admissible, 205 - granular material, 229
Strain rate tensor, 202, 208, 214, 218, 219 - linearized, 259, 260, 262
- generalized, 218, 219, 225, 230, 233 - plates, 286, 289, 296
Stream function, 236, 276, 283 - Tresca, 216, 229, 286
Stress tensor, 202 - Von Mises, 204, 216, 228, 240, 245, 256, 260,
- discrete, 238 274, 276, 286, 296
- space for, 216

Limit Analysis of Collapse States

Uploaded by

Copyright:

Available Formats

Limit Analysis of Collapse States

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Limit Analysis of Collapse States

Uploaded by

Copyright:

Available Formats

Limit Analysis of

HANDBOOK OF NUMERICAL ANALYSIS, VOL. IV

CHAPTER I. Introduction 199

CHAPTER II. The Duality Problem 209

4. Why so different from elasticity? 209

CHAPTER III. Discretization by Finite Elements 235

CHAPTER IV. Solution of the Discrete Problem 255

CHAPTER V. Limit Analysis for Plates 285

LIST OF SYMBOLS 307

SUBJECT INDEX 311

experience on such problems to help us choose between the possible combina-

1. The physical problem

FIG. 1.1. Test problem in plane strain.

teresting we make symmetric thin cuts at x = 0 half way to the mid-plane y = 0

numerical solution methods in limit analysis for a continuous medium known

2. The mathematical model

ij(u) =21 xj duxi+)+ (2.1)

a(r, ) = f(oa, E(U))

=- (V.r).u + f (v.r)-u. (2.5)

where N is the spatial dimension (2 or 3) of the problem. Clearly o satisfies the

3. Remarks on the development of limit analysis

The Duality Problem

4. Why so different from elasticity?

In this section we try to point out some of the mathematical characteristics of

a(o, u) = oru' dx = - u'u dx + o(1)u(1). (4.1)

a(a, u) = AF(u) = A (- u dx + au(1)) (4.2)

The solution is easily seen to be

1 = F(u*) = u* dx +au*(l) = a (-2uo +u*(1))

u*(1) = 2u0 + -, (4.6)

5. Existence of the collapse solution

ij I4uix + aixi) (5.1)

BD(2) = {u [L1(12)]N ] sij(u) M(Q), 1 < i,j < N}. (5.2)

BD(Q) is a Banach space with the norm

y(u) = ulI a1 Vu C BD(12) n [C ( f2 )]N . (5.4)

THEOREM 5.2 (Imbedding theorem for BD). Assume that 2 C RN is bounded

THEOREM 5.3. Assume that 12 C RN is bounded with Lipschitz continuous

5.1. The stress space X

aD(x) B(x) x E f2,

THEOREM 5.4. Let 12 C RN be bounded with Lipschitz continuous boundary. Let

or E [L1(12)]N2 , D [L (2)]N 2, Vor e [LP(f2)]N.

I tr(,o)llLq/R < Cq([ll0 IIL + IIVorlILP) -

PROOF. See Lemma 2.1 in CHRISTIANSEN [1986]. [

5.2. The flow space U

r-- - (V.or).u + I(v.o).uT (5.16)

(ij, 'ij)C(T) xC() = f(Vo)U2 + J(J.).UT. (5.17)

U is an incomplete space equipped with the norm

THEOREM 5.5. Assume that n2 C RN is bounded with Lipschitz continuous

y(ufl)=uT on T and y(u a )=OonS. (5.20)

BD(f2) is a subset of U through the imbedding u - (, Y(u)IT) for all u E

(ijij)(c())*xC(Q) = (ij, ij)M()xC) J(-O')y

E (AXij, 'i)c ()xC() = (m, v .)M(T)xC(T7)

for all symmetric or E [C(1)]N2 .

5.3. The duality of limit analysis

a(or, u) =- (V.o').ua + f(vor).u T (o, u) x U. (5.21)

By Theorem 5.2 this is a continuous bilinear form on the product space X x U.

F(u) = /f.u+ fgu. (5.23)

In the special case, where T = 0Q2, i.e., there is no boundary condition on u,

We shall need the following theorem.

Finally it is convenient to introduce a short notation for the set of admissible

A*= inf D(u)= inf supa(,u), (5.29)

1 = F(u) = u dx +au(l) = a (-2uo +u(1))

Green's formula for the collapse fields (,u). On X x U, a(o, u) is defined