class="bi x0 y0 w1 h1"
SOLUTIONS MANUAL
for
An Introduction to
The Finite Element M ethod
(Third Edition)
by
J. N. REDD Y
Department of Mechanical Engineering
Texas A & M University
College Station, Texas 77843-3123
PROPRIETARY A ND CONFIDENTIAL
This Manual is the proprietary property of The McGraw -Hill Companies, Inc.
(“McGraw-Hill”) and protected by copyright and other state and federal laws. By
opening and using this Man ual the user agrees to the following restrictions, and if
the recipient doe s not agree to these restrictions, the Manual should be promptly
returned unope ned to McGraw-Hill: This Manual is being provided only to
authorized professors and instructors for use in preparing for the classes
using the affiliated textbook. No other use or distribution of this Manual
is permitted. This Manual may not be sold and may not be distributed to
or used by any student or other third party. No part of this Manual
may be reproduced, displayed or distributed in any form or by any
means, electronic or otherwise, without the prior written permission of
the McGraw-Hill.
McGraw-Hill, New York, 2005
ii
iii
PREFACE
This solution manual is prepared to aid the instructor in discussing the solutions
to assigned problems in Chapters 1 through 14 from the book, An Introduction to
the Finite Element Method, Third Edition, McGraw—Hill, New York, 2006. Computer

J. N. Reddy
All t hat is not given is l ost.
iv
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
mgF
g
=
cvF
d
=
v
1
Chapter 1
INTRODUCTION
Problem 1.1: Newton’s second law can be expressed as
F = ma (1)
where F isthenetforceactingonthebody,m mass of the body, and a the
acceleration of the body in the direction of the net force. Use Eq. (1) to determine
the mathematical model, i.e., governing equation of a free-falling body. Consider
only the forces due to gravity and the air resistance. Assume that the air resistance
is linearly pr oportional to the velocity of the falling body.
Solution: From the free-body-diagram it follows that
m
dv
dt
= F
g
− F

(m
3
/day) and drained
at a r ate of q
0
(m
3
/day). Use the principle of conservation of mass to arrive at the
governing equation of the flow problem.
Solution: The conservation of mass requires
time rate of change in mass = mass inflow - mass outflow
The abov e equation for the problem at hand becomes
d
dt
(ρAh)=ρq
i
− ρq
0
or
d(Ah)
dt
= q
i
− q
0
where A is the area of cross section of the tank (A = πD
2
/4) and ρ is the mass density
of the liquid.
Problem 1.3: Consider the simple pendulum of Example 1.3.1. Write a computer

,v
i
). Table P1.3 contains
representative n umerical results.
Problem 1.4: An improvement of Euler’s method is provided by Heun’s method,
which uses the average of the derivatives at the tw o ends of the interval to estimate
the slope. Applied to the equation
du
dt
= f(t, u)(1)
Heun’s scheme has the form
u
i+1
= u
i
+
∆t
2
h
f(t
i
,u
i
)+f(t
i+1
,u
0
i+1
)
i

0.30 0.28185 0.37123 0.34225 -2.94148 -3.06403 -2.93371
0.35 0.13011 0.21803 0.19218 -3.10785 -3.35605 -3.17573
0.40 -0.02685 0.05023 0.03148 -3.14955 -3.53018 -3.29791
0.45 -0.18274 -0.12628 -0.13374 -3.06491 -3.57060 -3.29007
0.50 -0.33129 -0.30481 -0.29690 -2.85732 -3.46921 -3.15014
0.60 -0.58310 -0.63965 -0.59131 -2.11119 -2.85712 -2.50787
0.80 -0.78356 -1.05068 -0.91171 0.21536 -0.50399 -0.28356
1.00 -0.50591 -0.94062 -0.74672 2.41051 2.29398 2.19765
In books on numerical analysis, the second equation in (2) is called t he predictor
equation and the first equation is called the corrector equation. Apply Heun’s method
to Eqs. (1.3.4) and obtain the numerical solution for ∆t =0.05.
Solution: Heun’s method applied to the pair
dθ
dt
= v,
dv
dt
= −λ
2
sin θ
yields the following discrete equations:
θ
0
i+1
= θ
i
+ ∆tv
i
v
i+1

