78

Chapter 6
1.

Analytical Techniques and Solutions for Plastic Solids

6.1 Slip-line field theory

P
6.1.1. The figure shows the slip-line field for a rigid plastic
double-notched bar deforming under uniaxial tensile loading. a
The material has yield stress in shear k
6.1.1.1. Draw the Mohr’s circle representing the state of
stress at A. Write down (i) the value of f at this 
point, and (ii) the magnitude of the hydrostatic B
stress s at this point.
A 
C
6.1.1.2. Calculate the value of f at point B, and deduce the
magnitude of s . Draw the Mohr’s circle of stress  slip-line b slip-line
at point B, and calculate the horizontal and vertical
components of stress
6.1.1.3. Repeat 1.1 and 1.2, but trace the b slip-line from P
point C to point B.
6.1.1.4. Find an expression for the force P that causes
plastic collapse in the bar.

Q P
6.1.2. The figure shows a slip-line field for oblique indentation of a
rigid-plastic surface by a flat punch B A
6.1.2.1. Draw the Mohr’s circle representing the state of
stress at A. Write down (i) the value of f at this q 
point, and (ii) the magnitude of the hydrostatic stress
s at this point. a
6.1.2.2. f
Calculate the value of at point B, and deduce the
magnitude of s .
6.1.2.3. Draw the Mohr’s circle representation for the stress state at B, and hence calculate the
tractions acting on the contacting surface, as a function of k and q .
6.1.2.4. Calculate expressions for P and Q in terms of k, a and q , and find an expression for Q/P
6.1.2.5. What is the maximum possible value of friction coefficient Q/P? What does the slip-line
field look like in this limit?
79

6.1.3. The figure shows the slip-line field for a rigid plastic P
double-notched bar under uniaxial tension. The material
has yield stress in shear k. The slip-lines are logarithmic
spirals, as discussed in Section 6.1.3. b
6.1.3.1. Write down a relationship between the angle  a
y , the notch radius a and the bar width b. y
C
6.1.3.2. Draw the Mohr’s circle representing the r q B
state of stress at A. Write down (i) the value
of f at this point, and (ii) the magnitude of A b
the hydrostatic stress s at this point.
6.1.3.3. Determine the value of f and the P
hydrostatic stress at point B, and draw the
Mohr’s circle representing the stress state at this point.
6.1.3.4. Hence, deduce the distribution of vertical stress along the line BC, and calculate the force P
in terms of k, a and b.

6.1.4. The figure shows the slip-line field for a notched, rigid plastic
bar deforming under pure bending (the solution is valid for a y
y < p / 2 - 1 , for reasons discussed in Section 6.1.3). The solid A
has yield stress in shear k.
 b
6.1.4.1. Write down the distribution of stress in the triangular
region OBD M b M
6.1.4.2. Using the solution to problem 2, write down the O
stress distribution along the line OA
6.1.4.3. Calculate the resultant force exerted by tractions on d  b
the line AOC. Find the ratio of d/b for the resultant C
B D
force to vanish, in terms if y , and hence find an
equation relating a/(b+d) and y .
6.1.4.4. Finally, calculate the resultant moment of the tractions about O, and hence find a
relationship between M, a, b+d and y .
6.1.4.5. Show that the slip-line field is valid only for b+d less than a critical value, and determine an
expression for the maximum allowable value for b+d.

6.1.5. Consider the problem in 6.1.4. Propose a slip-line field solution that is valid for y > p / 2 - 1 , and
use it to calculate the collapse moment in terms of relevant material and geometric parameters.

6.1.6. The figure shows the slip-line field for a rigid plastic M
double-notched bar subjected to a bending moment. The
slip-lines are logarithmic spirals. b

6.1.6.1. Write down a relationship between the angle  a
y , the notch radius a and the bar width b. y
6.1.6.2. Draw the Mohr’s circle representing the D Cq
state of stress at A. Write down (i) the value B rA
of f at this point, and (ii) the magnitude of b
b
the hydrostatic stress s at this point.
M
 80

6.1.6.3. Determine the value of f and the hydrostatic stress just to the right of point B
6.1.6.4. Hence, deduce the distribution of vertical stress along the line BC
6.1.6.5. Without calculations, write down the variation of stress along the line BD. What happens
to the stress at point B?
6.1.6.6. Hence, calculate the value of the bending moment M in terms of b, a, and k.
6.1.6.7. Show that the slip-line field is valid only for b less than a critical value, and determine an
expression for the maximum allowable value for b.

6.1.7. Consider the problem in 6.1.6. Propose a slip-line field solution that is valid for y > p / 2 - 1 , and
use it to calculate the collapse moment in terms of relevant material and geometric parameters.

6.1.8. A rigid flat punch is pressed into the surface of an elastic-perfectly plastic half-space, with Young’s
modulus E , Poisson’s ratio n and shear yield stress k. The punch is then withdrawn.
6.1.8.1. At maximum load the stress state under the punch can be estimated using the rigid-plastic
slip-line field solution (the solution is accurate as long as plastic strains are much greater
than elastic strains). Calculate the stress state P
in this condition (i) just under the contact, and
(ii) at the surface just outside the contact. A B
6.1.8.2. The unloading process can be assumed to be
elastic – this means that the change in stress f
during unloading can be calculated using the
solution to an elastic half-space subjected to
b slip-line a  slip-line
uniform pressure on its surface. Calculate the
change in stress (i) just under the contact, and
(ii) just outside the contact, using the solution P/a
given in Section 5.2.8.
6.1.8.3. Calculate the residual stress (i.e. the state of A B
Unloading
stress that remains in the solid after unloading)
at points A and B on the surface.
81

6.2. Bounding Theorems in plasticity and their applications: Plastic Limit

Analysis

6.2.1. The figure shows a pressurized cylindrical cavity. The solid has yield
stress in shear k. The objective of this problem is to calculate an upper
bound to the pressure required to cause plastic collapse in the cylinder
6.2.1.1. Take a volume preserving radial distribution of velocity as the
r
collapse mechanism. Calculate the strain rate associated with a
the collapse mechanism p
6.2.1.2. Apply the upper bound theorem to estimate the internal
pressure p at collapse. Compare the result with the exact b
solution
e2 P
6.2.2. The figure shows a proposed collapse mechanism for e1
indentation of a rigid-plastic solid. Each triangle slides as a rigid w
block, with velocity discontinuities across the edges of the q A C
B
triangles.
6.2.2.1. Assume that triangle A moves vertically downwards.
Write down the velocity of triangles B and C
6.2.2.2. Hence, calculate the total internal plastic dissipation,
and obtain an upper bound to the force P
6.2.2.3. Select the angle q that minimizes the collapse load.

6.2.3. The figure shows a kinematically admissible e2

velocity field for an extrusion process. The velocity of e1 q D E
the solid is uniform in each sector, with velocity
discontinuities across each line. The solid has shear C
yield stress k. P
H O 2H
6.2.3.1. Assume the ram EF moves to the left at
constant speed V. Calculate the velocity of B
the solid in each of the three separate
regions of the solid, and deduce the q A F
magnitude of the velocity discontinuity
between neighboring regions
6.2.3.2. Hence, calculate the total plastic dissipation and obtain an upper bound to the extrusion
force P per unit out-of-plane distance
6.2.3.3. Select the angle q that gives the least upper bound.

e2
e1 D E
6.2.4. The figure shows a kinematically admissible r
velocity field for an extrusion process. Material
particles in the annular region ABCD move along
C q P
radial lines. There are velocity discontinuities across H 2H
the arcs BC and AD. B

A F
 82

6.2.4.1. Assume the ram EF moves to the left at constant speed V. Use flow continuity to write
down the radial velocity of material particles just inside the arc AD.
6.2.4.2. Use the fact that the solid is incompressible to calculate the velocity distribution in ABCD
6.2.4.3. Calculate the plastic dissipation, and hence obtain an upper bound to the force P.

6.2.5. The purpose of this problem is to extend the upper bound theorem to
pressure-dependent (frictional) materials. Consider, in particular, a solid t*
with a yield criterion and plastic flow rule given by
m b
f (s ij ) = s e + s kk - Y = 0
2
e&e �f e&e �3Sij m �
e&ijp = = � + d ij �
2
1+ m / 2 s
� ij 1 + m / 2 �2s e 2 �
2

Sij = s ij - (s kk d ij ) / 3 s e = 3Sij Sij / 2 e&e = 2e&ijpe&ijp / 3

where m < 1 is a friction coefficient like material parameter. The solid is subjected to a traction ti* on its
exterior boundary � R and a body force bi per unit volume in its interior. The solid collapses for loading
b ti* , b bi , where b is a scalar multiplier to be deterined.
p
6.2.5.1. Show that the rate of plastic work associated with a plastic strain rate e&ij can be computed

as Y e&e / 1 + m 2 / 2 vt vn
e2
6.2.5.2. We need to understand the nature of the plastic
dissipation associated with velocity discontinuities in e1 h
this material. We can develop the results for a velocity
discontinuity by considering shearing (and associated
dilatation) of a thin layer of material with uniform
thickness h as indicated in the figure. Assume that the strain rate in the layer is
homogeneous, and that the surface at x2 = h has a uniform tangential velocity vt and
normal velocity vn . Show that (i) the rate of plastic work per unit area of the layer can be
computed as s 22 vn + s12vt = Y 2vn2 + vt2 / 3(1 + m 2 / 2) , and (ii) to satisfy the plastic flow

rule the velocities must be related by vn / vt = 3m / 2 (1 - m 2 ) . Note that these results are
independent of the layer thickness, and therefore (by letting h � 0 ) also characterize the
dissipation and kinematic constraint associated with a velocity discontinuity in the solid.
6.2.5.3. To state the upper bound theorem for this material we introduce a kinematically admissible
velocity field v, which may have discontinuities across a set of surfaces Ŝ in the solid.
Define the strain rate distribution associated with v as
� vi � vj �
ˆ = 1 ��
e& + � e& ˆ = 2e&ˆ e&
ˆ
ij e ij ij / 3
2� ��x j �xi
�
�
ˆ / 2 1 + m 2 / 2 in the interior of the solid, and
The velocity field must satisfy eˆ&kk = 3me&e
must satisfy
83

�vn �= �vt �m 3 / 2 1- m2 �vn �= ( vi+ - vi- ) mi

�vt �= vi+ - vi- - �vn �mi
on Ŝ , where mi denotes a unit vector normal to Ŝ . Define the plastic dissipation function
as
Y eˆ&e Y
F( v ) =� 2
dV +� 2
�vt �dA - bi vi dV - ti*vi dA
� �
R 1+ m / 2 Sˆ 3(1 - m ) R �
R
Show that (i) F ( &
u ) = 0 &
u
, where denotes the actual velocity field in the solid at collapse,
and (ii) F( v ) �0
6.2.5.4. Hence, show that an upper bound to the load factor at collapse can be calculated as
Y eˆ&e Y
� 2 �
dV +
2
�vt �dA
R 1+ m / 2 ˆ 3(1 - m )
S
bL �
� bi vi dV + ti*vi dA
�
R �
R

6.2.6. As an application of the results derived in the preceding problem,

consider a soil embankment with vertical slope, as shown in the q
figure. The soil has mass density r and can be idealized as a h
frictional material with constitutive equation given in the preceding
problem. Using a collapse mechanism consisting of shearing and
dilatation along the line AB shown in the figure (the angle q for the
optimal mechanism must be determined), calculate an upper bound to
the admissible height h of the embankment.

6.2.7. The figure shows a statically indeterminate structure. All bars have
cross-sectional area A, Young’s modulus E and uniaxial tensile yield
stress Y. The solid is subjected to a cyclic load with mean value P and 450 450
amplitude DP as shown
6.2.7.1. Select an appropriate distribution of residual stress in the
structure, and hence obtain a lower bound to the shakedown
limit for the structure. Show the result as a graph of DP as a P
function of P
6.2.7.2. Select possible cycles of plastic strain in the structure, and hence obtain an upper bound to
the shakedown limit for the structure.
You should be able to find residual stresses and plastic strain cycles that make the lower and upper
bounds equal, and so demonstrate that you have found the exact shakedown limit.

6.2.8. Calculate upper and lower bounds to the shakedown limit for a L/2 P
beam subjected to three point bending as shown in the figure. Assume
the applied load varies cyclically with mean value P and amplitude
DP . L/2
 84

6.2.9. The stress state induced by stretching a large plate containing a

cylindrical hole of radius a at the origin is given by e2
� �3a 4 a 2 � 3a 2 �
s11 = s 0 �1+ � - �cos 4q - cos 2q � r
� �2r 4 r 2 � 2 r 2 �
� � � � q
a
�
�a 2 3a 4 � a 2 �
s 22 = s 0 �
�2 - 4 � cos 4q - cos 2q � e1
�
�r 2r � � 2r 2 �
�
� �
��3a 4 a 2 � a2 �
s12 = s 0 �� 4- 2� sin 4q - sin 2q �
��2r r � 2r 2 �
�� � �
Use these results to calculate lower and upper bounds to the shakedown limit for the solid (assume that
s 0 varies periodically between zero and its maximum value)
85

7. Chapter 7

Introduction to Finite Element Analysis in Solid Mechanics

