Stress Analysis by Finite Elements

R. A. WESTMANN, University of California, Los Angeles

•DETERMINATION of stress and displacement fields for realistic geometries often re-
quires use of numerical procedures. The applicability of any particular numerical
technique rests upon the speed of solution, accuracy of method, and capability of han-
dling complicated geometrices and different materials as well as mixed boundary conditions.
Two numerical methods are currently in vogue. The method of finite differences (1)
has been long established. In 1966, Schimming and Haas (5) described the application
of finite differences to the solution of a variety of problems in soil mechanics.
An alternate numerical technique, the so-called method of finite elements, was orig-
inally employed in the aerospace industry (7) and has been extensively developed recent -
ly. It has reached the stage of receiving textbook presentation (6); its use is becoming
widespread (2, 3, 4, 9). -
This presentatiOn-outlines the method of finite elements in one of its most elementary
forms. Attention is restricted to static two-dimensional problems involving only linear
materials. To illustrate the chief features and advantages of the technique several ex-
amples are solved and typical results presented.

BASIS OF PROCEDURE
This paper devotes attention to plane strain problems in classical elasto-statics. For
discussion purposes, the region of interest is denoted by A, the bounding surface by S,
and an x -y coordinate system utilized (Fig. 1). If u, v are the x, y components of the
displacement then the stress-displacement relations are given by

ax = (1 + v)~l - 2 v) ~ (1 - v) ~~ + v ~; f (la)

cry= (1 + v)E(i - 2 v) l v ~~ + (1 - v) ~; ~ (lb)

E ( au
r xy = 2 (1 + v) o Y +ox
av) (le)

where E, v are the usual engineering elastic constants. In addition, the stress com-
ponents must satisfy the differential equations of equilibrium:

acrx arxy
--+--+X=O (2a)
oX oY

ar aa
~ + - -y +Y= 0 (2b)
ax aY

where X, Y are the components of the body forces.

Paper sponsored by Committee on Mechanics of Earth Masses and Layered Systems .

46
 47
One suitable formulation of the problem is
s,
achieved by substituting from Eq. 1 into
Eq. 2 yielding the displacement equations of
equilibrium:
A
-T,
t
Ty

+ X (1 + r.i) (1 - 2 11) _ 0 ( 3a)

E -