4 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
Table P1.4: Numerical solutions of the nonlinear equation d
2
θ/dt
2
+ λ
2
sin θ =0
along with the exact solution of the linear equation d
2
θ/dt
2
+λ
2
θ =0.
Exact Approx. solution θ Exact Approx. solution v
t θ Euler’s Heun’s v Euler’s Heun’s
0.00 0.785398 0.785398 0.785398 -0.000000 -0.0000 00 -0.0 00000
0.05 0.769645 0.785398 0.771168 -0.628013 -0.5692 21 -0.5 69221
0.10 0.723017 0.756937 0.728680 -1.230833 -1.1384 42 -1.1 21957
0.20 0.545784 0.615453 0.564818 -2.266146 -2.2098 38 -1.1 21957
0.40 -0.026852 0.050228 0.015246 -3.149552 -3.5301 78 -3.073095
0.60 -0.583104 -0 .639652 -0.544352 -2.111190 -2.8571 21 -2.194398
0.80 -0.783562 -1 .050679 -0.787095 0.215362 -0.503993 -0.114453
1.00 -0.505912 -0 .940622 -0.587339 2.410506 2.293983 2.023807
PROPRIETARY AND CONFIDENTIAL
This Manual is the propri etary property of The McGraw-Hill Companies, Inc. (“McGraw-Hill”)
and protected by copyright and other state and federal laws. By opening and using this Manual the
user agrees to the following restrictions, and if the recipient does not agree to these restrictions, the
Manual should be promptly returned unopened to McGraw-Hill: This Manual is being provided

¶
¯
¯
¯
¯
x=0
=0 u(1) =
√
2
Solution: Following the three-step procedure, we write the w eak form:
0=
Z
1
0
v
∙
−
d
dx
(u
du
dx
)+f
¸
dx (1)
=
Z
1
0
∙

¸
dx (3)
Fo r this problem, the weak form does not c ontain an expression that is linear in both
u and v; the expression is linear in v but not linear in u. Therefore, a quadratic
functional does not exist for this case. The expressions for B(·, ·)and`(·)aregiven
by
B(v,u)=
Z
1
0
u
dv
dx
du
dx
dx (not linear in u and not symmetric in u and v)
`(v)=−
Z
1
0
vfdx (4)
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
6 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
♠ New Problem 2.1:
The instructor may assign the following problem:
−
d
dx

+ vu
¸
dx (symmetric)
`(v)=
Z
1
0
vx
2
dx +6v(1) (2)
I(u)=
1
2
B(u, u) − `(u)=
1
2
Z
1
0
"
(1 + 2x
2
)
µ
du
dx
¶
2
+ u
2