NOTE: The problems in the following section require a commercial finite element program. The problems
have been tested using the commercial version of ABAQUS/CAE Ver 6.6 (available from
http://www.simulia.com).

7.1. A guide to using finite element software

7.1.1. Please answer the following questions

7.1.1.1. What is the difference between a static and a dynamic FEA computation (please limit your
answer to a sentence!)
7.1.1.2. What is the difference between the displacement fields in 8 noded and 20 noded hexahedral
elements?
7.1.1.3. What is the key difference between the nodes on a beam element and the nodes on a 3D
solid element?
7.1.1.4. Which of the boundary conditions shown below properly constrain the solid for a plane
strain static analysis?
D

(i) (ii) (iii)

7.1.1.5. List three ways that loads can be applied to a finite element mesh
7.1.1.6. In a quasi-static analysis of a ceramic cutting tool machining steel, which surface would
you choose as the master surface, and which would you choose as the slave surface?
7.1.1.7. Give three reasons why a nonlinear static finite element analysis might not converge.

7.1.2. You conduct an FEA computation to calculate the natural frequency of vibration of a beam that is
pinned at both ends. You enter as parameters the Young’s modulus of the beam E, its area moment of
inertia I, its mass per unit length m and its length L. Work through the dimensional analysis to identify
a dimensionless functional relationship between the natural frequency and other parameters.
 86

7.1.3. Please answer the FEA related questions

7.1.3.1. What is the difference between a truss element and a solid element (please limit your
answer to a sentence!)
7.1.3.2. What is the difference between the displacement fields in 6 noded and 3 noded triangular
elements, and which are generally more accurate?
7.1.3.3. Which of the boundary conditions shown below properly constrain the solid for a static
analysis?

(a) (b) (c)

7.1.3.4. A linear elastic FE calculation predicts a maximum Mises stress of 100MPa in a

component. The solid is loaded only by prescribing tractions and displacements on its
boundary. If the applied loads and prescribed displacements are all doubled, what will be
the magnitude of the maximum Mises stress?
7.1.3.5. An FE calculation is conducted on a part. The solid is idealized as an elastic-perfectly
plastic solid, with Youngs modulus 210 GPa, Poisson ratio n = 0.3 . Its plastic properties
are idealized with Mises yield surface with yield stress 500MPa. The solid is loaded only
by prescribing tractions and displacements on its boundary. The analysis predicts a
maximum von-Mises stress of 400MPa in the component. If the applied loads and
prescribed displacements are all doubled, what will be the magnitude of the maximum
Mises stress?

7.1.4. The objective of this problem is to investigate

the influence of element size on the FEA 35cm
25cm
predictions of stresses near a stress concentration.
B C
Set up a finite element model of the 2D (plane R=2cm
strain) part shown below (select the 2D button 40cm
when creating the part, and make the little rounded C
20cm
radius by creating a fillet radius. Enter a radius of
2cm for the fillet radius). Use SI units in the
A D
computation (DO NOT USE cm!) 60cm

Use a linear elastic constitutive equation with E = 210GNm -2 , n = 0.3 . For boundary conditions, use
zero horizontal displacement on AB, zero vertical displacement on AD, and apply a uniform horizontal
displacement of 0.01cm on CD. Run a quasi-static computation.

Run computations with the following meshes:

7.1.4.1. Linear quadrilateral elements, with a mesh size 0.05 m
7.1.4.2. Linear quadrilateral elements, with mesh size 0.01 m
7.1.4.3. Linear quadrilateral elements with mesh size 0.005 m
87

7.1.4.4. Linear quadrilateral elements with a mesh size of 0.00125 m (this will have around 100000
elements and may take some time to run)
7.1.4.5. 8 noded (quadratic) quadrilateral elements with mesh size 0.005 m.
7.1.4.6. 8 noded (quadratic) quadrilateral with mesh size 0.0025 m.

For each mesh, calculate the maximum von Mises stress in the solid (you can just do a contour plot of
Mises stress and read off the maximum contour value to do this). Display your results in a table
showing the max. stress, element type and mesh size.

Clearly, proper mesh design is critical to get accurate numbers out of FEA computations. As a rough
rule of thumb the element size near a geometric feature should be about 1/5 of the characteristic
dimension associated with the feature – in this case the radius of the fillet. If there were a sharp corner
instead of a fillet radius, you would find that the stresses go on increasing indefinitely as the mesh size
is reduced (the stresses are theoretically infinite at a sharp corner in an elastic solid)

7.1.5. In this problem, you will run a series of tests to

compare the performance of various types of Built in Uniform pressure
element, and investigate the influence of mesh p=1 kN/m 2
design on the accuracy of a finite element
computation.
e2
e1 h=5.0
We will begin by comparing the behavior of
different element types. We will obtain a series of e3
finite element solutions to the problem shown
below. A beam of length L and with square cross L b=5.0
section b �b is subjected to a uniform distribution Elastic material
of pressure on its top face. E=210 GPa, n= Dimensions in cm

First, recall the beam theory solution to this problem. The vertical deflection of the neutral axis of the
beam is given by
u2 ( x3 ) = - {
w
24 EI }
( L - x3 ) 4 + 4 L3 x3 - L4
Here, w = bp is the load per unit length acting on the beam, and I = bh3 /12 is the area moment of
inertia of the beam about its neutral axis. Substituting and simplifying, we see that the deflection of the
neutral axis at the tip of the beam is
3 pL4
u2 ( L ) =
2 Eh3
Observe that this is independent of b, the thickness of the beam. A thick beam should behave the same
way as a thin beam. In fact, we can take b << h , in which case we should approach a state of plane
stress. We can therefore use this solution as a test case for both plane stress elements, and also plane
strain elements.
 88

7.1.5.1. First, compare the predictions of beam theory with a finite element solution. Set up a plane
stress analysis, with L=1.6m, h=5cm, E=210GPa, n = 0.33 and p=100 kN/m 2 . Constrain
both u1 and u2 at the left hand end of the beam. Generate a mesh of plane stress, 8 noded
(quadratic) square elements, with a mesh size of 1cm. Compare the FEM prediction with
the beam theory result. You should find excellent agreement.
7.1.5.2. Note that beam theory does not give an exact solution to the cantilever beam problem. It is
a clever approximate solution, which is valid only for long slender beams. We will check to
see where beam theory starts to break down next. Repeat the FEM calculation for L=0.8,
L=0.4, L=0.3, L=0.2, L=0.15, L=0.1. Keep all the remaining parameters fixed, including
the mesh size. Plot a graph of the ratio of the FEM deflection to beam theory deflection as
a function of L.

You should find that as the beam gets shorter, beam theory underestimates the deflection.
This is because of shear deformation in the beam, which is ignored by simple Euler-
Bernoulli beam theory (there is a more complex theory, called Timoshenko beam theory,
which works better for short beams. For very short beams, there is no accurate
approximate theory).
7.1.5.3. Now, we will use our beam problem to compare the performance of various other types of
element. Generate a plane stress mesh for a beam with L=0.8m, h=5cm, E=210GPa,
n = 0.33 , p=100 kN/m 2 , but this time use 4 noded linear elements instead of quadratic 8
noded elements. Keep the mesh size at 0.01m, as before. Compare the tip deflection
predicted by FEA with the beam theory result. You should find that the solution is
significantly less accurate. This is a general trend – quadratic elements give better results
than linear elements, but are slightly more expensive in computer time.
7.1.5.4. Run a similar test to investigate the
performance of 3D elements.
Generate the 3D meshes shown
above, using both 4 and 8 noded
hexahedral elements, and 4 and 10
noded tetrahedral elements. (Don’t
attempt to model 2 beams together as
shown in the picture; do the
computations one at a time otherwise
they will take forever).

Prepare a table showing tip deflection for 4 and 8 noded plane stress elements, 8 and 20
noded hexahedral elements, and 4 and 10 noded tetrahedral elements.

Your table should show that quadratic, square elements generally give the best
performance. Tetrahedra (and triangles, which we haven’t tried … feel free to do so if you
like…) generally give the worst performance. Unfortunately tetrahedral and triangular
elements are much easier to generate automatically than hexahedral elements.
89

7.1.6. We will continue our comparison of element types. Set up the beam problem again with L=1.6m,
h=5cm, E=210GPa, n = 0.33 , p=100 kN/m 2 , but this time use the plane strain mesh shown in the
figure.

There are 16 elements along the length of the beam and

5 through the thickness. Generate a mesh with fully
integrated 4 noded elements.

You should find that the results are highly inaccurate. Correct deformation
Similar problems occur in 3d computations if the
elements are severely distorted – you can check this out
too if you like. Undeformed
Integration points
This is due to a phenomenon know as `shear locking:’ FEM deformation
the elements interpolation functions are unable to
approximate the displacement field in the beam
accurately, and are therefore too stiff. To understand this, visualize the deformation of a material
element in pure bending. To approximate the deformation correctly, the sides of the finite elements
need to curve, but linear elements cannot do this. Instead, they are distorted as shown. The
material near the corners of the element is distorted in shear, so large shear stresses are generated in
these regions. These large, incorrect, internal forces make the elements appear too stiff.

There are several ways to avoid this problem. One approach is to use a Integration points
special type of element, known as a `reduced integration’ element.
Recall that the finite element program samples stresses at each
integration point within an element during its computation. Usually,
these points are located near the corners of the elements. In reduced
integration elements, fewer integration points are used, and they are
located nearer to the center of the element (for a plane stress 4 noded
quadrilateral, a single integration point, located at the element center, is FEM deformation
used, as shown). There is no shear deformation near the center of the
element, so the state of stress is interpreted correctly.

7.1.6.1. To test the performance of these reduced integration elements, change the element type in
your computation to 4 noded linear quadrilaterals, with reduced integration. You should
find much better results, although the linear elements are now a little too flexible.

7.1.6.2. Your finite element code may also contain more sophisticated element formulations
designed to circumvent shear locking. `Incompatible mode’ elements are one example. In
these elements, the shape functions are modified to better approximate the bending mode of
deformation. If your finite element code has these elements, try them, and compare the
finite element solution to the exact solution.
 90

7.1.7. This problem demonstrates a second type of element locking,

known as ‘Volumetric locking’. To produce it, set up the boundary
value problem illustrated in the figures. Model only one quarter of Symmetry
A
boundary
the plate, applying symmetry boundary conditions as illustrated in 10cm
the mesh. Assign an elastic material to the plate, with , E=210GPa, Symmetry
n = 0.33 . Run the following tests: e2 boundary

7.1.7.1. Run the problem with fully integrated 8 noded plane e1

strain quadrilaterals, and plot contours of horizontal, 100 MPa
10cm
vertical and Von-Mises equivalent stress.
7.1.7.2. Modify to increase Poisson’s ratio to 0.4999 (recall that
this makes the elastic material almost incompressible,
like a rubber). horizontal, vertical and Von-Mises
equivalent stress.. You should find that the mises stress
contours look OK, but the horizontal and vertical
stresses have weird fluctuations. This is an error – the
solution should be independent of Poisson’s ratio, so all
the contours should look the way they did in part 1.
The error you observed in part 2 is due to volumetric locking.
Suppose that an incompressible finite element is subjected to
hydrostatic compression. Because the element is incompressible, this loading causes no change in
shape. Consequently, the hydrostatic component of stress is independent of the nodal displacements,
and cannot be computed. If a material is nearly incompressible, then the hydrostatic component of
stress is only weakly dependent on displacements, and is difficult to compute accurately. The shear
stresses (Mises stress) can be computed without difficulty. This is why the horizontal and vertical
stresses in the example were incorrect, but the Mises stress was computed correctly. You can use two
approaches to avoid volumetric locking.
7.1.7.3. Use `reduced integration’ elements for nearly incompressible materials. Switch the element
type to reduced integration 8 noded quads, set Poisson’s ratio to 0.4999 and plot contours of
horizontal, vertical and Mises stress. Everything should be fine.
7.1.7.4. You can also use a special `Hybrid element,’ which computes the hydrostatic stress
independently. For fully incompressible materials, you must always use hybrid elements –
reduced integration elements will not work. Run the problem with hybrid elements using a
Poisson’s 0.4999, and plot the same stress contours. As before, everything should work
perfectly.
Clearly, elements must be selected with care to ensure accurate finite element computations. You should
consider the following guidelines for element selection:
 Avoid using 3 noded triangular elements and 4 noded tetrahedral elements, except for filling in regions
that may be difficult to mesh.
 6 noded triangular elements and 10 noded tetrahedral elements are acceptable, but quadrilateral and
brick elements give better performance.
 Fully integrated 4 noded quadrilateral elements and 8 noded bricks are usually specially coded to avoid
volumetric locking, but are susceptible to shear locking. They can be used for most problems, although
quadratic elements generally give a more accurate solution for the same amount of computer time.
 Use quadratic, reduced integration elements for general analysis work, except for problems involving
large strains or complex contact.
 Use quadratic, fully integrated elements in regions where stress concentrations exist. Elements of this
type give the best resolution of stress gradients.
 Use a fine mesh of linear, reduced integration elements or hybrid elements for simulations with very
large strains.
91

7.1.8. A bar of material with square cross section with base 0.05m and length 0.2m is made from an
isotropic, linear elastic solid with Young’s Modulus 207 GPa and Poisson’s ratio 0.3. Set up your
commercial finite element software to compute the deformation of the bar, and use it to plot one or
more stress-strain curves that can be compared with the exact solution. Apply a cycle of loading that
first loads the solid in tension, then unloads to zero, then loads in compression, and finally unloads to
zero again.

7.1.9. Repeat problem 7, but this time model the constitutive response of the bar as an elastic-plastic solid.
Use elastic properties listed in problem7, and for plastic properties enter the following data
Plastic Strain Stress/MPa
0 100
0.1 150
0.5 175
Subject the bar to a cycle of axial displacement that will cause it to yield in both tension and compression
(subjecting one end to a displacement of +/-0.075m should work). Plot the predicted uniaxial stress-strain
curve for the material. Run the following tests:
7.1.9.1. A small strain computation using an isotropically hardening solid with Von-Mises yield
surface
7.1.9.2. A large strain analysis using an isotropically hardening solid with Von-Mises yield surface.
7.1.9.3. A small strain computation using a kinematically hardening solid
7.1.9.4. A large strain analysis with kinematic hardening

7.1.10. Use your commercial software to set up a model of a

2D truss shown in the figure. Make each member of the
truss 2m long, with a 0.05m 2 steel cross section. Give
the forces a 1000N magnitude. Mesh the structure using
truss elements, and run a static, small-strain computation

Use the simulation to compute the elastic stress in all the

members. Compare the FEA solution with the analytical result.

7.1.11. This problem has several objectives: (i) To demonstrate FEA analysis with contact; (ii) To illustrate
nonlinear solution procedures and (iii) to demonstrate the effects of convergence problems that
frequently arise in nonlinear static FEA analysis.

Set up commercial software to solve the 2D (plane strain) R=40cm

contact problem illustrated in the figure. Use the Rigid surface
following procedure
 Create the part ABCD as a 2D deformable solid with a D C
40cm Plane strain
homogeneous section. Make the solid symmetrical
deformable solid
about the x2 axis A B
 Create the cylindrical indenter as a 2D rigid analytical 60cm
solid. Make the cylinder symmetrical about the x2
axis
 Make the block an elastic solid with E = 210GNm -2 , n = 0.3
 Make the rigid surface just touch the block at the start of the analysis.
 92

 Set the properties of the contact between the rigid surface and the block to specify a `hard,’
frictionless contact.
 To set up boundary conditions, (i) Set the vertical displacement of AB to zero; (ii) Set the horizontal
displacement of point A to zero; and (iii) Set the horizontal displacement and rotation of the
reference point on the cylinder to zero, and assign a vertical displacement of -2cm to the reference
point.
 Create a mesh with a mesh size of 1cm with plane strain quadrilateral reduced integration elements.

7.1.11.1. Begin by running the computation with a perfectly elastic analysis – this should run very
quickly. Plot a graph of the force applied to the indenter as a function of its displacement.
7.1.11.2. Next, try an elastic-plastic analysis with a solid with yield stress 800MPa. This will run
much more slowly. You will see that the nonlinear solution iterations constantly fail to
converge – as a result, your code should automatically reduce the time step to a very small
value. It will probably take somewhere between 50 and 100 increments to complete the
analysis.
7.1.11.3. Try the computation one more time with a yield stress of 500MPa. This time the
computation will only converge for a very small time-step: the analysis will take at least
150 increments or so.

7.1.12. Set up your commercial finite element software to conduct an explicit dynamic calculation of the
impact of two identical spheres, as shown below.
Use the following parameters:
 Sphere radius – 2 cm
 Mass density r = 1000 kg / m3
 Young’s modulus E = 210GNm -2 , Poisson ratio n = 0.3
 First, run an analysis with perfectly elastic spheres. Then
repeat the calculation for elastic-plastic spheres, with yield
stress s Y = 10000MNm -2 s Y = 5000 MNm -2
s Y = 1000 MNm-2 , s Y = 500 MNm -2 , and again with s Y = 50 MNm -2
 Contact formulation – hard contact, with no friction
 Give one sphere an initial velocity of v0 =100m/s towards the other sphere.

Estimate a suitable time period for the analysis and step size based on the wave speed.

Finally, please answer the following questions:

7.1.12.1. Suppose that the main objective of the analysis is to compute the restitution coefficient of
the spheres, defined as e = v A1 - v B1 / v A0 - v B 0 where v A0 , v A1 denote the initial and
final velocities of sphere A, and the same convention is used for sphere B. List all the
material and geometric parameters that appear in the problem.
7.1.12.2. Express the functional relationship governing the restitution coefficient in dimensionless
form. Show that for a perfectly elastic material, the restitution coefficient must be a
function of a single dimensionless group. Interpret this group physically (hint - it is the
ratio of two velocities). For an elastic-plastic material, you should find that the restitution
coefficient is a function of two groups.
7.1.12.3. If the sphere radius is doubled, what happens to the restitution coefficient? (DON’T DO
ANY FEA TO ANSWER THIS!)
7.1.12.4. Show that the kinetic energy lost during impact can be expressed in dimensionless form as
93

DK r v02 r v02
= f ( , )
R3 r v02 sY E
Note that the second term is very small for any practical application (including our
simulation), so in interpreting data we need only to focus on behavior in the limit
r v02 / E � 0 .
7.1.12.5. Use your plots of KE as a function of time to determine the change in KE for each analysis
case. Hence, plot a graph showing DK / R3 r v02 as a function of r v02 / s Y
7.1.12.6. What is the critical value of r v02 / s Y where no energy is lost? (you may find it helpful to
plot DK / R3 r v02 against log( r v02 / s Y ) to see this more clearly). If no energy is lost, the
impact is perfectly elastic.
7.1.12.7. Hence, calculate the critical impact velocity for a perfectly elastic collision between two
spheres of (a) Alumina; (b) Hardened steel; (c) Aluminum alloy and (d) lead
7.1.12.8. The usual assumption in classical mechanics is that restitution coefficient is a material
property. Comment briefly on this assumption in light of your simulation results.

7.2. A simple Finite Element program

7.2.1. Modify the simple FEA code in FEM_conststrain.mws to solve
problems involving plane stress deformation instead of plane strain (this 4 3
should require a change to only one line of the code). Check the modified
code by solving the problem shown in the figure. Assume that the block 2
has unit length in both horizontal and vertical directions, use Young’s 1
modulus 100 and Poisson’s ratio 0.3, and take the magnitude of the e2 2
distributed load to be 10 (all in arbitrary units). Compare the predictions of 1
the FEA analysis with the exact solution. e1

7.2.2. Modify the simple FEA code in FEM_conststrain.mws Axis of 6 5 4

to solve problems involving axially symmetric solids. The symmetry
figure shows a representative problem to be solved. It 2 4
1
represents a slice through an axially symmetric cylinder, r
which is prevented from stretching vertically, and ez 3

pressurized on its interior surface. The solid is meshed er 2

1 3
using triangular elements, and the displacements are
interpolated as
ui ( x1, x2 ) = ui( a ) N a ( x1 , x2 ) + ui(b) Nb ( x1 , x2 ) + ui(c) N c ( x1 , x2 )
where
 94

N a ( x1 , x2 ) =
(x 2 - x2(b) )(x( c)
1 ) (
- x1(b) - x1 - x1(b) )(x -x )
( c)
2
( b)
2

(x -x (a)
2
(b )
2 )(x(c )
1 - x1(b) ) -( x - x ) ( x - x )
(a)
1
(b )
1
(c)
2
(b)
2

N (x , x ) =
(x -x 2
(c )
2 )( x
( a)
1 - x1(c) ) - ( x - x )( x - x )
1
(c )
1
(a)
2
( c)
2

(x -x )(x ) -( x - x )( x - x )
b 1 2 (b ) (c) (a)
2 2 1 - x1(c) (b )
1
(c )
1
(a)
2
(c)
2

N (x , x ) =
(x -x 2
(a)
2 )( x
(b)
1 - x1( a ) ) - ( x - x )( x - x )
1
(a)
1
(b )
2
( a)
2

(x -x )(x ) -( x - x )( x - x )
c 1 2 (c ) (a) (b)
2 2 1 - x1( a) (c )
1
( a)
1
(b )
2
( a)
2
7.2.2.1. Show that the nonzero strain components in the element can be expressed as
�e rr �
�e � T
[ e ] = �e zz � [ u ] = �
�
ur( a ) u (za ) ur(b) u (zb ) ur(c ) u z(c ) ur( d ) u z( d ) �
�
�qq �
� �
�2e rz �

�� Na �Nb �Nc �
�� 0 0 0 �
r �r � r
� �
�0 �
Na �
Nb �Nc �
0 0
� � z � z �z �
[ B] = � �
�N a 0
Nb
0
Nc
0 �
�r r r �
�� Na � Na � Nb � Nb � Nc � Nc �
� �
�� z � r �z � r � z �r �
7.2.2.2. Let s = [ s rr s zz s qq s rz ] denote the stress in the element. Find a matrix [ D] that
satisfies [ s ] = [ D ] [ e ]
7.2.2.3. Write down an expression for the strain energy density U el of the element.
7.2.2.4. The total strain energy of each element must be computed. Note that each element
represents a cylindrical region of material around the axis of symmetry. The total strain
energy in this material follows as
W el = �2p rU el drdz
Ael
The energy can be computed with sufficient accuracy by evaluating the integrand at the
centroid of the element, and multiplying by the area of the element, with the result
W el = 2p Ael r U el
where r denotes the radial position of the element centroid, and U el is the strain energy
density at the element centroid. Use this result to deduce an
expression for the element stiffness, and modify the procedure (c)
elstif() in the MAPLE code to compute the element stiffness. t
7.2.2.5. The contribution to the potential energy from the pressure
acting on element faces must also be computed. Following
the procedure described in Chapter 7, the potential energy is e2 s L
(a) (b)
e1
95

�
P = - 2p r ti ui ds
0
where
s � s� s � s�
ui = ui( a ) + ui(c) �
1 - � r = r ( a ) + r (c) �
1- �
L � L� L � L�
and ui( a ) , ui(c) denote the displacements at the ends of the element face, and r ( a ) , r (c)
denote the radial position of the ends of the element face. Calculate an expression for P of
the form
P element = - [ t1 A t2 A t1B t2 B ] �
�u ( a ) u2( a) u1(c) u2(c) �
�1 �
where A and B are constants that you must determine. Modify the procedure elresid() to
implement modified element residual.
7.2.2.6. Test your routine by calculating the stress in a pressurized cylinder, which has inner radius
1, exterior radius 2, and is subjected to pressure p=1 on its internal bore (all in arbitrary
units), and deforms under plane strain conditions. Compare the FEA solution for
displacements and stresses with the exact solution. Run tests with different mesh densities,
and compare the results with the analytical solution.

7.2.3. Modify the simple FEA code in FEM_conststrain.mws to solve problems e2 (c)
which involve thermal expansion. To this end e1
7.2.3.1. Consider a generic element in the mesh. Assume that the material
inside the element has a uniform thermal expansion coefficient  ,
and its temperature is increased by DT . Let [B] and [D] denote the
matrices of shape function derivatives and material properties (a) (b)
defined in Sections 7.2.4, and let q = DT [ 1,1, 0]
T
denote a thermal
strain vector. Write down the strain energy density in the element, in terms of these
quantities and the element displacement vector u element .
7.2.3.2. Hence, devise a way to calculate the total potential energy of a finite element mesh,
accounting for the effects of thermal expansion.
7.2.3.3. Modify the FEA code to read the thermal expansion coefficient and the change in
temperature must be read from the input file, and store them as additional material
properties.
7.2.3.4. Modify the FEA code to add the terms associated with thermal expansion to the system of
equations. It is best to do this by writing a procedure that computes the contribution to the
equation system from one element, and then add a section to the main analysis procedure to
assemble the contributions from all elements into the global system of equations.
7.2.3.5. Test your code using the simple test problem

7.2.4. Modify the simple FEA code in FEM_conststrain.mws to solve plane

stress problems using rectangular elements. Use the following procedure. e2
To keep things simple, assume that the sides of each element are parallel to
(d) (c)
the e1 and e2 axes, as shown in the picture. Let (u1( a ) , u2( a ) ) , (u1(b) , u2(b) )
, (u1(c ) , u2(c ) ) , (u1( d ) , u2( d ) ) denote the components of displacement at H
nodes a, b, c, d. The displacement at an arbitrary point within the element (a) (b) e1
can be interpolated between values at the corners, as follows
B
 96

u1 = (1 - x )(1 - h )u1( a ) + x (1 - h )u1(b) + xh u1(c ) + (1 - x )h u1( d )

u2 = (1 - x )(1 - h )u2( a ) + x (1 - h )u2(b) + xh u2(c ) + (1 - x )h u2( d )
where
x = x1 / B, h = x2 / H
7.2.4.1. Show that the components of nonzero infinitesimal strain at an arbitrary point within the
element may be expressed as [ e ] = [ B ][ u ] , where
�e11 �
T
[ e ] = �e 22 �
�
� [ u] = �
u (a)
�1
u2( a ) u1(b ) u2(b ) u1(c ) u2(c ) u1( d ) u2( d ) �
�
�e12 �
�- 1 � x � 1 � x2 � x2 - x2 �
1- 2 � 1-
� B� � H �
0
B
�
� H
�
�
0
BH
0
BH
0
�
� �
� 1 � x1 � - x1 x1 1 � x1 � �
[ B] = � 0 - � 1- � 0 0 0 �1- �
�
H� B� BH BH H� B�
� �
� 1 � x1 � 1 � x2 � x1 1 � x2 � x1 x2 1 � x1 � - x2 �
- �1- � - �1- � - �1- � � 1- �
�
� 2H � B � 2B � H � 2 BH 2 B � H � 2 BH 2 BH 2H � B � 2 BH �
�
7.2.4.2. Modify the section of the code which reads the element connectivity, to read an extra node
for each element. To do this, you will need to increase the size of the array named connect
from connect(1..nelem,1..3) to connect(1..nelem,1..4), and read an extra integer node
number for each element.
7.2.4.3. In the procedure named elstif, which defines the element stiffness, you will need to make
the following changes. (a) You will need to modify the [B] matrix to look like the one in
Problem 11.1. Don’t forget to change the size of the [B] array from bmat:=array(1..3,1..6)
to bmat:=array(1..3,1..8). Also, note that as long as you calculate the lengths of the element
sides B and H correctly, you can use the [B] matrix given above even if node a does not
coincide with the origin. This is because the element stiffness only depends on the shape
of the element, not on its position. (b) To evaluate the element stiffness, you cannot assume
that [ B ]T [ D][ B] is constant within the element, so instead of multiplying [ B ]T [ D][ B] by
the element area, you will need to integrate over the area of the element
HB
K elem �=
�
� � �
�[ B ] T [ D ] [ B ]dx1dx2
0 0
Note that MAPLE will not automatically integrate each term in a matrix. There are various
ways to fix this. One approach is to integrate each term in the matrix separately. Let
[ A] = [ B ] [ D][ B ]
T

then for i=1..8, j=1..8 let

HB
kijelem = �
�aij dx1dx2
0 0
Use two nested int() statements to do the integrals. Note also that to correctly return a
matrix value for elstif, the last line of the procedure must read elstif=k, where k is the fully
assembled stiffness matrix.
7.2.4.4. Just before the call to the elstif procedure, you will need to change the dimensions of the
element stiffness matrix from k:=array(1..6,1..6) to k:=array(1..8,1..8).
7.2.4.5. You will need to modify the loop that assembles the global stiffness matrix to include the
fourth node in each element. To do this, you only need to change the lines that read
97

for i from 1 to 3 do
to
for i from 1 to 4 do
and the same for the j loop.

7.2.4.6. You will need to modify the part of the routine that calculates the residual forces. The only
change required is to replace the line reading
>pointer := array(1..3,[2,3,1]):
with
>pointer:= array(1..4,[2,3,4,1]):
7.2.4.7. You will need to modify the procedure that calculates element strains. Now that the strains
vary within the element, you need to decide where
to calculate the strains. The normal procedure t2=3.0
would be to calculate strains at each integration
point within the element, but we used MAPLE to E=10 1.0
evaluate the integrals when assembling the stiffness
matrix, so we didn’t define any numerical
integration points. So, in this case, just calculate the n= 1.0
strains at the center of the element.
e2
7.2.4.8. To test your routine, solve the problem shown in the
figure (dimensions and material properties are in
arbitrary units).
e1 2.0 1.0

7.2.5. In this problem you will develop and apply a finite element
method to calculate the shape of a tensioned, inextensible cable L
subjected to transverse loading (e.g. gravity or wind loading). x w(x) q(x)
The cable is pinned at A, and passes over a frictionless pulley at
B. A tension T is applied to the end of the cable as shown. A A B
(nonuniform) distributed load q(x) causes the cable to deflect by T
a distance w(x) as shown. For w<<L, the potential energy of the
system may be approximated as
L 2 L 1 2 3 4 5 6 7
T �dw �
V ( w) =� � �dx - qwdx
2 �dx � � B
0 0 A s
To develop a finite element scheme to calculate w, divide
w
a wb
the cable into a series of 1-D finite elements as shown. (a)
Consider a generic element of length l with nodes a, b at its l (b)
ends. Assume that the load q is uniform over the element,
and assume that w varies linearly between values wa , wb at the two nodes.

7.2.5.1. Write down an expression for w at an arbitrary distance s from node a, in terms of wa , wb ,
s and l.
7.2.5.2. Deduce an expression for dw / dx within the element, in terms of wa , wb and l
7.2.5.3. Hence, calculate an expression for the contribution to the potential energy arising from the
element shown, and show that element contribution to the potential energy may be
expressed as
1 �T / l -T / l ��
wa � ql / 2 �
�
V elem = [ wa , wb ] � �� � - [ wa , wb ] � �
2 �-T / l T / l ��
wb � ql / 2 �
�
 98

7.2.5.4. Write down expressions for the element stiffness matrix and residual vector.
7.2.5.5. Consider the finite element mesh shown in the figure. The loading q0 is uniform, and each
element has the same length. The cable tension is T. Calculate the global stiffness matrix
and residual vectors for the mesh, in terms of T, L, and q0 .
7.2.5.6. Show how the global stiffness matrix and q0
residual vectors must be modified to enforce
1 4
the constraints w1 = w4 = 0
7.2.5.7. Hence, calculate values of w at the two 2 3
intermediate nodes.
A B
L

Chapter 8
99

8. Theory and Implementation of the Finite Element Method

8.1. Generalized FEM for static linear elasticity

8.1.1. Consider a one-dimensional isoparametric quadratic element, illustrated in 3 2 1

the figure, and described in more detail in Section 8.1.5.
8.1.1.1. Suppose that the nodes have coordinates x11 = 0 , x12 = 1 , x13 = 2 . Using a parametric plot,
construct graphs showing the spatial variation of displacement in the element, assuming
that the nodal displacements are given by (a) u11 = 1, u12 = 0, u13 = 0 , (b)
u11 = 0, u12 = 1, u13 = 0 , (c) u11 = 0, u12 = 0, u13 = 1 , (d) u11 = 0, u12 = 0.5, u13 = 1
8.1.1.2. Suppose that the nodes have coordinates x11 = 0 , x12 = 1.75 , x13 = 2 . Plot graphs showing
the spatial variation of displacement in the element for each of the four sets of nodal
displacements given in 1.1.

8.1.2. Consider the finite element scheme to calculate displacements in an axially 1 3 2

loaded 1D bar, described in Section 8.15. Calculate an exact analytical
expression for the 3x3 stiffness matrix -1 x +1
x1
2 m (1 - n ) 2 N a ( x1 ) �
� N b ( x1 )
kab =
1 - 2n
h � �x1 �x1
dx1
x0
for a quadratic 1-D element illustrated in the figure.
Write a simple code to integrate the stiffness matrix x2 x2
numerically, using the procedure described in 8.1.5
and compare the result with the exact solution (for this 4 7
3
test, choose material and geometric parameters that
6
give 2 m (1 - n )h 2 /(1 - 2n ) = 1 ) Try integration
schemes with 1, 2 and 3 integration points.
4 x1 x1
1 5 2

8.1.3. Consider the mapped, 8 noded isoparametric element illustrated in the figure. Write a simple
program to plot a grid showing lines of x1 = constant and x 2 = constant in the mapped element, as
shown. Plot the grid with the following sets of nodal coordinates
 x1(1) = 0, x2(1) = 0 , x1(2) = 2.0, x2(2) = 2.0 ,
x1(3) = 2.0, x2(3) = 6.0 , x1(4) = 0.0, x2(4) = 4.0 ,
x1(5) = 1.0, x2(5) = 1.0 , x1(6) = 2.0, x2(6) = 4.0 ,
x1(7) = 1.0, x2(7) = 5.0 , x1(8) = 0.0, x2(8) = 2.0
 100

 x1(1) = 0, x2(1) = 0 , x1(2) = 2.0, x2(2) = 0. , x1(3) = 2.0, x2(3) = 2.0 , x1(4) = 0.0, x2(4) = 2.0
x1(5) = 1.5, x2(5) = 1.0 , x1(6) = 2.0, x2(6) = 1.0 , x1(7) = 1.0, x2(7) = 2.0 , x1(8) = 0.0, x2(8) = 1.0
Note that for the latter case, there is a region in the element where det � (
x j < 0 . This is
xi / � )
unphysical. Consequently, if elements with curved sides are used in a mesh, they must be designed
carefully to avoid this behavior. In addition, quadratic elements can perform poorly in large
displacement analyses.
x2
8.1.4. Set up an input file for the general 2D/3D linear elastic finite element code
6 5
7
provided to test the 6 noded triangular elements. Run the test shown in the
figure (dimensions and loading are in arbitrary units) and use a Young’s modulus 8 9 4 10 2
and Poissons ratio E = 100,n = 0.3 . Compare the FEA solution to the exact
solution. 1 2 3 x1
2
8.1.5. Set up an input file for the general 2D/3D linear elastic finite element
code provided to test the 4 noded tetrahedral elements. Run the test shown in 7 x3
the figure. Take the sides of the cube to have length 2 (arbitrary units), and 6
8
take E = 10,n = 0.3 . Run the following boundary conditions:
x1
8.1.5.1. u1 = u2 = u3 = 0 at node 1, u2 = u3 = 0 at node 2, and u3 = 0 at x2 5
4
nodes 3 and 4. The faces at x3 = 2 subjected to uniform traction 2
3
t3 = 2 (arbitrary units) 1

8.1.5.2. u1 = u2 = u3 = 0 at node 1, u2 = u1 = 0 at node 5, and u2 = 0 at nodes 3 and 6. The face at

x2 = 2 subjected to uniform traction t2 = 2 (arbitrary units)
In each case compare the finite element solution with the exact solution (they should be
equal)

8.1.6. Add lines to the 2D/3D linear elastic finite

element code to compute the determinant and 4 3 4 3 4 3
the eigenvalues of the global stiffness matrix.
Calculate the determinant and eigenvalues for a 1 2 1 2
mesh containing a single 4 noded quadrilateral
element, with each of the boundary conditions 1 2
shown in the figure. Briefly discuss the
implications of the results on the nature of solutions to the finite element equations.

8.1.7. Set up an input file for the mesh shown in the figure. Use material 1 1
properties E = 100,n = 0.3 and assume plane strain deformation. Run the
following tests 1
8.1.7.1. Calculate the determinant of the global stiffness matrix 1
8.1.7.2. Change the code so that the element stiffness matrix is
computed using only a single integration point. Calculate 1 the
determinant and eigenvalues of the global stiffness matrix.

8.1.8. Extend the 2D/3D linear elastic finite element code provided to solve problems involving anistropic
elastic solids with cubic symmetry (see Section 3.2.16 for the constitutive law). This will require the
following steps:
8.1.8.1. The elastic constants E ,n , m for the cubic crystal must be read from the input file
101

8.1.8.2. The orientation of the crystal must be read from the input
file. The orientation of the crystal can be specified by the x3
components of vectors parallel to the [100] and [010] 6 5
crystallographic directions. x1
4
8.1.8.3. The parts of the code that compute the element stiffness x2
matrix and stress (for post-processing) will need to be 2
modified to use elastic constants Cijkl for a cubic crystal. 3
1
The calculation is complicated by the fact that the
components of Cijkl must be expressed in the global coordinate system, instead of a
coordinate system aligned with the crystallographic directions. You will need to use the
unit vectors in 6.2 to calculate the transformation matrix Qij for the basis change, and use
the basis change formulas in Section 3.2.11 to calculate Cijkl .
8.1.8.4. Test your code by using it to compute the stresses and strains in a uniaxial tensile specimen
made from a cubic crystal, in which the unit vectors parallel to [100] and [010] directions
have components n[100] = cos q e1 + sin q e2 n[010] = - sin q e1 + cos q e 2 . Mesh the
specimen with a single cubic 8 noded brick element, with side length 0.01m. Apply
displacement boundary conditions to one face, and traction boundary conditions on another
corresponding to uniaxial tension of 100MPa parallel to the e1 axis. Run the following
tests: (i) Verify that if E ,n , m are given values that represent an istropic material, the
stresses and strains the element are independent of q . (ii) Use values for E ,n , m
representing copper (see Section 3.1.17). Define an apparent axial Young’s modulus for the
specimen as E (q ) = e11 / s 11 . Plot a graph showing
E (q ) as a function of q .
1
8.1.9. Extend the 2D/3D linear elastic finite element code to solve
problems involving body forces. This will require the following 5
steps
8.1.9.1. You will need to read a list of elements subjected to body forces, and the body force vector
for each element. The list can be added to the end of the input file.
8.1.9.2. You will need to write a procedure to calculate the contribution from an individual element
to the global system of finite element equations. This means evaluating the integral
bi N a ( x)dV
�
Ve(l )
over the volume of the element. The integral should be evaluated using numerical
quadrature – the procedure (other than the integrand) is essentially identical to computing
the stiffness matrix.
8.1.9.3. You will need to add the contribution from each element to the global force vector. It is
simplest to do this by modifying the procedure called globaltraction().
8.1.9.4. Test your code by using it to calculate the stress distribution in a 1D bar which is
constrained as shown in the figure, and subjected to a uniform body force. Use 4 noded
plane strain quadrilateral elements, and set E = 10,n = 0.25 (arbitrary units) and take the
body force to have magnitude 5. Compare the displacements and stresses predicted by the
finite element computation with the exact solution.

8.1.10. Implement the linear 3D wedge-shaped element shown in the figure the 2D/3D linear elastic finite
element code. To construct the shape functions for the element, use the shape functions for a linear
 102

triangle to write down the variation with (x1 , x 2 ) , and use a linear variation with x3 . Assume that
0 �xi �1 . You will need to add the shape functions and their derivatives to the appropriate procedures
in the code. In addition, you will need to modify the procedures that compute the traction vectors
associated with pressures acting on the element faces. Test your code by meshing a cube with two
wedge-shaped elements as shown in the figure. Take the sides of the cube to have length 2 (arbitrary
units), and take E = 10,n = 0.3 . Run the following boundary conditions:
8.1.10.1. u1 = u2 = u3 = 0 at node 1, u2 = u3 = 0 at node 2, and u3 = 0 at nodes 3 and 4. The faces
at x3 = 2 subjected to uniform traction t3 = 2 (arbitrary units) 7 x3
8.1.10.2. u1 = u2 = u3 = 0 at node 1, u2 = u1 = 0 at node 5, and u2 = 0 at 8
6
nodes 3 and 6. The face at x2 = 2 subjected to uniform traction 5 x1
4
t2 = 2 (arbitrary units) x2
In each case compare the finite element solution with the exact 2
3
solution (they should be equal). 1

8.2. The Finite Element Method for dynamic linear elasticity

8.2.1. Set up the general purpose 2D/3D dynamic finite element code
to solve the 1D wave propagation problem shown in the figure. x 1
Use material properties E = 10,n = 0.25, r = 1 (arbitrary units).
Mesh the bar with square 4 noded quadrilateral plane strain
20
elements. Assume that the left hand end of the bar is subjected to a
constant traction with magnitude p=5 at time t=0.
8.2.1.1. Calculate the maximum time step for a stable explicit computation.
8.2.1.2. Run an explicit dynamics computation (Newmark parameters b1 = 1/ 2, b 2 = 0 ) with a time
step equal to half the theoretical stable limit. Use the simulation to plot the time variation
of displacement, velocity and acceleration at x=0. Use at least 160 time steps. Compare
the numerical solution with the exact solution.
8.2.1.3. Repeat 1.2 with a time step equal to twice the theoretical stable limit.
8.2.1.4. Repeat 1.2 and 1.3 with Newmark parameters b1 = b 2 = 1/ 2 .
8.2.1.5. Repeat 1.2 and 1.3 with Newmark parameters b1 = b 2 = 1 . Try a time step equal to ten
times the theoretical stable limit.
8.2.2. The 2D/3D dynamic finite element code uses the row-sum method to compute lumped mass
matrices.
8.2.2.1. Use the row-sum method to compute lumped mass matrices for linear and quadratic
triangular elements with sides of unit length and mass density r = 1 (you can simply print
out the matrix from the code)
8.2.2.2. Use the row-sum method to compute lumped mass matrices for linear and quadratic
quadrilateral elements with sides of unit length and mass density r = 1
8.2.2.3. Modify the finite element code to compute the lumped mass matrix using the scaled
diagonal method instead. Repeat 2.1 and 2.2.
8.2.2.4. Repeat problem 1.2 (using Newmark parameters b1 = 1/ 2, b 2 = 0 ) using the two versions
of the lumped mass matrices, using quadratic quadrilateral elements to mesh the solid.
Compare the results with the two versions of the mass matrix.

L
x w(x)
A B
m
T
103

8.2.3. In this problem you will implement a simple 1D finite element method to calculate the motion of a
stretched vibrating string. The string has mass per unit length m, and is stretched by applying a tension
T at one end. Assume small deflections. The equation of motion for the string (see Section 10.3.1) is
d 2w d 2w
T =m
dx 2 dt 2
and the transverse deflection must satisfy w = 0 at x = 0, x = L .
8.2.3.1. Weak form of the equation of motion. Let d w be a kinematically admissible variation
of the deflection, satisfying d w = 0 at x = 0, x = L . Show that if
L L
d 2w dw dd w
�dt 2 d wdx + �
m T
dx dx
dx = 0
0 0
for all admissible d w , then w satisfies the equation of motion.
8.2.3.2. Finite element equations Introduce a linear finite
1 2 3 4 5 6 7
element interpolation as illustrated in the figure.
Calculate expressions for the element mass and
A B
stiffness matrices (you can calculate analytical s
expressions for the matrices, or evaluate the wa
integrals numerically)
wb
(a)
8.2.3.3. Implementation: Write a simple code to compute (b)
and integrate the equations of motion using the
l
Newmark time integration procedure described in
Section 8.2.5.
8.2.3.4. Testing: Test your code by computing the motion of the string with different initial
conditions. Try the following cases:
 w( x, t = 0) = ( L /10)sin( np x / L) , n=1,2…. These are the mode shapes for string, so the
vibration should be harmonic. You can calculate the corresponding natural frequency
by substituting w = cos wt sin( np x / L) into the equation of motion and solving for w
�x / 5 x < L/2
 w( x, t = 0) = �
�( L - x) / 5 x > L/2
Use geometric and material parameters L=10, m=1, T=1. Run a series of tests to investigate (i)
the effects of mesh size; (ii) the effects of using a lumped or consistent mass matrix; (ii) the
effects of time step size; and (iii) the effects of the Newmark parameters b1 , b 2 on the accuracy
of the solution.

8.2.4. Modify the 1D code described in the preceding problem to calculate the natural frequencies and
mode shapes for the stretched string. Calculate the first 5 natural frequencies and mode shapes, and
compare the numerical results with the exact solution.

8.2.5. The figure below shows a bar with thermal conductivity k and heat capacity c. At time t=0 the bar
is at uniform temperature, T = 0 . The sides of the bar
are insulated to prevent heat loss; the left hand end is L
held at fixed temperature T = 0 , while a flux of heat q *
q*
is applied into to the right hand end of the bar. In x qb
addition, heat is generated inside the bar at rate qb per T=0
unit length.
 104

The temperature distribution in the bar is governed by the 1-D heat equation
dT d 2T
c =k + qb
dt dx 2
and the boundary condition at x = L is
dT
k = q*
dx
In this problem, you will (i) set up a finite element method to solve the heat equation; (ii) show that the
finite element stiffness and mass matrix for this problem are identical to the 1-D elasticity problem
solved in class; (iii) implement a simple stepping procedure to integrate the temperature distribution
with respect to time.
8.2.5.1. Variational statement of the heat equation. Let d T be an admissible variation of
temperature, which must be (a) differentiable and (b) must satisfy d T = 0 on the left hand
end of the bar. Begin by showing that if the temperature distribution satisfies the heat
equation and boundary condition listed above, then
L L L
dT dT d d T
�
c dT + �
k qbd Tdx - q*d T ( x = 0) = 0
dx - �
dt dx dx
0 0 0
where d T ( x = 0) denotes the value of d T at the left hand end of the bar. To show this,
use the procedure that was used to derive the principle of virtual work. First, integrate the
first integral on the left hand side by parts to turn the d d T / dx into d T , then use the
boundary conditions and governing equation to show that the variational statement is true.
(Of course we actually rely on the converse – if the variational statement is satisfied for all
d T then the field equations are satisfied)
8.2.5.2. Finite element equations Now, introduce a linear 1D finite element interpolation scheme
as described in Section 8.2.5. Show that the diffusion equation can be expressed in the
form
dT b
Cab + K abT b (t = 0) = f a (t = 0) 1 < a �N
dt
T1 = 0
where T b , b=1,2…N denotes the unknown nodal values of temperature, Cab is a heat
capacity matrix analogous to the mass matrix defined in 8.2.5, and K ab is a stiffness
matrix. Derive expressions for the element heat capacity matrix, and the element stiffness
matrix, in terms of relevant geometric and material properties.
8.2.5.3. Time integration scheme Finally, we must devise a way to integrate the temperature
distribution in the bar with respect to time. Following the usual FEA procedure, we will
use a simple Euler-type time stepping scheme. To start the time stepping, we note that T a
is known at t=0, and note that we can also compute the rate of change of temperature at
time t=0 by solving the FEM equations
dT b
Cab + K abT b (t = 0) = f a (t = 0) 1 < a �N
dt
T1 = 0
For a generic time step, our objective is to compute the temperature T a and its time
derivative T&a at time t + Dt . To compute the temperature at the end of the step, we write
105

dT a
T a (t + Dt ) �T a (t ) + Dt
dt
and estimate dT / dt based on values at the start and end of the time step as
dT dT (t ) dT (t + Dt )
= (1 - b ) +b
dt dt dt
where 0 �b �1 is an adjustable numerical parameter. The time derivative of temperature
at t + Dt must satisfy the FEA equations
dT b (t + Dt )
Cab + K abT b (t + Dt ) = f a (t + Dt ) 1 < a �N
dt
T 1 (t + Dt ) = 0
Hence, show that the time derivative of temperature at time t + Dt can be calculated by solving
dT b (t + Dt ) �b dT b (t ) �
[ ab
C + D t b K ab ] = - K T
�
ab � (t ) + (1 - b ) Dt �+ f a (t + Dt )
�
1 < a �N
dt � dt �
T 1 (t + Dt ) = 0
whereupon the temperature at t + Dt follows as
dT a (t ) dT a (t + Dt )
T a (t + Dt ) �T a (t ) + Dt (1 - b ) + Dt (1 - b )
dt dt
8.2.5.4. Implementation. Modify the 1D dynamic FEA code to calculate the temperature variation
in a 1D bar.
8.2.5.5. Testing: Test the code by setting the bar length to 5; x-sect area to 1; heat capacity=100,
thermal conductivity=50; set the heat generation in the interior to zero (bodyforce=0) and
set the heat input to the right hand end of the bar to 10 (traction=10), then try the following
tests:
 To check the code, set the number of elements to L=15; set the time step interval to dt=0.1and
the number of steps to nstep=1000. Also, set up the code to use a lumped mass matrix by
setting lumpedmass=true and explicit time stepping b = 0 . You should find that this case runs
quickly, and that the temperature in the bar gradually rises until it reaches the expected linear
distribution.
 Show that the explicit time stepping scheme is conditionally stable. Try running with dt=0.2
for just 50 or 100 steps.
 Show that the critical stable time step size reduces with element size. Try running with dt=0.1
with 20 elements in the bar.
 Try running with a fully implicit time stepping scheme with b = 1 , 15 elements and dt=0.2.
Use a full consistent mass matrix for this calculation. Show that you can take extremely large
time steps with b = 1 without instability (eg try dt=10 for 10 or 100 steps) – in fact the
algorithm is unconditionally stable for b = 1

8.3. Finite element method for hypoelastic materials

8.3.1. Set up an input file for the example hypoelastic finite element code
described in Section 8.3.9 to calculate the deformation and stress in a
hypoelastic pressurized cylinder deforming under plane strain
conditions. Use the mesh shown in the figure, with appropriate
 106

symmetry boundary conditions on x1 = 0 and x2 = 0 . Apply a pressure of 50 (arbitrary units) to the

internal bore of the cylinder and leave the exterior surface free of traction. Use the following material
properties: s 0 = 10, e 0 = 0.001, n = 2,n = 0.3 . Plot a graph showing the variation of the radial
displacement of the inner bore of the cylinder as a function of the internal pressure. HEALTH
WARNING: The example code uses fully integrated elements and will give very poor results for large
values of n.

8.3.2. Solve the following coupled nonlinear equations for x and y using Newton-Raphson iteration.

( ) ( )
1/ 3 1/ 3
x x2 + y 2 + 1 -5= 0 y x2 + y 2 + 1 +3=0

8.3.3. Consider the simple 1D bar shown in the figure. Assume

that the boundary conditions imposed on the bar are identical L
to those for the 1D elastic bar discussed in Section 8.1.5, so
that u1 ( x1 ) is the only nonzero displacement component in the t * b
1
bar. Assume that the material has a hypoelastic constitutive
equation, which has a nonlinear volumetric and deviatoric x1
response, so that the stress-strain relation has the form
2 eij e
s ij = Sij + s kk d ij / 3 Sij = f (e e ) s kk = f (e ) kk
3 ee e
where
1 2
eij = e ij - e kk d ij ee = eij eij e = e iie jj
3 3
� � 2 �
� � 1 + n2 � n x � 1 �
s0 - - � - x �e 0
� � 2 � n - 1�
� � ( n - 1) �n - 1 e 0 �
f ( x) = � �
� 1/ n
�s �x � x �e 0
�0� � e
�
�
� 0
8.3.3.1. Show that the virtual work principle can be reduced to
L L
�
� d v1
u1 �� �u � �f (�
� u /� x) � x1 > 0
u1 / �
� ��x1 �
s 11 � bd v1dx1 - t *d v1 (0) = 0
dx1 - � s11 � 1 �= � 1 1
��x1 ��x1 � �- f (�u1 / �
x1 ) �
u1 / �x1 < 0
0 0
8.3.3.2. Introduce a 1D finite element interpolation for u1 and d v1 following the procedure
outlined in Section 8.1.5 to obtain a nonlinear system of equations for nodal values of u1a .
Show that the nonlinear equations can be solved using a Newton-Raphson procedure, by
repeatedly solving a system of linear equations
K ab dw1b + R a - F a = 0
where dw1b denotes a correction to the current approximation to u1a . Give expressions for
K ab and R a .
8.3.3.3. Extend the simple 1D linear elastic finite element code described in Section 8.1.5 to solve
the hypoelasticity problem described here.
107

8.3.3.4. Test your code by (i) calculating a numerical solution with material properties s 0 = 5 , n=2,
e 0 = 0.1 and loading b=0, t1* = 10 . and (ii) calculating a numerical solution with material
properties s 0 = 5 , n=10, e 0 = 0.1 and loading b=10, t1* = 0 . Compare the numerical
values of stress and displacement in the bar with the exact solution.

s ij / �
8.3.4. Calculate the tangent moduli � e kl for the hypoelastic material described in the preceding
problem.

8.3.5. In this problem you will develop a finite element code to solve dynamic problems involving a
hypoelastic material with the constitutive model given in Section . Dynamic problems for nonlinear
materials are nearly always solved using explicit Newmark time integration, which is very
straightforward to implement. As usual, the method is based on the virtual work principle
�2ui d vi
�
��t 2 d vi dV + �
r s ij [e kl ] dV - � bid vi dV - � ti*d vi dA = 0
�xj
R R R �R2

ui = ui* on �
1R
8.3.5.1. By introducing a finite element interpolation, show that the virtual work principle can be
reduced to a system of equations of the form
( M abu&&ib + Ria - Fia ) = 0
and give expressions for M ab , Ria , Fia .
8.3.5.2. These equations of motion can be integrated using an explicit Newmark method, using the
following expressions for the acceleration, velocity and displacement at the end of a generic
time-step
Dt 2 a
uia (t + Dt ) �uia (t ) + Dtu&ia (t ) + u&
&i (t )
2
a -1 � b a
u&
&i (t + Dt ) = M ab � - Ri [ui (t + Dt )] + Fib (t ) �
�
a a a a
u&i (t + Dt ) �u&i (t ) + Dt � (1 - b )u&& (t ) + b1u&
& �
i (t + Dt ) �
� 1 i
A lumped mass matrix should be used to speed up computations. Note that the residual
force vector Ria is a function of the displacement field in the solid. It therefore varies with
time, and must be re-computed at each time step. Note that this also means that you must
apply appropriate constraints to nodes with prescribed accelerations at each step.
Implement this algorithm by combining appropriate routines from the static hypoelastic
code and the Newmark elastodynamic code
provided. x2
8.3.5.3. Test your code by simulating the behavior of a 1D
(plane strain) shown in the figure. Assume that the
x1
bar is at rest and stress free at t=0, and is then
subjected to a constant horizontal traction at x1 = 10 for t>0. Fix the displacements for the
node at x1 = x2 = 0 and apply u1 = 0 at x1 = 0, x2 = 1 . Take the magnitude of the traction to
be 2 (arbitrary units) and use material properties s 0 = 1, e 0 = 0.01, n = 2,n = 0.3, r = 10 .
Take b1 = 0.5 in the Newmark integration, and use 240 time steps with step size 0.01 units.
Plot a graph showing the displacement of the bar at x1 = 10 as a function of time. Would
 108

you expect the vibration frequency of the end of the bar to increase or decrease with n?
Test your intuition by running a few simulations.

8.4. Finite element method for large deformations: hyperelastic materials

8.4.1. Write an input file for the demonstration hyperelastic finite element code described in Section 8.4.7
to calculate the stress in a Neo-Hookean tensile specimen subjected to uniaxial tensile stress. It is
sufficient to model the specimen using a single 8 noded brick element. Use the code to plot a graph
showing the nominal uniaxial stress as a function of the stretch ratio l = l / l0 . Compare the finite
element solution with the exact solution.

8.4.2. Extend the hyperelastic finite element code described in Section 8.4.7 to solve problems involving a
Mooney-Rivlin material. This will require the following steps:
e
8.4.2.1. Calculate the tangent stiffness Cijkl for the Mooney-Rivlin material
8.4.2.2. Modify the procedure called Kirchoffstress(…), which calculates the Kirchoff stress in the
material in terms of Bij , and modify the procedure called materialstiffness(…), which
computes the corresponding tangent stiffness. Run the following tests on the code, using
an 8 noded brick element to mesh the specimen:
 Subject the specimen to a prescribed change in volume, and calculate the corresponding stress
in the element. Compare the FE solution with the exact solution.
 Subject the specimen to a prescribed uniaxial tensile stress, and compare the FE solution with
the exact solution
 Subject the specimen to a prescribed biaxial tensile stress, and compare the FE solution to the
exact solution.

8.4.3. Modify the hyperelastic finite element code described in Section 8.4.7 to apply a prescribed true
traction to the element faces. To do this, you will need to modify the procedure that calculates the
element distributed load vector, and you will need to write a new routine to compute the additional term
in the stiffness discussed in Section 8.4.6. Test your code by repeating problem 1, but plot a graph
showing true stress as a function of stretch ratio.

8.5. Finite element method for inelastic materials

8.5.1. Set up the demonstration viscoplastic finite element code
e3
described in Section 8.5.7 to calculate the stress-strain relation for
the viscoplastic material under uniaxial tension. Mesh the specimen 5
8
with a single 8 noded brick element, using the mesh shown in the
6
7
1
4
2 e2
e1 3
109

figure. Apply the following boundary constraints to the specimen: u1 = u2 = u3 = 0 at node 1;

u2 = u3 = 0 at node 2, u3 = 0 at nodes 3 and 4. Apply a uniform traction whose magnitude increases
from 0 to 20 (arbitrary units) in time of 2 units on face 2 of the element.
8.5.1.1. Run a simulation with the following material parameters:
E = 10000, n = 0.3
Y = 15, e 0 = 0.5 n = 10, e&0 = 0.1 m = 10
and plot a graph showing the variation of traction with displacement on the element
8.5.1.2. Modify the boundary conditions so that only the constraint u3 = 0 is enforced at node 2.
Run the code to attempt to find a solution (you will have to abort the calculation). Explain
why the Newton iterations do not converge.
8.5.1.3. Repeat 1.1 with m = 200 , the correct boundary conditions, and with a maximum traction of
18 units. In this limit the material is essentially rate independent. Compare the predicted
traction-displacement curve with the rate independent limit.

8.5.2. Modify the viscoplastic finite element program to apply a constant (nominal) uniaxial strain rate to
the specimen described in the preceding problem, by imposing an appropriate history of displacement
on the nodes of the mesh. Test the code by plotting a graph showing uniaxial stress-v-strain in the
specimen, for material parameters E = 10000, n = 0.3 Y = 15, e 0 = 100. n = 100, e&0 = 0.1 m = 4
(in this limit the material is essentially a power-law creeping solid with constant flow stress) with an
applied strain rate of e&33 = 0.1 . Compare the numerical solution with the exact solution.

8.5.3. Modify the viscoplastic finite element program so that instead of using a power-law function to
represent the variation of flow-stress s 0 (e e ) with accumulated plastic
strain e e , the flow stress is computed by interpolating between a user- s
defined series of points, as indicated in the figure (the flow stress is
constant if the plastic strain exceeds the last point). Test your code by
using it to calculate the stress-strain relation for the viscoplastic
material under uniaxial tension. Use the mesh described in Problem 1, ee
and use material parameters E = 10000, n = 0.3 e&0 = 0.1 m = 150
(this makes the material essentially rigid and rate independent, so the stress-strain curve should follow
the user-supplied data points).

8.5.4. Modify the viscoplastic finite element code described in Section 8.5.7 to solve problems involving a
rate independent, power-law isotropic hardening elastic-plastic solid, with incremental stress-strain
relations
De ij = De ije + De ijp
1 +n � n � 3 Sij
De ije = �Ds ij - Ds kk d ij � D e ijp = De e
E � 1 + n � 2 se
3 2 p p
Sij = s ij - s kk d ij / 3 s e = Sij Sij De e = De De
2 3 ij ij
and a yield criterion
1/ n
� ee �
s e - Y0 �1+ � = 0
� e0 �
Your solution should include the following steps:
 110

( n +1)
8.5.4.1. Devise a method for calculating the stress s ij at the end of a load increment. Use a
fully implicit computation, in which the yield criterion is exactly satisfied at the end of the
load increment. Your derivation should follow closely the procedure described in Section
8.5.4, except that the relationship between s e( n+1) and De e must be calculated using the
yield criterion, and you need to add a step to check for elastic unloading.
8.5.4.2. Calculate the tangent stiffness � s ij( n+1) / �
De kl for the rate independent solid, by
differentiating the result of 4.1.
8.5.4.3. Implement the results of 4.1 and 4.2 in the viscoplastic finite element code.
8.5.4.4. Test your code by using it to calculate the stress-strain relation for the viscoplastic material
under uniaxial tension. Use the mesh, loading and boundary conditions described in
Problem 1, and use material properties E = 10000, n = 0.3 Y = 18, e 0 = 0.5 n = 10

8.5.5. Modify the viscoplastic finite element code described in Section 8.5.7 to solve problems involving a
rate independent, linear kinematic hardening elastic-plastic solid, with incremental stress-strain
relations
De&ij = De&ije + De&ijp
1 +n � n � 3 Sij -  ij
De ije = �Ds ij - 1 + n Ds kk d ij � D e ijp = De e
E � � 2 Y
3 2 p p
Sij = s ij - s kk d ij / 3 s e = ( Sij -  ij )( Sij -  ij ) De e = De De
2 3 ij ij
and a yield criterion and hardening law

se -Y = 0 D ij = cDe e
( Sij - ij )
Y
Your solution should include the following steps:
( n +1)
8.5.5.1. Devise a method for calculating the stress s ij at the end of a load increment. Use a
fully implicit computation, in which the yield criterion is exactly satisfied at the end of the
load increment. Your derivation should follow closely the procedure described in Section
8.5.4, except that the relationship between s e( n+1) and De e must be calculated using the
yield criterion and hardening law, and you need to add a step to check for elastic unloading.
8.5.5.2. Calculate the tangent stiffness � s ij( n+1) / �
De kl for the rate independent solid, by
differentiating the result of 8.5.5.1.
8.5.5.3. Implement the results of 8.5.5.1 and 8.5.5.2 in the viscoplastic finite element code.
8.5.5.4. Test your code by using it to calculate the stress-strain relation for the viscoplastic material
under uniaxial tension. Use the mesh, loading and boundary conditions described in
Problem 1, and use material properties E = 10000, n = 0.3 Y = 18, c = 100

8.5.6. In this problem you will develop a finite element code to solve dynamic problems involving
viscoplastic materials. Dynamic problems for nonlinear materials are nearly always solved using
explicit Newmark time integration, which is very straightforward to implement. As usual, the method is
based on the virtual work principle
111

�2ui d vi
�
��t 2 d vi dV + �
r s ij [e kl ] dV - � bid vi dV - � ti*d vi dA = 0
�xj
R R R �R 2

ui = ui* on �
1R
8.5.6.1. By introducing a finite element interpolation, show that the virtual work principle can be
reduced to a system of equations of the form
( M abu&&ib + Ria - Fia ) = 0
and give expressions for M ab , Ria , Fia .
8.5.6.2. To implement the finite element method, it is necessary to calculate the stress s ij in the
solid. Idealize the solid as a viscoplastic material with constitutive equations described in
Section 8.5.1. Since very small time-steps must be used in an explicit dynamic calculation,
it is sufficient to integrate the constitutive equations with respect to time using an explicit
method, in which the plastic strain rate is computed based on the stress at the start of a time
( n +1)
increment. Show that the stress s ij at time t + Dt can be expressed in terms of the
( n)
stress s ij at time t, the increment in total strain De ij during the time interval Dt and
material properties as
� m ( n) �
E � �s e( n ) � 3 Sij � E
( n +1) (n)
s ij = s ij + �De ij - Dte&0 exp(-Q / kT ) � �
�+ De kk
1 +n � �s ( n ) � 2 ( n)
s � 3(1 - 2n )
� �0 � e
�
8.5.6.3. The equations of motion can be integrated using an explicit Newmark method using the
following expressions for the acceleration, velocity and displacement at the end of a generic
time-step
Dt 2 a
uia (t + Dt ) �uia (t ) + Dtu&ia (t ) + u&
&i (t )
2
a -1 � b a
u&
&i (t + Dt ) = M ab � - Ri [ui (t + Dt )] + Fib (t ) �
�
u&ia (t + Dt ) �u&ia (t ) + Dt � &a (t ) + b1u&
(1 - b )u& & a �
i (t + Dt ) �
� 1 i
A lumped mass matrix should be used to speed up computations. Note that the residual
force vector Ria is a function of the displacement field in the solid. It therefore varies with
time, and must be re-computed at each time step. Note that this also means that you must
apply appropriate constraints to nodes with prescribed accelerations at each step.
Implement this algorithm by combining appropriate routines from the static viscoplastic
code and the Newmark elastodynamic code
provided. x2
8.5.6.4. Test your code by simulating the behavior of a 1D
(plane strain) shown in the figure. Assume that the
x1
bar is at rest and stress free at t=0, and is then
subjected to a constant horizontal traction at x1 = 10 for t>0. Fix the displacements for the
node at x1 = x2 = 0 and apply u1 = 0 at x1 = 0, x2 = 1 . Take the magnitude of the traction to
be 2 (arbitrary units) and use material properties s 0 = 1, e 0 = 0.01, n = 2,n = 0.3, r = 10 .
Take b1 = 0.5 in the Newmark integration, and use 240 time steps with step size 0.01 units.
Plot a graph showing the displacement of the bar at x1 = 10 as a function of time.
 112

8.6. Advanced element formulations – incompatible modes; reduced

integration; and hybrid elements
8.6.1. Volumetric locking can be a serious problem in computations involving
nonlinear materials. In this problem, you will demonstrate, and correct,
locking in a finite element simulation of a pressurized hypoelastic cylinder.
8.6.1.1. Set up an input file for the example hypoelastic finite element
code described in Section 8.3.9 to calculate the deformation and
stress in a hypoelastic pressurized cylinder deforming under
plane strain conditions. Use the mesh shown in the figure,
with appropriate symmetry boundary conditions on x1 = 0 and
x2 = 0 . Apply a pressure of 20 (arbitrary units) to the internal
bore of the cylinder and leave the exterior surface free of
traction. Use the following material properties: s 0 = 10, e 0 = 0.001, n = 2,n = 0.3 . Plot a
graph of the variation of the radial displacement of the inner bore of the cylinder as a
function of the applied pressure. Make a note of the displacement at the maximum
pressure.
8.6.1.2. Edit the code to reduce the number of integration points used to compute the element
stiffness matrix from 9 to 4. (modify the procedure called `numberofintegrationpoints’).
Repeat the calculation in 8.6.1.1. Note the substantial discrepancy between the results of
8.6.1.1 and 8.6.1.2 – this is caused by locking. The solution in 5.2, which uses reduced
integration, is the more accurate of the two. Note also that using reduced integration
improves the rate of convergence of the Newton-Raphson iterations.

8.6.2. Modify the hypoelastic finite element code described in Section 8.3.9
to use selective reduced integration. Check your code by (a) repeating the
calculation described in problem 8.6.1; and (b) running a computation
with a mesh consisting of 4 noded quadrilateral elements, as shown in the
figure. In each case, calculate the variation of the internal radius of the
cylinder with the applied pressure, and plot the deformed mesh at
maximum pressure to check for hourglassing. Compare the solution
obtained using selective reduced integration with the

8.6.3. Run the simple demonstration of the B-bar method described in Section 8.6.2 to verify that the
method can be used to solve problems involving near-incompressible materials. Check the code with
both linear and quadratic quadrilateral elements.

8.6.4. Extend the B-bar method described in Section 8.6.2 to solve problems involving hypoelastic
materials subjected to small strains. This will require the following steps:
8.6.4.1. The virtual work principle for the nonlinear material must be expressed in terms of the
modified strain measures e ij and de ij defined in Section 8.6.2. This results in a system of
nonlinear equations of the form
113

*
�
s ij [e ij (uk )]de ij dV - �
bid vi dV - �ti d vi dA = 0
R R �2 R

ui = ui* on �1R
which must be solved using the Newton-Raphson method. Show that the Newton-Raphson
procedure involves repeatedly solving the following system of linear equations for
corrections to the displacement field dwkb
K aibk dwkb = - Ria + Fia
s pq
� a b
K aibk = ��e kl B pji Bqlj dV Ria = s kj �
� e (wb ) �B a dV Fia = bi N a dV +
� �
ti N* a
dA
�kl i � kji
R R R �
R
a
where Bijk is defined in Section 8.6.2.
8.6.4.2. Modify the hypoelastic code provided to compute the new form of the stiffness matrix and
element residual. You will find that much of the new code can simply be copied from the
small strain linear elastic code with the B-bar method
8.6.4.3. Test the code by solving the problems described in 8.6.1 and 8.6.2.

8.6.5. Extend the B-bar method described in Section 8.6.2 to solve problems involving hypoelastic
materials subjected to small strains, following the procedure outlined in the preceding problem.

8.6.6. In this problem you will extend the B-bar method to solve problems involving finite deformations,
using the hyperelasticity problem described in Section 8.4 as a representative example. The first step is
to compute new expressions for the residual vector and the stiffness matrix in the finite element
approximation to the field equations. To this end
 New variables are introduced to characterize the volume change, and the rate of volume
change in the element. Define
1 1 1
h= det(F ) dV � h&= JF ji-1F&
ij dV = � JLkk dV �
Vel Vel Vel
Vel Vel Vel
Here, the integral is taken over the volume of the element in the reference configuration.
1/ n
 The deformation gradient is replaced by an approximation Fij = Fij ( h / J ) , where n=2
for a 2D problem and n=3 for a 3D problem, while J=det(F).
 The virtual velocity gradient is replaced by the approximation
d Lij = d Lij + d ij (dh&/h - d Lkk ) / n
The virtual work equation is replaced by
*
� �
t ij [ Fkl ]d Lij dV0 - r 0bid vi dV0 - �ti d vih dA0 = 0
V0 V0 �2V0

8.6.6.1. Verify that Lij = Fik F&kj-1

8.6.6.2. In calculations to follow it will be necessary to calculate fij = �
h /�
Fij . Find an expression
for fij .
8.6.6.3. The virtual work equation must be solved for the unknown nodal displacements by
Newton-Raphson iteration. Show that, as usual, the Newton-Raphson procedure involves
repeatedly solving the following system of linear equations for corrections to the
displacement field dwkb
 114

K aibk dwkb = - Ria + Fia

and derive expressions for Ria and Kiakb .

Chapter 9

Modeling Material Failure

9.
115

9.1. Summary of mechanisms of fracture and fatigue under static and cyclic
loading

9.1.1. Summarize the main differences between a ductile and a brittle material. List a few examples of
each.

9.1.2. What is the difference between a static fatigue failure and a cyclic fatigue failure?

9.1.3. Explain what is meant by a plastic instability, and explain the role of plastic instability in causing
failure

9.1.4. Explain the difference between `High cycle fatigue’ and `Low cycle fatigue’

9.1.5. Summarize the main features and mechanisms of material failure under cyclic loading. List
variables that may influence fatigue life.

9.2. Stress and strain based fracture and fatigue criteria

9.2.1. A flat specimen of glass with fracture strength s TS , Young’s modulus E and Poisson’s ratio n is
indented by a hard metal sphere with radius R, Young’s modulus E2 and Poisson’s ratio n 2 . Using
solutions for contact stress fields given in Chapter 5, calculate a formula for the load P that will cause
the glass to fracture, in terms of geometric and material parameters. You can assume that the critical
stress occurs on the surface of the glass.

9.2.2. The figure shows a fiber reinforced composite laminate.

e2 e1
(i) When loaded in uniaxial tension parallel to the fibers, it fails at a stress
of 500MPa.
(ii) When loaded in uniaxial tension transverse to the fibers, it fails at a
stress of 250 MPa.
(iii) When loaded at 45 degrees to the fibers, it fails at a stress of 223.6
MPa
The laminate is then loaded in uniaxial tension at 30 degrees to the fibers.
Calculate the expected failure stress under this loading, assuming that the
material can be characterized using the Tsai-Hill failure criterion.

9.2.3. A number of cylindrical specimens of a brittle material with a 1cm radius and length 4cm are tested
in uniaxial tension. It is found 60% of the specimens withstand a 150MPa stress without failure; while
30% withstand a 170 MPa stress without failure.
9.2.3.1. Calculate values for the Weibull parameters s 0 and m for the specimens
9.2.3.2. Suppose that a second set of specimens is made from the same material, with length 8cm
and radius 1cm. Calculate the stress level that will cause 50% of these specimens to fail.
 116

9.2.4. A beam with length L, and rectangular cross-section b �h is made from a brittle material with
Young’s modulus E, Poisson’s ratio n , and the failure probability distribution of a volume V0 is
characterized by Weibull parameters s 0 and m .
9.2.4.1. Suppose that the beam is loaded in uniaxial tension parallel to its length. Calculate the
stress level s 0 corresponding to 63% failure, in terms of geometric and material
parameters.
9.2.4.2. Suppose that the beam is loaded in 3 point bending. Let s R = 3PL /(bh 2 ) denote the
maximum value of stress in the beam (predicted by beam theory). Find an expression for
the stress distribution in the beam in terms of s R
9.2.4.3. Hence, find an expression for the value of s R that corresponds to 63% probability of
failure in the beam. Calculate the ratio s R / s 0 .

9.2.5. A glass shelf with length L and rectangular cross-section b �h is used to display cakes in a bakery.
As a result, it subjected to a daily cycle of load (which may be approximated as a uniform pressure
acting on it surface) of the form p(t ) = p0 (1 - t / T ) where 0 < t < T is the time the store has been open,
and T is the total time the store is open each day. As received, the shelf has a tensile strength s TS 0 ,
and the glass can be characterized by static fatigue parameters  and m . Find an expression for the
life of the shelf, in terms of relevant parameters.

9.2.6. A cylindrical concrete column with radius R, cross-sectional radius R, and length L is subjected to a
monotonically increasing compressive axial load P. Assume that the material can be idealized using the
constitutive law given in Section 9.2.4, with the compressive yield stress-v-plastic strain of the form
1/ m
p
�e p �
Y (e ) = s 0 � �
�e 0 �
� �
where s 0 , e 0 and m are material properties. Assume small strains, and a homogeneous state of stress
and strain in the column. Neglect elastic deformation, for simplicity.
9.2.6.1. Calculate the relationship between the axial stress P / A and strain d / L , in terms of the
plastic properties c, s 0 , e 0 and m
9.2.6.2. Calculate the volume change of the column, in terms of P / A , c, s 0 , e 0 and m
9.2.6.3. Suppose that the sides of the column are subjected to a uniform traction q. Repeat the
calculations in parts 9.2.6.1 and 9.2.6.2.

9.2.7. Suppose that the column described in the previous problem is encased in a steel tube, with (small)
wall thickness t. The steel can be idealized as a rigid perfectly plastic material with yield stress Y .
Calculate the relationship between the axial stress P / A and strain d / L , in terms of geometric and
material properties.

9.2.8. Extend the viscoplastic finite element program described in Section 8.5 to model the behavior of a
porous plastic material with constitutive equations given in Section 9.2.5. This will involve the
following steps:
( n +1) ( n +1)
9.2.8.1. Develop a procedure to calculate the stress s ij , the void volume fraction V f , the
effective strain measures e m( n +1) , e e( n+1) at the end of the time increment, given their values
at the start of the increment and given an increment in plastic strain De ij . You should use a
117

fully implicit update, as discussed in Section 8.5. The simplest approach is to set up, and
( n +1)
solve, three simultaneous nonlinear equations for ( V f e m( n+1) , e e(n+1) ) using Newton-
Raphson iteration, and subsequently compute the stress distribution.
9.2.8.2. s ij( n+1) / �
Calculate the tangent stiffness � De ij for the material
9.2.8.3. Implement the new constitutive equations in the viscoplastic finite element program
9.2.8.4. Test your code by simulating the behavior of a uniaxial tensile specimen subjected to
monotonic loading.

9.2.9. A specimen of steel has a yield stress of 500MPa. Under cyclic loading at a stress amplitude of 200
MPa it is found to fail after 104 cycles, while at a stress amplitude of 100MPa it fails after 105 cycles.
This material is to be used to fabricate a plate, with thickness h, containing circular holes with radius
a<<h. The plate will be subjected to constant amplitude cyclic uniaxial stress far from the holes, and
must have a life of at least 105 cycles. What is the maximum stress amplitude that the plate can
withstand?

9.2.10. A spherical pressure vessel with internal radius a and external radius b=1.5a is repeatedly
pressurized from zero internal pressure to a maximum value p . The sphere has yield stress Y, and its
fatigue behavior of the vessel (under fully reversed uniaxial tension) can be characterized by Basquin’s
law s a N b = C .
9.2.10.1. Find an expression for the fatigue life of the vessel in terms of p , and relevant geometric
and material properties. Assume that the effects of mean stress can be approximated using
Goodman’s rule. Assume that p / Y < 2(1 - a3 / b3 ) / 3
9.2.10.2. Suppose that the vessel is first pressurized to its collapse load and then unloaded, so as to
induce a distribution of residual stress in the cylinder. It is subsequently subjected to a
cyclic pressure with magnitude p . Calculate the mean stress and stress amplitude as a
function of position in the vessel wall, and hence deduce an expression for its fatigue life.

9.2.11. The figure shows a solder joint on a printed circuit board. The
printed circuit board can be idealized as a pinned-pinned beam with e3 L/2
thickness h, length L, Young’s modulus E and mass density r . The
board vibrates in its fundamental mode with a frequency
e1 h
w = (p h / L) E /12 r and mode shape u3 = A sin(p x1 / L) . The yield
stress of solder is so low it can be neglected, and it is firmly bonded L
to the printed circuit board. As a result, it is subjected to a cyclic
plastic strain equal to the strain at the surface of the beam. The fatigue life of solder can be
characterized by a Coffin-Manson law De p N b = C . Find an expression for the time to failure of the
solder joint, in terms of relevant geometric and material parameters.

9.2.12. A specimen of steel is tested under cyclic loading. It is found to have a fatigue threshold
s 0 = 75MPa , and fails after 103 cycles when tested at a stress amplitude s a = 1.5s 0 . Suppose that,
in service, the material spends 80% of its life subjected to stress amplitudes s a < s 0 , 10% of its life at
s a = 1.1s 0 , and the remainder at s a = 1.2s 0 . Calculate the life of the component during service
(assume that the mean stress s m = 0 during both testing and service).
 118

9.3. Modeling failure by crack growth – linear elastic fracture mechanics

9.3.1. Using the equations for the crack tip fields find an expression for the maximum shear stress
distribution around the tip of a plane-strain Mode I crack. Hence, plot approximate contours of
successive yield zones

9.3.2. Briefly describe the way in which the concept of stress intensity factor can be used as a fracture
criterion.

9.3.3. A welded plate with fracture toughness K IC contains a residual stress

distribution
� � 2b
a a �b �
r 22 = s R � - arcsinh � �� 2a
� x2 + a2 b �
�a �
� 1 �
along the line of the weld. A crack with length 2a lies on the weld line. The r22
solid is subjected to a uniaxial tensile stress s 22 = s 0 . Find an expression for
the critical value of s 0 that will cause the weld to fracture, in terms of K IC ,
s R and a.

9.3.4. Hard, polycrystalline materials such as ceramics often contain a distribution of inter-granular
residual stress. The objective of this problem is to estimate the influence of this stress distribution on
crack propagation through the material. Assume that
 The solid has mode I fracture toughness K IC
 As a rough estimate, the residual stress distribution can be idealized as s 22 = s R sin(p x1 / L) ,
where L is of the order of the grain size of the solid and s R is the magnitude of the stress.
 A long (semi-infinite) crack propagates through the solid – at some time t,the crack tip is located at
x1 = c
 The solid is subjected to a remote stress, which induces a mode I stress intensity factor K I� at the
crack tip
9.3.4.1. If the solid is free of residual stress, what value of K I� that causes fracture.
9.3.4.2. Calculate the stress intensity factor induced by the residual stress distribution, as a function
of c.
9.3.4.3. What value of K I� is necessary to cause crack propagation through the residual stress
field? What is the maximum value of K I�, and for what crack tip position c does it occur?

e2
d

e1
119

9.3.5. A dislocation, with burgers vector b = b1e1 + b2 e2 and line direction e3 lies a distance d ahead of a
semi-infinite crack. Calculate the crack tip stress intensity factors.

9.3.6. The figure shows a simple model that is used to

estimate the size of the plastic zone at a crack tip. The s
crack, with length 2a, together with the plastic zones with
x2
length rp , are considered together to be a crack with Y Y
x1
length 2( a + rp ) . The solid is loaded by uniform stress
s at infinity. The region with length rp near each crack
rp a a rp
tip is subjected to traction Y acting to close the crack.
Using the solutions in Section 9.3.3, calculate an s
expression for the Mode I crack tip stress intensity factor.
rp �ps �
Show that K I = 0 if = sec � �- 1
a �2Y �

9.3.7. Suppose that an ASTM compact tension specimen is used to measure

the fracture toughness of a steel. The specimen has dimensions W = 40mm
and B = 20 mm. The crack length was 18.5 mm, and the fracture load was 2rp
15kN.
a
B
9.3.7.1. Calculate the fracture toughness of the steel.
9.3.7.2. If the steel has yield stress 800MPa, was this a valid
measurement?
W

9.3.8. Find expressions for the Mode I and II stress intensity factors for the angled
crack shown. If  = 450 , what is the initial direction of crack propagation? 2a

Confirm your prediction experimentally, using a center-cracked specimen of paper.

9.3.9. A large solid contains a crack with initial length 2c0 . The solid has plane-strain
fracture toughness K IC , and under cyclic loading the crack growth rate obeys Paris
Ds
law
da n
= C ( DK I )
dN
9.3.9.1. Suppose that the material is subjected to a cyclic uniaxial stress with
amplitude Ds and mean stress Ds (so the stress varies between 0 and 2c0
 120

2 Ds ). Calculate the critical crack length that will cause fracture, in terms of K IC and
Ds
9.3.9.2. Calculate an expression for the number of cycles of loading that are necessary to cause a
crack to grow from an initial length 2c0 to fracture under the loading described in 7.1
9.3.9.3. Show that the number of cycles to failure can be expressed in the form of Basquin’s law
(discussed in Section 9.2.7) as N ( Ds ) b = D , where b and D are constants. Give
expressions for b and D in terms of the initial crack length, the fracture toughness, and the
material properties in Paris law

9.4. Energy methods in fracture mechanics

M
9.4.1. The figure shows a double-cantilever beam fracture specimen
that is loaded by applying moments to the ends of the beams. h
Define the compliance of the solid as C = 2q / M 2q
9.4.1.1. Derive an expression relating the crack tip energy hB
a
release rate to compliance and dimensions of the
specimen. M
9.4.1.2. Use the result of 9.4.1.1 to calculate the crack tip stress intensity factors for the specimen in
terms of M and relevant geometric and material properties

h P
9.4.2. The figure shows a thin film with thickness h, Film
n
Youngs modulus E and Poisson’s ratio . To test the
interface between film and substrate, a delaminated 2a
region of the film with width 2a is subjected to a
vertical force (per unit out of plane distance) P. By Substrate
idealizing the delaminated region as a beam, and using
the method of compliance, estimate the crack tip energy release rate in terms of P and relevant material
and geometric properties.

9.4.3. The figure shows a thin rectangular strip of material

with height 2h and out-of-plane thickness B. The D
material can be idealized as a linear elastic solid with L h
Young’s modulus E and Poisson’s ratio n . The strip is
damaged by an array of widely-spaced cracks with
2a h
length 2a and spacing L>>a. It is loaded by a uniform
tensile traction t acting on the top and bottom surface,
which induces a displacement D .
9.4.3.1. Write down an expression for the compliance D / t of the undamaged strip (i.e. with no
cracks)
9.4.3.2. Write down a relationship between the compliance of the strip and the crack tip energy
release rate.
9.4.3.3. Estimate the crack tip energy release rate using the energy release rate for an isolated crack
in an infinite solid. Hence, find an expression for the compliance of the cracked strip, in
terms of relevant geometric and material parameters.
121

9.4.4. The figure shows a double cantilever beam specimen that is P

loaded by forces applied to the ends of the beams. Evaluate the
J integral around the path shown to calculate the crack tip
energy release rate. You can use elementary beam theory to
b h
estimate the strain energy density, stress, and displacement in
a hB
the two cantilevers. The solution must, of course, be
independent of b. P

9.4.5. Use the J integral, together with the solution for the stress and displacement field near the tip of a
crack given in Section 9.3.1, to calculate the relationship between the crack tip energy release rate and
the stress intensity factor for a Mode I crack.

9.4.6. The figure shows a thin film with thickness h, thermal

expansion coefficient  f , Young’s modulus E and Poisson’s Film h
ratio n on a large substrate with thermal expansion coefficient
 s . The film is initially perfectly bonded to the substrate and
stress free. The system is then heated, inducing a thermal Substrate
stress in the film. As a result, the film delaminates from the
substrate, as shown in the figure.
9.4.6.1. Calculate the state of stress a distance d>>h ahead
of the advancing crack tip
9.4.6.2. Assume that the film is stress free a distance d>>h behind the crack tip. By directly
calculating the change in energy of the system as the crack advances, find an expression for
the crack tip energy release rate
9.4.6.3. Check your answer to 9.4.5.2 by evaluating the J integral around the path indicated in the
figure.

9.5. Plastic fracture mechanics

9.5.1. Explain briefly the main concepts underlying the use of the J integral as a fracture criterion in
components experiencing large-scale plastic deformation.

9.5.2. The figure shows the tip of a semi-infinite crack in an

elastic-plastic material with a bi-linear uniaxial stress-
strain curve, as indicated in the figure. To provide some s e2
H e1
E
e
 122

insight into the nature of the crack tip fields, the constitutive behavior can be approximated as a
hypoelastic material, characterized by a strain energy density W such that s ij = � W /� e ij . Suppose that
the solid is subjected to remote mode I loading (so that the shear stresses s12 = s13 = 0 on x2 = 0 ).
9.5.2.1. Construct the full stress-strain equations for the hypoelastic material, using the approach
described in Section 3.3
9.5.2.2. Consider a material point that is very far from the crack tip, and so is subjected to a very
low stress. Write down the asymptotic stress field in this region, in terms of an arbitrary
constant K I� that characterizes the magnitude of the remote mode I loading
9.5.2.3. Consider a material point that is very close to the crack tip, and so is subjected to a very
large stress. Write down the asymptotic stress field in this region, in terms of an arbitrary
constant K Itip that characterizes the magnitude of the near tip stresses.
9.5.2.4. Using the path independence of the J integral, find a relationship between K I�, K Itip , and
the slopes E , H of the uniaxial stress-strain curve
9.5.2.5. Suppose that the material fractures when the stress at a small distance d ahead of the crack
tip reaches a critical magnitude s 0 . Assume that the critical distance is much smaller than
the region of high stress considered in 2.4. Calculate the critical value of K I� that will
cause the crack to grow, in terms of relevant material parameters.
9.5.2.6. Consider a finite sized crack with length a in the hypoelastic material. Assume that the
solid is subjected to a remote uniaxial stress far from the crack. Discuss qualitatively how
the stress field around the crack evolves as the remote stress is increased. Discuss the
implications of this behavior on the validity of the fracture criterion derived in 9.5.2.5.

9.6. Linear Elastic Fracture mechanics of interfaces

9.6.1. Calculate values for the elastic constants  , b and the crack tip singularity parameter e for the
following bi-material interfaces:
(a) Aluminum on glass (b) A glass fiber in a PVC matrix
(c) Nickel on titanium carbide (d) Copper on Silicon
123

9.6.2. A center-cracked bi-material specimen is made by bonding Al to

SiO 2 . It contains a crack with length 10mm, and is loaded to failure in x2
uniaxial tension. It is found to fail at a stress level s 22 = 200 MPa.
9.6.2.1. Calculate the fracture toughness and the corresponding x1
phase angle of loading, using a characteristic length
L = 100 m m a a
9.6.2.2. Calculate the fracture phase angle if L = 1mm is chosen as
the characteristic length.
9.6.2.3. Calculate the distance behind the crack tip where (according
to the asymptotic solution) the crack faces first overlap

9.6.3. A bi-material interface is made by bonding two materials together. The material above the
interface has shear modulus and Poisson’s ratio m1 ,n1 ; the material below the crack has shear modulus
and Poisson’s ratio m2 ,n 2 Due to roughness, a residual stress distribution
s 22 ( x1 ) = s 0 sin(2p x1 / l )
acts on the bi-material interface. Suppose that the interface contains a long (semi-infinite) crack, with
crack tip located at x1 = a . Calculate the crack tip stress intensity factors as a function of the elastic
mismatch parameter e and other relevant parameters. Deduce expressions for the energy release rate
and the phase angle.

Chapter 10

10.
Approximate theories for solids with special shapes:
rods, beams, membranes, plates and shells
 124

10.1. Preliminaries: Dyadic notation for vectors and tensors

10.1.1. Let {e1 , e2 , e3} be a Cartesian basis. Express the identity tensor as a dyadic product of the basis
vectors

10.1.2. {e1 , e2 , e3} and {m1, m 2 , m 3} be two Cartesian bases. Show that the tensor R = mi �ei can be
visualized as a rigid rotation (you can show that R is an orthogonal tensor, for example, or calculate the
change in length of a vector that is multiplied by R).

10.1.3. Let a and b be two distinct vectors (satisfying a‫ ׹‬b 0 ). Let S = a �b . Find an expression for all
u=0
the vectors u that satisfy S �

10.1.4. Find the eigenvalues and eigenvectors of the tensor S = a �b in terms of a, b, and their magnitudes
(Don’t forget to find three independent eigenvectors).

10.1.5. Let S = a �b . Find the condition on a and b necessary to ensure that S is orthogonal.

10.1.6. Let ai be three linearly independent vectors. Define ai to be three vectors that satisfy
�1 i= j
ai �
ai = �
�0 i �j
a j , g ij = ai �
and let gij = ai � a j and g ij = ai �
a j denote the 27 possible dot products of these vectors.
10.1.6.1. Find expressions for ai in terms of vector and scalar products of ai
�
10.1.6.2. Let S = Sij ai �a j be a general second order tensor. Find expressions for S ij , S�
i i
j, Sj
satisfying
�
S = S ij ai �a j = S�
i j i j
j ai �a = S j a �ai

in terms of Sij and gij , g ij and g ij

10.1.6.3. ( i
)( j
) i
(
Calculate ai �a �a j �a . What does the tensor ai �a represent? )
10.1.6.4. ( i
) ( )
Express ai �a in terms of ai �a j and appropriate combinations of gij , g ij and g ij

Express ( ai �a ) in terms of ( a �a ) and appropriate combinations of gij , g ij and g ij

i i i
10.1.6.5.
10.1.6.6. Let F denote a homogeneous deformation gradient, satisfying bi = F � ai . Express F in
terms of dyadic products of ai and bi .
10.1.6.7. Find an expression for FT F in terms of scalar products of bi and dyadic products of ai ,
i.e. find components U ij satisfying FT F = U ij ai �a j
10.1.6.8. Find an expression for the Lagrange strain tensor E = (FT F - I ) / 2 in terms of dyadic
products of ai , i.e. find Eij satisfying (FT F - I ) / 2 = Eij ai �a j , in terms of scalar
products of bi and appropriate combinations of gij , g ij and g ij
125

10.2. Motion and Deformation of slender rods

10.2.1. The figure shows an inextensible rod that is bent into a helical shape.
The shape of the helix can be characterized by the radius r of the
generating cylinder, and the number of turns n in the helix per unit axial
length. Consider a point on the axis of the rod specified by the polar
coordinates (r ,qˆ, z ) .
10.2.1.1. Write down an expression for z in terms of r, n and q .
10.2.1.2. Write down the position vector of the point as components in
the {e1 , e2 , e3} basis.
10.2.1.3. Calculate an expression for the unit vector m3 that is tangent
to the rod, in terms of the basis vectors {e1 , e 2 , e3} and
appropriate coordinates.
10.2.1.4. Assume that m 2 is perpendicular to the axis of the cylinder.
Use this and the solution to 1.3 to find expressions for the basis vectors m1 , m 2 in terms of
{e1 , e2 , e3} .
10.2.1.5. Calculate the normal and binormal vectors to the curve and hence deduce an expression for
the torsion of the curve.
10.2.1.6. Deduce an expression for the curvature vector of the rod.
10.2.1.7. Suppose that the stress state in the deformed rod is a simple axial distribution s m3 �m3 .
Calculate the stress components in the {e1 , e2 , e3} basis.

10.2.2. An initially straight, inextensible slender bar with length L and

circular cross-section with radius a is bent into a circle with radius R by e3 m3
terminal couples, as shown in the figure. Assume that cross-sections of e1
the rod remain circles with radius a and remain transverse to the axis of m1
the rod after deformation.
10.2.2.1. Write down an expression for the position vector r ( x3 ) of a
R
position vector on the axis of the deformed rod, expressing x3
your answer as components in the { e1, e 2 , e3 } basis
10.2.2.2. Find an expression for the basis vectors { m1, m 2 , m 3 } as a
function of arc-length s, expressing each unit vector as components in { e1, e2 , e3 } . Hence
find an expression for the orthogonal tensor R that maps { e1 , e2 , e3 } onto { m1, m 2 , m3 } .
10.2.2.3. Write down the curvature vector for the deformed rod, and verify that
dm1 dm 2 dm3
= -k 2m3 + k 3m 2 = k1m3 - k 3m1 = -k1m 2 + k 2m1
ds ds ds
10.2.2.4. Write down expressions for the deformation gradient in the rod, expressing your answer as
both components in { e1, e2 , e3 } and { m1, m 2 , m3 }
10.2.2.5. Find an expression for the Lagrange strain tensor in the rod, expressing your answer as both
components in { e1, e 2 , e3 } and { m1, m 2 , m3 } . Neglect second-order terms.
10.2.2.6. Hence deduce expressions for the Material stress and Cauchy stress in the rod.
10.2.2.7. Calculate the resultant internal moment and force acting on a generic internal cross-section
of the rod.
 126

10.2.2.8. Show that the internal moment satisfies the equations of equilibrium.

10.2.3. Consider a deformable rod, as shown in the figure. Let x3

denote the arc-length of a point on the axis un-deformed rod from
some arbitrary origin, and let y denote the twist of the rod, as
defined in Section 10.2.2. In addition, let r ( x3 ) and
v = dr / dt = vi mi denote the position vector and velocity of this
point on the deformed rod, let s denote its arc-length after
deformation, and let κ = k i mi denote the curvature vector of the
rod, where mi are basis vectors aligned with the deformed rod as
discussed in Section 10.2.2 Show that
10.2.3.1. The time derivative of a unit vector tangent to the rod
can be computed as
dt dm3 d �dr � dr � dr d �dr � � �dv1 � �dv �
� = � �- � � � � � = � - v2k 3 + v3k 2 � m1 + � 2 + v1k 3 - v3k1 �
m2
dt dt ds �dt � ds �ds ds �dt � � �ds � �ds �
10.2.3.2. The time derivative of the the rate of stretching of the rod’s centerline is related to its
velocity by
d �ds � ds d �dr � ds �dv3 �
� �= m3 � � �= � - v1k 2 + v2k1 �
dt �dx3 � dx3 ds �dt � dx3 �ds �
10.2.3.3. The angular velocity of the basis vectors can be calculated as
dv �dv � �dv �
ω = m3 � + y&m3 = - � 2 + v1k 3 - v3k1 � m1 + � 1 - v2k3 + v3k 2 �
m 2 + y&m3
ds �ds � �ds �
10.2.3.4. The angular acceleration of the basis vectors can be calculated as
2
dω da �dv � dv dy �dv �dv � � d y
α +�-+‫״‬m -�3 == 2 � m m
3� 3 � � m m
3 � 3� m3
dt ds �ds � ds dt �ds �ds � dt 2

10.3. Simplified versions of the general theory of deformable rods

10.3.1. Consider a flexible, inextensible cable subjected to

b e1
transverse loading p = - pe1 (e.g. due to gravity) as illustrated in
the figure. e3
-y1
10.3.1.1. Express the basis vector e1 in terms of m1 , m 3 and m1
s y3
(derivatives of) y1 , y3 .
10.3.1.2. Find an expression for the curvature vector for the
q
cable in terms of (derivatives of) y1 , y3 m3
127

10.3.1.3. Find the two equilibrium equations relating the axial tension T3 to the external loading and
geometry of the cable, by substituting M1 = M 2 = M 3 = T1 = T2 = 0 into the general
equations of motion for a rod.

10.3.2. Consider a long, straight rod, with axis parallel to e3 , which is subjected to pure
twisting moments Q = Qe3 acting at its ends. The rod may be idealized as a linear
elastic solid with shear modulus m . The deformation of the rod can be characterized
by the twist y ( x3 ) and the transverse displacement of the cross-section u3 ( x ) .
Assume that the only nonzero internal moment component is M 3 , and the nonzero
internal stress components are s 3 . Simplify the general governing equations for a
deformable rod to obtain:
10.3.2.1. A simplified expression for the curvature tensor κ for the deformed rod,
in terms of y ( x3 )
10.3.2.2. Equations of equilibrium and boundary conditions for M 3 and s 3
10.3.2.3. Expressions relating s 3 to y ( x3 ) and the warping function w. Show that the equilibrium
equation for s 3 reduces to the governing equation for the warping function given in
Section 10.2.10.
10.3.2.4. Expressions relating M 3 to y ( x3 ) .

10.3.3. An initially straight beam, with axis parallel to m1 p1(x3)

the e3 direction and principal axes of inertia (0 ) m3 p3(x3) P (L )
P3 3
parallel to e1 , e2 is subjected to a force per unit
u1(x3) e
length p = p1e1 + p3e3 . The beam has Young’s 1
r
modulus E and mass density , and its cross- x3
e3
section has area A and principal moments of inertia
I1 , I 2 , I 3 . Assume that a large axial internal force Nm3 is developed in the beam, either by a horizontal
force per unit length p3 or horizontal forces P3(0) e3 , P3( L )e3 acting at the ends of the beam. Suppose
that the beam experiences a finite transverse displacement u = u1e1 , so that the stretch of the beam and
its curvature must be approximated by
2 2 2
ds � � u3 � �� u1 � �u3 1 �� u1 � κ �d u1 e
�� 1+ � 1 + � � �1 + + � � 2
dx3 � � x3 � �� x3 � �x3 2 �� x3 � dx32
Show that the static equilibrium equations for the displacement components can be reduced to
�� 2 2
d 4u1 u 1 ��u ��d 2u du � 1 ��u ��
EI 2 + EA � 3 + � 1 �� 1 - p3 1 + p1 � 1 - � 1 ��= 0
dx34 �� x 2 ��x3 ��dx32 dx3 � 2 ��x3 ��
� 3 � � �
2 2
d ��� u3 1 �� u �� � 1 �� u ��
EA + � 1 ��+ p3 � 1 - � 1 ��= 0
dx3 �� x 2 �� x3 �� � 2 �� x3 ��
� 3 � � �
and list the boundary conditions on the ends of
the beam. e1
e3
p1(x3)
A B
u1(x3)
T0 T0
x3
 128

10.3.4. The goal of this problem is to derive the equation of motion for an inextensible stretched string
subjected to small displacements by a direct application of the principle of virtual work. Assume that at
some instant the string has transverse deflection u1 ( x3 ) and velocity v1 ( x3 ) as indicated in the figure.
10.3.4.1. Write down an equation for the curvature of the string, accurate to first order in u1 ( x3 )
10.3.4.2. Write down an expression for the relative velocity of the end B of the string relative to A, in
terms of u1 ( x3 ) and v1 ( x3 ) .
10.3.4.3. Write down the rate of virtual work done by the transverse forces p1 ( x3 ) in terms of a
virtual velocity d v1 ( x3 )
10.3.4.4. Write down the rate of virtual done by the applied tension T0 .
10.3.4.5. Hence use the principle of virtual work to derive the equation of motion for the string.

10.4. Exact solutions to problems involving slender rods

10.4.1. A slender, linear elastic rod has shear modulus m and an elliptical cross-section,
as illustrated in the figure. It is subjected to equal and opposite axial couples with
magnitude Q on its ends. Using the general theory of slender rods, and assuming that
the rod remains straight:
10.4.1.1. Write down the internal force and moment distribution in the rod
10.4.1.2. Calculate the twist per unit length of the shaft y
10.4.1.3. Find an expression for the displacement field in the shaft
10.4.1.4. Find an expression for the stress distribution in the shaft
10.4.1.5. Find an expression for the critical couple Q that will cause the shaft to
yield

10.4.2. A slender, linear elastic rod has shear modulus m and an equilateral triangular
cross-section, as illustrated in the figure. It is subjected to equal and opposite axial
couples with magnitude Q on its ends. Using the general theory of slender rods, and
assuming that the rod remains straight:
10.4.2.1. Write down the internal force and moment distribution in the rod
10.4.2.2. Calculate the twist per unit length of the shaft y
10.4.2.3. Find an expression for the displacement field in the shaft
10.4.2.4. Find an expression for the stress distribution in the shaft
10.4.2.5. Find an expression for the critical couple Q that will cause the shaft to
yield

10.4.3. The figure shows a flexible cable, subjected to a

transverse force per unit length p = pi mi and forces
p
m3
P (0) = Pi(0) mi and P ( L ) = Pi( L )mi acting at its ends. In a
s m1 e1
e3
P(0) P(L)
129

flexible cable, the area moments of inertia can be neglected, so that the internal moments M i �0 . In
addition, the axial tension
�
T3 = s 33dA
is the only nonzero internal force. Show that under these conditions the virtual work
A
equation reduces to
L0 L L
d d s&
( pid vi ) ds - �
r Aaid vi ds - �
�dx3 T3dx3 + � P (0)d vi � - �
�i �
Pi( L )d vi � = 0
x3 =0 � �
x3 = L
0 0 0
where d vi is the virtual velocity of the cable, and d s&is the corresponding rate of change of arc-length
along the cable. Show that if the virtual work equation is satisfied for all d vi and compatible d s&, the
internal force and curvature of the cable must satisfy

10.4.4. The figure shows a flexible string, which is e1

supported at both ends and subjected to a tensile L
e3
force T0 . The string is subjected to a uniform p
transverse for p per unit length. Calculate the
u1(x3)
deflection u1 ( x3 ) of the string, assuming small T0
T0
deflections.
x3
e1
10.4.5. The figure shows a flexible string with b
length L, which is pinned at both ends. The e3
string is subjected to a uniform transverse for p p
per unit length. Calculate the deflection u1 ( x3 ) u1(x3)
of the string, assuming small deflections. x3

10.4.6. The figure shows a flexible cable with length L and weight m
per unit length hanging between two supports under uniform vertical b e1
gravitational loading. In a flexible cable, the area moments of e3
inertia can be neglected, so that the internal moments M i �0 -y1
m1
10.4.6.1. Write down the curvature vector of the cable in terms of s y3
the angle q ( s) shown in the figure q
10.4.6.2. Hence, show that the equations of equilibrium for the m3
cable reduce to
dT3 dq
+ m sin q = 0 T3 + m cos q = 0
ds ds
10.4.6.3. Hence, show that T3 cosq = H , where H is a constant. Interpret the equation physically.
10.4.6.4. Deduce that w( y3 ) = tan q = -m( s - L / 2) / H � dw / dy3 = -m 1 + w2 / H
10.4.6.5. Hence, deduce that w( y3 ) = dy1 / dy3 = - sinh( my3 / H ) and calculate y1 as a function of
y3
10.4.6.6. Finally, calculate the internal forces in the cable.
 130

10.4.6.7. Show that, as L � b the full solution approaches the

small deflection solution calculated in Problem 5.
Find the value of L / b for which the discrepancy
between the full solution and the small deflection
solution is 10%.
b e1 b
10.4.7. The figure shows an e3
inextensible cable with h2
m1 h
weight per unit length m, and s 1

length 2L that is pinned at q

both ends. The cable is
m3
supported by a frictionless
pulley midway between the
two ends. Find all the possible equilibrium values of the sags h1 , h2 of the cable. Display your results
by plotting a graph showing the equilibrium values of h1 / b as a function of L / b . You will need to
solve problem 6 before attempting this one.

e1
10.4.8. The figure shows a flexible string, which is e3
supported at both ends and subjected to a tensile L
force T0 . The string has mass per unit length m and
can be approximated as inextensible. Calculate the u1(x3)
natural frequencies of vibration and the T0 T0
corresponding mode shapes, assuming small x3
transverse deflections.

10.4.9. Estimate the fundamental frequency of vibration for the stretched string described in the preceding
problem, using the Rayleigh-Ritz method. Use the approximation u1 ( x3 ) = Ax3 ( L - x3 ) for the mode
shape. Compare the estimate with the exact solution derived in problem 10.4.8.

10.4.10.The figure shows an Euler-Bernoulli beam with Young’s modulus E, area moments of inertia I1 , I 2
and length L, which is clamped at x3 = L and pinned at x3 = 0 . It is subjected to a uniform load p per
unit length. Calculate the internal moment and shear force in the beam, and calculate the transverse
deflection.

10.4.11. The figure shows an initially straight, inextensible elastic rod, with
Young’s modulus E, length L and principal in-plane moments of area
I1 = I 2 = I , which is subjected to end thrust. The ends of the rod are
constrained to travel along a line that is parallel to the undeformed rod, but
the ends are free to rotate. Use the small-deflection solution for beams
subjected to significant axial force given in Section 10.3.3 to calculate the
value of P required to hold the rod in equilibrium with a small nonzero
131

deflection, and find an expression for the deflected shape. Compare the predicted deflection with the
exact post-buckling solution given in Section 10.4.3.

10.4.12.The theory describing small-deflections of beams

subjected to significant axial force given in Section
10.3.3 can be extended to obtain an approximate large
deflection solution. Consider the beam shown in the
figure. The beam has Young’s modulus E, cross-
sectional area A, and principal transverse moments of
inertia I1 , I 2 . The bar is subjected to load per unit
length p = p1e1 + p3e3 , and axial forces P (0) = P3(0)e3 , P ( L ) = P3( L )e3 at its two ends. Assume that the
displacement field can be described as u = u1 ( x3 )e1 + u3 ( x3 )e3 . Deformation measures are to be
expanded up to second order in transverse deflection, so that
 The axial stretch can be approximated as
2 2
ds � � u � �� u � �u 1 ��
u �
1 + 3 � 1 + � 1 � �1 + 3 + � 1 �
��
dx3 � � x3 � �� x3 � �x3 2 ��
x3 �
d 2u1
 The curvature can be approximated as κ � e2
dx32
10.4.12.1. Show that the static equilibrium equations for the displacement components can be reduced
to
�� 2 2
d 4u1 � u3 1 �� u1 �� � d 2u1 du1 � 1 ��u1 ��
EI 2 + EA + � � - p3 + p1 1 - � ��= 0
�
dx34 �� x3 2 �� x3 ��dx32 dx3 � 2 ��x3 ��
� � � �
2 2
d ��� u3 1 �� u �� � 1 �� u ��
EA + � 1 ��+ p3 � 1 - � 1 ��= 0
dx3 �� x 2 �� x3 �� � 2 �� x3 ��
� 3 � � �
and list the boundary conditions on the ends of the beam.
10.4.12.2. Solve the governing equations for the beam problem described in Problem 14. Compare
the predicted deflection with the exact post-buckling solution given in Section 10.4.3, for
the limiting case of an inextensible beam.

10.4.13.An initially straight tent-pole with Young’s modulus E and 3h/2

hollow circular cross-section with external radius a, moment of
inertia I is to be bent into an arc with height h and base 3h / 2 as
shown in the figure. Calculate expressions for h
10.4.13.1. The force P required to bend the pole into shape
10.4.13.2. The total length of the pole P
10.4.13.3. The maximum stress in the pole
P
 132

10.4.14.An initially straight, elastic rod with Young’s modulus E,

area moments of inertia I1 = I 2 = I and axial effective inertia
J 3 is subjected to an axial couple Q = Qe3 , which remains P P
q
fixed in direction as the rod deforms.
10.4.14.1. Show that the straight rod, with an appropriate R  
twist is a possible equilibrium configuration for
all values of Q, and calculate the value of twist
10.4.14.2. Show that, for a critical value of Q, the rod may adopt a helical shape, with one complete
turn and arbitrary height h and radius r. Calculate the critical value of Q
10.4.14.3. What can you infer about the stability of a straight rod subjected to end couples?

10.4.15.The figure shows a rod, which is a circular arc with radius R in its stress free configuration, and is
subjected to load per unit length p = p1m1 + p3m 3 and forces P (0) , P ( L ) on its ends that cause a small
change in its shape. In this problem, we shall neglect out-of-plane deformation and twisting of the rod,
for simplicity. Let s = Rq denote the arc length measured along the undeformed rod, and let
u = u1 (q )m1 + u3 (q )m3 the displacement of the rod’s centerline.
10.4.15.1. Note that approximate expressions for the resulting (small) change in arc length and
curvature of the rod can be calculated using the time derivatives given in Section 10.2.3.
Hence, show that
d d s 1 �du3 �
 The derivative of the change in arc-length of the deformed rod is = � + u1 (q ) �
ds R �dq �
1 �d 2u1 du3 �
 The change in curvature vector is d κ = 2 � - �
m2
R � �dq
2 dq �
�
10.4.15.2. The geometric terms in the equilibrium equations listed in Section 10.2.9 can be
approximated using the geometry of the undeformed rod. Show that internal forces
T = T1m1 + T3m3 and internal moment M = M 2m 2 must satisfy the following static
equilibrium equations
dT1 dT3 dM 2
- T3 + Rp1 = 0 + T1 + Rp3 = 0 + RT1 = 0
dq dq dq
10.4.15.3. Assume that the rod is elastic, with Young’s modulus E and area moment of inertia
I1 = I 2 = I , and can be idealized as inextensible. Show that under these conditions the
axial displacement u3 must satisfy
EI �d 6u3 d 4u3 d 2u3 � dp1
� +2 + �+ + p3 = 0
R4 ��dq
6
dq 4 dq 2 �
� dq m3
and write down expressions for the boundary conditions at m1 p1
the ends of the rod. m 2
p 3
10.4.15.4. As a particular example, consider a rod which is a s
semicircular arc between q = �  , subjected to R
equal and opposite forces acting on its ends, as
shown in the figure. Assume that the displacement q
and rotation of the rod vanish at q = 0 , for
simplicity. Calculate u1 (q ) and u3 (q ) for the rod.
133

10.5. Motion and Deformation of thin shells

e3 eR
10.5.1. A spherical-polar coordinate system is to be used to
describe the deformation of a spherical shell with radius R. The P ef
two angles (f ,q ) illustrated in the figure are to be used as the
coordinate system (x1, x 2 ) for this geometry. eq
10.5.1.1. Write down the position vector r in terms of (f ,q ) O R
, expressing your answer as components in the basis
f
{e1 , e2 , e3} shown in the figure. e2
10.5.1.2. Calculate the covariant basis vectors ( m1 , m 2 , m3 ) e1
in terms of {e1 , e2 , e3} and (f ,q )
10.5.1.3. Calculate the contravariant basis vectors
( m1, m2 , m3 ) in terms of {e1 , e 2 , e3} and (f ,q )
10.5.1.4. Calculate the covariant, contravariant and mixed components of the metric tensor g
10.5.1.5. Calculate the covariant, contravariant and mixed components of the curvature tensor κ for
the shell
i
10.5.1.6. Find the components of the Christoffel symbol Gb for the coordinate system
10.5.1.7. Suppose that under loading the shell simply expands radially to a new radius r. Find the
components of the mid-plane Lagrange strain tensor
1
( )
γ = g b m �m b = gb - gb m �m b and the components of the curvature change
2
tensor Dκ = (k b - k b )gl m l �m b
 

10.5.1.8. Suppose that the shell is elastic, with Young’s modulus E and Poisson’s ratio n . Calculate
the contravariant components of the internal force T b and internal moment M b
10.5.1.9. Find the physical components of the internal force and moment, expressing your answer as
components in the spherical-polar basis of unit vectors {e R , eq , ef }

10.5.2. A cylindrical-polar coordinate system is to be used to describe e3 ez

the deformation of a cylindrical shell with radius r. The angles and
axial distance (q , z ) illustrated in the figure are to be used as the P eq
z r
coordinate system (x1, x 2 ) for this geometry. er
10.5.2.1. Write down the position vector r in terms of (q , z ) , O
expressing your answer as components in the basis
q
{e1 , e2 , e3} shown in the figure. e2
10.5.2.2. Calculate the covariant basis vectors ( m1 , m 2 , m3 ) in e1
terms of {e1 , e2 , e3} and (q , z )
10.5.2.3.
1 2
(
Calculate the contravariant basis vectors m , m , m
3
)
in terms of {e1 , e2 , e3} and (q , z )
10.5.2.4. Calculate the covariant, contravariant and mixed components of the metric tensor g
 134

10.5.2.5. Calculate the covariant, contravariant and mixed components of the curvature tensor κ for
the shell
i
10.5.2.6. Find the components of the Christoffel symbol Gb for the undeformed shell
10.5.2.7. Suppose that under loading the shell simply expands radially to a new radius r , without
axial stretch. Find the covariant components of the mid-plane Lagrange strain tensor γ
and the covariant components of the curvature change tensor Dκ
10.5.2.8. Suppose that the shell is elastic, with Young’s modulus E and Poisson’s ratio n . Calculate
the contravariant components of the internal force T b and internal moment M b .
10.5.2.9. Find the physical components of the internal force and moment, expressing your answer as
components in the spherical-polar basis of unit vectors {e r , eq , e z }

10.5.3. The figure illustrates a triangular plate, whose geometry can be described
by two vectors a and b parallel to two sides of the triangle. The position b
vector of a point is to be characterized using a coordinate system (x1, x 2 ) by
setting r = x1a + x 2b where 0 �x1 �1 , 0 �x 2 �1 .
10.5.3.1. Calculate the covariant basis vectors ( m1 , m 2 , m3 ) in terms of a a
and b
10.5.3.2.
1
( 2 3
)
Calculate the contravariant basis vectors m , m , m in terms of a and b
10.5.3.3. Calculate the covariant, contravariant and mixed components of the metric tensor g
10.5.3.4. Suppose that the plate is subjected to a homogeneous deformation, so that after deformation
its sides lie parallel to vectors c and d . Find the mid-plane Lagrange strain tensor
1
( )
γ = g b m �m b = gb - gb m �m b , in terms of a, b, c and d
2
10.5.3.5. Suppose that the plate is elastic, with Young’s modulus E and Poisson’s ratio n . Calculate
the contravariant components of the internal force T b

10.5.4. The figure illustrates a triangular plate. The position points in the plate e2
is to be characterized using the height z and angle q as the coordinate
system (x1, x 2 ) .
10.5.4.1. Calculate the covariant basis vectors ( m1 , m 2 , m3 ) q
expressing your answer as components in the basis z e1
{e1 , e2 , e3} shown in the figure.
10.5.4.2.
1
( 2
Calculate the contravariant basis vectors m , m , m
3
)
as components in {e1 , e2 , e3}
10.5.4.3. Calculate the covariant, contravariant and mixed components of the metric tensor g
i
10.5.4.4. Find the components of the Christoffel symbol Gb for the undeformed plate
10.5.4.5. Suppose that the plate is subjected to a deformation such that the position vector of a point
that lies at ( z ,q ) in the undeformed shell has position vector r = z tan(q +  )e1 + ze2 after
deformation Find the mid-plane Lagrange strain tensor γ = g b m �m b , in terms of
z ,q , 
10.5.4.6. Suppose that the plate is elastic, with Young’s modulus E and Poisson’s ratio n . Calculate
the contravariant components of the internal force T b
135

10.5.4.7. Find the physical components of the internal force T as components in the {e1 , e2 , e3}
basis.

10.6. Simplified versions of general shell theory – flat plates and membranes

10.6.1. Consider a shell that is so thin that the internal moments

M  all vanish. Find the simplified equations of motion for
the internal forces T b and the transverse force V  in terms
of relevant geometric parameters.

10.6.2. The figure shows a thin circular plate with thickness h, mass
density r , Young’s modulus E and Poisson’s ratio n that is
simply supported at its edge and is subjected to a pressure
distribution acting perpendicular to its surface. The goal of this
problem is to derive the equations governing the transverse
deflection of the plate in terms of the cylindrical-polar
coordinate (r ,q ) system shown in the figure.
10.6.2.1. Write down the position vector of a point on the
mid-plane of the undeformed plate in terms of (r ,q ) , expressing your answer as
components in the {e1 , e2 , e3} basis.
10.6.2.2. ( 1 2 3
)
Calculate the basis vectors ( m1 , m 2 , m3 ) and m , m , m , expressing your answer as
components in the basis {e1 , e2 , e3} shown in the figure.
i
10.6.2.3. Find the components of the Christoffel symbol Gb for the undeformed plate;
10.6.2.4. Calculate the contravariant components of the metric tensor g
10.6.2.5. Find the basis vectors ( m1 , m 2 , m 3 ) for the deformed plate, neglecting terms of order
( �u3 / �r ) 2 , ( �u3 / �q ) 2 , etc
10.6.2.6. Show that the curvature tensor has components
�2u3 1� u3 �2u3 �u �2u3
k11 = - k12 = k 21 =
- k 22 = -r 3 -
� r2 r � q �� r q � r �q2
10.6.2.7. Express the internal moments M b in the plate in terms of E ,n and u3 and its derivatives.
10.6.2.8. Write down the equations of motion for the plate in terms of M b and V 
10.6.2.9. Hence, show that the transverse displacement must satisfy the following governing
equation
 136

Eh3 ��2 1 � 1 �2 � ��2u3 1 �u3 1 �2u3 � �2u3

� + + �
� + + �+ rh = p3
12(1 - n 2 ) �
��r2 r � q2 �
r r2 � ��
�r
�
2 r � q2 �
r r2 � � �t2

10.7. Solutions to Problems Involving Membranes, Plates and Shells

10.7.1. A thin circular membrane with radius R is subjected to in-

plane boundary loading that induces a uniform biaxial
membrane tension with magnitude T0 . Vertical displacement
of the membrane is prevented at r = R . The membrane is
subjected to a uniform out-of-plane pressure with magnitude p
on its surface. Calculate the displacement field in the
membrane, assuming small deflections. This problem can be
solved quite easily using Cartesian coordinates.

10.7.2. The figure shows a thin circular plate with thickness h,

Young’s modulus E and Poisson’s ratio n . The edge of the
plate is clamped, and its surface is subjected to a uniform out-
of-plane pressure with magnitude p. Calculate the
displacement field and internal moment and shear force in the
plate, assuming small deflections. This problem can easily be
solved using Cartesian coordinates.

10.7.3. The figure shows a circular elastic plate with Young’s

modulus E, Poissions ratio n . The plate has thickness h and
radius R, and is clamped at its edge. The goal of this problem
is to calculate the mode shapes and natural frequencies of
vibration of the plate. The solution to Problem 10.6.2 may be
used.
10.7.3.1. Show that the general solution to the equation of
motion for the freely vibrating plate is given by
(
u3 (r ,q ) = AJ n (k( m,n) r )sin(nq + q0 ) + BYn (k( m,n) r )sin(nq + q1 )

)
+ CI n ( k( m,n ) r )sin(nq + q 2 ) + BK n (k(m,n) r )sin( nq + q3 ) cos(w( m,n) t + f )
where A,B,C,D, f q0 ...q3 are arbitrary constants, J n , Yn are Bessel functions of the first
and second kinds, and I n , K n are modified Bessel functions of the first and second kinds,
with order n, while k( m,n) and w( m,n) are a wave number and vibration frequency that are
related by
k(2m,n) = w(m,n ) 12(1 -n ) r / Eh 2
137

10.7.3.2. Show that most general solution with bounded displacements at r=0 has the form
u3 (r ,q ) = � ( sin(nq )
A1J n (k( m,n ) r ) + B1I n (k( m,n) r ) �
� �
+�
A2 J n (k( m,n) r ) + B2 I n (k( m,n) r ) �
� � )
cos(nq ) cos(w( m,n)t + f )
10.7.3.3. Write down the boundary conditions for u3 at r=R, and hence show that the wave numbers
k( m,n) are roots of the equation
J n (k( m,n) R ) I n +1 (k( m,n ) R) + I n (k( m,n) R ) J n +1 (k( m,n) R ) = 0
10.7.3.4. Show that the corresponding mode shapes are given by
U ( m,n ) ( r ,q ) = A � sin(nq + q1 )
I n ( k ( m , n ) R ) J n ( k( m , n ) r ) - J n ( k( m , n ) R ) I n ( k ( m , n ) r ) �
� �
10.7.3.5. Calculate k( m,n) for 0<n<3, 1<m<3 and tabulate your results. Compare the solution to the
corresponding membrane problem solved in Section 10.7.2.
10.7.3.6. Plot contours showing the mode shapes for the modes with the lowest four frequencies.

10.7.4. Repeat the preceding problem for a plate with a simply supported edge. You should find that the
wave numbers k( m,n) are the roots of the equation
(1 - n ) �
k( m,n) RJ n ( k( m,n ) R) I n +1 ( k(m,n) R ) + I n (k( m,n ) R ) J n +1 (k( m,n) R) �
� �

( )
2
- 2 k( m , n ) R I n ( k ( m, n ) R ) J n ( k ( m , n ) R ) = 0
Calculate the natural frequencies for n = 0.3 .

10.7.5. Use Rayleigh’s method, with a suitable guess for the mode shape, to estimate the fundamental
frequency of vibration of the circular plate described in Problem 10.7.3. Compare the approximate
solution with the exact solution.

10.7.6. Use Rayleigh’s method, with a suitable guess for the mode shape, to estimate the fundamental
frequency of vibration of the circular plate described in Problem 10.7.4. Compare the approximate
solution with the exact solution.

10.7.7. A thin film with thickness h f is deposited on the surface

of a circular wafer with radius R and thickness h. Both film
and substrate have Young’s modulus E, and Poisson’s ratio n .
p p
An inelastic strain e11 = e 22 = e 0 is introduced into the film by
some external process (e.g. thermal expansion), which
generates stresses in the film, and also causes the substrate to
bend. In Section 10.7.3, expressions were derived relating the substrate curvature to the mismatch
strain in the film. These expressions are only
valid if the substrate curvature is small. For
mismatch strains, the wafer buckles, as shown
 138

in the figure. The goal of this problem is to estimate the critical value of mismatch strain that will cause
the wafer to buckle.
10.7.7.1. Assume that the displacement of the mid-plane of the wafer-film system is
u3 = k1x12 / 2 + k 2 x22 / 2 u1 = A1x1 + A2 x13 + A3 x1x22 u2 = B1x1 + B2 x23 + B3 x2 x12
Calculate the distribution of strain in the film and substrate, and hence deduce the total
strain energy density of the system. It is best to do this calculation using a symbolic
manipulation program.
10.7.7.2. Calculate the values of Ai , Bi that minimize the potential energy of the plate, and hence
show that the two curvatures satisfy
k1k 22 (1 - n 2 ) R 4 + 16h 2k1 + 16 h 2nk 2 - Ce 0 = 0
k 2k12 (1 - n 2 ) R 4 + 16h2k 2 + 16h 2nk1 - Ce 0 = 0
and find a formula for the constant C.
10.7.7.3. Hence, plot a graph showing the equilibrium values of the normalized curvature
k1R 2 1 + n /(4h) as a function of a suitably normalized measure of mismatch strain, for a
Poisson’s ratio n = 0.3 .

10.7.8. The figure shows an elastic plate with Young’s modulus E,

Poisson’s ratio n and thermal expansion coefficient  . The plate is
circular, with radius R and its edge is clamped. The plate is initially
stress free and is then heated to raise its temperature by DT , inducing
a uniform internal force Tb = - EhDT db /(1 - n ) in the plate. The
goal of this problem is to calculate the critical temperature that will
cause the plate to buckle, using the governing equations listed in
Section 10.6.2.
10.7.8.1. Assume an axially symmetric buckling mode, so that u3 = w(r ) with r = x x . Show
that w satisfies the governing equation
Eh3 �d 4 w 2 d 3 w 1 d 2 w 1 dw � EhDT ��2 w 1 � w�
� + - + � + � + �= 0
12(1 -n 2 ) �
�dr 4 r dr 3
r 2
dr 2
r 3 dr � (1 - n ) � 2
� � � r r �
r �
�
10.7.8.2. Show that the general solution to this equation is
w( r ) = A + B log( r ) + CJ 0 ( kn r ) + DY0 (kn r )
where A,B,C,D, are arbitrary constants, J 0 , Y0 are Bessel functions of the first and second
kind of order zero, and kn is the wave number for the nth buckling mode. Find an
expression for kn in terms of DT and relevant geometric and material properties.
10.7.8.3. Show that the wave number kn satisfies J1 (kn R ) = 0 , and deduce an expression for the
critical temperature change at which the plate can first buckle.
10.7.8.4. Plot the buckling mode associated with this temperature.

10.7.9. Repeat the preceding problem for a plate with a simply supported edge. You should find that the
wave numbers kn are the roots of the equation
(1 - n ) J1 (kn R) - kn RJ 0 (kn R ) = 0
Calculate the critical buckling temperature and plot the corresponding buckling mode for n = 0.3 .
139

10.7.10.The interface between a thin film and its substrate contains

a tunnel crack with length 2a. The film is initially stress free,
then heated to raise its temperature by DT , inducing a
uniform biaxial stress field s b = - EDT db /(1 - n ) in the
film. At a critical temperature, the film buckles as shown in
the figure. For R>>h the buckled film can be modeled as a
plate with clamped edge. The goal of this problem is to
calculate the critical temperature required to cause the film to buckle, and to calculate the crack tip
energy release rate as the buckled film delaminates from the substrate. Assume that the displacement of
the mid-plane of the film has the form u = u ( x2 )e 2 + w( x2 )e3 , and use the Von-Karman plate bending
theory of Section 10.6.3.
10.7.10.1. Write down expressions for the mid-plane strain g b and the curvature change tensor
Dkb for the film in terms of u and w , and hence find a formula for the internal force
Tb and moment M b in terms of u , w , DT and material properties.
10.7.10.2. Hence, show that the equations of equilibrium in terms of u and w reduce to
� 2
dT22 d 4 w 12(1 -n 2 ) d 2w Eh � du 1 �dw ��
=0 - T22 =0 T22 = -(1 + n )DT + + � ��
dx2 dx24 Eh3 dx22 12(1 -n 2 ) �
�
dx2 2 �dx2 ��
�
10.7.10.3. Assume that w = dw / dx2 = 0 at x2 = ma . Hence, calculate the smallest value of T22 for
which a solution with nonzero w exists.
10.7.10.4. Deduce an expression for the critical temperature required to cause the film to buckle
10.7.10.5. Assume that DT exceeds the critical value. Calculate the resulting displacement field
10.7.10.6. Hence calculate the decrease in energy of the film during the buckling
10.7.10.7. Use the preceding result to deduce the crack tip energy release rate during delamination.

10.7.11. Consider a thin circular plate with thickness h, Young’s modulus

E, Poisson’s ratio n and thermal expansion coefficient  . Let (r ,q )
be a cylindrical-polar coordinate system and let { e r , eq , e z } be a triad
of unit vectors parallel to the natural (covariant) basis vectors for the
coordinate system. Assume that the deformation of the plate is axially
symmetric, so that the displacement field can be expressed as
u = u (r )e r + w(r )e z . Show that the Von-Karman equations
governing the motion of the plate can be reduced to the following
form
2
du 1 �dw � u
10.7.11.1. Mid-plane strain-deflection relation g rr = + � � g qq =
dr 2 �dr � r
d 2w 1 dw
10.7.11.2. Displacement-curvature relation Dk rr = - Dkqq = -
2 r dr
dr
Eh Eh
10.7.11.3. Membrane stress-strain relation Trr = 2
( g rr + ng qq ) Tqq = ( gqq + ng rr )
(1 - n ) (1 - n 2 )
10.7.11.4. Moment-curvature relations
 140

Eh3 Eh3
M rr = ( Dk rr + nDkqq ) M qq = ( Dkqq + nDk rr )
12(1 - n 2 ) 12(1 -n 2 )
10.7.11.5. Equations of motion
dTrr 1 d 2u dM rr 1
+ (Trr - Tqq ) = r h + ( M rr - M qq ) - Vr = 0
dr r dt 2 dr r
dVr Vr d 2w
+ - Trr Dk rr - Tqq Dkqq + p3 = r h
dr r dt 2

10.7.12.The Von-Karman equations for a circular plate cannot be

solved analytically, even for the case of axially symmetric
deformations. However, it is very straightforward to set up a
simple spreadsheet to calculate the potential energy of the
plate, and find a numerical solution by using the spreadsheet’s
solver function to minimize the potential energy. To this end,
consider a circular plate with radius R, thickness h, Young’s
modulus E, Poisson’s ratio n and thermal expansion
coefficient  . Assume that the plate is initially stress free, and
is then heated to increase its temperature uniformly by DT . At the same time, the surface of the plate is
subjected to a uniform pressure p acting perpendicular to its surface. Assume that the deformation of the
plate is axially symmetric, so that the displacement field can be expressed as u = u (r )er + w(r )e z .
10.7.12.1. Show that the potential energy of the plate can be expressed as
R
p Eh �
U=
2 �
(1 - n ) �
�
(
� )
(g rr - DT )2 + (g qq - DT ) 2 + 2n (g qq - DT )(g rr - DT ) rdr
�0
R � R
h2
+
12 �
( D k 2
rr + Dk 2
qq + 2n D k rr )
Dkqq rdr
�
�
�
- 2p p � w(r )rdr
0 � 0
g , g
where rr qq and D k rr , Dk qq are related to the displacements by the equations listed in
the preceding problem.
10.7.12.2. To obtain a numerical solution, let ri = iDr ui , wi ,
with i=0,1,2…N and Dr = R / N , denote the
ui
position and displacements of a set of N+1 equally
spaced points along a radius of the disk. The 0 1
derivatives of u and w can be approximated using wi
simple difference formulas. For example, define ri N
ui��du / dr = ( ui - ui -1 ) / Dr , Dr
wi��dw / dr = ( wi - wi -1 ) / Dr , i=1,2…N to be the derivatives of u and w at
r = (ri + ri -1 ) / 2 ; similarly, approximate the second derivative of w as
wi�
��d 2 w / dr 2 = ( wi� ) / Dr at each i=1,2…N-1. The value of d 2 w / dr 2 at r=R can be
+1 - wi�
estimated by extrapolation as w� ��( 2 w�
N N �) The integrals in the expression for the
�-1 - w�
N
potential energy can be evaluated using a simple piecewise-constant approximation (or if
you prefer piecewise linear) to the integrand. Use this approach to set up an EXCEL
spreadsheet to calculate the total potential energy of the plate. You should find N=40
sufficient for most practical purposes.
141

10.7.12.3. Check your spreadsheet by calculating the potential energy of a heated flat plate with
w = u = 0 , and compare the prediction with the exact solution.
10.7.12.4. The deflection of the plate under loading can be calculated by using the solver feature of
EXCEL to determine the values of ui and wi that minimize the potential energy. To test
your spreadsheet, use it to calculate the deflection of a circular plate with thickness
h/R=0.02, with a clamped boundary, which is subjected to a pressure
pR3 (1 - n 2 ) / Eh3 = 0.0025 . You will need to enforce the clamped boundary condition
using a constraint. You can set w ' N = 0 but you will find that this makes the plate slightly
too stiff. A better result is obtained by extrapolating the slope of the plate to r=R based on
the slope of the last two segments, and enforcing the constraint ( 3w ' N - w ' N -1 ) / 2 = 0 .
Show that, for this pressure, the numerical solution agrees with the predictions of the exact
small-deflection solution.
10.7.12.5. Show that if the deflection of the center of the plate is comparable to the plate thickness, the
small deflection solution predicts deflections that are substantially greater than the
numerical solution of the Von-Karman equations. This is because the in-plane strains
begin to significantly stiffen the plate. Plot a graph showing the deflection of the center of
the plate (normalized by plate thickness) as a function of normalized pressure
pR3 (1 - n 2 ) / Eh3 .

10.7.13. The interface between a thin film and its substrate contains a
circular crack with radius R. The film has Young’s modulus E,
Poisson’s ratio n and thermal expansion coefficient  . The
substrate can be idealized as rigid. The film is initially stress free,
then heated to raise its temperature by DT , inducing a uniform
biaxial stress field s b = - EDT db /(1 - n ) in the film. At a
critical temperature, the film buckles as shown in the figure. For R>>h the buckled film can be
modeled as a plate with clamped edge, so that the critical buckling temperature is given by the solution
to problem 10.7.8. When the film buckles, some of the strain energy in the film is relaxed. This
relaxation in energy can cause the film to delaminate from the substrate.
10.7.13.1. Show that the change in strain energy of the system during the formation of the buckle can
be expressed in dimensionless form by defining dimensionless measures of displacement,
position and strain as
uˆ = u / R DT wˆ = w / R rˆ = r DT / R
2
duˆ 1 �dwˆ � uˆ d 2 wˆ 1 dwˆ
gˆrr = + � � gˆqq = Dkˆrr = - Dkˆqq = -
drˆ 2 �drˆ � rˆ drˆ2 rˆ drˆ
so that
1� 1
b2 � 2 �
(1 + n ) DU
U0
= -� �
�rr
2 2
gˆ + gˆqq + 2ngˆrr gˆqq + 2 �
Dkˆ rr + Dkˆqq + 2nDkˆ rr Dkˆqq �
12 �
x dx + 2�
�� ( )
gˆrr + gˆqq x dx
0� � 0
where b = h / R DT and U 0 = p R 2 hE ( DT ) 2 /(1 - n ) is the total strain energy of the
circular portion of the film before buckling.
10.7.13.2. The implication of 10.7.13.1 is that the change in strain energy (and hence the normalized
displacement field which minimizes the potential energy) is a function of material and
geometric parameters only through Poisson’s ratio n and the dimensionless parameter b .
 142

Use the spreadsheet developed in problem 10.7.12 to plot a graph showing (1 + n )DU / U 0
as a function of b for a film with n = 0.3 . Verify that the critical value of b
corresponding to DU = 0 is consistent with the solution to problem 10.7.8.
10.7.13.3. Find an expression for the crack tip energy release rate G = - ( � DU / �R ) / ( 2p R ) and the
dimensionless function (1 + n )DU / U 0 = f ( b ,n ) . Hence, use the results of 10.7.13.2 to
plot a graph showing the crack tip energy release rate (suitably normalized) as a function of
b , for a film with n = 0.3 .

10.7.14.The figure shows a thin-walled, spherical dome with

radius R, thickness h and mass density r . The dome is open
at its top, so that the shell is bounded by spherical polar angles
q0 < q < q1 . Calculate the internal forces induced by
gravitational loading of the structure, using the membrane
theory of shells in Section 10.7.7.

10.7.15.The figure shows a thin-walled conical shell with thickness h and

mass density r . Calculate the internal forces induced by gravitational
loading of the structure, using the membrane theory of shells in Section
10.7.7. Use the cylindrical-polar coordinates ( z ,q ) as the coordinate
system.
143

Appendix A

11.
Vectors and Matrices

A.
A.1 Calculate the magnitudes of each of the vectors shown below
A.1.1 r = 3i + 6 j + 2k
A.1.2 r = 16i + 6 j
A.1.3 r = -9.6 j + 2.4i - 4.6k

A.2 A vector has magnitude 3, and i and j components of 1 and 2, respectively. Calculate its k component.

A.3 Let {i,j,k} be a Cartesian basis. A vector a has magnitude 4 and subtends angles of 30 degrees and 100
degrees to the i and k directions, respectively. Calculate the components of a in the basis {i,j,k}

A.4 Find the dot products of the vectors listed below

A.4.1 a = i+ j+k , b = i- j+k
A.4.2 a = 3i + 2 j , b = 4i + 6k
A.4.3 a = xi + yj + zk , b = xi - hj + k

A.5 Calculate the angle between each pair of vectors listed in Problem A.4. – i.e. find the angle q ( a , b)
between a and b in each case
b
a 1500
c
 144

A.6 The vectors a and b shown in the figure have magnitudes a = 3 , b = 5 . Calculate a  b .

A.7 Find the cross products of the vectors listed below

A.7.1 a = i+ j+k , b = i- j+k
A.7.2 a = 3i + 2 j , b = 4i + 6k
A.7.3 a = xi + yj + zk , b = x i - h j +  k

A.8 The vectors a and b shown in the figure below have magnitudes a = 3 , b
b = 5 . Calculate a  b . What is the direction of a  b ? a 1500
c

A.9 Let a=5i-6k; b=4i+2k Let r=5k. Express r as components parallel to a and b, i.e. find two scalars 
and b such that r =  a + b b

A.10 Find the direction cosines of the following vectors

A.10.1 a = i + j+k
A.10.2 a = 3i + 2 j

A.11Let a = 5i + 6 j - 3k , b = 3 j + 6k , c = 45i - 30 j + 15k .

A.11.1 Verify that a, b and c are mutually perpendicular and that c = a  b
A.11.2 In view of (1), three unit vectors {e1 , e 2 , e 3} parallel to a b and c can form a Cartesian
basis. Calculate the components of {e1 , e 2 , e 3} in the {i,j,k} basis.
A.11.3 Let r=4i+6k. Calculate the components of r in {e1 , e 2 , e 3} .

A.12 Let a and b be two unit vectors. Let  and b be two scalars, and let c be a vector such
that
c = a + bb
A.12.1 Prove that
(b � b) - ( a �
c)(a � c) (a � b) - (b �
c)(a � c)
= b=
1- ( a �
b) 1- ( a �
b)
2 2

A.12.2 Let {i,j,k} be a Cartesian basis. Consider the three vectors a = (3 / 5)i + (4 / 5) j ,
b = k , c = i + k . For this set of vectors, calculate values for
(b � b ) - (a �
c)(a � c) (a � b) - (b �
c)(a � c)
= b=
1- ( a �
b) 1- ( a �
b)
2 2

A.12.3 By substituting values, show that if a, b, c,  and b have the values given in
(3.2), then c � a + b b
A.12.4 In 30 words or less, explain why the vectors used in (3.2) and (3.3) cannot satisfy
c =  a + b b for any values of  and b . (You may use as many mathematical
equations and symbols as you like!).

A.13 Let a = 3i - 4 j , b = i + 2 j - 3k be two vectors, and let r = xi + yj + zk be a third vector with

unknown components
145

A.13.1 Solve the equation a + r = b

A.13.2 r =1
Find a solution (any solution will do) to the equation a �
A.13.3 Show that the equation a � r = 1 has more than one solution (e.g. by finding another
one!)
A.13.4 Suppose that the vector r satisfies
r�a=0 r�b=0
r
Show that must be parallel to a �b
A.13.5 Give all the solutions to the equations
a = 0 r�
r� b = 0 r�
r =1

A.14 Suppose that the three vectors {a, b, c} are used to define a (non-Cartesian) basis. This
means that a general vector r is expressed as r =  a + b b + g c , where ( , b , g ) denote the
components of r in the basis {a, b, c} . Let r1 = 1a + b1b + g 1c and r2 =  2a + b 2b + g 2c denote
two vectors, and let r3 = r1 �r2 . Calculate formulas for the components of r3 in {a, b, c} , i.e.
find formulas in terms of (1, b1 , g 1 ) and ( 2 , b 2 , g 2 ) for ( 3 , b 3 , g 3 ) such that
r3 =  3a + b3b + g 3c .

A.15 Here is a nice matrix.

�1 3 0�
[ A] = �3 2 0�
�
�
�
�0 0 5 �
�
A.15.1 Find det([ A] )
A.15.2 Find [ A] -1 (Don’t try to use the general expression for the inverse of a matrix – this matrix
can be inverted trivially)
A.15.3 Find the eigenvalues and eigenvectors of [ A] . (You can write down one of the eigenvalues
and eigenvectors by inspection. The other two can be found using the formulae for a 2 �2
matrix)

A.16 Consider a square, symmetric matrix

�1 4�
[ A] = �
�4 3�
�
A.16.1 Find the spectral decomposition of [ A] , i.e. find a diagonal matrix [ L ] and an orthogonal
matrix [ Q ] such that [ A] = [ Q ] [ L ] [ Q ] T
A.16.2 Hence, calculate [ A]1/2

A.17 The exponential of a matrix is defined as

�
[ A]k
exp([ A]) = �
k =1 k !
A.17.1 Show that the exponential of a diagonal matrix [ L ] = diag (l1, l2 ...ln ) is given by
 146

exp ( [ L ] ) = diag [ exp(l1 ),exp(l2 )...exp(ln ) ]

A.17.2 Let [ A] be a square, symmetric n �n matrix, and let [ L ] = diag (l1, l2 ...ln ) and [ Q ] be a
diagonal and orthogonal matrix, respectively, such that [ A] = [ Q ] [ L ] [ Q ] T . Find an
expression for exp ( [ A]) in terms of [ Q ] and [ L ]
A.17.3 Calculate the exponential of the matrix given in A.16.

A.18 Find a way to compute the log of a square, symmetric matrix.

A.19 A vector displacement field u(x) has components

u1 = -l x1 +  x2 x3
u2 = -l x2 -  x1 x3
u3 = b x3
where ( x1, x2 , x3 ) are the components of the position vector, and  , b , l are constants. Find expressions for
the matrix of components of grad(u) and find an expression for the divergence of u.

A.20 The figure shows a three noded, triangular finite element. (3,3,2)
In the basis { e1 , e2 , e3 } , the nodes have coordinates (2,1,3);
(4,2,3)
(4,2,3), (3,3,2). All dimensions are in cm.
A.20.1 Find the area of the element. Use vector algebra to do e2
m1
this – you do not need to calculate the lengths of the (2,1,3)
sides of the triangle or any angles between the sides. e1 m2 m3
A.20.2 Find the angles subtended by the sides of the e3
triangle
A.20.3 Find two unit vectors normal to the plane of the
element, expressing your answer as components in the basis { e1, e2 , e3 }
A.20.4 Let { m1, m 2 , m 3 } be a basis with m1 parallel to the base of the element, m 2 normal to the
element, and m3 chosen so as to ensure that { m1, m 2 , m3 } are right handed basis vectors.
Find the components of m1 , m 2 and m3 in the basis { e1, e 2 , e3 } . (To do this, first write
down m1 , choose one of the solutions to the preceding problem for m 2 , and then form m3
by taking the cross product of m 2 and m3 )
A.20.5 Set up the transformation matrix [ Q ] that relates vector components in { m1, m 2 , m 3 } to
those in { e1 , e 2 , e3 }
A.20.6 The displacement vector at the center of the element has coordinates (4,3,2) mm in the
basis { e1, e 2 , e3 } . Find its components in the basis { m1, m 2 , m3 }
147

Appendix B

12.
Brief Introduction to Tensors

B.
B.1 Let { e1 , e2 , e3 } be a Cartesian basis. Vector u has components (1, 2,0) in this basis, while tensors S and
T have components
�1 6 0� �1 -1 3 �
� � �
T �� 6 2 0 � S �� 2 4 2� �
�0 � � �
� 0 5 � �1 2 6 �
� �
B.1.1 Calculate the components of the following vectors and tensors
v = Tu v = u � T V =S +T V =S� T V = ST
B.1.2 Find the eigenvalues and the components of the eigenvectors of T.
B.1.3 Denote the three (unit) eigenvectors of T by m1 , m 2 , m3 (It doesn’t matter which
eigenvector is which, but be sure to state your choice clearly).
B.1.4 Let { m1, m 2 , m 3 } be a new Cartesian basis. Write down the components of T in
{ m1, m 2 , m3} . (Don’t make this hard: in the new basis, T must be diagonal, and the
diagonal elements must be the eigenvalues. Do you see why this is the case? You just need
to get them in the right order!)
B.1.5 Calculate the components of S in the basis { m1, m 2 , m 3 } .

B.2 Let S and T be tensors with components

 148

�1 2 -2 � 1 -1 3 �
�
T ��� -1 4 3 �
�
�
2 4 2�
S �� �
�
�1 2 6 � � �
1 2 1�
� �
B.2.1 Calculate S : T and S ��T
B.2.2 Calculate trace(S) and trace(T)

B.3 Let { e1, e2 , e3}

be a Cartesian basis, and let m1 = (5e1 + 6e 2 - 3e3 ) / 70 , m 2 = (3e2 + 6e3 ) / 45 ,
m3 = (9e1 - 8e2 + 3e3 ) / 154 be three unit vectors. The components of a tensor T in { e1 , e2 , e3 } are
�1 2 0�
T ���2 4 0� �
�
�0 0 5 �
�
B.3.1 Verify that { m1, m 2 , m3 } is also Cartesian basis.
B.3.2 Calculate the components of T in { m1, m 2 , m3 }

B.4 Let { e1, e 2 , e3 } be a Cartesian basis, and let m1 = (5e1 + 6e2 - 3e3 ) , m 2 = (3e 2 + 6e3 ) ,
m3 = (9e1 - 8e2 + 3e3 ) be three vectors.
B.4.1 Calculate the components in { e1 , e2 , e3 } of the tensor T that satisfies mi = T �
ei .
B.4.2 Calculate the eigenvalues and eigenvectors of T.
B.4.3 Calculate the components of T in a basis of unit vectors parallel to m1 , m 2 , m3 .

B.5 Let S be a tensor, and let

�S (e) S12( e) ( e) �
S13 �S (m ) S12
(m ) (m) �
S13
�11 � �11 �
(e ) ( e) ( e) � (m ) (m ) (m) �
[ S (e ) ] = �S21 S 22 S 23 [ S (m ) ] = �S21 S22 S23
� � � �
� (e ) ( e) ( e) � � (m ) (m ) (m) �
S31 S32 S33 S31 S32 S33
� � � �
l q
denote the components of S in Cartesian bases e1 , e 2 , e 3 and { m1, m 2 , m 3} , respectively. Show that

{
(e )
} { }
(m )
the trace of S is invariant to a change of basis, i.e. show that trace [ S ] = trace [ S ] .

B.6 Show that the inner product of two tensors is invariant to a change of basis.

B.7 Show that the outer product of two tensors is invariant to a change of basis.

B.8 Show that the eigenvalues of a tensor are invariant to a change of basis. Are the eigenvectors similarly
invariant?

B.9 Let S be a real symmetric tensor with three distinct eigenvalues li and corresponding eigenvectors
m ( j ) = 0 for i �j .
m (i ) . Show that m (i ) �
149

B.10 Let S be a real symmetric tensor with three distinct eigenvalues li and corresponding normalized
eigenvectors m (i ) satisfying m (i ) �
m (i ) = 1 . Use the results of B.9 to show that
3
b=
S� �li (m(i) �b) m(i)
i =1
for any arbitrary vector b, and hence deduce that
3
S= �li m(i) �m(i)
i =1

B.11 Use the results of B.10. to find a way to calculate the square root of a real, symmetric tensor.

B.12Let
�0 a12 a13 �
T ��
- a12
� 0 a23 �
�
- a13
�
� - a23 0 ��
Find expressions for the eigenvalues and eigenvectors of T in terms of its components aij

Appendix C

13.
Index Notation

C.
C.1 Which of the following equations are valid expressions using index notation? If you decide an
expression is invalid, state which rule is violated.
� s ij �2ui
(a) s ij = Cklij e kl (b) �kkk = 0 (c) + bi = r (d) �ijk �ijk = 6
�xj �t2

C.2 Match the meaning of each index notation expression shown below with an option from the list
(a) l = Tij Sij (b) Eij = Tik Skj (c) Eij = S kiTkj (d) ai =�kij b j ck (e) l = aibi
(f) d ij (g) Tij m j = l mi (h) ai = Sij b j (i) Aki Akj = d ij (j) Aij = A ji
(1) Product of two tensors
(2) Product of the transpose of a tensor with another tensor
(3) Cross product of two vectors
(4) Product of a vector and a tensor
(5) Components of the identity tensor
(6) Equation for the eigenvalues and eigenvectors of a tensor
(7) Contraction of a tensor
(8) Dot product of two vectors
(9) The definition of an orthogonal tensor
(10) Definition of a symmetric tensor
 150

C.3 Write out in full the three equations expressed by

1 �2uk �2ui b r �2ui
+ + i=
1 - 2n �
xk �
xi � xk m m �
xk � t2

C.4 Let a, b, c be three vectors. Use index notation to show that

a �(b �c) = ( a �
c) b - ( a �
b) c

C.5 Let A and B be tensors with components Aij and Bij . Use index notation to show that
( A �B ) T = BT �
AT

C.6 Let A , B and C be tensors with components Aij , Bij and Cij . Use index notation to show that
B ) : C = ( AT �
( A� C) : B = ( C �
BT ) : A

C.7 Let
En
Cijkl =
E
2 ( 1 +n )
(d ild jk + d ikd jl +) d d
( 1 +n ) ( 1 - 2n ) ij kl
1 +n n
Sijkl =
2E
( E
)
d ild jk + d ik d jl - d ijd kl
T
be two tensors. Calculate ijkl = S C
ijpq pqkl

�R �2 R
C.8 Let R = xk xk . Calculate and
�
xi �
xi �
xj

C.9 The stress-strain relations for an isotropic, linear elastic material are
1 +n n
e ij = s ij - s kk d ij + DT d ij
E E
Calculate the inverse relation giving stresses in terms of strains.

C.10Let s ij denote a symmetric second order tensor, and let

s e = Sij Sij Sij = s ij - s kk d ij / 3
Show that
�f 3 Sij
=
s ij 2 s e
�

C.11The strain energy density for a hypoelastic material is given by

151

( n +1) / 2 n
1 2ns 0e 0 �I 2 �
U = KI12 + �2 �
6 n +1 � �
�e 0 �
where
I1 = e kk I2 =
1
2
(
e ij e ij - e kk e pp / 3 )
Show that the stresses follow as
(1- n ) / 2 n
�U K �I � �e ij - e kk d ij / 3 �
s ij = = e kk d ij + s 0 � 2 � � �
e ij 3
� �e 2 � � e0 �
�0 �

C.12Let Fij denote the components of a second order tensor and let J = det(F) denote the determinant of
F. Show that
� J
= JF ji-1
�Fij

C.13 The strain energy density of a hyperelastic material with a Neo-Hookean constitutive relation is given
by
m K
U = 1 ( I1 - 3) + 1 ( J - 1) 2
2 2
where
F F
I1 = ki ki J = det(F)
J 2/3
Show that
1 �W m � 1 �
s ij = Fik = 1 �Fik F jk - Fkl Fkl d ij �+ K1 ( J - 1) d ij
J � Fkj J �5/3 3 �
You may use the solution to problem C.9

C.14 A hypoelastic material has a stress-strain relation given by

2 eij E 1
s ij = Sij + s kk d ij / 3 Sij = s e s kk = e kk
3 ee 1 - 2n 3
where
� 2 2
� 1 + n - � n - e e � - 1 e �e
2 �n - 1 e � e 0
1 2 se � � ( n - 1) � 0� n -1
eij = e ij - e kk dij ee = eij eij =�
3 3 s0 � 1/ n
�e e �
�� � e e �e 0
��e 0 �
�
and E = ns 0 / e 0 is the slope of the uniaxial stress-strain curve at e e = 0 . Show that

ds ij 4 e e 2 d d +d d d d E
= ( Et - Es ) ij 2kl + Es ( ik jl il jk - ij kl ) + d kl dij
d e kl 9 ee 3 2 3 9(1 - 2n )
where Es = s e / e e is the secant modulus, and Et = ds e / d e e
 152

Appendix D

14.
Vectors and Tensors in Polar Coordinates

D.
D.1 A geo-stationary satellite orbits the earth at radius 41000 km in the equatorial
plane, and is positioned at 0 o longitude. A satellite dish located in Providence,
Rhode Island (Longitude 410 40.3' N , Latitude 71034.6' W ) is to be pointed at
the satellite. In this problem, you will calculate the angles  and b to
position the satellite. Let {i,j,k} be a Cartesian basis with origin at the center of
the earth, k pointing to the North Pole and i pointing towards the intersection
of the equator (0o latitude) and the Greenwich meridian (0 o longitude). Define a
spherical-polar coordinate system ( R, f ,q ) with basis vectors {e R , ef , eq } in
the usual way. Take the earth’s radius as 6000 km.
D.1.1 Write down the values of ( R, f ,q ) for Providence, Rhode Island
D.1.2 Write down the position vector of the satellite in the Cartesian {i,j,k}coordinate system
D.1.3 Hence, find the position vector of the satellite relative to the center of the earth in the
{e R , ef , eq } basis located at Providence, Rhode Island.
D.1.4 Find the position vector SP of the satellite relative to Providence, Rhode Island, in terms of
basis vectors {e R , ef , eq } located at Providence, Rhode Island.
D.1.5 Find the components of a unit vector parallel to SP, in terms of basis vectors {e R , ef , eq }
located at Providence, Rhode Island.
153

D.1.6 Hence, calculate the angles  and b

D.2 Calculate the gradient and divergence of the following vector fields
D.2.1 v = R 2e R
D.2.2 v = R sin QeF
D.2.3 v = eR / R2

D.3 Show that the components of the gradient of a vector field in spherical-polar coordinates is
� �vR 1 � vR vq 1 � vR vf �
� - - �
�� R R � q R R sin q � f R �
�� v 1� vq vR 1 � vq vf �
v �Ѻ+-� q cot q �
�� R R � q R R sin q � f R �
�� �
�f
v 1� vf 1 � vf v v
+ cot q q + R �
�
�� R R � q R sin q � f R R�
�

D.4 In this problem you will derive the expression given in k eR

Appendix D for the gradient operator associated with polar
coordinates. P ef
D.4.1 Consider a scalar field f ( R,q , f ) . Write down an
expression for the change df in f due to an eq
infinitesimal change in the three coordinates O R
( R,q , f ) , to first order in (dR, dq , d f ) . f j
D.4.2 Write down an expression for the change in position
vector dr due to an infinitesimal change in the
three coordinates ( R,q , f ) , to first order in
i
(dR, dq , df ) , expressing your answer as
components in the {e R , eq , ef } basis.
D.4.3 Hence, find expressions for (dR, dq , d f ) in terms of (e R � dr, eq �
dr , ef �
dr )
D.4.4 Finally, substitute the result of (c) into the result of (a) to obtain an expression relating df to
dr . Rearrange the result into the form df = [ �f ] �r d , and hence deduce the expression for
the gradient operator.

D.5 In this problem you will derive the expression given in Appendix D k ez
for the gradient operator associated with polar coordinates.
D.5.1 Consider a scalar field f (r ,q , z ) . Write down an P eq
z r
expression for the change df in f due to an infinitesimal er
change in the three coordinates (r ,q , z ) , to first order in O
(dr , dq , dz ) .
D.5.2 Write down an expression for the change in position q j
vector dr due to an infinitesimal change in the three i
 154

coordinates (r ,q , z ) , to first order in (dr , dq , dz ) , expressing your answer as components in

the {e r , eq , e z } basis.
D.5.3 Hence, find expressions for (dr , dq , dz ) in terms of (er � dr , eq �
dr , e z �
dr )
D.5.4 Finally, substitute the result of (c) into the result of (a) to obtain an expression relating df to
dr . Rearrange the result into the form df = [ �f ] �r d , and hence deduce the expression for
the gradient operator.

D.6 Show that the components of the divergence of a tensor field S in spherical-polar coordinates are
��S RR S 1� Sq R S 1 � Sf R 1 �
� + 2 RR +
q
+ cot q q R +
R sin q � f
- (
Sqq + Sff � )
� � R R R � R R �
��S Rq S Rq 1 � Sqq Sqq 1 � Sfq Sq R Sff �
-+++++‫׹‬
S � � 2 cot q cot q �
�� R R R � q R R sin q � f R R �
��
S S �S S �Sff 1 �
� Rf + 2 Rf + sin q qf + cosq qf + 1 + (
Sf R + Sfq � )
��
� R R R � q R R sin q � f R �
�

D.7 Show that the components of the gradient of a vector field in cylindrical-polar coordinates are
�� vr 1 � vr vq � vr �
�� -
r r � q r �z �
� �
�� vq 1 � vq vr � vq �
v �Ѻ+ �
� r r � q r �z �
� �
�� vz 1� vz �vz �
�
�� r r � q �z �
�

D.8 Consider a rigid sphere, as k eR k er

shown in the figure. The
sphere is rotated through an P eF P ef
angle  about k and has
instantaneous angular eQ eq
velocity & about the k O R O r
direction. Let ( R, Q, F ) F j
f
denote the spherical-polar j
coordinates of a point in the
sphere prior to deformation, i i
and let {e , e , e }
R Q F denote
the spherical-polar basis Undeformed Deformed
vectors associated with this
point. Let (r ,q , f ) denote the spherical-polar coordinates of the same point after deformation, and let
{e r , eq , ef } denote the corresponding basis vectors.
D.8.1 Write down the deformation mapping relating (r ,q , f ) to ( R, Q, F)
155

D.8.2 Write down the velocity field in the sphere in terms of & and (r ,q , f ) , expressing your
answer as components in {e r , eq , ef }
D.8.3 Find the spatial velocity gradient L = v ��y as a function of (r ,q , f ) , expressing your
answer as components in {e r , eq , ef } .
D.8.4 Show that the deformation gradient can be expressed as F = er �e R + eq �eQ + ef �eF
and find a similar expression for F -1 . (its enough just to show that material fibers parallel
to the basis vectors in the undeformed solid are parallel to the basis vectors in the deformed
solid)
D.8.5 & -1
Use the result of (d) and (e) to verify that L = FF

Appendix E

15.
Miscellaneous Derivations

E.
E.1 Let Aij be a time varying tensor. Show that
dA-pq1 dAij
= - A-pi1 A-jq1
dt dt

E.2 Consider a deformable solid. Let V0

denote a closed region within the
undeformed solid, and let V be the n
same region of the solid in the S0 R0 S
V0 u(x) b
deformed configuration. Let A0 V T(n)
denote the area of the surface e2
surrounding V0 , and let A denote the y R
x t
area of the same surface after
Deformed
deformation. Let x denote the e1 Original Configuration
position of a material particle in the Configuration
y e3
solid before deformation and let be
the position of the same point after deformation. Define
 156

�yi
Fij = J = det(F)
�
xj
Show that
dA d � -1 dFkl dFkj �
= �ni ni dA = ��Fkl d ij - Fik-1 ni n j dA
�
dt dt � dt dt �
�
V �
V