Figure l.
_ 11 ) 0 v +.!.~
2 2
(1 - 211) o v + (l
2 2 2 2 oxay
ox oY

(1 + 11) (1 - 2 11) =0 (3b)

+Y E .

which must be satisfied throughout the region A. Statement of the problem is completed
upon observing that on the boundary S either the displacements or certain components of
the stress Eqs. 1 are specified.
While Eqs. 3 and associated boundary conditions are a suitable starting point for the
method of finite differences, an alternate approach is used for finite elements. An
equivalent formulation of the problem may be stated in terms of the minimization of an
integral instead of the solution of a system of partial differential equations.

r
The integral for the potential energy n of a body in a state of plane strain (1) is given
by -

n = 2 (1 + 11~( 1 - 211) f
A
11 ( 1 - II) [ ( ~~ +( ~ ; )2 J+ 2 II ~ ~ ~ ;
2
+ (1-211) (au
- +av)
- dx dy
2 aY ox

(4)

The double integration is carried out over the entire area A, whereas the line integral
is evaluated on only the put of the boundary, 81, on which the surface tractions Tx,
Ty are prescribed (Fig. 1).
The theorem of minimum potential energy then states: among all the displacement
fields which satisfy the geometric boundary conditions (u = U, v = V on 82, Fig. 1), that
displacement field which makes the potential energy nan absolute minimum is the solu-
tion to Eqs. 3 and associated boundary conditions.
A well known technique for obtaining approximate analytical solutions is as follows.
First, displacement fields which are functions of arbitrary parameters are selected.
These displacements (which must satisfy the geometric boundary conditions) are sub-
stituted into Eq. 4 and the indicated integration completed. Minimization of the po-
tential energy n with respect to the arbitrary parameters in the displacement field then
yields the approximate solution. The accuracy of this solution is consistent with the
approximations made in the assumed displacements.
The minimization process of Eq. 4 may be illustrated as follows. Let u, v be given
functions of x, y and parameters Ci which satisfy the geometric boundary conditions on
82;
N
u= I: 1 (x, y, c.)
i =
f.
1 1
(5a)
48

L.
Figure 2. Figure 3 .

v= j = ~N + 1 fj (x, y, cj) (5b)

Substitution from Eqs. 5 into Eq. 4 and completion of the spatial integration eliminates
the dependence upon x, y . The potential energy n is then a known function of only the
parameters Ci.

i = 1, 2, •.. , M

In this way the minimization of n is reduced to a standard maxima-minima problem

in the differential calculus of several variables. The minimization is assured upon re-
quiring that

i = 1, 2, . . . , M (6)

This leads to M simultaneous equations in the M unknown parameters Ci . Provided u,

v (Eqs. 5) are linear functions of the parameters, Eq. 6 yields simultaneous linear
equations.

FINITE ELEMENTS
In principle, the method of finite elements is identical to the preceding approach. The
technique is simply a numerical procedure which systemizes the selection of the dis-
placement fields and at the same time enables a greater degree of freedom to be intro-
duced into the assumed displacements.
To begin, the region A is divided into triangles , as s hown in F i gure 2 (other approac he s
and element shapes are possible; the method a s outline d here i s in its si mplest form) .
The nth triangle is illustrated in Figure 3. T he i, j, k nodes have coor dinates Xi, Yi ;
Xj , Yj; Xk, Yk defining the shape, size and location of the triangular element. The dis-
placement field in the nth element is assumed to be linear;

un = un (x, y) = a nx + b ny + c n (7a)

v = v (x, y) = d x + e y + f (7b)
n n n n n

where the constants an' bn' en, dn' en' fn are to be determined.
 49

Due to the linearity of the field, compatibility between elements is assured provided
adjacent elements have the same displacement components at common nodes. It is this
factor which determines the constants an, bn, en, dn. en, fn. Letting ui, Vi; Uj, Vj ;
uiv Vk denote the nodal displacement at the i, j, k nodes leads to the following relations
which the constants must satisfy:

u. == a x. + b y. + c (Ba)
1 n1 n1 n

(Bb)

(Be)

v.1 = d x. + e n1
n1
y. + f
n
(Bd)

V. == d X. + e y. + f (Be)
J n1 nJ n

(Bf)

Solution of Eqs. B permits the expression of an, bn, en, etc., in terms of the unknown
nodal displacements and the coordinates of the nodes:

an ={(yj - yk) ui + (yk - Yi)uj +(yi - Yj )uk} ~ (9a)

bn ={(xk-xj) ui +(xi - xk)uj +( xj - xi)uk} ~ (9b)

~
en = { ( xjyk - yjxk) ui + ( xky i - xiyk) uj + ( xiyj - y ixj) uk } 2 (9c)

dn ={(yj -yk)vi +(yk - Yi)vj +(yi - Yj )vk} 2 1~ (9d)

en={ (xk - xj )vi +(xi - xk )vj +(xj - xi) vk} ~ (9e)

where ~n is the area of the nth element.

(9g)

Substitution from Eqs. 9 into Eqs. 7 then gives the displacement field of the element
in terms of the unknown nodal displacements and the known geometry of the triangle.
The displacements obtained in this way for the entire region are continuous as all the
elements have the same nodal displacements at common nodes.
50

The stress state corresponding to the selected displacement field is readily calcu-
lated. From Eqs. 7, it follows that

oUn
-ax = a n (lOa)

av
n
--=e (lOb)
ay n

au av
____E+--.E:=b +d (lOc)
ay ax n n

so the components of stress in the nth element are

(lla)

(llb)

(llc)

The state of stress in each element is constant but, of course, varies from element to
element.
It remains to substitute from Eqs. 7 into the expression for the potential energy. The
potential energy of the nth element then becomes

+ (l -l v) ( bn + dn r} .An - { { x ( anx + bny +en)

n

+Y(dn' +•ny +'n) }dxdy {s, )!Tx(•nx + bny +en)

+ T y (d nx + e ny + fn ) } d L (12)

The boundary integral on (s1) n occurs only if an edge of the nth element forms part of
the 81 boundary of the region A. The total potential energy of the system is just the sum
of the potential energy of all the elements. If there are N elements altogether· then
 51

-J ~/{x(v 1
n
+ bny +en)+ Y(dn' + V + 'n)} dxdy

- ~ (s,).
n
J 0 1
{T (a x +
x n
bny + c n ) + Ty (ctnx + e ny + fn )}ctL (13)

In Eq. 13, the first sum is quadratic in the coefficients an, bn, dn, en and therefore
quadratic in the nodal displacements Ui, Vi. The second and third sums are linear in
the coefficients and therefore linear in the nodal displacements.
The minimization of the potential energy is performed with respect to the 2 M un-
known nodal displacements:

on E o 1 - v ( a 2 + e 2 ) + 2 va e
--=
aui 2(1+v)(l-2v)oui (
-
n~l
N {
( ) n n nn

~=
ovi
E
2(1+v)(l-2v) ovi (
_o_
n~l
N {
(
1- v
)( n
a2 + e2
n) + 2 va e
nn

+(I -,' v) ( bn + dn )'} "n )- a'v 1 ( J 1

{[H•nx+bn" en)

+Y (ctnx + e ny+ fn ) } d x d y)- -av0 i ( n ~~ 1 j ( ) l~ Tx (anx + bny+ c n )

81

+Ty ( dn' + eny + fn)} dL) o 0 i o I, 2, , ~, , M (14b)

52
The summation indicated in Eqs. 14 need only be carried out for the n elements im-
mediately surrounding the i th node . This leads to a system of simultaneous linear
algebraic equations which may be stated in the matrix fo rm

U2

K21 . . . . . . . . . . = (15)

uM B2M

K2Ml' ' ' ' ' ~M2M vM

The stiffness matrix, Kij, arises from the first sums in Eqs. 14; the nonhomogeneous
terms Bi arise from the second and third sums in Eqs. 14 and depend only on the body
forces and boundary conditions. One very important feature is that the stiffness matrix
Kij is symmetric. This follows from the quadratic nature of the pertinent terms in the
expressions for the potential energy.
To illustrate the preceding formulation a simple example is presented in detail. Con-
sider an isosceles triangle with sides t, t, ,_f2' t and loaded and supporled as shown in
Figure 4. Each side is subjected to a uniform normal pressure, p, so that the triangle
is in overall equilibrium. The supports are only included to prevent rigid body motion.
The triangular area is considered as one element (n = 1) with nodes 1, 2, 3 (Fig. 4).
From the support conditions it is clear that

(16a)

(16b)

Noting that ~1 = t 2/2, Eqs. 9 yield

(17a)

(17b)

In this particular example, the body forces X, Y

vanish but there is a contrib:ution to the potential
energy from the boundary forces. Since all three
sides of triangle are on the 81 boundary of the re-
gion A, they must be included in the boundary inte -
gral. Completion of the details of this integration
gives

Figure 4.
 53

Substituting from Eqs. 17 and Eq. 18 into Eq. 13 gives the expression for the po-
tential energy of the system:

Il -_ 2 (1 + v)E (l (1 - v)
_ 2 v) ~ I ( U22 +Vs2) + VU2Vs + (l - 2 V ) Us+
2 2pt (u2 +Vs)
(19)
Minimization of Eq. 19 requires that

leading to the matrix equation

2 (1 + v) (1 - 2 v)
E [lo:v
v
1 - v
0
0
0
2 (1 - 2 v)
] U2
Vs
us
l -pt
-i-
0
1
1 (20)

The solution of Eq. 20 is

v - pt (1 + v) (1 - 2 v) , us= 0 (21)
u2= s= E

Substitution of these values into Eqs. 11 gives the element stresses

(22)

a result which could have been anticipated from the start.

'l
e e e e
",10.i~o.~ A
A A
f>
"
0 ti., I
0. 2
I
0 .3
I
0 .4
I
o_y 0 .6
I
0 .7
I I
0.8 0,9
I
•.Io
R 1v

3 .0

2.0
0

2.0 e -x

1.0

0.5

0 0.1 0.2 0.3 0 .4 0-' 0 .6 0.7 08 0 .9 >O

x
Fi 1 - - - - 1 - p- 2' ---l
Figure 5. Figure 6.
54

v=0.35

30
"v
%.
2.0

1.0 -- - --- - - - -

l--2.__j
0 0.2 0.4 0 .6 0.8 10

X/o
Figure 7.

o.~ ~-----------------,
·-
2=
c-
A
-~
-0
Jl:Q , 35

T.y
•J...o,3
0.2
... -~~ ,

..
YI
-- .
01 '
', \\
o l...~~::;::===;;i:;==~;===:;o;:.iis=~'~1.o
l--20 - - I
x;.

Figure 8.

10 ~-----------------,

-~
A
II =0.35 ~-
D- -D
o.a

~
YI
P/20
-t--- ·
0 .6
A

z . ~I

Figure 9.
 55

In general Eqs. 15 are not as simple as in

the preceding example, Eq. 20. Invariably the
solution must be obtained numerically on the
digital computer. Several standard routines
are available for the efficient solution of large
systems of banded, symmetric, linear equa-
tions. Once the nodal displacements are deter-
t. mined from Eqs. 15, the approximate displace-
t500,000 ~ t 613 NODES , 560 ELEMENTS
5 1/2 MINUTES
ment field is fixed and the stress state in each
~ element may be evaluated from Eqs. 9 and 11.
Accuracy obtained depends heavily on the

I
number, size, and distribution of elements as
well as the character of the stress fields. [ It
should be noted that the formulation is improper
for incompressible materials, 11 = 1/2; for
v
9 (
I
E , v:Q. 3
1
' /
nearly incompressible materials, say v >0.490,
an alternate formulation must be adopted (see
r"\ _/
Herrmann, 10).] Only the experience of solv-
ing exampleproblems and comparing with exact
solutions permits one to assess the errors in-
herent in the approximate technique. The first
of the following examples helps to serve this
- purpose. The other three examples have been
I
9
I selected to illustrate other essential features
of the method of finite elements.
Figure 10.
Example 1
The first example is concerned with the linear
elastic analysis of a test configuration for con-
crete specimens. The specimen is a long circular cylinder subjected to diametrically
opposed line loads (Fig. 5). A closed form solution for this problem may be found in
Muskhelishvili (8); the numerical results are compared with this exact solution in
Figure 5. -
These numerical results are indicated by circles and triangles. The different sym-
bols are used to distinguish between elements of different size. Since the stress states
are averaged over the element, it might be expected that the numerical stresses best
correspond to the centroid of the element;
therefore, the elements closer to the di-
ameter s ho u l d give results in closer
agreement with the exact solution.
The exact solution shows that the nor-
mal stress across the vertical diameter 0.4

is constant but that a lateral concentrated 0 .6

force exists directly under the load. The
o.e
approximate solution deviates from this
close to the boundary. This deviation is 1.0

partly due to the fact that the triangular

elements are not centered on the vertical
~V(X,Ol
diameter. p

In addition to assessing the accuracy

of the method, this example also demon-
strates its application to a problem in-
volving a curved boundary.
II =0, 3
Example 2
The s e c o n d example deals with the
linear elastic analysis of a plane strain Figure 11.
56
version of the triaxial test for soils. Fig-
ure 6 illustrates the geometry and distribu-
tion of elements. The specimen is loaded
by rigid platens, and it is assumed that the
displacements across the platen-soil inter-
face are continuous. The components of
stress at different cross-sections are shown
in Figures 7, 8, and 9; it is clear that in
the upper one-quarter of the specimen the
elastic stress state is far from uniform.
It should be mentioned that the nonuni-
formity of the stress field depends on the
value of Poisson's ratio of the material.
Figure 12 . In the event the material has a Poisson's
ratio of zero, the stress state is uniform,
cry = P/ 2a, crx = Tx.y = 0. In this example
Poisson's ratio is taken to be v = 0. 35.
This example was selected both for its
interest in soil mechanics and to illustrate
how mixed boundary conditions are handled.

Example 3
The third example is concerned with the
problem of an elastic footing resting on an
elastic foundation. There is complete con-
tinuity between footing and foundation and
the footing is centrally loaded by a normal
force of P = 10 6 •
750 NODES
In Figure 10, the geometry is indicated
696 ELEMENTS as well as the distribution of elements; each
st MINUTE rectangle is composed of four triangles. Due
to the symmetry of the problem only one-
half of the geometry is shown. As the soil-
structure interaction is of prime interest,
this problem was solved for a variety of
ratios of the elastic moduli of the footing
and foundation materials.
Figure 11 shows the normal interface
displacement; the constant C has been left
undetermined on purpose, because the dis-
placements for this type of problem in-
Figure 13. v o 1v in g a ha If plane are not uniquely
determined (11).
In Figure 12, the interface normal stress
is presented for different ratios of material
moduli. For comparison purposes the an-
0 .0
1.0 2 .0 3 .0 •.O 6 .0 alytical result for a rigid footing, E2/ E1 =
oo, is included as well. The assumption of
-0.70
a uniform pressure dist r ibution under the
5
footing leads to a constant value of cry/10 =
-- o,eo
2. 5.
~ -0.90
2
10

E •10 8
Example 4
1
"1 & ,...2 =0 . 25 The final example deals with the stress
-1 .20
analysis of a series of linear elastic in-
clusions in a linear elastic layer loaded by
Figure 14. anormalsurface pressure p= lOOO(Fig.13).
 57
3 .0 Quantities of chief interest are the normal
2. 5 surface displacements and the stresses in-
SECTION A-A

2.0 i\ duced in the inclusion and surrounding ma-

' terial. Figures 14 and 15 illustrate these
~ quantities for several different ratios of the
10'
material moduli. In Figure 14, the stresses
are not plotted for those values of x for which
~ " IO O
EI the numerical solution did not give satisfactory
0 .0 3.0 4.0 5. 0 resolution.
x
This exam p 1 e illustrates how parameter
studies of soil-structure problems involving
several materials can be made.
3. 0

~- ~
SECTION 8-8
CONCLUSIONS
2.0 The preceding examples illustrate some of
.::..:i_
10' ..~
the chief advantages of the finite element
method. The ease with which problems in-
10
volving complex geometrices (two-dimensional),
~::1100
I nonhomogeneous mate r i a 1 properties, and
mixed boundary condiUons are handled make
30 4 ,0
this technique particularly attractive. In ad-
dition, it is not necessary to introduce ficti-
Figure 15. tious points at the boundaries or interfaces
and the formulation automatically leads to a
symmetric banded matrix. This latter point
is an important simplification as such systems are amenable to efficient numerical
solution.
One disadvantage is that the method of finite elements is not as widely applicable as
finite differences. Only those problems which have an equivalent variational formula-
tion may be attacked by finite elements while any partial differential equation may be
approximated by finite difference operators. Another disadvantage of the finite element
method, as presented herein, is that the stresses are averaged across any one element. It
is never clear what point in the triangular element corresponds to the averaged stress state.
The running times have been indicated for each example. All problems were run on
an IBM 7094 digital computer; the commercial costs are about $10 per minute. The
computer costs generally are small compared to the personnel time necessary for prop-
er evaluation of all output.
Many universities produce graduate students well versed in the details of finite ele-
ments. In addition, short courses and extension programs are being offered so those
interested may learn the techniques involved. Because of this, it is anticipated that the
method of finite elements will become a standard tool for elastic analysis within the
next few years.
Currently, industrial and academic groups are actively engaged in further develop-
ment and extension of the method. Problems in three dimensions, nonlinear material
properties, large deformations, viscoelasticity, and thermoelasticity are being re-
searched. With the continued improvement in computer technology it is apparent that
our analysis capabilities are significantly improved. Shortcomings in analysis need no
longer be a stumbling block in the solution of realistic problems.

ACKNOWLEDGMENTS
The finite element analysis discussed here was undertaken under N. S. F. grant GK-
626 in the Division of Engineering and Applied Science at the California Institute of
Technology. This work forms part of a general investigation into constitutive relations
for soils and their use in practical problems. The work is under the general super-
vision of R. F. Scott. The numerical details for Example 1 were done by A. Levine;
the second example is due to T. Y. Chang. The last two examples were done by the
author under an internal grant in the Department of Engineering at the University of
California, Los Angeles.
58

REFERENCES
1. Wang, C. T. Applied Elasticity. McGraw-Hill, 1953.
2. Wilson, E. L. Structural Analysis of Axisymmetric Solids. AIAA Jour., Vol. 3,
No. 12, p. 2269-2274, Dec. 1965.
3. Herrmann, L. R. Elastic Torsional Analysis of Irregular Shapes. Jour. of En -
gineering Mechanics Division, ASCE, Dec. 1965.
4. Zienkiewicz, 0., Mayer, P., and Cheung, Y. K. Solution of Anisotropic Seepage
by Finite Elements. Jour. of Engineering Mechanics Division, ASCE, Feb. 1966.
5. Schimming, B. B., and Haas, H. J. Numerical Determination of Stresses in Earth
Masses. Highway Research Record 145, p. 109-126, 1966.
6. Rubinstein, M. F. Matrix Computer Analysis of Structures. Prentice-Hall, 1966.
7. Turner, M. J., Clough, R. W., Martin, H. C., and Topp, L. J. Stiffness and
Deflection Analyses of Complex Structures. Jour. of Aeronautical Sciences,
Vol. 23, No. 9, p. 805, Se!Jt. 1956.
8. Muskhelishvili, N. I. Some Basic Problems of the Mathematical Theory of Elas-
ticity. P. Noordhoff Ltd., p. 330-334, 1963.
9. Argyris, J. H., and Patton, P. C. Computer Oriented Research in a University
Milieu. Applied Mechanics Rev., Vol. 19, No. 12, p. 1029-1039.
10. Herrmann, L. R. Elasticity Equations for Incompressible and Nearly Incompress-
ible Materials by a Variational Theorem. AIAA Jour., Vol. 3, No. 10,
p. 1896-1900, Oct. 1965.
11. Timoshenko, S., and Goodier, J. N. Theory of Elasticity, McGraw-Hill, p. 89-
92, 1951.