dx
2
Ã
EI
d
2
w
dx
2
!
−
d
dx
(
EA
dw
dx
"
du
dx
+
1
2
µ
dw
dx
¶
2
#)
= q

Solution: The first step of the formulation is to multiply each equation with a weight
function, say v
1
for the first equation and v
2
for the second equation, and integrate
ov er the interval (0,L). In the second step, carry out the integration-by-parts once
in the firstequation,twiceinthefirst term of the second equation, and once in the
second part of the second equation. Then use the fact that v
1
(0) = v
1
(L)=0(because
u is specified there), v
2
(0) = v
2
(L) = 0 (because w is specified), and (dv
2
/dx)(0) = 0
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
SOLUTIONS MANUAL 7
(because dw/dx is specified at x = 0). In addition, we have EI(d
2
w/dx
2
)=M
0

L
0
(
EI
d
2
v
2
dx
2
d
2
w
dx
2
+ EA
dv
2
dx
dw
dx
"
du
dx
+
1
2
µ
dw
dx

0
(
EA
2
"
µ
du
dx
¶
2
+
du
dx
µ
dw
dx
¶
2
+
1
2
µ
dw
dx
¶
4
#
+
EI
2

+ a
12
∂u
∂y
¶
−
∂
∂y
µ
a
21
∂u
∂x
+ a
22
∂u
∂y
¶
+ f =0 in Ω
u = u
0
on Γ
1
,
µ
a
11
∂u
∂x
+ a

and t
0
are known functions on portions Γ
1
and Γ
2
of
the boundary Γ: Γ
1
+ Γ
2
= Γ.
Solution: Multiplying with the weight function v and integrating by parts, we obtain
the weak
0=
Z
Ω
∙
∂v
∂x
µ
a
11
∂u
∂x
+ a
12
∂u
∂y
¶

n
x
+
µ
a
21
∂u
∂x
+ a
22
∂u
∂y
¶
n
y
¸
ds
=
Z
Ω
∙
∂v
∂x
µ
a
11
∂u
∂x
+ a
12

where v =0onΓ
1
. The bilinear form (symmetric only if a
12
= a
21
) and li near form
are:
B(v,u)=
Z
Ω
µ
a
11
∂v
∂x
∂u
∂x
+ a
12
∂v
∂x
∂u
∂y
+ a
21
∂v
∂y
∂u
∂x

11
µ
∂u
∂x
¶
2
+2a
12
∂u
∂x
∂u
∂y
+ a
22
µ
∂u
∂y
¶
2
#
dxdy
−
Z
Ω
uf dxdy +
Z
Γ
2
ut
0

∂x
+ v
∂v
∂y
= −
1
ρ
∂P
∂y
+ ν
Ã
∂
2
v
∂x
2
+
∂
2
v
∂y
2
!
∂u
∂x
+
∂v
∂y
=0
⎫

∂u
∂x
n
x
+
∂u
∂y
n
y
¶
−
1
ρ
Pn
x
=
ˆ
t
x
ν
µ
∂v
∂x
n
x
+
∂v
∂y
n
y

∙
w
1
µ
u
∂u
∂x
+ v
∂u
∂y
¶
−
1
ρ
∂w
1
∂x
P + ν
µ
∂w
1
∂x
∂u
∂x
+
∂w
1
∂y
∂u
∂y

2
∂y
P + ν
µ
∂w
2
∂x
∂v
∂x
+
∂w
2
∂y
∂v
∂y
¶¸
dxdy
−
Z
Γ
2
w
2
ˆ
t
y
ds
0=
Z
Ω

∂ψ
∂y
∂ζ
∂x
=0
⎫
⎪
⎬
⎪
⎭
in Ω
Assume that all essential boundary conditions are specified to be zero.
Solution: First, we note the the identity
−w∇
2
ψ = −w∇· ∇ψ = −∇· (w∇ψ)+∇w · ∇ψ
and then use the Green—Gauss theorem to obtain
−
Z
Ω
w∇
2
ψ dxdy =
Z
Ω
[−∇ · (w∇ψ)+∇w · ∇ψ] dxdy
= −
I
Γ
w

· ∇ζ + w
2
µ
∂ψ
∂x
∂ζ
∂y
−
∂ψ
∂y
∂ζ
∂x
¶¸
dxdy (2)
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
10 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
Problem 2.6: Compute the coefficientmatrixandtheright-handsideoftheN-
parameter Ritz approximation of the equation
−
d
dx
∙
(1 + x)
du
dx
¸
=0 for 0<x<1
u(0) = 0,u(1) = 1

dx
dφ
j
dx
dx (1a)
F
i
= −B(φ
i
, φ
0
)=−
Z
1
0
(1 + x)
dφ
i
dx
dφ
0
dx
dx (1b)
The approximation functions φ
0
and φ
i
should be chosen such that
φ
0

(4a)
F
i
=
1
(1 + i)(2 + i)
(4b)
For the two-parameter (N =2)case,wehave
B
11
=
1
2
,B
12
= B
21
=
17
60
,B
22
=
7
30
,F
1
=
1
6

1
+ c
2
φ
2
= x +
55
131
(x − x
2
) −
20
131
(x
2
− x
3
)
=
1
131
(186x − 75x
2
+20x
3
)
The exact solution is given by
u
exact
=

Z
1
0
(1 + x)cosπx cos πxdx
B
12
=
Z
1
0
(1 + x)
dφ
1
dx
dφ
2
dx
dx =2π
2
Z
1
0
(1 + x)cosπx cos 2 πxdx= B
21
B
22
=
Z
1
0

Z
1
0
(1 + x)cosπx cos
πx
2
dx
F
2
= −
Z
1
0
(1 + x)
dφ
2
dx
dφ
0
dx
dx = −π
2
Z
1
0
(1 + x)cos2πx cos
πx
2
dx
Using the following trigonometric identities,

=
½
−
1
9
(6π −10)
68
225
+
4π
15
¾
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
12 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
and the solution is
U
2
(x)=c
1
sin πx + c
2
sin 2πx +sin
πx
2
= −0.12407 sin πx +0.02919 sin 2πx +sin
πx
2
Problem 2.8 A steel rod of diameter d = 2 cm, length L = 25 cm, and t hermal

Ak
=
βπD
1
4
πD
2
k
=
4β
kD
= 256 m
2
P being the perimeter and A the cross sectional area of the rod. The boundary
conditions are
θ(0) = T (0) − T
∞
= 100
◦
C,
µ
k
dθ
dx
+ βθ
¶
¯
¯
¯
¯

dφ
i
dx
dφ
j
dx
+ cφ
i
φ
j
¶
dx +ˆcφ
i
(L)φ
j
(L)(2a)
F
i
= −B(φ
i
, φ
0
)=−
Z
L
0
µ
dφ
i
dx

11
=
499
300
,B
12
= B
21
=
133
400
,B
22
=
91
1200
,F
1
= −832 ,F
2
= −
424
3
or
⎡
⎣
499
300
133
400

= −1, 033.3859 ,c
2
=2, 667.2635
The two-parameter Ritz solution is given by
θ(x)=100− 1033.3859x + 2667.2635x
2
θ(0.125) = 12.503
◦
C , θ(0.25) = 8.3575
◦
C
Problem 2.9: Set up the equations for the N-parameter Ritz approximation of
the following equations associated with a simply supported beam and subjected to a
uniform transverse load q = q
0
:
d
2
dx
2
Ã
EI
d
2
w
dx
2
!
= q
0

(i +1)(j +1)
i + j − 1
¸
F
i
=
q
0
(L)
i+2
(1 + i)(2 + i)
Note that the expre ssion given above for B
ij
is not valid when i =1andj =
1, 2, ···,N;wehave,
B
11
=4EIL, B
1j
= B
j1
=2EIL
j
, (j>1)
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
14 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
For N =1theRitzcoefficien t is given by c
1

L
2
24EI
x(L − x)=
q
0
L
4
24EI
x
L
(1 −
x
L
)
(b) Choose φ
0
=0andφ
i
=sin
iπx
L
.Thecoefficients are given by
B
ij
=
EIL
2
µ
iπ

iπ
¶
5
=
4q
0
L
4
EI
µ
1
iπ
¶
5
Hence, the solution becomes
w
2
(x)=c
1
φ
1
+ c
3
φ
3
=
4q
0
L
4

= q
0
L
"
L
i
a
+
i
a
Z
L
0
x
i−1
cos ax dx
#
− q
0
"
−
L
i+1
a
+
i +1
a
Z
L
0

3
and the solution is c
1
= c
2
L =2q
0
L
2
/(3EIπ
3
).
(b) Choose φ
0
=0andφ
i
=sin
iπx
L
.Thecoefficients B
ij
are the same as in Problem
2.9(b). The coefficients F
i
are given by F
1
= f
0
L/2andF
i

0
δ(x −
1
2
L), where δ(x)istheDirac
delta function (i.e., a point load Q
0
is applied at the center of the beam).
Solution: The coefficien ts F
i
are given by
(a) F
i
= Q
0
µ
L
2
¶
i+1
(b) F
i
= Q
0
(−1)
i−1
for i odd, and F
i
=0 for i even
Note that c

X
j=1
c
j
ψ
j
≡
N
X
j=1
c
j
cos
jπx
L
(1)
Substitution of Eq. (1) into the weak forms (S = GAK and D = EI)
0=
Z
L
0
∙
GAK
dv
1
dx
µ
dw
dx
+ Ψ

11
][K
12
]
[K
21
][K
22
]
¸½
{b}
{c}
¾
=
½
{F
1
}
{F
2
}
¾
(3)
where
K
11
ij
=
Z
L

ij
=
Z
L
0
GAKψ
i
dφ
j
dx
dx , K
22
ij
=
Z
L
0
µ
EI
dψ
i
dx
dψ
j
dx
+ GAK ψ
i
ψ
j
¶

µ
iπ
L
¶µ
jπ
L
¶
+
kL
2
,K
12
ij
= GAK
L
2
µ
iπ
L
¶
= K
21
ji
,
K
22
ij
=
L
2

clamped at the left end (x = 0) and subjected to an end mom ent, M
0
at x = L.
The latter can be assigned with (a) algebraic or (b) trigonometric approximation
functions. For example, for Problem 2a, we have the following (M,N)-parameter
Ritz solution with algebraic polynomials,
w
M
=
M
X
j=1
b
j
φ
j
≡
M
X
j=1
b
j
x
j
, Ψ
N
=
N
X
j=1

= −M
0
ψ
i
(L)(2)
Fo r the choice of approximation functions, φ
i
= ψ
i
= x
i
,thecoefficients can be
evaluated as,
K
11
ij
= GAK
ij
i + j − 1
(L)
i+j−1
,K
12
ij
= GAK
i
i + j
(L)
i+j
K

(L)
i+j−1
+ GAK
1
i + j +1
(L)
i+j+1
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
SOLUTIONS MANUAL 17
For M = N =1,wehave
b
1
=
q
0
L
3
6CEI
µ
3EI
GAKL
2
+1
¶
+
M
0
L

q
0
L
GAK
,c
1
= −
1
CEI
Ã
q
0
L
2
6
+ M
0
!
b
2
= −
q
0
L
2
12EI
µ
1 −
6EI
GAKL

2
2EI
Ψ(x)=
q
0
x
6EI
(−3L
2
+3Lx − x
2
) −
M
0
x
EI
(6)
Problem 2.13: Solve the P oisson equation governing heat conduction in a square
region:
−k∇
2
T = g
0
T =0 onsides x =1 and y =1 (1)
∂T
∂n
= 0 (insulated) on sides x =0 and y =0 (2)
using a one-parameter Ritz approximation of the form
T
1

¸
dxdy (4)
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
18 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
The coefficients B
11
and F
1
are given by
B
11
=
Z
1
0
Z
1
0
k
µ
∂φ
1
∂x
∂φ
1
∂x
+
∂φ

64
45
k (5a)
F
1
=
Z
1
0
Z
1
0
g
0
φ
1
dxdy
=
Z
1
0
Z
1
0
g
0
(1 − x
2
)(1 − y
2

(0) = θ(0) ,
∙
dφ
0
dx
+ˆcφ
0
¸
x=L
=0 (1)
and φ
i
must be selected such that it satisfies the homogeneous form of all specified
boundary conditions:
φ
i
(0) = 0 ,
∙
dφ
i
dx
+ˆcφ
i
¸
x=L
=0 (2)
To construct these functions, we begin with φ
0
= a+ bx, and determine the constants
a and b such that φ

The next function should be higher order than φ
1
; and there are two choices:
φ
2
= a + bx + cx
3
and φ
2
= a + bx
2
+ cx
3
.Forthefirst choice, we obtain,
φ
2
= x
∙
1 −
1+ˆcL
3+ˆcL
(
x
L
)
2
¸
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.

x=L
=0
Find a two-parameter Galerkin approximation of the problem using trigonometric
approximation functions, when (a) f = f
0
cos(πx/L)and(b) f = f
0
.
Solution: For this problem, we c an choo se φ
0
= 0 or a constan t (i.e., the solution
can be determined only within a constant) and φ
i
=cosiπx/L. The residual is given
by
R = −
N
X
i=1
c
j
d
2
φ
j
dx
2
− f
The weighted-residual statements are given by
0=

Rdx=(
2π
L
)
2
L
2
c
2
−
Z
L
0
f cos
2πx
L
dx
For (a) f = f
0
cos
πx
L
,weobtainc
1
=
f
0
L
2
π

the exact solution u
0
=1− x
2
. Use (a) the Galerkin method, (b) the least-squares
method, and (c) the Petrov—Galerkin method with weight function w =1.
PROPRIETARY MATERIAL.
c
°The M cGraw-Hill Companies, Inc. All rights reserved.
20 AN INTRODUCTION TO THE FINITE ELEMENT METHOD
Solution: We must cho ose φ
0
such that it satisfies all specified boundary conditions:
φ
0
(0) = 1 , φ
0
(1) = 0 (1)
and φ
i
must be selected such that it satisfies the homogeneous form of all specified
boundary conditions:
φ
i
(0) = 0 , φ
i
(1) = 0 (2)
Obviously, the following choice w ould meet the requirements,
φ
0

− 4
= −2
h
(1 − x)+c
1
(x − x
2
)
i
(−2c
1
)+[−1+c
1
(1 − 2x)]
2
− 4
= −3+2c
1
+(c
1
)
2
(4)
(a) The weighted-residual statement for the Galerkin method is given by
0=
Z
1
0
(x − x
2

2
.
(b) The least-squares statement is given by
0=
Z
1
0
dR
dc
1
Rdx=
Z
1
0
2(1 + c
1
)
h
−3+2c
1
+(c
1
)
2
i
dx
which giv es three solutions, (c
1
)
1