ELECTRONIC NOTES

CEE598 - Finite Elements

for Engineers: Module 1
S. D. Rajan
Department of Civil Engineering
Arizona State University
Tempe, AZ 85287-5306

Initial Version: Fall 1998
Last Revision: Fall 2010


Finit e Element s for Engineers
No part of this text may be reproduced or transmitted in any form or by any means, electronic or mechanical,
including photocopying, recording, or information storage or retrieval systems, for any purpose other than the user's
personal use, without the express written permission of the S. D. Rajan, Department of Civil Engineering, Arizona
State University, Tempe, AZ 85287-5306.

The software and manuals are provided "as is" without any warranty of any kind. No warranty is made regarding the
correctness, accuracy, reliability, or otherwise concerning the software and the manuals.

1998-2010 Dr. S. D. Rajan
Department of Civil Engineering
Arizona State University
Tempe, AZ 85287-5306
Phone 480.965.1712 Fax 480.965.0557
e-mail s.rajan@asu.edu





Table of Cont ent s
I nt r oduc t i on 1
Syl l abus 3
I nst r uc t or 4
Not at i on 5
Lesson Pl an 6
R E V I E W
Sol i d Mec hani c s 8
Heat Tr ansf er 9
Fl ui d Mec hani c s 10
El ec t r omagnet i c s & El ec t r ost at i c s 11
Li near Al gebr a & Numer i c al Anal ysi s 12
Di f f er ent i al Equat i ons & Cal c ul us 13
Revi ew Ex er c i ses 14
T H E S I X M A J O R S T E P S
Over vi ew of Fi ni t e El ement s 22
The Si x St eps 31
Revi ew Ex er c i ses 48
T H E G A L E R K I N S M E T H O D
Met hod of Wei ght ed Resi dual s 50
Appl yi ng t he Gal er k i ns Met hod 55
Revi ew Ex er c i ses 63
1 D P R O B L E M S
The El ement Conc ept 64
Revi ew Ex er c i ses 76
Sol i d Mec hani c s Ex ampl es 80
Heat Tr ansf er and El ec t r ost at i c Ex ampl es91
Hi gher -Or der El ement s 102
Mesh Ref i nement and Conver genc e 110
Revi ew Ex er c i ses 115

Usi ng t he 1DBVP

Pr ogr am 120

Sol vi ng Li near Al gebr ai c Equat i ons 123

I ndex 140

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
S. D. Raj an, 1998-2005 1
Int roduct ion
Change is inevitable, growth is optional. Anon.
his course is the first in a three-part series or modules titled Finite Elements for
Engineers that meets the mathematics requirements for the Master of Engineering
(M. Eng.) degree in the College of Engineering at Arizona State University.
Who should take this course?
Finite elements has become the defacto industry standard for solving multi-disciplinary
engineering problems that can be described by equations of calculus. Applications cut across
several industries by virtue of the applications solid mechanics (civil, aerospace, automotive,
mechanical, biomedical, electronic), fluid mechanics (geotechnical, aerospace, electronic,
environmental, hydraulics, biomedical, chemical), heat transfer (automotive, aerospace,
electronic, chemical), acoustics (automotive, mechanical, aerospace), electromagnetics
(electronic, aerospace) and many, many more.
Course Objectives
- To understand the basic ideas behind the finite element method.
- To understand and apply the Galerkins Method in solving (multi-disciplinary) one-
dimensional engineering problems.

Prerequisites
- Mathematics Linear Algebra, Numerical Analysis, Partial and/or ordinary differential
equations.
- Knowledge of undergraduate core material from (at least two of) solid mechanics, fluid
mechanics, heat transfer and electromagnetics. Knowledge in one area should be at an
advanced level (aligned with the undergraduate major).
- Knowledge of a high level programming language and use of computer-based tools.

Instructor-Student Interaction
To successfully meet the course objectives it is necessary that the students avail themselves of all
the resources discussion forums, e-mail, chat rooms, libraries. Keep the instructor and
teaching assistant informed of all your concerns. The web pages connected with this course will
T
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 2
contain instructions on how to communicate with the instructor regarding the questions you
may have or turning in the assignments etc.

Computer Programs and Computer Aids
There is one computer program 1DBVP

that is included on the CD-ROM. This is a

program that works under Windows 95/98, Windows NT as a console application. You
can either create the data using an editor or interactively in the program. The program
creates a text output file that can be viewed in the editor or printed.

The usage of other computer programs (including commercial programs that you may
have access to) that are more sophisticated and practical will follow in Modules 2 and 3.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 3
Syllabus
"A man with one watch knows what time it is. A man with two watches is
never sure. Segal's Law.
Suggested Texts
1

T1: Chandrupatla and Belegundu, Introduction to Finite Elements in Engineering, 2
nd

Edition, Prentice-Hall.
T2: Bickford, A First Course in the Finite Element Method, Irwin/McGraw-Hill.
T3: Burnett, Finite Element Analysis From Concepts to Applications, Addison-Wesley.

Other Texts
T4: Lewis, Morgan, Thomas and Seetharamu, The Finite Element Method in Heat Transfer
Analysis, Wiley.
T5: Reddy and Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics,
CRC Press.
T6: Jin, The Finite Element Method in Electromagnetics, Wiley.
T7: Cook, Malkus and Plesha, Concepts and Applications of Finite Element Analysis, 3
rd
Edition, Wiley.
T8: Reddy, An Introduction to the Finite Element Method, 2
nd
Edition, McGraw-Hill.
T9: Silvester and Ferrari, Finite Elements for Electrical Engineers, 3
rd
Edition, Cambridge
Press.
T10: Huebner and Thornton, The Finite Element Method for Engineers, 2
nd
Edition, Wiley.
T11: Zienkiewicz and Taylor, The Finite Element Method Volume 1, 4
th
Edition, McGraw-
Hill.
T12: Bathe, Finite Element Procedures, Prentice-Hall.

Outline (Lesson Plan)
- Review of Background Material
- Overview of Finite Elements. The Six Major Steps
- The Galerkins Method
- Applying Galerkins Method to One-Dimensional Problems
- Solving problems using the 1DBVP

program

1
Texts T1, T2 and/or T3 are strongly recommended.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 4

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 5
Inst ruct or
"Minds are like parachutes. They only function when they are open. Sir James
Dewar.
Subramaniam (Subby) D. Rajan
Professor of Civil and Environmental Engineering
I received my doctoral degree in Structural Mechanics from the University of Iowa in 1983. I
joined the Department of Civil Engineering at Arizona State University during Fall 1983.
Over the years I have taught a variety of undergraduate and graduate courses
Introduction to Engineering, Statics, Deformable Solids, Computer Programming for Civil
Engineers, Structural Analysis and Design, Matrix and Computer Applications in
Structural Engineering, Stress Analysis, Finite Elements for Civil Engineers, Advanced
Finite Element Analysis, Applied Optimal Design, and CAD Tools for Engineers. On and
off I have conducted workshops on Finite Elements, Design Optimization and Use of
Engineering Computer Tools.

My research interests are primarily in the areas of finite element analysis, design
optimization and engineering software development. The research efforts have been aimed
at solving engineering design problems in biomedical, civil structures, automotive,
mechanical, aerospace and electronic packaging areas. I have had the pleasure of
collaborating with my colleagues on research projects sponsored by the government
agencies such as National Science Foundation, NASA, Arizona DOT, and FHWA, and
private sector companies such as Motorola, Allied Signal, Allied American. My consulting
clients have included Special Devices Inc., York International, U-Haul, APS, Interstate
Assessment Technology etc.

You can reach me
- By phone at (480)965-1712
- By e-mail at s.rajan@asu.edu
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 6
Not at ion
As complexity rises, precise statements lose meaning, and meaningful statements
lose precision. Lotfi Zadeh
Vectors

1 n
a column vector with n rows

i
a element i of vector a

m 1
b row vector with m columns

Matrices

n m
A matrix with m rows and n columns

ij
A element row i and column j of matrix

Others
' y Derivative of y (or,
dx
dy
)
L (Units of) length
F (Units of) force
M (Units of) mass
t (Units of) time
T (Units of) temperature
E (Units of) energy

Greek Alphabets
Lower
Case
Upper
Case
English
Word
Lower
Case
Upper
Case
English
Word
Lower
Case
Upper
Case
English
Word
o
A
alpha |
B
Beta

I
Gamma
o
A
delta c
E
epsilon ,
Z
zeta
q
H
eta u O theta i
I
lota
k K
kappa A
lambda

M
mu
v N nu

xi o O omicron
t
H
pi

P
rho o
E
sigma
t
T
tau u
T
upsilon |
u
phi
_
X
chi

+
psi e
O
omega
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 7

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 8
Lesson Plan
"Problems cannot be solved by the same level of thinking that created them. A.
Einstein.
Module 1 is divided into four topics. Each topic has several lessons designed to focus on the
critical issues. With each lesson there is a set of objectives. I have also listed the relevant pages
from the list of textbooks that appear in the syllabus. There are several review problems at the
end of every topic. Solutions to most problems are also provided. Note that the set of problems
represents the minimal set needed to understand the material. You should solve more problems
from some of the referenced texts.
Topic 1 covers most of the prerequisites. Hence I have not listed the books to refer to nor
added my comments. Your undergraduate books or similar books are adequate. The
review problems should be very handy and bring you up to speed.

Topic 1: Review of Background Material.
Lesson 1: Solid Mechanics.
Lesson 2: Heat Transfer.
Lesson 3: Fluid Mechanics.
Lesson 4: Electrostatics & Electromagnetics
Lesson 5: Linear Algebra and Numerical Analysis.
Lesson 6: Differential Equations and Calculus
Review Exercises

Topic 2 is an introduction to finite elements. I have provided almost all of the material
text, examples and exercises.

Topic 2: The Six Major Steps.
Lesson 1: Overview of Finite Elements.
Lesson 2: The Six Steps.
Review Exercises

Topics 3 and 4 form the backbone of Module 1. The solution technique that is very
commonly used, the Galerkins Method, forms Topic 3. You should spend an adequate
amount of time here to understand the basics.

Topic 3: The Galerkins Method.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 9
Lesson 1: Method of Weighted Residuals.
Lesson 2: Applying the Galerkins Method.
Review Exercises

In Topic 4 we will look at the engineering problems described as one-dimensional
boundary-value (BVP) problems. These problems will be solved using the Galerkin-Based
Finite Element Method.

Topic 4: One-Dimensional Problems.
Lesson 1: The Element Concept.
Lesson 2: Solid Mechanics Examples.
Lesson 3: Heat Transfer and Fluid Flow Examples.
Lesson 4: Higher-Order Elements.
Lesson 5: Mesh Refinement and Convergence.
Review Exercises

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 10
Topic 1: Review
I'll play with it first and tell you what it is later. Miles Davis.
Lesson 1: Solid Mechanics
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To recognize and know the properties of simple structural systems
Truss
Beam
Frame

- To understand the concepts associated with
Equilibrium
Stress and Strain
Strain-Displacement Relations
Stress-Strain Relations
Simple theories of failure
Plane Elasticity Plane stress and strain

- To be able to compute
Strain energy
Work potential
Total Potential Energy
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 11
Lesson 2: Heat Transfer
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To understand the concepts associated with
Conservation of energy
Newtons Law of Cooling
Fouriers Law
Conduction, convection and radiation

- To recognize and know the properties of simple heat transfer problems involving
Steady-state conduction
Steady-state conduction and convection

- To be able to compute
Thermal energy
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 12
Lesson 3: Fluid Mechanics
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To understand the concepts associated with
Mass balance
Energy balance
Momentum balance
Pipe flow

- To recognize and know the properties of simple fluid mechanics problems involving
Inviscid compressible and incompressible flow
Viscous flow (newtonian/non-newtonian, laminar/transition/turbulent,
compressible/incompressible)
Flow through porous media
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 13
Lesson 4: Electromagnetics & Electrostatics
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To understand the concepts associated with
Electrostatic Fields
Magnetostatic Fields

- To recognize and know the numerical solution to
Maxwells Equations
Wave Equations
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 14
Lesson 5: Linear Algebra & Numerical Analysis
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To understand the concepts associated with
Vectors
Matrices
Polynomial approximation and interpolation

- To recognize and know the numerical solution to
Linear algebraic equations
Linear eigenvalue problem
Numerical integration
Numerical differentiation

- To be able to write a program in a high-level language to generate
A library of matrix operations (addition, transpose, multiplication etc.)
The solution to a set of linear algebraic equations
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 15
Lesson 6: Differential Equations and Calculus
Objectives: Knowledge of the following ideas and topics is crucial to understanding the
subsequent lessons. The major objectives are listed below.

- To understand the concepts associated with
Ordinary Differential Equations
Partial Differential Equation
Divergence Theorem

- To recognize and know the properties of
Initial Value Problems
Boundary Value Problems
Elliptic and Parabolic Partial Differential Equations

- To be able to
Integrate by parts
Derive classical solution to simple ODEs
Derive classical solution to simple PDEs
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 16
Review Exercises
SOLID MECHANICS
Problem T1L1-1
In a plane strain problem for a material with modulus of elasticity psi E ) 10 ( 20
6
= and
Poissons ratio 3 . 0 = v , the state of stress at a point is given by
{ } ksi ksi
xy y x
7 , 0 , 20 = = = t o o . Compute (a)
z
o , and (b)
x
c .

Problem T1L1-2
In a two-dimensional (plane stress) body, the y x displacements designated as ) , ( v u are
given as (in mm)
y y x y x u 8 3 ) , (
2 2
+ =
xy y x y x v 8 6 3 ) , ( + =
Determine the state of strain at the point ) 4 , 2 ( . Using the material properties from
Problem 1, determine the state of stress at the point.

Problem T1L1-3
In a solid body, the state of stress at a point is given by
{ }
zx yz xy z y x
t t t o o o , , , , , ={ }MPa 5 , 10 , 20 , 15 , 5 , 10 . Find the state of stress on a plane
whose normal has the direction cosines |
.
|

\
|

3
1
,
3
1
,
3
1
. Find the principal stresses at the
point and the orientations of the principal planes.

Problem T1L1-4
Consider the slender rod as shown in the Fig. T1LR-1. The rod has a circular cross-section
of radius cm 0.5 . The strain at any point is given as x x
x
8 2
2
= c . What is the strain
energy in the bar if the bar is made of aluminum?

L= 20cm
u,x

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 17

Fig. T1LR-1
Problem T1L1-5
Consider the slender rod as shown in Fig. T1LR-2. The rod has a circular cross-section of
radius cm 0.5 and is made of aluminum. What is the potential energy of the bar?

L= 20cm
u,x
1 kN

Fig. T1LR-2

Problem T1L1-6
A simply supported beam is shown in the Fig. T1LR-3. Compute the largest transverse
displacement and rotation in the beam. EI is a constant.

5 kN/m
10 m
A B


Fig. T1LR-3

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 18
HEAT TRANSFER
Problem T1L2-1
The inside and outside surfaces of a window glass are at C

40 and C

10 , respectively. The
glass has a thermal conductivity of
C m
W

8 . 0 . The glass is 75 cm by 50 cm and 1.25 cm

thick. Determine the heat loss through the glass over a period of 45 minutes.

Problem T1L2-2
An insulating board is to be used to limit the heat loss to
2
50
m
W
for a temperature
difference of C

125 across it. The thermal conductivity of the material is

C m
W

04 . 0 .
Determine the thickness of the insulating board.

Problem T1L2-3
The inside surface of an insulation layer is maintained at C T

200
1
= . The outside surface is
dissipating heat by convection into air at C T

20 =

. The insulating layer has a thickness of

cm 5 and thermal conductivity of
C m
W

5 . 1 . What is the minimum value of the heat

transfer coefficient at the outside surface, if the temperature
2
T at the outside surface
should not exceed C

100 .

Problem T1L2-4
Determine the interface temperature
1
T and the surface temperature
3
T of the composite
wall shown in the Fig. T1LR-4.
C m
W
k

= 1 . 0
1
,
C m
W
k

= 0 . 1
2
and
C m
W
k

= 0 . 2
3
.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 19
2 cm
150 C
k k k
T T T T
3 2 1
1 2 3 0
50 C
4 cm 3 cm

Fig. T1LR-4

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 20
FLUID MECHANICS
Problem T1L3-1
A pipeline with a 28 cm inside diameter is carrying liquid at a flow rate of
s
m
3
03 . 0 . A
reducer is placed in the line, and the outlet diameter is 15 cm. Determine the velocity at the
beginning and end of the reducer.

Problem T1L3-2
For the piping system shown below, the distance from A to B is 4500 ft, and the main line
is made of 10-nominal, schedule 40 wrought iron pipe. The attached loop is 8-nominal,
schedule 40 wrought iron pipe. The flow
1
Q is
s
ft
3
3 . 0 of water. Determine the flow rate
in both branches.

A
A
A
A
,V ,Q
Q
,Q
,Q
1
2
3
1 1
1
2
3
B


Fig. T1LR-5

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 21
LINEAR ALGEBRA & NUMERICAL ANALYSIS
Problem T1L5-1
Define the following terms.
(a) Symmetric matrix
(b) Singular matrix
(c) Rank of a matrix and rank deficiency
(d) Positive definite matrix

Problem T1L5-2
Given the matrix

(
(
(

=
10 8 6
8 3 5
6 5 4
A
Compute its determinant. Is the matrix positive definite?

Problem T1L5-3
Solve the set of linear equations given below.
(
(
(

10 2 1
2 8 3
1 3 7

3
2
1
x
x
x
=

2
0
4

Problem T1L5-4
Numerically integrate the following expressions and compare the answers with the
analytical solution.
(a)
}
+
5
3
2 3
) 10 3 12 ( dx x x x
(b)
}

+

1
1
2
13 6
1 2
dx
x x
x


Problem T1L5-5
Numerically differentiate the following expressions and compare the answers with the
analytical solution.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 22
(a) 10 8 12 ) (
3
+ = x x x f at 3 = x .
(b)
10 4
4 3
) (
+

=
x
x
x f at 1 = x .

Problem T1L5-6
Given a set of four points ) 8 , 1 ( , ) 12 , 2 ( , ) 10 , 3 ( , and ) 8 , 4 ( , do a least-squares fit of a
quadratic polynomial.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 23
DIFFERENTIAL EQUATIONS & CALCULUS
Problem T1L6-1
Differentiate between an ordinary differential equation (ODE) and a partial differential
equation (PDE)?

The boundary conditions to a differential equation can be classified as essential
(Dirichlet), natural (Neumann) or mixed (Robin). What are the traits of each of these types
of boundary conditions?

What is meant by a problem being well-posed?

Problem T1L6-2
Explain the terms - Initial Value Problems, Boundary Value Problems, Elliptic and
Parabolic Partial Differential Equations.

Problem T1L6-3
Integrate the following expressions
(a)
}
dx x x cos
2

(b)
}

dx xe
x


Problem T1L6-4
Solve the differential equations given below.
(a) 0 1 2
2
= + + x xy
dx
dy
x 0 ) 1 ( = y
(b) ( ) 0
) (
1 =
|
.
|

\
|
+
dx
x dy
x
dx
d
2 1 < < x
1 ) 1 ( = = x y 1 ) 1 (
2
=
|
.
|

\
|
+
= x
dx
dy
x
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 24
Topic 2: The Six Major St eps
A common mistake that people make when trying to design
something completely foolproof is to underestimate the ingenuity of
complete fools. Douglas Adams
Lesson 1: Overview of Finite Elements
Objectives: In this lesson we will look at the what and when of finite element analysis.
The major objectives are listed below.

- To understand what is meant by finite element method.
- To understand when the finite element method can be used.

Suggested Reading
Reference Pages
T1 1-2
T2 4-7
T3 3-32

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 25
What is the Finite Element Method?

Definition
The Finite Element Method (FEM) is a numerical technique to obtain approximate
solutions to a wide variety of engineering problems where the variables are related by
means of algebraic, differential and integral equations.

Ponder over the keywords that are underlined.

Approach
Briefly, the engineering problem that is described by governing differential or integral
equations is transformed to a set of algebraic equations or eigenequations. These equations
are solved numerically.

Flavors of FEM
FEM is not one technique or methodology. There are several different approaches to
formulating the problem solution. Three approaches are presented below.

Direct Stiffness Method: Historically this was the first approach used in the area of structural
mechanics. We will see the details of this method in the next lesson.

Variational Approach: The original problem is converted to a functional
2
that is
extremized. The extremum gives the solution. We will look at this approach very briefly
with respect to solid mechanics problems.

Weighted Residuals Approach: An approximate solution is assumed symbolically. The error
in the approximate solution is weighted over the domain of the problem and minimized.
We will use this approach as the major approach throughout the course.

Let us revisit the definition. What are the characteristics of the governing equations and the
solution? The characteristics will determine whether we have a
- Time-independent problem static (solid mechanics), steady-state (heat transfer)
- Time-dependent problem dynamic and wave propagation (solid mechanics), diffusion
(heat transfer)
- Linear problem Solution obtained in one pass since the problem parameters are
constants
- Nonlinear problem - Solution obtained in several passes since the problem parameters
are not constants but interdependent

2
A functional is a function of a function.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 26
- Eigenproblem modal analysis (solid mechanics), acoustics, heat transfer

There are further differences in problem formulation p-finite elements, stress hybrid
finite elements (solid mechanics), edge elements (electromagnetics).

Usage
There are three major stages in any finite element solution.
Pre-Processing: The physical system is converted to a finite element model. This is the most
labor intensive and time-consuming step.
Solution: Once the pre-processor has been used to build the model, the finite element
analysis can take place. This is the most compute intensive step.
Post-Processing: A post-processor is a program or part of a program that is used to examine
the results from the FE analysis. Post-processing is usually graphical but can also be query
or report based.

Applications of the Finite Element Method
Finite elements have been used to solve problems in a variety of industries such as
Automotive stress analysis, heat transfer, crashworthiness.
Civil structures stress analysis, structural dynamics.
Mechanical stress analysis, dynamics.
Aerospace structural analysis, fracture mechanics, fluid flow.
Biomedical stress analysis, fluid flow.
Electromagnetics parasitic, harmonic.
Semiconductor (Electronic) Packaging stress analysis, contact, impact.

Commercial Codes
Some of the commercial programs that are finite-element based are listed below.
ABAQUS: http://www.hks.com/
ADINA: http://www.adina.com/
ALGOR: http://www.algor.com/
ANSYS: http://www.ansys.com/
COSMOS: http://www.cosmosm.com/
MSC-NASTRAN: http://www.msc.com/

Specialized Preprocesors
Parametric Technology: http://www.ptc.com/
PATRAN: http://www.msc.com
Solidworks: http://www.solidworks.com/
Unigraphics: http://www.ugs.com

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 27
Finite Element Terminology (the basics)

Fig. T2L1-1 shows a model of a roof truss that is commonly used in residential buildings.
The truss members are usually made of wood. Note that the model contains numbers (1
through 14) that denote locations of certain key points on the truss.

MODEL
UNDO : V4.4 T0A-112
1
2
4
5
7
8
3
2 5
6
2
11 11
12 12
8
1 99 7
8
13 13
14 14
5
7 10 10 4
11
9
12
7
13
7
14
10
X
Y
Z

Fig. T2L1-1 Model of roof truss

Fig. T2L1-2 shows the finite element model of a barricade that is used to test air bag
explosive charges. The barricade is made of two chambers - the supply chamber and the main
chamber. The top of the chambers is vented (opening as shown in the figure). The isometric
view shows a collection of triangles used to model the walls.
Using these two examples, we will look at some of the commonly used terms with finite
element models.
Elements
Elements can be thought of as the basic building blocks. The collection of the finite
number elements forms the finite element mesh. In Fig. T2L1-1, the elements are the truss
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 28
members between two key points. In Fig. T2L1-2, the triangles are the elements. The finite
element mesh
3
is an approximation to the physical problem.



Fig. T2L1-2 Finite element model of a barricade


Nodes
Nodes are the key points. Nodes are used to define elements. In Fig. T2L1-1, two nodes are
needed to define an element a straight line. In Fig. T2L1-2, three nodes are needed to
define an element the triangle.

Element Properties
Since the elements have a physical sense, they are described in terms of element properties -
geometry, dimensions and material. To describe the members in the truss, we need to
know the shape of the members (typically rectangular), the dimensions (sides of the
rectangle) and the material properties (properties of the wood). To describe the triangles,
we need to know the thickness of the triangle (the thickness of the walls etc. of the
chambers) and the material properties of the material making the chamber walls or sides.

Loads

3
The collection of finite elements forms the finite element mesh.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 29
In both the above examples, loads are applied to the structure. The responses that we are
looking for are quantities such as displacements, stresses, strains etc. that result from the
loads applied on the structure.

Boundary Conditions
Finally, in order to compute the response quantities, we need to know how the truss or
the barricade is supported. Since the supports are on the boundary of the structure, they
are called boundary conditions. As we will see later, there are other examples of boundary
conditions apart from supports.


An Illustrative Example T2L1-1

Problem Statement
Consider the following problem. We need to estimate the area of a circle of radius R. We
will assume that we do not know the exact answer. However, we know some of the
rudiments of trigonometry and geometry.

Solution
Step 1: We will assume that we know the formula for the area of a triangle. We will split
the circle into a collection of triangles (that can be used as an approximation to the circle)
as shown as Mesh A. Note that the mesh of triangles will provide a lower bound to the
exact solution.

Mesh A
Fig. T2L1-3

We can see the concept of elements (and nodes) with this example. The continuous region
(circle) is represented in each mesh by a finite number of triangles. Each triangle is a finite
element. Note that in the mesh every triangle is the same as all the others. When we have
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 30
this situation, the mesh is known as a uniform mesh. While it is not always necessary to
work with a uniform mesh, it is appropriate for this problem.

Step 2: For a typical triangle that subtends an angle u at the center of the circle as shown
below

R
h
u/2
u/2
b

Fig. T2L1-4

) 2 / sin(u R b = (T2L1-1)
) 2 / cos(u R h = (T2L1-2)

If there are n triangles in the mesh,
n
t
u
2
= . Hence, the area of one triangle is

|
.
|

\
|
=
n
R
a
e
t 2
sin
2
2
(T2L1-3)

Step 3: Now the estimate the area of the circle we use the obvious concept the sum of the
areas of all the triangles is the area of the circle. In other words

=
|
.
|

\
|
= =
n
e
e
n
A
n
nR
a A
1
2
) (
2
sin
2
t
(T2L1-4)

Step 4: Given the numerical value of R , we can estimate the area of the circle by
substituting the numerical values in the above equation. It is apparent that the accuracy is a
function of n . In other words, increasing n decreases the error. However note that one
has to contend with the numerical problems by increasing n .

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 31
100
150
200
250
300
350
10 20 30 40
n
A
r
e
a
Area of a circle of radius 10 units

Fig. T2L1-5

Now, consider the collection of triangles that can be used as an approximation to the circle
as shown as Mesh B. The mesh of triangles will provide an upper bound to the exact
solution.

Mesh B

Fig. T2L1-6


Conclusions
Lets summarize the steps. First, we split the circle (problem domain) into simple shape
whose properties are known. The process of splitting the problem domain into elements
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 32
(or, creating the finite element mesh) is known as discretization. Second, we computed the
properties of a typical element. Third, we estimated the area of the circle used the concept
of summing the areas of the individual triangles. This is known as assembly where the
properties of the elements are assembled to form the properties of the system. In this
example, the property of the system was a scalar (a single unknown). Last, the equation
was solved numerically to obtain the solution.

In the next lesson we will see these steps in greater detail.


When should finite elements be used?

Finite elements while extremely powerful should be used with caution. FEM should NOT
be used as a black box.

Let us look at some of the sources of error in a finite element solution.

Translating the Physical System to the FE Model
There are inherent approximations when a physical system is transformed into a finite
element model. The approximations include not only those involving the geometry,
boundary conditions, and like, but also the approximations inherent in the FE formulation
(arising out of the assumptions).

Errors in Parameter Values
The material properties, the loads etc. are best estimates when input into a model.

Numerical Errors
The solution is obtained numerically and is susceptible to numerical errors truncation
and round-off.

Of these three, the first two play a very major role for most of the problems.

Finally, note that there are some classes of problems that can be solved more appropriately
by other techniques such as finite differences, finite strips, boundary elements, etc.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 33
Lesson 2: The Six Steps
Objectives: In this lesson we will look at the six steps in a typical finite element solution.
- To understand the basics of the direct stiffness method.
- To understand the concepts associated with primary unknowns, elements, nodes,
element equations, assembly, boundary conditions and solution.
- To be able to apply the method to solving a variety of one-dimensional problems.

Suggested Reading
Reference Pages
T1 The direct stiffness method is not discussed here. However you
may read the pages 9-13 that discusses the Rayleigh-Ritz Method
that is closely tied to the Variational approach that we will see
later.
T2 While the term direct stiffness method is not used, pages 7-37
illustrate the basic ideas.
T3


You can also refer to Huebner and Thornton, The Finite Element Method for Engineers, 2
nd

Edition, Wiley (Chapter 2).

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 34
The Six Steps

We will look at the six major steps in a typical finite element problem by using examples
from different engineering areas. These examples will deal with discrete systems. These are
systems similar to the roof truss example shown earlier. The elements are essentially one-
dimensional. Their geometry and properties are described by a single spatial variable x .
Refer to the Illustrative Example T2L1-1 from the previous lesson to reflect on some of the
terminology.

Step 1: Discretization
The process of breaking the system or structure into elements is quite simple (unlike
creating the triangles in the barricade example) for systems involving discrete elements.
The discrete elements are the finite elements.

Step 2: Element Equations
To develop the equations that capture the behavior of a typical element, we need to know
the physics of the problems. In the next few sections we will illustrate how this is done for
discrete problems in a variety of engineering problems.

We will look at some engineering applications in the areas of solid mechanics, heat
transfer, fluid flow and electrical networks. The motivation is to show the similarities
between these different specialty areas and point out that the differences are purely based
on the physics of the problem.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 35
Linear Springs (Primary unknowns: Displacements)
The physical problem is one involving a collection of springs that are acted upon by
external forces. The intent is to calculate the deformation and force in every spring. We
will use the following terminology - the primary unknown is displacement (this is what we
wish to compute first) and the force in the spring is a derived variable (one that can be
computed provided we have the values of the primary variables). Consider the free-body
diagram of a typical linear spring i as shown in the figure below.

x
a
b
k
u
u
f
f
i
a
b
a
b

Fig. T2L2-1

The spring element is defined in terms of two nodes that are labeled a and b .
a
u
represents the displacement of node a and
a
f is the force acting at node a . Note that the
displacements and forces are all shown acting in the positive coordinate direction. In a
general derivation we always assume that all the unknown quantities are positive (the
solution to an actual problem will tell us the sign associated with a specific unknown). The
spring constant of the element is
i
k . In any linear spring

k
f
u = (T2L2-1)
Now let us look at the free-body diagrams of the two nodes.

x
a
b
f
f ) )
k k (u (u - u - u
a
b
i i
a b b a


Using the concept of equilibrium

b i a i b a i a
u k u k u u k f = = ) ( (T2L2-2)

a i b i a b i b
u k u k u u k f = = ) ( (T2L2-3)

These two equations can be written in the form of matrix equations as

)
`

=
)
`

b
a
b
a
i
f
f
u
u
k
1 1
1 1
(T2L2-4)

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 36
or,
1 2 1 2 2 2
= f u k (T2L2-5)

where
2 2
k is the element stiffness matrix

1 2
u is the element nodal displacements vector

1 2
f is the element nodal forces vector

These are the element equations that capture the concepts of equilibrium-compatibility for
a typical spring.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 37
1D Heat Flow (Primary unknowns: Temperature)
The physical problem is one involving a collection of thermal elements that are essentially
in a one-dimensional heat flow. The intent is to calculate the temperature distribution and
the heat flux. Consider the diagram of a typical thermal element i as

x
k
T
q
q
T
a b
t
i
a
a
b
b

Fig. T2L2-2

Using Fouriers Law

dx
dT
kA q = (T2L2-6)
where q is the heat flux or energy flow in the x direction(positive if entering the
element), k is the thermal conductivity that is assumed to be a constant, T is the
temperature and A is the area. With reference to the typical element, at the two nodes we
have

t
T T
A k q
b a
i a
) (
= (T2L2-7)

t
T T
A k q
a b
i b
) (
= (T2L2-8)
These equations can be written in the matrix form as
)
`

=
)
`

b
a
b
a
i
q
q
T
T
t
A k
1 1
1 1
(T2L2-9)

or,
1 2 1 2 2 2
= q T k (T2L2-10)

where
2 2
k is the element thermal conductivity matrix
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 38

1 2
T is the element nodal temperature vector

1 2
q is the element nodal flux vector
These are the element equations that capture the concepts of conservation of energy or
energy balance for a typical thermal element.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 39
Pipe Flow (Primary unknowns: Pressure)
The physical problem is one involving a collection of pipes. The intent is to calculate the
pressure (at the ends of the pipe) and the flow. Let us assume that the pipes are circular and
the flow is fully developed laminar flow. Consider a typical fluid pipe element i as

p
p
Q
Q
a
L
b
b
a
b
a
Fig. T2L2-3

Using Darcys Law, the pressure drop between the ends of a pipe

4 2 1
128
D
QL
p p
t

= (T2L2-11)
where L is the length of the pipe, Q is the flow (positive if entering the element), is the
dynamic viscosity and D is the pipe diameter. The flows entering the ends of the typical
pipe are given as
) (
128
4
b a a
p p
L
D
Q =

t
(T2L2-12)
) (
128
4
a b b
p p
L
D
Q =

t
(T2L2-13)
These equations can be written in the matrix form as
)
`

=
)
`

b
a
b
a
Q
Q
p
p
L
D
1 1
1 1
128
4

t
(T2L2-14)

or,
1 2 1 2 2 2
= Q p k (T2L2-15)

where
2 2
k is the element fluidity matrix

1 2
p is the element nodal pressure vector

1 2
Q is the element nodal flow vector

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 40
These are the element equations that capture the concepts of conservation of mass or mass
balance for a typical fluid pipe element.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 41
Electrical Network (Primary unknowns: Voltage)
The physical problem is one involving an electrical network. The intent is to calculate the
voltage (at the ends of the resistors or elements). Consider a typical element i as

x
a
b
r
V
V
I
I
i
a
b
a
b

Fig. T2L2-4

Using Ohmss Law, the voltage drop between the ends of an element
rI V V =
2 1
(T2L2-16)
where r is the resistance, I is the current (positive if entering the element), and V is the
voltage. The current entering the ends of the typical element are given as
) (
1
b a
i
a
V V
r
I = (T2L2-17)
) (
1
a b
i
b
V V
r
I = (T2L2-18)
These equations can be written in the matrix form as
)
`

=
)
`

b
a
b
a
i
I
I
V
V
r 1 1
1 1
1
(T2L2-19)

or,
1 2 1 2 2 2
= I V k (T2L2-20)

where
2 2
k is the element resistivity matrix

1 2
V is the element nodal voltage vector

1 2
I is the element nodal current vector

These are the element equations that capture the concepts of electrical flow and continuity
in a typical resistor.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 42
Step 3: Assembly
The obvious question is What do we do with the element equations? The answer is provided
with respect to an example shown below.

Illustrative Example T2L2-1
The figure shows a collection of 4 springs that is loaded with two concentrated forces.

1
1
2
3 4
2
4
3
X
5
k
k
k
P
P
k
2
1
3
3
4
4

Fig. T2L2-5

The first task is to number the nodes. We have labeled the nodes starting at 1,
consecutively all the way to 5. Designating a node as 1 or 2 is arbitrary. Next, the
elements are labeled starting at 1 consecutively all the way to 4. Hence, the finite element
mesh has 5 nodes and 4 elements. Each node has one unknown variable associated with it
the displacement at the node. In FE terminology, we state that there is one degree of
freedom (DOF) per node. The problem is described in terms of a total of 5 DOF -
{ }
5 4 3 2 1
, , , , D D D D D . What we have done here is what needs to be done in Step 1 with
respect to discretizing a problem.

Lets go ahead and implement Step 2 for this problem.

Element 1: 2 = a and 3 = b
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 43
)
`

=
)
`

1
3
1
2
3
2
1
1 1
1 1
f
f
D
D
k (T2L2-21)

Element 2: 1 = a and 2 = b
)
`

=
)
`

2
2
2
1
2
1
2
1 1
1 1
f
f
D
D
k (T2L2-22)



Element 3: 2 = a and 4 = b
)
`

=
)
`

3
4
3
2
4
2
3
1 1
1 1
f
f
D
D
k (T2L2-23)

Element 4: 4 = a and 5 = b
)
`

=
)
`

4
5
4
4
5
4
4
1 1
1 1
f
f
D
D
k (T2L2-24)

A slightly different notation is used above with respect to applying Eqn. (T2L2-4). The
element nodal displacements are mapped from the generic
a
u and
b
u to the appropriate
global designation
n
D . The forces acting at the ends of the element are denoted
j
i
f
meaning that the force is acting at end i of element j (the superscript j is necessary to
differentiate the contributions made by different elements that share a node).

The assembly process is now to take these 8 equilibrium-compatibility equations and
construct 5 equilibrium-compatibility equations for the entire system. The system
equations look like

(
(
(
(
(
(

5
4
3
2
1
5
4
3
2
1
55 54 53 52 51
45 44 43 42 41
35 34 33 32 31
25 24 23 22 21
15 14 13 12 11
F
F
F
F
F
D
D
D
D
D
K K K K K
K K K K K
K K K K K
K K K K K
K K K K K
(T2L2-25a)

or,
1 5 1 5 5 5
= F D K (T2L2-25b)
where K is the global (or, structural, or, system) stiffness matrix, D and F are the vector
of global nodal displacements and global nodal forces respectively.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 44
But how do we construct Eqns. (T2L2-25)? Now back to the element equations. The first
equation for element 1 is the equilibrium-compatibility of element 1 along
2
D . Similarly,
the second equation for element 2 and the first equation for element 3 are the equilibrium-
compatibility along
2
D . Hence these three equations combine to generate the second
global equation. In a similar manner we can collect the appropriate equations for the other
four DOFs.

However, one recommended way to assemble is to go through sequentially one element at
a time and update Eqn. (T2L2-25). Hence after using the Eqns. (T2L2-21) for element 1 we
have

(
(
(
(
(
(

0
0
0
0 0 0 0 0
0 0 0 0 0
0 0 0
0 0 0
0 0 0 0 0
1
3
1
2
5
4
3
2
1
1 1
1 1
f
f
D
D
D
D
D
k k
k k
(T2L2-26)

After using Eqns. (T2L2-22) for element 2 we have

+
=

(
(
(
(
(
(

+

0
0
0 0 0 0 0
0 0 0 0 0
0 0 0
0 0
0 0 0
1
3
2
2
1
2
2
1
5
4
3
2
1
1 1
1 2 1 2
2 2
f
f f
f
D
D
D
D
D
k k
k k k k
k k
(T2L2-27)

After using Eqns. (T2L2-23) for element 3 we have

+ +
=

(
(
(
(
(
(

+ +

0 0 0 0 0 0
0 0 0
0 0 0
0
0 0 0
3
4
1
3
3
2
2
2
1
2
2
1
5
4
3
2
1
3 3
1 1
3 1 3 2 1 2
2 2
f
f
f f f
f
D
D
D
D
D
k k
k k
k k k k k k
k k
(T2L2-28)

Finally, after all the four elements have been assembled

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 45

+
+ +
=

(
(
(
(
(
(

+ +

4
5
4
4
3
4
1
3
3
2
2
2
1
2
2
1
5
4
3
2
1
4 4
4 4 3 3
1 1
3 1 3 2 1 2
2 2
0 0 0
0 0
0 0 0
0
0 0 0
f
f f
f
f f f
f
D
D
D
D
D
k k
k k k k
k k
k k k k k k
k k
(T2L2-29)

Now comparing the system load vector to the Fig. T2L2-5, 0
3
2
2
2
1
2
= + + f f f ,
3
1
3
P f = and
4
4
4
3
4
P f f = + .
1
F and
5
F are not known and are the support reactions that develop at the
two supports. Interestingly enough 0
1
= D and 0
5
= D , a fact that we will use in the next
step.
Hence the equations reduce to

(
(
(
(
(
(

+ +

5
4
3
1
5
4
3
2
1
4 4
4 4 3 3
1 1
3 1 3 2 1 2
2 2
0
0 0 0
0 0
0 0 0
0
0 0 0
F
P
P
F
D
D
D
D
D
k k
k k k k
k k
k k k k k k
k k
(T2L2-30)

The rest of the problem will be solved using numerical data. Let

4 1
10 k
in
lb
k = = ,
in
lb
k 15
3
= ,
in
lb
k 20
2
= , lb F 100
3
= , lb F 50
4
=


Step 4: Imposition of Boundary Conditions

Eqns. (T2L2-29) cannot be solved for two reasons first we do not know the RHS vector
completely and second, the K matrix is singular! The boundary conditions for the given
problem refer to the essential boundary conditions (EBC), i.e. 0
1
= D and 0
5
= D .
Imposing the EBCs will make it possible to solve the system equations.

Since 0
1
= D , the first equation is not necessary to solve for the other unknowns. Similarly,
since 0
5
= D , the last equation is not necessary to solve for the other unknowns. These
two conditions can be simply imposed by modifying the first and the last equations and
the columns corresponding to
1
D and
5
D , i.e. the first and the last columns
4
. Hence,


4
The LHS is KD and involves matrix multiplication.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 46

(
(
(
(
(
(

+ +
0
0
0
1 0 0 0 0
0 0 0
0 0 0
0 0
0 0 0 0 1
4
3
5
4
3
2
1
4 3 3
1 1
3 1 3 2 1
P
P
D
D
D
D
D
k k k
k k
k k k k k
(T2L2-31a)

or, 1 5
1 5
5 5

= F D K (T2L2-31b)

where the parameters are the modified global stiffness, global displacement and modified
global load vector respectively.

This method of imposing the EBC is known as the direct or elimination approach.
Examine the above equations the first and the last equations are decoupled from the rest
and solving them yields 0
1
= D and 0
5
= D which is what we set out to do.

Now substituting the numerical values, we have

(
(
(
(
(
(

0
50
100
0
0
1 0 0 0 0
0 25 0 15 0
0 0 10 10 0
0 15 10 45 0
0 0 0 0 1
5
4
3
2
1
D
D
D
D
D
(T2L2-32)


Step 5: Solution of Primary Unknowns

The system equations can be solved using a variety of techniques. Some of the more important
techniques used in FE programs are
(a) Direct Method: Gaussian Elimination Technique (examples include
T
LDL or Cholesky
Decomposition for symmetric, positive definite coefficient matrix), and
(b) Iterative Method: Preconditioned Conjugate Gradient Method.
While the theory is somewhat standard for these techniques, there exist tens of ways of
implementing the solution techniques.

Solving Eqns. (T2L2-32), we have { } { } 0 , 5 , 15 , 5 , 0 , , , ,
5 4 3 2 1
= D D D D D in .
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 47


Step 6: Obtaining Derived Variables

Using the nodal displacements we can now compute the forces in the springs and the
support reactions.

Using Eqns. (T2L2-4) for each element (units are in and lb ).

Element 1:
)
`

=
)
`

=
)
`

b
a
f
f
100
100
15
5
1 1
1 1
10 (T2L2-33)



Element 2:
)
`

=
)
`

=
)
`

b
a
f
f
100
100
5
0
1 1
1 1
20 (T2L2-34)

Element 3:
)
`

=
)
`

=
)
`

b
a
f
f
0
0
5
5
1 1
1 1
15 (T2L2-35)

Element 4:
)
`

=
)
`

=
)
`

b
a
f
f
50
50
0
5
1 1
1 1
10 (T2L2-36)

Lets interpret one of these results by selecting Element 2. The FBD for the element is

x
1
100 100
2
k
2


Fig. T2L2-6

indicating that the spring is in equilibrium (and in tension). Draw the FBDs for all the
other elements.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 48
Finally, lets draw the FBD for nodes 1 and 2.

X
1
2
100
100
0
100
F
1


Fig. T2L2-7

Again note that each node should be in equilibrium. Using the FBD of node 1

+
= + = = 100 100 0
1 1
F F F
x

and the reaction at node 1 (a support) is lb 100 . An examination of Node 2 FBD shows
that it is in equilibrium.


Some Important Observations

- The element stiffness matrix k is symmetric and singular.
- The global stiffness matrix K is symmetric, banded and singular.
- The modified global stiffness matrix K is symmetric, banded and nonsingular (positive
definite).
- It will be necessary to look at the imposition of the EBC again in order to handle
situations where the EBCs are non-homogeneous.

Consider the set of equations shown below.

(
(
(

3
2
1
3
2
1
33 32 31
23 22 21
13 12 11
b
b
b
x
x
x
A A A
A A A
A A A
(T2L2-38)

We will first rewrite the equations as
1 3 13 2 12 1 11
b x A x A x A = + +
2 3 23 2 22 1 21
b x A x A x A = + + (T2L2-39)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 49
3 3 33 2 32 1 31
b x A x A x A = + +

Let the EBC be such that c x =
2
, a known constant. Imposing the EBC in the first and last
equations, we have
c A b x A b x A x A
12 1 2 12 1 3 13 1 11
= = +
c A b x A b x A x A
32 3 2 32 3 3 33 1 31
= = +

Since c x =
2
, the second equation can be modified to
c x x x = + +
3 2 1
) 0 ( ) 1 ( ) 0 (
Collecting the three modified equations, we now have

(
(
(

c A b
c
c A b
x
x
x
A A
A A
32 3
12 1
3
2
1
33 31
13 11
0
0 1 0
0
(T2L2-40)
which can be solved for the unknown vector x . Note that c can be any known number
including zero.



Summary

The six major steps in a typical finite element solution are as follows.

Step 1: Discretization
The problem domain is discretized into a collection of simple shapes (or, elements).

Step 2: Element Equations
Knowing the physics of the problem once can develop the element equations for a typical
element. These equations are symbolic, and can be evaluated numerically by substituting
the numerical values for the symbols (problem parameters). The primary unknowns are
evaluated at the nodes of the elements.

Step 3: Assembly
The element equations for each element in the FE mesh needs to be assembled into a set of
global equations that represent the properties of the entire system.

Step 4: Imposition of Boundary Conditions
The global equations cannot be solved unless the boundary conditions are imposed. The
boundary conditions reflect which of the primary unknowns have known values.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 50
Step 5: Solution of Primary Unknowns
The modified global equations from Step 4 are solved for the primary unknowns. One has
to be careful about issues like condition number of the stiffness matrix and round-off
errors.

Step 6: Obtaining Derived Variables
The evaluated values of the primary unknowns are used to obtain derived values using the
relationships between the primary unknowns and the derived variables.

Finally, one must use the computed answers to check whether the physical solution makes
sense. In the example, we checked for the equilibrium of the nodes and elements.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 51
Review Exercises
Problem T2L1-1
Estimate the area of the circle of unit radius using Mesh A and Mesh B. Compare the two
solutions.
Problem T2L2-1
A uniform bar is made up of two sections is loaded by concentrated nodal forces as shown.
Use a two-element finite element model to compute the nodal displacements, element
forces and stresses. Hint: The stiffness k of a bar of length L , uniform cross-sectional area
A and modulus of elasticity E is
L
AE
k = .

0.5 k
3
2
A = 2 in
60 in
1 k
120 in
E = 10 psi

E = 5(10 ) psi
2
1
A= 1.25 in
2
2
7
6


Problem T2L2-2
A composite wall consists of layers of aluminum, copper, and steel. The surface
temperatures are K

264 and K

311 . Use a three-element model to determine the

interfacial temperatures and the heat flow per unit area.

3
4
Al
k
Cu
Steel
2 in
4 in 1 in
264 K
2 Al 202
Cu 388
Steel 45
W
m- K
}

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 52
Problem T2L2-3
Laminar incompressible flow occurs in the branched network of circular pipes. If
s
m
3
085 . 0 fluid enters and leaves the piping network, use a three-node, four-element model
to compute the fluid nodal pressures and the volume flow rate in each pipe.

Pipe Diameter ( cm) Length ( m )
1 7.6 30.5
2 5.1 45.8
3 5.1 45.8
4 10.2 91.5
The dynamic viscosity is
2
3
) 10 ( 96 . 0
m
s N

.

1
2
3
1
3
3
4
3
P = 0 0.085m /s

Problem T2L2-4
The voltages at the output terminals of the DC circuit have the values shown. Use a four-
node finite element model to compute the voltages at each node and the current in each
resistor.
1
2
3
4
1
1
4
4
4 5
3
2
V =20 V
V =0
I
I
O O
O
O
1
2
3
4

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 53
Topic 3: The Galerkins Met hod
"Engineering problems are under-defined, there are many solutions, good, bad
and indifferent. The art is to arrive at a good solution. This is a creative activity,
involving imagination, intuition and deliberate choice. Ove Arup
Lesson 1: Method of Weighted Residuals.
Objectives: In this lesson we will look at Method of Weighted Residuals to solve one-
dimensional differential equations and the Galerkins method in particular.

- To understand the Method of Weighted Residuals.
- To understand the Galerkins Method and concepts associated with trial solution,
optimizing criteria, accuracy.

Suggested Reading
Reference Pages
T1 13-16
T2 Section 2.3
T3 Chapter 3


F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 54
Background

One-dimensional boundary-value problems have several applications some of which we
saw in the context of the direct stiffness method in the previous topic. We will list some
additional problems below.

Solid Mechanics Transverse deflection of a cable
Hydrodynamics One-dimensional flow in an
inviscid, incompressible fluid
Magnetostatics One-dimensional magnetic
potential distribution
Heat Conduction One-dimensional heat flow in a
solid medium
Electrostatics One-dimensional electric potential
distribution

Consider a one-dimensional boundary value problem (also known as equilibrium problem)
given by
) ( ) ( ) (
2
2
x F y x Q
dx
dy
x P
dx
y d
= + + (T3L1-1)
If P and Q are constants, then
) (
2
2
x F Qy
dx
dy
P
dx
y d
= + + (T3L1-2)
The total solution is
..... ) ( ' ) ( ) (
1 0
2 1
+ + + + = x F C x F C Be Ae x y
x x
(T3L1-3)
where the constants of integration A and B can be found by substituting the two
boundary conditions into Eqn. (T3L1-2). Note that we need two boundary conditions for
the problem to be well-posed.

The boundary conditions are of three types.
(a) The function y may be specified. This is known as the Dirichlet or Essential boundary
condition.
(b) The derivative ' y may be specified. This is known as the Neumann or Natural
boundary condition.
(c) The function y and the derivative ' y may be specified. This is known as the Mixed or
Robin boundary condition.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 55
In the FE solution, an approximate or trial solution
~
( ) y x is constructed and solved for
5
.
The FE approach has three distinct operations.
(a) A trial solution
~
( ) y x is constructed.
(b) An optimizing criterion is applied to
~
( ) y x .
(c) An estimation of the accuracy of
~
( ) y x is made.

Trial Solution
The trial solution is constructed with a finite number of terms as
) ( ...... ) ( ) ( ) ( ) ; (
2 2 1 1 0
~
x a x a x a x a x y
n n
| | | | + + + + = (T3L1-4)
where ) (x
i
| are known trial or basis functions and the coefficients
i
a are undetermined
parameters known as degrees of freedom (DOF). The purpose of the trial function ) (
0
x | is
to satisfy some or all of the boundary conditions. The most common form of trial
solutions is to use polynomials. We will see more about this later.

Optimizing Criterion
The optimizing criterion is used to generate the appropriate equations so that we can solve
for the numerical values of the coefficients
i
a . As you can guess, the optimizing criterion is
not unique and different approaches define what is meant by the best possible
approximation to the exact solution. The two most common forms are
(a) The Method of Weighted Residuals (MWR) applicable when the problem is described
by differential equations, and
(b) The Ritz Variational Method (RVM) - applicable when the problem is described by integral
(or, variational) equations.

In this lesson, the focus is on the former method. In the MWR, the criteria minimize an
expression of error in the differential equation. In the RVM, an attempt is made to
extremize (typically, minimize) a physical quantity. We will see this approach in the
second module of this course.

Accuracy Estimate
Since the FE solution is usually approximate and since for practical problems it is virtually
impossible to generate or estimate the exact solution, we need a measure to tell us how

5
Those familiar with the finite difference solution will note that there are two approaches to solving the boundary-value
ODE the shooting (initial-value) method and the equilibrium (boundary-value) method.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 56
close the FE solution is to the exact solution. We will illustrate the ideas and estimates
later.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 57
Method of Weighted Residuals

There are at least four different optimizing criteria and we will see them in this lesson.
Consider the differential equation (T3L1-1) rewritten as
0 ) ( ) ( ) (
2
2
= + + x F y x Q
dx
dy
x P
dx
y d
(T3L1-5a)
or, 0 ) ( ) ( = x F y (T3L1-5b)

Substituting the trial solution, we have, in general
0 ) ( ) ( ) ; (
~
= = x F y a x R (T3L1-6)

) ; ( a x R is known as the residual or error in the solution. In the MWR, the optimizing
criterion is to find the numerical values for
i
a which will make ) ; ( a x R as close to zero as
possible for all values of x throughout the domain of the problem. Once a specific
criterion is applied, a set of algebraic equations is produced. As you can see, the process is
to transform the original (linear) ODE to a set of linear algebraic equations.

Collocation Method
In this method, for each of the parameter
i
a a point
i
x is chosen and the residual is set to
zero at that point.
0 ) ; ( = a x R
i
n i ,.., 1 = (T3L1-7)
The points
i
x are known as the collocation points. Note that we are setting the error in the
residual, not the error in the solution, to zero.

Subdomain Method
In this method, for each of the parameter
i
a an interval
i
x A is chosen and the average of
the residual is set to zero.
0 ) ; (
1
=
A
}
A
i
x i
dx a x R
x
n i ,.., 1 = (T3L1-8)
The intervals
i
x A are called the subdomains.

Least-squares Method
In this method, with respect to each
i
a we minimize the integral over the entire domain,
the square of the residual.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 58
0 ) ; (
2
=
c
c
}
O
dx a x R
a
i
n i ,.., 1 = (T3L1-9)

The Galerkins Method
In this method, for each
i
a we require that a weighted average of the residual over the
entire domain be zero. The weighting functions are the trial functions ) (x
i
| associated
with each
i
a .
0 ) ( ) ; ( =
}
O
dx x a x R
i
| n i ,.., 1 = (T3L1-10)

The natural question is Which of these techniques is the most appropriate? A detailed
answer is outside the scope of this course. However, experience has shown that the
Galerkins Method is the most suitable for the type of finite element applications to be
seen in this course.

In the next lesson, we will apply the Galerkins Method in the classical sense. In the next
topic we will overcome the drawbacks of the classical approach and use the finite element
concepts.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 59
Lesson 2: Galerkins Method
Objectives: In this lesson we will look at applying the Galerkins Method (in the
theoretical or classical sense) and understanding its traits.

- To apply the Galerkins Method.
- To understand the concepts of trial solution, accuracy, classical solution, numerical
solution and convergence.

Suggested Reading
Reference Pages
T1
T2 56-60
T3 Chapter 3



F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 60
An Illustrative Example T3L2-1

Let us look at an example to see how the Galerkins Method works. Consider the problem

DE:
2
2 ) (
x dx
x dy
x
dx
d
=
|
.
|

\
|
2 1 s s x (T3L2-1)
BC: 2 ) 1 ( = = x y (T3L2-2)

2
1
2
=
|
.
|

\
|

= x
dx
dy
x (T3L2-3)

Classical (or, Theoretical) Solution
Let us assume the trial solution as

3
4
2
3 2 1
~
x a x a x a a y + + + = (T3L2-4a)
Hence,

2
4 3 2
~
3 2 x a x a a
dx
y d
+ + = (T3L2-4b)
In order to satisfy the BC, we need

3
4
2
3 2 1
~
) 1 ( ) 1 ( ) 1 ( 2 ) 1 ( a a a a x y + + + = = =
or, 2
4 3 2 1
= + + + a a a a (T3L2-5)
And,

4 3 2
2
~
24 8 2
2
1
a a a
dx
y d
x
x
= =
|
|
|
.
|

\
|

=

or,
4
1
12 4
4 3 2
= + + a a a (T3L2-6)
Eqns. (T3L2-5) and (T3L2-6) are known as constraint equations since the 4
i
a s in Eqn.
(T3L2-4a) are no longer independent. The two constraint equations render only 2 out of
the 4
i
a s as independent parameters (or, DOF). Using the two constraint equations to
write
1
a and
2
a in terms of the other two, we have

4 3 2 1
2 a a a a = (T3L2-7)

4 3 2
12 4
4
1
a a a = (T3L2-8)
Substituting these two equations in (T3L2-4a)
) 11 )( 1 ( ) 3 )( 1 ( ) 1 (
4
1
2
2
4 3
~
+ + + = x x x a x x a x y (T3L2-9)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 61
This equation is now in the familiar form
) ( ) ( ) ( ) ; (
2 2 1 1 0
~
x a x a x a x y | | | + + = (T3L2-10)
where ) 1 (
4
1
2 ) (
0
= x x |
) 3 )( 1 ( ) (
1
= x x x |
) 11 )( 1 ( ) (
2
2
+ = x x x x |
and we will assume for the sake of convenience that
1
a and
2
a in (T3L2-10) are the same as
3
a and
4
a in (T3L2-9).

We will now write the residual as
0
2 ) (
) ; (
2
~
=
|
|
|
.
|

\
|
=
x dx
x y d
x
dx
d
a x R (T3L2-11)
Substituting (T3L2-10) in the above equation, we have

2 2
2
1
2
) 4 3 ( 3 ) 1 ( 4
4
1
) ; (
x
a x a x a x R + + = (T3L2-12)
Now using the Galerkins Method (Eqn. (T3L1-10))
0 ) 3 )( 1 (
2
) 4 3 ( 3 ) 1 ( 4
4
1
2
1
2 2
2
1
=
|
.
|

\
|
+ +
}
dx x x
x
a x a x
0 ) 11 )( 1 (
2
) 4 3 ( 3 ) 1 ( 4
4
1
2
2
1
2 2
2
1
= +
|
.
|

\
|
+ +
}
dx x x x
x
a x a x (T3L2-13a)
Integrating yields two equations

+
+
=
)
`

(
(
(

2 ln 24
16
211
2 ln 8
6
29
2
81
5
41
5
41
3
5
2
1
a
a
(T3L2-13b)
And solving
138 . 2
1
= a and 348 . 0
2
= a
Hence, substituting in Eqn. (T3L2-10) yields

839 . 4 629 . 4 138 . 2 348 . 0
) 11 )( 1 ( 348 . 0 ) 3 )( 1 ( 138 . 2 ) 1 (
4
1
2 ) ; (
2 3
2
~
+ + =
+ + =
x x x
x x x x x x a x y

(T3L2-14)

There is an important concept associated with a derived term - flux that we have not seen
before and should have. Flux is defined as
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 62

dx
dy
x x = ) ( t (T3L2-15)
Flux has a physical interpretation and we will discuss the concepts in the next topic.
Substituting in Eqn. (T3L2-15), we have

x x x
x x x x x x a x
629 . 4 276 . 4 043 . 1
) 2 )( 2 ( 043 . 1 ) 2 ( 276 . 4 ) 2 (
4
1
2
1
) ; (
2 3
~
+ =
+ + + = t
(T3L2-16)
Finally, note that this solution is called the theoretical (or, classical) solution because of the
manner in which the two boundary conditions were imposed. In this methodology, the
BCs were imposed on the trial solution itself so that the trial solution would satisfy the
BCs exactly. However, this in no way implies that the solution is correct in the interior of
the problem domain.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 63
Some Important Observations

- Why did we start with a trial solution that is a cubic polynomial? Since there are two
BCs that must be imposed, the lowest order polynomial that we could have used would
have been a quadratic polynomial. Imposing the BCs then would have left one free
parameter (or, one DOF). The choice of using a cubic polynomial was an arbitrary
choice and we were left with two free parameters.
- The classical way of applying the BCs can become extremely cumbersome with more
complex problems.
- How do we know whether the solution is good? One way to answer this question is to
look at the concept of convergence. The exact solution to the problem is
x
x
x y ln
2
1 2
) ( + = (T3L2-17a)
and

2
1 2
) ( =
x
x t (T3L2-17b)
Convergence can be checked by starting with a low-order polynomial and increasing
the order of the polynomial gradually. The different trial functions should converge to
a solution. If we had started with a quadratic polynomial, the Galerkins solution
would have been

531 . 3 958 . 1 427 . 0
) 3 )( 1 ( 427 . 0 ) 1 (
4
1
2 ) ; (
2
~
+ =
+ =
x x
x x x a x y
(T3L2-18)
and
x x x x x 958 . 1 854 . 0 ) 2 ( 854 . 0 ) 2 (
4
1
2
1
2
~
+ = + = t (T3L2-19)
In the previous section, with the trial solution as a cubic polynomial the solution was

839 . 4 629 . 4 138 . 2 348 . 0
) 11 )( 1 ( 348 . 0 ) 3 )( 1 ( 138 . 2 ) 1 (
4
1
2 ) ; (
2 3
2
~
+ + =
+ + =
x x x
x x x x x x a x y

and

x x x
x x x x x x a x
629 . 4 276 . 4 043 . 1
) 2 )( 2 ( 043 . 1 ) 2 ( 276 . 4 ) 2 (
4
1
2
1
) ; (
2 3
~
+ =
+ + + = t

What if we had started with a quartic polynomial? The solution is then
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 64
277 . 5 848 . 5 3725 . 3 888 . 0 0864 . 0
) 31 )( 1 ( 0864 . 0
) 11 )( 1 ( 8881 . 0 ) 3 )( 1 ( 3725 . 3 ) 1 (
4
1
2 ) ; (
2 3 4
2 3
2
~
+ + =
+ + +
+ + + =
x x x x
x x x x
x x x x x x a x y
(T3L2-20)
and
x x x x
x x x x x x x x x x a x
848 . 5 745 . 6 664 . 2 346 . 0
) 4 2 )( 2 ( 346 . 0 ) 2 )( 2 ( 664 . 2 ) 2 ( 745 . 6 ) 2 (
4
1
2
1
) ; (
2 3 4
2
~
+ + =
+ + + + + = t

(T3L2-20)
We will now plot the three trial functions and the exact solution both for the function and
the flux.

1.25
1.50
1.75
2.00
1.0 1.2 1.4 1.6 1.8 2.0
0.0864*x*x*x*x-0.888*x*x*x+3.3725*x*x-5.848*x+5.277
-0.348*x*x*x+2.138*x*x-4.629*x+4.839
0.427x*x-1.958*x+3.531
2/x+0.5*ln(x)
x
y
Comparison of Solutions

Fig. T3L2-1 Comparison of the function

While some difference can be noted between the quadratic and the exact solution, there is
very little difference between the cubic, quartic and the exact solutions. Note that all the
solutions have the same function value at the left boundary point 1 = x .
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 65
0.25
0.50
0.75
1.00
1.25
1.50
1.0 1.2 1.4 1.6 1.8 2.0
-0.346*x*x*x*x+2.664*x*x*x-6.745*x*x+5.848*x
1.043*x*x*x-4.276*x*x+4.626*x
-0.854*x*x+1.958*x
(2/x)-0.5
x
t
Comparison of Solutions

Fig. T3L2-2 Comparison of the flux

The differences between the exact solution and the three Galerkins solutions with respect
to the flux is a different matter altogether. The error in the flux from the quadratic
solution is large. However, the error decreases with the cubic and the quartic functions.
Note that all the solutions have the same flux value at the right boundary point 2 = x . The
message here is quite clear (a) small errors in the function do not translate to small errors
in the flux that involve the derivatives of the function, and (b) increasing order
(polynomials) trial functions yield better and converging solutions
6
.

6
Later we will see that the trial functions must have certain properties for this to be true. If we keep on increasing order
of the trial solution will we converge to the exact solution for this problem?
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 66
Summary

The residual function is the function corresponding to the original differential equation
(with all the terms on the LHS) in which an approximate solution is substituted. It
measures how close the approximation is to satisfying the DE but does not tell us how
close the approximation is to the exact solution. The Method of Weighted Residuals
converts the original DE into a set of algebraic equations that are much easier to solve.
There are different approaches used in terms of weighting the residual. Of all the
techniques discussed, the Galerkins Method is by far the most suitable for the type of
problems encountered in engineering analysis.

In the lessons in this topic, we saw how to assume the approximate solution called the trial
solution and use the Galerkins Method to generate the algebraic equations. The trial
solution is usually a polynomial because of the properties that they possess (continuous,
differential, easy to handle etc.). We saw how to enforce the boundary conditions on the
trial solution. We also looked at the flux term. Finally, we also saw how to obtain more
accurate solutions by increasing the order of the polynomial in the trial solution.

There are two problems with this classical (or, theoretical approach). First, the trial
solution is valid for the entire problem domain. Hence it is not possible to accurately
model problems in which there are known discontinuities in the solution and the flux.
Second, the manner in which the boundary conditions are enforced can become
cumbersome for more complex problems. In this next topic, we will see how to overcome
both these drawbacks.


F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 67
Review Exercises
Problem T3L2-1
Consider the differential equation
0
4
2
2
= +
y
dx
y d
t < < x 0
with 1 ) 0 ( = y and 0 ) ( = t y . Find the solution using the Galerkins Method by assuming
the trial solution as
3
4
2
3 2 1
) ; ( x a x a x a a a x y + + + = . Compare with the exact solution.

Problem T3L2-2
Consider the differential equation
( ) x x x x
dx
dy
x
dx
d
110 351 204 30
12
1
2 3 4 2
+ + =
|
.
|

\
|
4 0 < < x
with 1 ) 0 ( = y and 0 ) 4 ( = y . Find the solution using the Galerkins Method by assuming
the trial solution as (i) a quadratic polynomial, and (ii) a cubic polynomial. Compare to the
exact solution.

Problem T3L2-3
Consider the differential equation
( ) 0
) (
1 =
|
.
|

\
|
+
dx
x dy
x
dx
d
2 1 < < x
with 1 ) 1 ( = = x y 1 ) 1 (
2
=
|
.
|

\
|
+
= x
dx
dy
x
(a) Using the Galerkins method with a quadratic polynomial for the trial solution obtain an
approximate solution for both the function and the flux.
(b) Obtain a second approximate solution using a cubic polynomial for the trial solution.
(c) Compare the two solutions with each other and the exact solution.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 68
Topic 4: 1D Problems
""Mach 2 travel feels no different." a passenger commented on an early Concorde
flight. "Yes," Sir George replied. "That was the difficult bit."- Sir George
Edwards, co-director of Concorde development Quoted in Kenneth
Owen, "Concorde, New Shape in the Sky"
Lesson 1: The Element Concept
Objectives: In this lesson we will look at Galerkins Method to solve one-dimensional
boundary-value problems.

- To understand the 1D Boundary Value Problems (BVP).
- To apply the Galerkins Method to solving the 1D BVP.

Suggested Reading
Reference Pages
T1 Chapter 3
T2 Chapter 2
T3 Chapters 4-9

Integration by parts
Evaluate ( ) ( ) I f x g x dx =
}
.
Define ( ) u f x =
( ) dv g x dx =
Compute ( ) du f x dx ' =
( ) v g x dx =
}

Then ( ) ( ) ( ) ( ) I u x dv u x v x v x du = =
} }

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 69
Overcoming Drawbacks
In the previous topic we looked at applying the Galerkins Method in the classical (or,
theoretical) sense. To make the approach general and useful we need to deviate to a new
format as shown below.

Step 1: Let the trial solution be assumed as (Eqn. (T3L1-4))
~
0 1 1 2 2 0
1
( ; ) ( ) ( ) ( ) ...... ( ) ( ) ( )
n
n n j j
j
y x a x a x a x a x x a x
=
= + + + + = +

(T4L1-1)

We will stick with the example from the previous topic. Using the definition of the
residual, we have
2
~
2 ) (
) ; (
x dx
x y d
x
dx
d
a x R
|
|
|
.
|

\
|
= (T4L1-2)
Since we have n parameter trial function, we need n residual equations (see Eqn. (T3L1-
10))
0 ) ( ) ; ( =
}
b
a
x
x
i
dx x a x R | n i ,..., 1 = (T4L1-3)
Notice that we have changed the limits of integration (instead of using 1 and 2 as the
limits) just to make the derivation general and more useful. Substituting, we have
0 ) (
2 ) (
2
~
=
(
(

|
|
|
.
|

\
|
}
b
a
x
x
i
dx x
x dx
x y d
x
dx
d
| n i ,..., 1 = (T4L1-4)

Step 2: Integrate by parts the highest derivative term in the residual equations (which
would be first term containing the second-order derivative). Rearranging the terms, we
have

b
a
b
a
b
a
x
x
i i
x
x
i
x
x
dx
y d
x dx
x
dx
dx
d
dx
y d
x
(
(

|
|
|
.
|

\
|
=
} }
| |
|
~
2
~
2
n i ,..., 1 = (T4L1-5)
Three observations here (a) the highest order derivative in Eqn. (T4L2-5) is now lower
(only first order derivatives compared to second-order in Eqn. (T4L1-4)), and (b) the
stiffness term is on the LHS and the loading terms are on the RHS, and (c) the loading
term contains the boundary flux term.

Step 3: Substituting Eqn. (T4L1-1) in Eqn. (T4L1-5) and noting that

=
+ =
n
i
i
i
dx
x d
a
dx
x d
dx
a x y d
1
0
~
) ( ) ( ) ; ( | |
(T4L1-6)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 70
we have

b
a
b
a
b
a
z
x
i
x
x
i j
n
j
j
x
x
i
dx
y d
x dx
x
a dx
dx
d
x
dx
d
(
(

|
|
|
.
|

\
|
=
|
|
.
|

\
|
}

}
=
| |
|
|
~
2
1
2

dx
dx
d
x
dx
d
b
a
x
x
i 0
| |
}
n i ,.., 1 = (T4L1-7)

Let
dx
dx
d
x
dx
d
K
j
x
x
i
ij
b
a
|
|
}
=
and dx
dx
d
x
dx
d
dx
y d
x dx
x
F
b
a
b
a
x
x
i
i
x
x
i i
0
~
2
2 | |
| |
} }

(
(

|
|
|
.
|

\
|
= (T4L1-8)
Then we can rewrite Eqn. (T4L1-7) in the matrix notation as

(
(
(
(
(
(

n n nn n n
n
n
F
F
F
a
a
a
K K K
K K K
K K K
.
.
.
.
. .
. . . . .
. . . . .
. .
. .
2
1
2
1
2 1
2 22 21
1 12 11
(T4L1-9)
or,

1 1
=
n n n n
F a K (T4L1-10)
The above equations are remarkably similar to Eqns. (T2L2-25a) encountered in the Direct
Stiffness Method. The stiffness matrix K is symmetric since

ij
i
x
x
j
ji
K dx
dx
d
x
dx
d
K
b
a
= =
}
|
|
(T4L1-8b)

Step 4: Use a specific form of the trial solution
In the previous topic we used a polynomial as the trial solution. Polynomials are easy to
deal with and we will stick with that approach here. Let us assume the solution as

=
= + + =
n
j
j j
x a x a x a a y
1
2
3 2 1
~
) ( | (T4L1-11)
implying that
1 ) (
1
= x | x x = ) (
2
|
2
3
) ( x x = | (T4L1-12)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 71
The only difference between Eqn. (T4L1-1) and the above equation is the
0
| term. In this
modified approach we will be applying the BCs numerically and hence there is no need for
the
0
| term. Now we are ready to compute the terms in Eqn. (T4L1-8). Using (T4L1-12)
0
1
=
dx
d|
1
2
=
dx
d|
x
dx
d
2
3
=
|
(T4L1-13)

We will compute a typical term in the stiffness matrix
( )
3 3
23
3
2
) 2 )( )( 1 (
a b
x
x
x x dx x x K
b
a
= =
}
(T4L1-14)
and a force term (in two parts)

a
b
x
x
x
x
dx x
x
F
b
a
ln 2
2
2
int
2
= =
}
(T4L1-15)
and

b
x
a
x
bnd
x
dx
y d
x x
dx
y d
x F
b a
|
|
|
.
|

\
|

|
|
|
.
|

\
|
=
~ ~
2
(T4L1-16)
where
int
2
F is the interior load term and
bnd
F
2
is the boundary load term. Similarly, the flux
can be expressed as

2
3 2
~
~
2 x a x a
dx
y d
x = = t (T4L1-17)

Step 5: Generate the algebraic equations as
( ) ( )
( ) ( )

|
|
.
|

\
|

(
(
(
(
(
(

2
~
2
~
~ ~
~ ~
3
2
1
4 4 3 3
3 3 2 2
| |
| |
| |
) ( 2
ln 2
1 1
2
3
2
0
3
2
2
1
0
0 0 0
b x a x
b x a x
x x
a b
a
b
a b
a b a b
a b a b
x x
x x
x x
x
x
x x
a
a
a
x x x x
x x x x
b a
b a
b a
t t
t t
t t
(T4L1-18)
We can now proceed further by substituting the problem data.

Step 6: Substituting the numerical data.
1 =
a
x 2 =
b
x

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 72

(
(
(
(
(
(

= =
= =
= =
) 4 ( | |
) 2 ( | |
| |
2
2 ln 2
1
15
3
14
0
3
14
2
3
0
0 0 0
2
~
1
~
2
~
1
~
2
~
1
~
3
2
1
x x
x x
x x
a
a
a
t t
t t
t t
(T4L1-19)

Step 7: Applying the BCs
First we will apply the natural boundary condition
2
1
2
=
|
.
|

\
|

= x
dx
dy
x

(
(
(
(
(
(

=
=
=
2 |
1 |
2
1
|
2
2 ln 2
1
15
3
14
0
3
14
2
3
0
0 0 0
1
~
1
~
1
~
3
2
1
x
x
x
a
a
a
t
t
t
(T4L1-20)
Since the NBC is applied to the system equations it does not guarantee that the final
solution will satisfy the NBC. In the classical approach, we enforced the NBC up front on
the trial solution itself so that the final solution did satisfy the NBC.

Applying the EBC 2 ) 1 ( = = x y is different than what we saw in the Direct Stiffness
Method (Topic 2). This is because there is no equation directly associated with ) (x y .
Imposing the EBC on the trial solution itself (Eqn. (T4L1-11))
2
3 2 1
= + + a a a (T4L1-21)
We can write
3
a in terms of the other two parameters as

2 1 3
2 a a a = (T4L1-22)
Eliminating
3
a from Eqns. (T4L1-20) using the above equation, we have

=
)
`

(
(
(
(
(
(

=
=
=
34 |
3
31
2 ln 2 |
2
3
|
15
3
14
15
3
14
2
3
3
14
0 0
1
~
1
~
1
~
2
1
x
x
x
a
a
t
t
t
(T4L1-23)
Eliminating the last equation (from above, Eqn. (1)-Eqn. (3), Eqn.(2)-Eqn. (3))

=
)
`

(
(
(

2 ln 2
3
71
2
65
6
43
3
31
3
31
15
2
1
a
a
(T4L1-24)

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 73
The final equations are symmetric, do not contain any boundary terms and can be solved
numerically.

Step 8: Solving the system equations
719 . 3
1
= a 254 . 2
2
= a (T4L1-25)
Substituting in Eqn. (T4L1-22),
535 . 0
3
= a . (T4L1-26)
Hence

2
~
535 . 0 254 . 2 719 . 3 x x y + = (T4L1-27)
and

2
~
070 . 1 254 . 2 x x = t (T4L1-28)

While we have overcome some obstacles with this procedure, the process has still a major
drawback. We again revisit the issue of imposing the boundary conditions.

The Element Concept (Finally!)
In the preceding section, the trial solution was assumed to be applicable for the entire
problem domain. We will depart from this approach and go through the process of
discretization (or, creating a mesh with elements). This will enable us to develop a very
general methodology making it convenient to do a variety of things including applying the
boundary conditions.

One-Element Solution
Step 4: We will once again stick with the same problem as the last section. Let us assume
that the domain is discretized into a single element. But just as in the Direct Stiffness
Approach, let us now assume a typical element as shown below.

x
y
y
y
x
x
1
2
1
2
1
L
2


Fig. T4L1-1

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 74
The element is described by two nodes that are labeled 1 and 2 and are located at
1
x and
2
x respectively, with the length of the element as L . Let us also assume that the solution
varies linearly over the element and is
1
y at node 1 and
2
y at node 2 noting that at this
stage these values (called nodal values) are unknowns. We can describe the variation of the
solution over the element as a linear interpolation using the two nodal values. In other
words

2 2 1 1 2 1
~
) ( ) ( ) ( y x y x x a a x y | | + = + = (T4L1-29)
This is indeed the trial solution that we have used as the starting point in the preceding
sections. Using the end conditions

1 1
) ( y x x y = =
2 2
) ( y x x y = = (T4L1-30)
and substituting in the first part of Eqn. (T4L1-29) we obtain

1 2 1 1
x a a y + =

2 2 1 2
x a a y + = (T4L1-31)
Solving for the
i
a s, we have

1 2
2 1 1 2
1
x x
y x y x
a

=
1 2
1 2
2
x x
y y
a

= (T4L1-32)
Therefore,
x
x x
y y
x x
y x y x
x y
1 2
1 2
1 2
2 1 1 2
~
) (

= (T4L1-33)
Rearranging

2
1
1
2
2
1 2
1
1
1 2
2
~
) ( y
L
x x
y
L
x x
y
x x
x x
y
x x
x x
x y

+

= (T4L1-34)
The trial functions

L
x x
=
2
1
|
L
x x
1
2

= | (T4L1-35)
are known as the shape functions. They have special properties as shown below.

x
|
| |
1
1
x
x
1
1 2
2
1
L
2


Fig. T4L1-2
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 75

Going back to Step 4 in the preceding section, we have no ) (
0
x |
7
and the
i
| s are given by
Eqn. (T4L1-35). Hence for a typical element

)
`

=
)
`

2
1
2
1
22 21
12 11
F
F
y
y
k k
k k
(T4L1-36)
Let us evaluate couple of typical terms.

L
x x
dx
L
x
L
dx
dx
d
x
dx
d
k
x
x
x
x
2 1 2 1
12
2
1 1
) (
1
2
1
2
1
+
=
|
.
|

\
|
|
.
|

\
|

= =
} }
| |
(T4L1-37)

1
2
1 2 1
1 2
int
1
ln
2 2 2
2
1
x
x
x x x
dx
x
F
x
x

+ = =
}
| (T4L1-38)

b
x
bnd
dx
y d
x F
|
|
|
.
|

\
|
=
~
2
(T4L1-39)
Step 5: With all the other terms evaluated similarly,

( ) ( )
( ) ( )

)

+
=
)
`

+ +
+ +
=
=
2
~
1
~
1
2
2
1
2
1
2
1
2 1 2 1
2 1 2 1
|
|
ln
2 2
ln
2 2
2
1
x
x
x
x
L x
x
x
L x
y
y
x x x x
x x x x
L
t
t
(T4L1-40)

These are the element equations similar to Eqn. (T2L2-4) etc. The flux expression is of the
form
) (
2 1
1 2
~
~
y y
x x
x
dx
y d
x

= = t (T4L1-41)

Step 6: Substituting the numerical values
1
1
= x
2
2 x =
we have

+
)
`

+
=
)
`

=
=
2
~
1
~
2
1
|
|
2 ln 2 1
2 ln 2 2
3 3
3 3
2
1
x
x
y
y
t
t
(T4L1-42)

Step 7: Applying the boundary conditions
First we will apply the natural boundary condition
2
1
2
=
= x
t . This results in

7
It is not necessary to have this term since the BCs will be imposed numerically.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 76

+
)
`

+
=
)
`

=
2
1
|
2 ln 2 1
2 ln 2 2
3 3
3 3
2
1 1
~
2
1
x
y
y
t
(T4L1-43)
Now applying the EBC 2 ) 1 (
1
= = = y x y is done in a manner described in the section
containing Eqn. (T2L2-39). The equations reduce to

=
)
`

(
(
(

2 ln 2
2
7
2
2
3
0
0 1
2
1
y
y
(T4L1-43)

Step 8: Solving
2
1
= y 409 . 1
2
= y (T4L1-44)
Hence, the approximate solution over the element is found by substituting these values in
Eqns. (T4L1-34) and (T4L1-41) yielding
x y 591 . 0 591 . 2
~
= (T4L1-45)
and
x 591 . 0
~
= t (T4L1-46)

Comparing this solution to the exact solution shows that the results are quite in error. This
is because of the linear trial solution that was assumed and because of the fact that we used
only one element.

Two-Element Solution
The finite element mesh for the two-element solution is shown below.

x
1
1
2
3
y
y
y
1 2
x = 1 x = 1.5
x = 2
3 2


Fig. T4L1-3

This is a uniform mesh since both the elements are geometrically identical. The unknowns
at the three nodes are labeled
3 2 1
, , y y y . Just as in the Direct Stiffness Method, we will
generate the element equations.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 77
Element 1: 1
1
= x and 5 . 1
2
= x

( ) ( )
( ) ( )

|
.
|

\
|

|
.
|

\
|
+

+
=
)
`

+ +
+ +
=
=
1
5 . 1
~
1
1
~
2
1
|
|
1
5 . 1
ln
5 . 0
2
5 . 1
2
1
5 . 1
ln
5 . 0
2
1
2
5 . 1 1 5 . 1 1
5 . 1 1 5 . 1 1
) 5 . 0 ( 2
1
x
x
y
y
t
t
(T4L1-47)

Element 2: 5 . 1
1
= x and
2
2 x =

( ) ( )
( ) ( )

|
.
|

\
|

|
.
|

\
|
+

+
=
)
`

+ +
+ +
=
=
2
2
~
2
5 . 1
~
3
2
|
|
5 . 1
2
ln
5 . 0
2
2
2
5 . 1
2
ln
5 . 0
2
5 . 1
2
2 5 . 1 2 5 . 1
2 5 . 1 2 5 . 1
) 5 . 0 ( 2
1
x
x
y
y
t
t
(T4L1-48)

We need to expand the notation to differentiate between the two elements, i.e.
1
1
~
|
|
.
|

\
|
= x
t
represents the flux at 1 = x for element 1. Assembling the two equations to create the
system equations, we obtain

|
.
|

\
|

|
.
|

\
|
+

+
=

(
(
(

=
=
2
2
~
1
1
~
3
2
1
|
0
|
3
4
ln 4 1
9
8
ln 4
2
3
ln 4 2
5 . 3 5 . 3 0
5 . 3 0 . 6 5 . 2
0 5 . 2 5 . 2
x
x
y
y
y
t
t
(T4L1-49)
Note that the second term in the boundary load vector is zero. This is because when we
assemble the system equations we assume that

2
5 . 1
~
1
5 . 1
~
| |
|
.
|

\
|
=
|
.
|

\
|
= = x x
t t (T4L1-50)
In other words the flux is assumed to be continuous across the two elements. As before,
since this constraint is not enforced at the system level (but a mere substitution is made for
the flux on the RHS), the final solution will not show this flux continuity across the
elements
8
.

Now we need to impose the boundary conditions. Imposing the NBC first, we have

8
This is similar to the imposition of the natural boundary condition at the system level.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 78

+
=

(
(
(

=
2
1
0
|
3
4
ln 4 1
9
8
ln 4
2
3
ln 4 2
5 . 3 5 . 3 0
5 . 3 0 . 6 5 . 2
0 5 . 2 5 . 2
1
~
3
2
1
x
y
y
y
t
(T4L1-51)
Now imposing the EBC yields

1
2
3
2
1 0 0
8
0 6.0 3.5 4ln 5
9
0 3.5 3.5
1 4
4ln
2 3
y
y
y


(


(
= +
` `
(

(
)

)
(T4L1-52)

Solving the equations yields
2
1
= y 551 . 1
2
= y 365 . 1
3
= y (T4L1-53)

Hence, for element 1 (using Eqns. (T4L1-34) and (T4L1-41))
898 . 2 898 . 0
5 . 0
1
551 . 1
5 . 0
5 . 1
2 ) (
1
~
+ =
|
.
|

\
|

+
|
.
|

\
|

=
|
.
|

\
|
x
x x
x y (T4L1-54)
x x 898 . 0 ) (
1
~
=
|
.
|

\
|
t (T4L1-55)
and for element 2
109 . 2 372 . 0
5 . 0
5 . 1
365 . 1
5 . 0
0 . 2
551 . 1 ) (
2
~
+ =
|
.
|

\
|

+
|
.
|

\
|

=
|
.
|

\
|
x
x x
x y (T4L1-54)
x x 372 . 0 ) (
2
~
=
|
.
|

\
|
t (T4L1-55)

Lets compare the one-element and the two-element solutions.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 79
0.5
1.0
1.5
2.0
1.0 1.2 1.4 1.6 1.8 2.0
0.898*x
0.372*x
-0.372*x+2.109
-0.898*x+2.898
0.591*x
2.591-0.591*x
x
y
(
x
)

a
n
d

t
(
x
)
Comparison of the Solutions

Fig. T4L1-4

(1) Both the solutions satisfy the EBC at . 1 = x They should, as we have imposed the EBC
on the final system equations.
(2) There is an improvement in the two-element solution with respect to ) (x y (An error
of 4.6% and 1.3% for the one-element and two-element solutions, respectively at
2 = x ).
(3) The error in the one-element flux is large throughout the domain. The error in the flux
from the two-element solution is much less (An error of 136% and 49% for the one-
element and two-element solutions, respectively at 2 = x ).
(4) While the function is continuous throughout the domain, the flux is discontinuous for
the two-element model at the element interface at 5 . 1 = x .

Review of the Steps
Step 1: Assume the trial solution of the form

=
= + + + =
n
i
i i
n
n
x y x a x a a a x y
1
1 0
~
) ( .... ) ; ( |
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 80
where n is equal to the number of DOF in the element,
i
y are the nodal values and ) (x
i
|
are the shape functions. Construct the residual and substitute the trial solution in the
residual equations. Note that there are as many residual equations as there are DOF in the
element.

Step 2: Integrate by parts the highest derivative term. Integrating by parts will not only
lower the highest derivative term but also generate a boundary (flux) term.

Step 3: Rewrite the equations so that the stiffness related terms are on the left and the load
terms (interior and boundary) are on the right.

Step 4: To generate the equations completely we need to assume the exact form of the trial
solution, i.e. fix the value of n . This will enable us to generate the terms in the stiffness
matrix and the load vector. We can now generate the element equations in the form

1 1
=
n n n n
f u k
We also need to generate the expression for the flux.

Step 5: Using the problem data, we can generate the element equations for all the elements
in the model. These equations are then assembled into the system or global equations of
the form

1 1
=
N N N N
F U K
for a model with N system DOF.

Step 6: Using the problem data, we impose the NBCs first. Then we impose the EBCs
resulting in equations of the form

1 1
=
N N N N
F U K

Step 7: Solve the system equations for the primary nodal unknowns U.

Step 8: Using the primary unknowns we can compute the flux in each element.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 81
Review Exercises
In the following problems the quadratic element is to be used. To obtain the shape functions
for the quadratic element refer to Lesson 4.
Problem T4L1-1
Consider the differential equation
( ) 0
) (
1 =
|
.
|

\
|
+
dx
x dy
x
dx
d
2 1 < < x
with 1 ) 1 ( = = x y 1 ) 1 (
2
=
|
.
|

\
|
+
= x
dx
dy
x
(a) Using the element concept, use the linear polynomial as the trial solution. Apply the
steps outlined in this section to obtain an approximate solution for both the function
and the flux.
(b) Repeat the problem but now use a quadratic polynomial. Compare the two solutions.
(c) Are these solutions different than those obtained in Problem T3L2-3?

Problem T4L1-2
Consider the differential equation
( ) x x x x
dx
dy
x
dx
d
110 351 204 30
12
1
2 3 4 2
+ + =
|
.
|

\
|
4 0 < < x
with 1 ) 0 ( = y and 0 ) 4 ( = y . Using the element concept, use a quadratic polynomial as the
trial solution. Obtain an approximate solution for both the function and the flux.
Compare this to the solution obtained for Problem T3L2-2. (Hint: While a linear element
goes with a linear polynomial and requires an element with 2 nodes, a quadratic element
requires an element with 3 nodes. See Lesson 4 for details.)

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 82
One-Dimensional Boundary Value Problem

The one-dimensional BVP is described by (with the positive x direction left to right)
DE: ) ( ) ( ) (
) (
) ( x f x y x
dx
x dy
x
dx
d
= +
|
.
|

\
|
| o
b a
x x x < < (T4L1-56)
BCs: At
a
x x = either
a
y y = or
a a
c y d t = +
At
b
x x = either
b
y y = or
b b
c y d t = + (T4L1-57)

where ) (x o , ) (x | , and ) (x f are known functions,
a
y ,
b
y ,
a
c ,
b
c ,
a
d and
b
d are
constants, and
dx
dy
o t = is the flux. The BCs are EBC, NBC or mixed. Note that the
mixed BC reduces to a NBC by setting 0
a
c = and 0
b
c = . We will look at specific
engineering problems that are governed by this differential equation in the next lesson.

We will use the Galerkins Method to generate the element equations. The following steps
pertain to a typical element in the mesh.

Step 1: Assume the trial solution as

=
=
n
j
j j
x y a x y
1
~
) ( ) ; ( | (T4L1-58)
As in the previous section, we do not have the ) (
0
x | term since the boundary conditions
will be imposed numerically. We will drop the ~(tilde) notation for convenience sake.
Substituting in the residual equations and integrating over the domain O of the element,
we have
0 ) ( ) ( ) ( ) (
) (
) ( =
(

+
|
.
|

\
|

}
O
dx x x f x y x
dx
x dy
x
dx
d
i
| | o n i ,.., 2 , 1 = (T4L1-59)

Step 2: Integrating by parts the highest-order derivative, we have
| |
I
O O
=
(

+
} }
i i i
i
dx x f dx x y x
dx
d
dx
dy
x | t | | |
|
o ) ( ) ( ) ( ) ( n i ,.., 2 , 1 = (T4L1-60)
where I represents the boundary of the element.
Step 3: Using Eqn. (T4L1-58) in (T4L1-60)
=
(

} }
=
O O
j
n
j
j i
j
i
y dx x x x dx
dx
d
x
dx
d
1
) ( ) ( ) ( ) ( | | |
|
o
|

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 83
| |
I
O
}

i i
dx x x f t| | ) ( ) ( n i , ., 2 , 1 = (T4L1-61)
Step 4: Let us use the linear interpolation two-node element from the earlier section ) 2 ( = n
L
x x
=
2
1
|
L
x x
1
2

= | (T4L1-62)
Substituting Eqn. (T4L1-62) in Eqn. (T4L1-61), we have the element equations as

)
`

+
)
`

=
)
`

bnd
bnd
F
F
F
F
y
y
k k
k k
2
1
int
2
int
1
2
1
22 21
12 11
(T4L1-63)
Let us look at a typical stiffness term first.
dx
x x
x x
x
x x
x x
dx
x x
x
x x
k
x
x
x
x
|
|
.
|

\
|

|
|
.
|

\
|

+
|
|
.
|

\
|

|
|
.
|

\
|

=
} }
1 2
2
1 2
2
1 2 1 2
11
) (
1
) (
1
2
1
2
1
| o (T4L1-64)
To evaluate the above equation we need to know ) (x o and ) (x | . The functions can then be
substituted and the integral evaluated. Instead, we will assume that we know the numerical
value of these two functions at the centroid of the element, i.e. at
2
2 1
x x
x
+
= . For this element
in which the solution is assumed to vary linearly, this is the most accurate point to evaluate the
integral. Let us denote the centroidal values (constants) as o and | . Now, integrating

3
11
L
L
k
| o
+ = (T4L1-65)
We can use a similar strategy with the load terms. When all the terms are evaluated we have,

| |
| |

)

(
(
(
(

+ +
+ +
I
I
2
1
2
1
2
2
3 6
6 3
t|
t|
| o | o
| o | o
L f
L f
y
y
L
L
L
L
L
L
L
L
(T4L1-66)
The term that requires special treatment here is the last term. In general
| | | | | |
1 2
x i x i i
t| t| t| =
I
(T4L1-67)
From Eqn. (T4L1-57), we have
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 84
( ) cy d t = + for
a
x x = and ( ) cy d t = + for
b
x x = .
Substituting this in Eqn. (T4L1-67) for 1 = i
| | | | | | | | ( )
2 1 1 1
1 1 1 1 1 1
0
x x x x
cy d c y d t t t t
I
= = = + = (T4L1-68)
Similarly, for 2 = i
| | | | | | | | ( )
2 1 2 2
2 2 2 2 2 2
0
x x x x
cy d c y d t t t t
I
= = = + = + (T4L1-69)
Note that the unknown y appears in the two equations and must be moved to the LHS.
Therefore the modified element equations are
1
1 2
2
1 0 0 0
3 6
0 0 0 1
6 3
L L
y
L L
c c
y
L L
L L
o | o |
o | o |
( (
+ +
( (
( (

( (
+ =
` ( (
( (

)
+ +
( (



1
2
2
2
f L
d
d
f L




+
` `

)

)
(T4L1-70)
The
1
c term exists provided
a
x x =
1
. Similarly,
2
c term exists provided
b
x x =
2
. The flux at
the center of the element is given by
( )
1 2
y y
L dx
dy
= =
o
o t (T4L1-71)
The element equations are ready and we now need to look at engineering problems that are
described as 1D BVP to illustrate these steps and the rest of the solution.
Finally a warning boundary conditions can be dangerous to your health! Applying the BCs
incorrectly is one of the most common forms of errors.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 85
Lesson 2: Solid Mechanics Examples
Objectives: In this lesson we will look at Galerkins Method to solve one-dimensional
solid mechanics problems.

- To understand the relationship between some of the solid mechanics problems and the
general 1D BVP discussed in the previous lesson.
- To solve an example so as to understand the process.

Suggested Reading
Reference Pages
T1 Chapter 3
T2 Chapter 2
T3 Chapter 6

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 86
Axial Deformation of an Elastic Rod

Consider the elastic rod as shown in figure below. The axial displacement is denoted
) (x u u = and the axial loading (body force) on the rod is ) (x w w = . A sample set of
boundary conditions is shown where 0 = u (EBC) on the left end and the force F (NBC)
on the right end is zero.

u, x
w
u=0
F=0


Fig. T4L2-1

The governing differential equation is given as
) ( ) (
) (
) ( ) ( x A x w
dx
x du
x E x A
dx
d
=
|
.
|

\
|
(T4L2-1)
with the boundary conditions as
At
a
x x = ,
a a
u x x u = = ) ( or NBC (T4L2-2)
At
b
x x = ,
b b
u x x u = = ) ( or NBC (T4L2-3)
The NBC is of the form

x x
n X o = (T4L2-3a)
where
x
n is the direction cosine of the outward normal, X is the force per unit area and is
positive if it acts in the positive x direction. Let us look at some possibilities with respect
to the boundary conditions.
Rod has a known displacement
a
u (incl. zero) at the left end

a
u u =
There is a concentrated force
a
F applied at the left end in the positive x direction ( ) 1
x
n =
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 87

a
x x
a
dx
du
AE F
=
= or
a
a
x x
du
E
dx
o
=
=
There is a concentrated force
b
F applied at the right end in the positive x direction ( ) 1
x
n =

b
x x
b
dx
du
AE F
=
= or
b
b
x x
du
E
dx
o
=
=

Using the process discussed in the previous lesson, Eqn. (T4L1-70) reduces to

(
(
(

2
1
2
1
2
2
1 1
1 1
F
F
L wA
L wA
u
u
L
AE
(T4L2-4)
1 2 1 2 2 2
= f u k (T4L2-5)
A dimensional analysis will show that -
2
L A ,
2
F
L
E , L u and
L
F
wA . The LHS and
the RHS have units of F . While we will look formally at systems modeled with discrete
structural elements (truss and frame) in the second module, it should be noted that the
stiffness matrix (in the local coordinate system aligned with the axis of the element) is the
stiffness matrix for a truss element provided
(a) the end of the element are pin connections,
(b) the cross-sectional area and the modulus of elasticity are constant within the element,
and
(c) 0 ) ( = x f since no element loads are allowed on a truss member.

An Illustrative Example T4L2-1

Figure shows a bar made of a material with modulus of elasticity of 200 GPa. Compute the
nodal displacements, element stresses and support reactions.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 88
150 mm 150 mm
X
2
2
P= 300 kN
250 mm
400 mm
300 mm


Fig. T4L2-2a
Solution: Let us use m and N as the problem units. Let us use a three-element model
placing a node where the concentrated force acts
9
. The FE model is shown below. We have
EBCs at the left and the right ends.

2
2 3
4 3 1
1
0.3 m 0.15 m 0.15 m
X


Fig. T4L2-2b

Element 1:
2
9
) 10 ( 200
m
N
E = ,
2 6
) 10 ( 250 m A

= , m L 15 . 0 = , 0 = w . Hence the element
equations are

(
(
(

1
2
1
1
2
1
8 8
8 8
) 10 ( 333 . 3 ) 10 ( 333 . 3
) 10 ( 333 . 3 ) 10 ( 333 . 3
F
F
U
U

Element 2:
2
9
) 10 ( 200
m
N
E = ,
2 6
) 10 ( 250 m A

= , m L 15 . 0 = , 0 = w . Hence the element
equations are

(
(
(

2
2
2
1
3
2
8 8
8 8
) 10 ( 333 . 3 ) 10 ( 333 . 3
) 10 ( 333 . 3 ) 10 ( 333 . 3
F
F
U
U


9
This is not necessary since we can use a Dirac Delta function ) (x w to define the concentrated force and then compute
the equivalent forces acting at the nodes of the element.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 89
Element 3:
2
9
) 10 ( 200
m
N
E = ,
2 6
) 10 ( 400 m A

= , m L 30 . 0 = , 0 = w . Hence the element
equations are

(
(
(

3
2
3
1
4
3
8 8
8 8
) 10 ( 667 . 2 ) 10 ( 667 . 2
) 10 ( 667 . 2 ) 10 ( 667 . 2
F
F
U
U


Assembling the equations, we have

(
(
(
(
(
(
(
(
(

3
2
3
1
1
4
3
2
1
8
0
) 10 ( 300
667 . 6 667 . 6 0 0
667 . 6 6 333 . 3 0
0 333 . 3 667 . 6 333 . 3
0 0 333 . 3 333 . 3
10
F
F
U
U
U
U

Imposing the boundary conditions, we have

(
(
(
(
(
(
(
(
(

0
0
) 10 ( 300
0
1 0 0 0
0 6 333 . 3 0
0 333 . 3 667 . 6 0
0 0 0 1
10
3
4
3
2
1
8
U
U
U
U

Solving, the nodal displacements are
{ } { } m U U U U ) 10 ( 0 , 46 . 3 , 23 . 6 , 0 , , ,
4
4 3 2 1

=

The force in each element is computed as follows. Note that flux is negative of the member
force.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 90
Element 1
N U U
L
E A
F ) 10 ( 077 . 2 ) 10 )( 23 . 6 (
15 . 0
) 10 ( 5
) (
5 4
7
1 2
1
1 1
1
= = =

(Tension)
This should also be equal and opposite to the support reaction at the left end, i.e.
= N R
left
) 10 ( 077 . 2
5
.
Element 2
N U U
L
E A
F ) 10 ( 2333 . 9 ) 10 )( 23 . 6 46 . 3 (
15 . 0
) 10 ( 5
) (
4 4
7
2 3
2
2 2
2
= = =

(Compression)
Element 3
N U U
L
E A
F ) 10 ( 2333 . 9 ) 10 )( 46 . 3 0 (
30 . 0
) 10 ( 8
) (
4 4
7
3 4
3
3 3
3
= = =

(Compression)
This should also be equal and opposite to the support reaction at the right end, i.e.
= N R
right
) 10 ( 2333 . 9
4
.
We can now examine the results. The displacements satisfy the EBC at the left and right ends.
The force equilibrium is satisfied since the sum of the two support reactions is equal to the
applied force. Also the forces in elements 2 and 3 are equal.
Thermal Loads Example T4L2-2

Consider a bar of length L , constant cross-section A, modulus of elasticity E and
coefficient of thermal expansion, o . Let the temperature change in the bar be T A . The
initial strain is given by


0
T c o = A (ET4L2-2-1)

The equivalent nodal loads due to the temperature change in the element for the 1D-C
0

linear element is given by


1
2 1
2
1
1
f
AE T
f
o

= = A
` `
) )
f (ET4L2-2-2)
The resulting load vector is shown in Fig. ET4L2-2-1.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 91

f
1
f
2

Fig. ET4L2-2-1

Example (a)
Fig. ET4L2-2-2 shows a bar made of steel. The length of the bar is 2 m. The cross-sectional
area is
2
0.001m . The bar is initially at 25 C

. The temperature is increased to 100 C

.
Compute the stress in the bar.


Fig. ET4L2-2-2

Solution: We will solve this problem using a one element mesh with the linear element.
The element equations are as follows.


( ) ( )
( ) ( )( ) ( )
9
1 9 6
2
0.001 200 10 1 1 1
0.001 200 10 11.7 10 75
1 1 1 2
D
D

(
=
` `
(

) )

or
8 8
1
8 8
2
175500 10 10
175500 10 10
D
D
(
=
` `
(

) )


Since these are the system equations, applying the essential BC at the left end, we have


8
2 2
10 175500 0.001755m D D = =

This indicates that the bar is expanding, as it should! Now

2 1
0.001755
0.0008775
2
D D
L
c

= = =

( )( )
6
0
11.7 10 75 0.0008775 c

= =
( ) ( )
9
0
200 10 0.0008775 0.0008775 0 E o c c = = =
Since the bar is unstrained (free to expand), there is no stress in the bar.

Example (b)
Fig. ET4L2-2-3 shows a bar made of steel. The length of the bar is 2 m. The cross-sectional
area is
2
0.001m . The bar is initially at 25 C

. The temperature is increased to 100 C

.
Compute the stress in the bar.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 92


Fig. ET4L2-2-3

Solution: We will solve this problem using a two element mesh with the linear elements.
The element equations are as follows.

Element 1


( ) ( )
( ) ( )( ) ( )
9
1 9 6
2
0.001 200 10 1 1 1
0.001 200 10 11.7 10 75
1 1 1 1
D
D

(
=
` `
(

) )

or
( )
( )
8 8
1
8 8
2
175500 2 10 (2)10
175500 2 10 (2)10
D
D
(
=
` ` (

) )


Element 2


( ) ( )
( ) ( )( ) ( )
9
2 9 6
3
0.001 200 10 1 1 1
0.001 200 10 11.7 10 75
1 1 1 1
D
D

(
=
` `
(

) )

or
( )
( )
8 8
2
8 8
3
175500 2 10 (2)10
175500 2 10 (2)10
D
D
(
=
` ` (

) )



Assembling the equations, applying the essential BC at the left and right ends, we have


8
2 2
(4)10 0 0 D D = =

This indicates that the bar is neither expanding nor contracting, as it should!

Element 1

2 1
0
D D
L
c

= =

( )( )
6
0
11.7 10 75 0.0008775 c

= =
( ) ( )
9 8
0
200 10 0 0.0008775 1.755 10 175 ( ) E Pa MPa C o c c o = = = =
( )( )
8
1.755 10 0.001 175500 F A N o = = =

Element 2
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 93

3 2
0
D D
L
c

= =

( )( )
6
0
11.7 10 75 0.0008775 c

= =
( ) ( )
9 8
0
200 10 0 0.0008775 1.755 10 175 ( ) E Pa MPa C o c c o = = = =
( )( )
8
1.755 10 0.001 175500 F A N o = = =

Since the bar is completely restrained (prevented from expanding), the bar is in
compression as illustrated by the stress in both bars. The FBDs of the elements and node 2
are shown in Fig. ET4L2-2-4.

175500 175500
1
175500 175500
2
175500
175500
2


Fig. ET4L2-2-4 Element and nodal FBDs

Example (c)
Fig. ET4L2-2-5 shows a bar made of steel and aluminum of equal lengths. The total length
of the bar is 2 m. The cross-sectional area is
2
0.001m . The entire bar is initially at 25 C

.
The temperature is uniformly increased to 100 C

. Compute the stress distribution in the

bar.

Steel
Aluminum

Fig. ET4L2-2-5

Solution: We will solve this problem using a two element mesh with the linear elements.
The element equations are as follows.

Element 1 (Steel)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 94


( ) ( )
( ) ( )( ) ( )
9
1 9 6
2
0.001 200 10 1 1 1
0.001 200 10 11.7 10 75
1 1 1 1
D
D

(
=
` `
(

) )

or
( )
( )
8 8
1
8 8
2
175500 2 10 (2)10
175500 2 10 (2)10
D
D
(
=
` ` (

) )



Element 2 (Aluminum)

( ) ( )
( ) ( )( )( )
9
2 9 6
3
0.001 70 10 1 1 1
0.001 70 10 23.0 10 75
1 1 1 1
D
D

(
=
` `
(

) )

or
( )
( )
8 8
2
8 8
3
120750 0.7 10 (0.7)10
120750 0.7 10 (0.7)10
D
D
(
=
` ` (

) )



Assembling the equations, applying the essential BC at the left and right ends, we have


8
2 2
(2.7)10 175500 120750 54750 0.000202778 D D m = = =

This indicates that the steel part of the bar is expanding and the aluminum part of the bar
is contracting.


Element 1 (Steel)

2 1
0.000202778
D D
L
c

= =

( )( )
6
0
11.7 10 75 0.0008775 c

= =
( ) ( )
9 8
0
200 10 0.000202778 0.0008775 1.349 10 135 ( ) E Pa MPa C o c c o = = = =

( )( )
8
1.349 10 0.001 134900 F A N o = = =

Element 2 (Aluminum)

3 2
0.000202778
D D
L
c

= =

( ) ( )
6
0
23.0 10 75 0.001725 c

= =
( ) ( )
9 8
0
70 10 0.000202778 0.001725 1.349 10 135 ( ) E Pa MPa C o c c o = = = =

( )( )
8
1.349 10 0.001 134900 F A N o = = =
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 95

Since the bar is completely restrained (prevented from expanding), the bar is in
compression as illustrated by the stress in both bars. The FBDs of the elements and node 2
are shown in Fig. ET4L2-2-6.

134900 134900
1
134900 134900
2
134900
134900
2


Fig. ET4L2-2-6 Element and nodal FBDs

Note that the support reaction at the left end is ( ) 134900N and the support reaction at
the right end is ( ) 134900N .



F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 96
Lesson 3: Heat Transfer and Electrostatic Examples
Objectives: In this lesson we will look at Galerkins Method to solve one-dimensional heat
transfer and fluid flow problems.

- To understand the relationship between some of the heat transfer and fluid flow
problems and the general 1D BVP discussed in the previous lesson.
- To solve an example so as to understand the process.

Suggested Reading
Reference Pages
T1 Sections 10.1 and 10.2
T2 Chapter 2
T3 Chapter 6


F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 97
One-Dimensional Heat Conduction and Convection Problem

Consider a long, thin rod as shown in figure below. The temperature ) (x T T = , ) (x q q =
is the heat flux, ) (x k k = is the thermal conductivity, ) (x Q Q = is the interior volume heat
source, ) (x A A = is the cross-sectional area, ) (x h h = is the convective heat transfer
coefficient, ) (x l l = is the circumference,

T is the ambient temperature. A sample set of

boundary conditions is shown with an NBC on the left end and EBC on the right end. We
could also have a case involving a mixed BC on the left or right ends.

Convective heat loss
dx
x
L
T

q
Q
A
A+ dA T= T
0
perimeter, l
q
a
h
q+ dq

Fig. T4L3-1


The governing differential equation is given as

+ = +
|
.
|

\
|
T x l x h x A x Q x T x l x h
dx
x dT
x k x A
dx
d
) ( ) ( ) ( ) ( ) ( ) ( ) (
) (
) ( ) ( (T4L3-1a)
or, for a cylindrical rod with a constant cross-sectional area

+ = +
|
.
|

\
|
T
A
hl
x Q x T
A
hl
dx
x dT
x k
dx
d
) ( ) (
) (
) ( (T4L3-1a)
with the boundary conditions as
At
a
x x = ,
a
T T = or
a
c q = or ) (

=
a a
T T h q (T4L3-2a)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 98
At
b
x x = ,
b
T T = or
b
c q = or ) (

=
b b
T T h q (T4L3-2b)
Looking ahead to two and three-dimensional problems, these boundary conditions are
special cases of the general form (for specified heat flow)

S z z y y x x
q n q n q n q = + + (T4L3-3a)
if heat
S
q is flowing into the surface S , and ( )
z y x
n n n , , are the direction cosines of the
outward normal from the surface. Similarly, for free convection from surface S , we have
) (

= + + T T h n q n q n q
S z z y y x x
(T4L3-3b)

Possible boundary conditions at the left end ( ) 1
x
n =
Left end is at a known temperature,
a
T T =

a
T T =
Heat (
a
q ) is flowing into the left end

a
q q =
Left end is insulated
0 = q
Free convection is taking place at the left end (ambient temperature is
a
T

and the convective

coef is
a
h ,
a
T T

> )

=
a a a
T h T h q or

+ =
a a a
T h T h q (T3L3-3c)

Possible boundary conditions at the right end ( ) 1
x
n =
Right end is at a known temperature,
b
T T =

b
T T =
Heat (
b
q ) is flowing out of the right end

b
q q =
Right end is insulated
0 = q
Free convection is taking place at the right end (ambient temperature is
b
T

and the convective

coef is
b
h ,
b
T T

> )

=
b b b
T h T h q (T3L3-3d)

Comparing these equations to the general form of the 1D BVP
10
,
) ( ) ( x T x y = ) ( ) ( x k x = o
A
hl
x = ) ( | (T4L3-4a)

10
The sign convention is as follows energy (or, heat flow) into the surface or boundary is positive.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 99

+ T
A
hl
x Q x f ) ( ) ( q
dx
dT
k = = t (T4L3-4b)
Compare
a a
c y d t = + and

+ =
a a a
T h T h q
and
b b
c y d t = + and

=
b b b
T h T h q , and we have
1 a a
c h c = =
2 b b
c h c = =
1 a a a
d h T d

= =
2 b b b
d h T d

= = (T4L3-4c)

Using the natural and mixed boundary conditions as discussed above, and substituting the
above in to Eqn. (T4L1-70), we have
1
1 2
2
1 0 0 0
3 6
0 0 0 1
6 3
k hlL k hlL
T
L A L A
h h
T
k hlL k hlL
L A L A
(
(
+ +
(
(
( (

(
(
+ + =
` ( (
( (

)
+ +
(
(
(

+
+

2
2
1
1
2
1
2
T h
T h
q
q
T
A
hl
Q
T
A
hl
Q
L
(T4L3-5a)
1 2 1 2 2 2
= f u k (T4L3-5b)
Note that on the LHS and the RHS some of the components will be zero depending on
whether the boundary condition is NBC or mixed. A dimensional analysis will show that -
T T ,
2
L A ,
tLT
E
k ,
T tL
E
2
h , L = l ,
3
tL
E
Q ,
2
tL
E
c , the LHS and the RHS have
the units
2
tL
E
.
Illustrative Example T4L3-1
The figure shows a composite wall made of three materials. The outer temperature is C

20 .
Convection heat transfer place on the inner surface of the wall with C T

800 =

and
C m
W
h

=
2
25 . The thermal conductivities are
C m
W
k
C m
W
k
C m
W
k

= 50 , 30 , 20
3 2 1
.
We need to determine the temperature distribution in the wall.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 100
0.3 m 0.15 m
0.15 m
T
k k k
h, T
0
0 0
1 2 3
0
= 20 C

Fig. T4L3-2a
Solution: We will use a three-element model as shown. Note the following:
(a) This is a problem where the convective heat exchange (gain) takes place from the left end.
The mixed boundary condition can be expressed as
( )
x x a a a
q n q h T T

= = . The leading
minus sign on the right is because heat is flowing into the surface. Or,
a a a a
q h T h T

= + (same
as before).
(b) There is no convective heat exchange from the top and bottom of the wall that are assumed
to be very tall compared to the thickness of the wall. Hence, l h = = 0 . We will assume a unit
area for all computations, i.e
2
1m A = . There is no internal heat generation, i.e. 0 = Q .
(c) The right end is at a specified temperature (EBC). For the sake of convenience, we will not
include the boundary flux load terms by assuming inter-element flux continuity.

2
2 3
4 3 1
1
0.3 m 0.15 m 0.15 m
X

Fig. T4L3-2b
Element 1:
C m
W
k

= 20 , m L 3 . 0 = ,
2
1m A = , 0 = h , 0 = l ,
1 2
25
W
h
m C
=

, C T

800
1
=

,
0
2
= h , 0 = Q , 0
1
= c and 0
2
= c . The element equations are as follows.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 101

(
(
(
(

+
0
) 800 ( 25
3 . 0
20
3 . 0
20
3 . 0
20
25
3 . 0
20
2
1
T
T

Element 2:
C m
W
k

= 30 , m L 15 . 0 = ,
2
1m A = , 0 = h , 0 = l ,
1
0 h = , 0
2
= h , 0 = Q , 0
1
= c
and 0
2
= c . The element equations are as follows.

(
(
(
(

0
0
15 . 0
30
15 . 0
30
15 . 0
30
15 . 0
30
3
2
T
T

Element 3:
C m
W
k

= 50 , m L 15 . 0 = ,
2
1m A = , 0 = h , 0 = l ,
1
0 h = , 0
2
= h , 0 = Q , 0
1
= c
and 0
2
= c . The element equations are as follows.

(
(
(
(

0
0
15 . 0
50
15 . 0
50
15 . 0
50
15 . 0
50
4
3
T
T

Assembling the equations

(
(
(
(
(
(
(
(
(

0
0
0
000 , 20
3333 . 333 3333 . 333 0 0
3333 . 333 3333 . 533 0 . 200 0
0 0 . 200 6667 . 266 6667 . 66
0 0 6667 . 66 6667 . 91
4
3
2
1
T
T
T
T

There is no natural boundary condition (the mixed BC was taken care of when the element
equations were generated). To enforce the EBC 20
4
= T , we use the elimination approach. The
modified equations are
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 102
1
2
3
4
91.6667 66.6667 0 0 20, 000
66.6667 266.6667 200.0 0 0
0 200.0 533.3333 0 6666.6667
0 0 0 1 20
T
T
T
T
(
(


(
=
` `
(

(

) )

Solving the equations we have (note decreasing temperature from node 1 to 4),
{ } { } C T T T T

20 , 14 . 57 , 05 . 119 , 76 . 304 , , ,
4 3 2 1
=
We can now compute the flux in each element.
Element 1:

2 1 2
1
1 1
12380.95 ) 76 . 304 05 . 119 (
3 . 0
20
) (
m
W
T T
L
k
= = = t
Element 2:

2 2 3
2
2 2
12380.95 ) 05 . 119 14 . 57 (
15 . 0
30
) (
m
W
T T
L
k
= = = t
Element 3:

2 3 4
3
3 3
12380.95 ) 1 . 57 0 . 20 (
15 . 0
50
) (
m
W
T T
L
k
= = = t
We will now examine the results. The solution satisfies the EBC at the right end. The flux is
constant throughout the domain as it should be.

Another Illustrative Example T4L3-2
Fig. T4L3-3 shows a circular cross-section pin fin. It has a diameter of 0.3125 and a length of
5. At the root, the temperature is F T

150
0
= . The ambient temperature is F T

80 =

, the
convective coefficient is
F ft h
BTU
h

=
2
6 , and the thermal conductivity is
F ft h
BTU
k

= 8 . 24 .
Determine the temperature distribution in the fin.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 103
d
0
L
T

Fig. T4L3-3
Solution: We will use a two-element model to solve the problem. The FE model is shown
below.
2
2
3 1
1
2.5 2.5
X

Fig. T4L3-4
Note the following - (a) In this problem the left end is tied to an EBC. (b) There is no
convective heat exchange from the right end (NBC with 0 = q ). (c) There is no internal heat
generation, i.e. 0 = Q . We will select ft as the units for length.
Element 1:
F ft h
BTU
k

= 8 . 24 , ft L 208 . 0 = ,
2 4
2
) 10 ( 326 . 5
4
ft
d
A

= = t , F T

80 =

,
F ft h
BTU
h

=
2
6 , ft d l 0818 . 0 = = t , 0
1
= h , 0
2
= h , 0 = Q , 0
1
= c and 0
2
= c . The element
equations are as follows.

=
)
`

(
(
(
(

+

+

+
) 80 (
4 326 . 5
) 0818 . 0 ( 6
) 80 (
4 326 . 5
) 0818 . 0 ( 6
2
208 . 0
) 4 326 . 5 ( 3
) 208 . 0 )( 0818 . 0 ( 6
208 . 0
8 . 24
) 4 326 . 5 ( 6
) 208 . 0 )( 0818 . 0 ( 6
208 . 0
8 . 24
) 4 326 . 5 ( 6
) 208 . 0 )( 0818 . 0 ( 6
208 . 0
8 . 24
) 4 326 . 5 ( 3
) 208 . 0 )( 0818 . 0 ( 6
208 . 0
8 . 24
2
1
e
e
T
T
e e
e e

or,

(
(
(

7667
7667
12 . 183 285 . 87
285 . 87 12 . 183
2
1
T
T

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 104
Element 2: Same as element 1 except we use
2
T and
3
T .
Assembling, we have

(
(
(
(
(
(

7667
15334
7667
12 . 183 285 . 87 0
285 . 87 25 . 366 285 . 87
0 285 . 87 12 . 183
3
2
1
T
T
T

The NBC term is zero. To apply the EBC we use the elimination technique resulting in

(
(
(
(
(
(

7667
28427
150
12 . 183 285 . 87 0
285 . 87 245 . 366 0
0 0 1
3
2
1
T
T
T

Solving, we have
{ } { } F T T T

97 . 88 , 82 . 98 , 150 , ,
3 2 1
=
The temperature decreases from left to right as expected. The element flux is computed as
follows.
Element 1:
2 1 2
1
1 1
6102 ) 150 82 . 98 (
208 . 0
8 . 24
) (
ft h
BTU
T T
L
k

= = = t
Element 2:
2 2 3
2
2 2
1174 ) 82 . 98 97 . 88 (
208 . 0
8 . 24
) (
ft h
BTU
T T
L
k

= = = t
Let us carry out a few checks. The solution satisfies the EBC of F

150 at 0 = x . The flux at

the right end must be zero. However, with this crude two-element model the error is large. To
obtain a better accuracy we must refine the mesh, i.e. add more elements to the mesh. We will
look at this aspect in the next two sections.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 105
Electrostatic Field
Fig. T4L3-4(a) and (b) show an essentially one-dimensional electrostatic field.
dx
x
L
D
A
A+ dA
D+ dD

Fig. T4L3-4 (a) A dielectric rod with its lateral surface electrically insulated
D

D+ dD
w
A

Fig. T4L3-4 (b) A wall or panel of dielectric material
The conservation of electric charge as per Gauss Law states that
( )( ) DA Adx D dD A dA + = + + (T4L3-6)
where
( ) D x displacement of electric charge,
2
Q L
( ) x charge density,
3
Q L
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 106
Neglecting second-order terms we have
( ) ( ) ( ) ( ) ( )
d
D x A x x A x
dx
= (T4L3-7)
The constitutive equation is given as
( ) ( ) ( ) D x x E x c = (T4L3-8)
where
( ) E x electrostatic field, V L
( ) x c permittivity, C L
c is a material property that measures how easily charged particles are displaced in a dielectric.
It is related to another property, the dielectric constant K , by the relation
0
K c c = where
0
c
is the permittivity of a vacuum. In addition we also need the relation

( )
( )
d x
E x
dx
u
= (T4L3-9)
where ( ) x u is the electrostatic potential (units: V). Eqn. (T4L3-9) is the corollary of
Coulombs Law which is a physical law describing the nature of electrostatic forces between
charged particles. Combining Eqn. (T4L3-8) and (T4L3-9) yields the constitutive equation

( )
( ) ( )
d x
D x x
dx
c
u
= (T4L3-10)
Substituting Eqn. (T3L3-10) into (T4L3-7) yields the differential equation

( )
( ) ( ) ( ) ( )
d d x
x A x x A x
dx dx
c
u | |
=
|
\ .
(T4L3-11)
The potential ( ) x u is the unknown function. The displacement ( ) D x plays the role of a flux.
The boundary conditions are either essential ( ( ) x u is known) or natural ( ( ) D x is known).
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 107
Lesson 4: Higher-Order Elements
Objectives: In this lesson we will look at the characteristics of elements that use higher-
order interpolation.

- To understand the properties of higher-order elements.
- To derive and use in an example a higher-order element.

Suggested Reading
Reference Pages
T1 Section 3.9
T2 Section 2.9
T3 Chapter 8


F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 108
Higher-Order Interpolation

So far we have seen the element concept illustrated using linear interpolation on an
element described by two nodes. We also saw in the classical Galerkins solution
methodology that superior results were obtained using higher order interpolation. Can we
not tie this concept to the element concept? The answer is a resounding Yes.

Fig. T4L4-1 shows the linear element (a) on which the solution is defined by a linear
polynomial, (b) is described by two nodes with nodal values
1
y and
2
y , and (c) the
resulting shape functions that use the nodal values to interpolate the solution.
2 2 1 1 2 1
~
) ( ) ( ) ( y x y x x a a x y | | + = + = (T4L1-29)

x
x
y
y
y
x
x
1
2
1
2
1
L
2
|
| |
1
1
1 2


Fig. T4L4-1

Fig. T4L4-2 shows the quadratic element (a) on which the solution is defined by a
quadratic polynomial, (b) is described by three nodes with nodal values
1
y ,
2
y and
3
y , and
(c) the resulting shape functions that use the nodal values to interpolate the solution.
3 3 2 2 1 1
2
3 2 1
~
) ( ) ( ) ( ) ( y x y x y x x a x a a x y | | | + + = + + = (T4L4-1)

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 109
x
x
y
y
y
y
x
x x
1
3
2
1
3 2
1
L/2 L/2
2 3
|
|
|
|
1
1
1
1
3
2

Fig. T4L4-2
A few comments are in order.
(a) For a one-dimensional element, we need to have nodes at the ends of the elements so that
one element can be tied to the element next to it. So we need a minimum of two nodes.
(b) When we have one degree of freedom per node, the number of nodes in the element is
equal to the number of coefficients in the trial solution. In other words, we need a total of
two nodes for the linear element and a total of three nodes for the quadratic element. Only
then will we have sufficient equations to write the coefficients of the polynomials in terms
of the nodal values.
(c) The number of shape functions is also equal to the number of nodal values. We should now
look at the properties of the shape functions. First, the trial solution must be complete a
quadratic polynomial must have the constant, linear and quadratic terms (e.g. having the
constant, quadratic and cubic terms is not allowed). This is done to ensure that the element
can assume all possible solution modes. Second, the shape functions satisfy the following
conditions
ij j i
x o | = ) ( (T4L4-2)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 110
where
ij
o is the Kroneckers Delta. In other words, the shape function will have a unit value at
the node it is associate with and zeros at the other nodes (see the plots shown above). This also
ensures that the shape functions are linearly independent.
Now on with the quadratic element. Using Eqn. (T4L4-1) and the nodal conditions, we have

(
(
(
(
(
(

3
2
1
2
1
0
2
3 3
2
2 2
2
1 1
1
1
1
y
y
y
a
a
a
x x
x x
x x
(T4L4-3)
These equations are similar to Eqn. (T4L1-31). Solving for the
i
a s and collecting like terms, we
have

3
2 3 1 3
2 1
2
3 2 1 2
3 1
1
3 1 2 1
3 2
~
) )( (
) )( (
) )( (
) )( (
) )( (
) )( (
) ( y
x x x x
x x x x
y
x x x x
x x x x
y
x x x x
x x x x
x y


+


+


= (T4L4-4)
Hence,

) )( (
) )( (
) (
3 1 2 1
3 2
1
x x x x
x x x x
x


= | (T4L4-5a)

) )( (
) )( (
) (
3 2 1 2
3 1
2
x x x x
x x x x
x


= | (T4L4-5b)

) )( (
) )( (
) (
2 3 1 3
2 1
3
x x x x
x x x x
x


= | (T4L4-5c)
There is an easier way to compute the shape functions that we shall see in the next module. We
still have two questions to answer before we can generate the element equations for the
quadratic element. First, Where is node 2 located within the element? Again, a detailed
answer will be generated in Module 2. For the time being let us assume that it is located at the
center of the element. Hence,

2
) ( ) (
2 3 1 2
L
x x x x = = L x x = ) (
1 3
(T4L4-6)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 111
Second, How do we handle the ) (x o , ) (x | and ) (x f terms? We will develop a
sophisticated technique in Module 2. For the time being let us assume that they are constants
within the element. Hence,

(
(
(
(
(
(

(
(
(
(
(
(
(

+
(
(
(
(
(
(
(

0 0 0
0 0 0
0 0 1
30
4
30
2
30
30
2
30
16
30
2
30 30
2
30
4
3
7
3
8
3
3
8
3
16
3
8
3 3
8
3
7
1
g
L L L
L L L
L L L
L L L
L L L
L L L
| | |
| | |
| | |
o o o
o o o
o o o


1
3 2
3
0 0 0
0 0 0
0 0 1
y
h y
y
( (
( (
+ =
`
(
(

(
(
)

3
1
0
6
6
4
6
c
c
fL
fL
fL
(T4L4-7)
The flux in the element is given by

~
~
2 3 1 3
1 2 2 2
2(2 ) ( 4)(2 )
( )
d y x x x x x x
x y y
dx L L
t o o

|
= = + +

\
1 2
3 2
2(2 ) x x x
y
L

|
|
.

( ) ( ) ( ) ( )
1 2 3 1 2 3 2 1 3 3 1 2 2
4 8 4 4 2 2 2 2 x y y y x y y x y y x y y
L
o
( = + + +

(T4L4-8)
In a similar manner we can generate higher order elements involving cubic, quartic, quintic etc.
trial functions with four, five, six etc. nodes per element. The manner in which we identify
these elements is usually as follows. These elements are designated as 1D-C
m
interpolation order
element where
m
C denotes that the problem variable and its derivatives up to order m are
continuous across element boundaries. In the element formulation that we have discussed so
far, only the problem variable is continuous across the element boundaries. Hence the elements
are designated 1D-C
0
linear element, or 1D-C
0
quadratic element etc.



F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 112
Illustrative Example T4L4-1
Let us resolve the last problem from the previous section. This time we will use the quadratic
element. Let us assume that the number of nodes is the same. The FE model is shown in Fig.
T4L4-3.
2 3 1
1
2.5
2.5
X

Fig. T4L4-3
Element 1:
F ft h
BTU
k

= 8 . 24 , ft L 416 . 0 = ,
2 4
2
) 10 ( 326 . 5
4
ft
d
A

= = t , F T

80 =

,
F ft h
BTU
h

=
2
6 , ft d l 0818 . 0 = = t , 0
1
= h , 0
2
= h , 0 = Q , 0
1
= c and 0
2
= c . The element
equations are as follows ( k = o ,
A
hl
= | ,
A
hlT
f

= )

\
|
(
(
(
(
(
(
(

) 416 . 0 ( 3
) 8 . 24 ( 7
) 416 . 0 ( 3
) 8 . 24 ( 8
) 416 . 0 ( 3
) 8 . 24 (
) 416 . 0 ( 3
) 8 . 24 ( 8
) 416 . 0 ( 3
) 8 . 24 ( 16
) 416 . 0 ( 3
) 8 . 24 ( 8
) 416 . 0 ( 3
) 8 . 24 (
) 416 . 0 ( 3
) 8 . 24 ( 8
) 416 . 0 ( 3
) 8 . 24 ( 7
+
=

|
|
|
|
|
|
|
|
.
|
(
(
(
(
(
(
(

3
2
1
30
) 416 . 0 )( 52 . 921 ( 4
30
) 416 . 0 )( 52 . 921 ( 2
30
) 416 . 0 )( 52 . 921 (
30
) 416 . 0 )( 52 . 921 ( 2
30
) 416 . 0 )( 52 . 921 ( 16
30
) 416 . 0 )( 52 . 921 ( 2
30
) 416 . 0 )( 52 . 921 (
30
) 416 . 0 )( 52 . 921 ( 2
30
) 416 . 0 )( 52 . 921 ( 4
T
T
T

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 113

6
) 416 . 0 )( 4 . 73721 (
6
) 416 . 0 )( 4 . 73721 ( 4
6
) 416 . 0 )( 4 . 73721 (

Or,

(
(
(
(
(
(

3 . 5111
20445
3 . 5111
7 . 181 38 . 116 4256 . 1
38 . 116 33 . 488 38 . 116
4256 . 1 38 . 116 7 . 181
3
2
1
T
T
T

Applying the boundary conditions (EBC for
1
T ), we have

(
(
(
(
(
(

2 . 5325
37902
150
7 . 181 38 . 116 0
38 . 116 33 . 488 0
0 0 1
3
2
1
T
T
T

Solving,
{ } { } F T T T

26 . 93 , 84 . 99 , 150 , ,
3 2 1
=
The flux in the element is computed using Eqn (T4L4-8) and are given by

2
( 0.0879') 6382
BTU
x
h ft
t = =

2
( 0.3281') 383.1
BTU
x
h ft
t = =

The reason for selecting the two locations will be explained in the next module. Bear in mind
that the flux distribution in this element is, as Eqn. (T4L4-8) shows, linear unlike the linear
element in which the flux is a constant. Lets compare the two results.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 114

Temperature ( F

)
Flux
|
|
.
|

\
|

2
ft h
BTU

Node Linear
Elements
Quadratic
Element
Linear
Elements
Quadratic
Element
1 150 150
2 98.92 99.84 6102 6382
3 88.97 93.26 1174 383.1

The nodal temperatures are close but the flux values are quite different.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 115
Lesson 5: Mesh Refinement and Convergence
Objectives: In this lesson we will look at the approach to applying the FE method to solve
a problem.

- To understand the concepts associated with mesh refinement and convergence.
- To resolve one or more previous examples with different meshes and use the concept
of convergence to obtain the solution.

Suggested Reading
Reference Pages
T1
T2 Section 2.12
T3 Chapter 9

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 116
Mesh Refinement

In the previous lessons we learnt two important facts. First, increasing the number of
degrees of freedom increases the accuracy of the solution. In other words, if we keep on
increasing the number of elements (and nodes) in the FE model, we should obtain better
solutions that converge to the exact result. This is known as h-convergence. The h notation
refers to the size of the elements. Second, increasing the order of the polynomial in the
trial solution also increases the accuracy of the solution. In other words, if we keep on
increasing the element order (and hold the number of elements constant), we should obtain
better solutions that converge to the exact result. This is known as p-convergence. The p
notation refers to the polynomial order. We could also combine the two and obtain what
is known as hp-convergence.

The advantages of low-order finite elements are (a) the size of the element matrices is small,
(b) the computational ease with which the elements can be generated, and (c) the size of the
hand-band width of the structural stiffness matrix is small. The major disadvantage is that
the convergence is slow. The comments pertaining to higher-order elements are just the
opposite.

The FE mesh need not be uniform nor contain just one type of element. It may be
advantageous to use the lower-order elements where the solution does not change rapidly
(flux has a low value) and use the higher-order elements in regions where higher accuracy is
required. By the same token, the mesh can be finer where higher accuracy is required.

The biggest disadvantage of any FE solution is the computational expense in solving the
system equations. Hand calculations (with help from equation solvers in calculators) are
tedious if the number of unknowns is greater than about 10. However, the finite element
method is a numerical technique that is ideal for computer-based solutions. As we
mentioned in the first topic, the amount of human time taken to create the data and
examine the results is far greater than the time taken to solve most FE problems.

We will illustrate these ideas using the results from the computer program. If you wish you
may jump to the section Using the 1DBVP

Program to familiarize yourself with the

program before reading the next section.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 117
An Illustrative Example T4L5-1

Let solve Example T4L3-2 using the concepts discussed earlier. We will monitor three
response quantities temperature at the right end, the flux at the left end that is the highest
flux in the model and the flux at the right end that should converge to zero. The results are
summarized below.

Model
ID
Number of
elements
Element
Type
Temperatur
e at right
end
Flux at
right end
Flux at left
end
1 2 Linear 89 1 174 6 102
2 4 - do - 90.5 542 7 804
3 8 - do - 90.9 266 8 970
4 16 - do - 91 132 9 664
5 32 - do - 91 66 10 044
6 64 - do - 91 33 10 244
7 128 - do - 91 17 10 346

1 2 Quadratic 91.1 417 8 039
2 4 - do - 91 220 9 137
3 8 - do - 91 111 9 767
4 16 - do - 91 55 10 102
5 32 - do - 91 28 10 275
6 64 - do - 91 14 10 362

The temperature at the left end converges very rapidly to F

91 . The flux at the right end

appears to converge linearly. The flux at the left end converges much more slowly for the
linear element compared to the quadratic element.

Figs. T4L5-1 and T4L5-2 show the plots for the flux at the left and right ends for the linear
and quadratic elements.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 118
0
250
500
750
1000
1250
0 15 30 45 60 75
Quadratic Element
Linear Element
Number of Elements
F
l
u
x
Flux at Right End

Fig. T4L5-1

6000
7000
8000
9000
10000
11000
0 15 30 45 60 75
Quadratic Element
Linear Element
Number of Elements
F
l
u
x
Flux at the Left End

Fig. T4L5-2
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 119
Summary

We achieved two primary goals with the lessons in this topic. First, we generalized the
Galerkins Method by tying it to the element concept. This made it easier to generate the
trial solutions that are valid over an arbitrary problem subdomain (the element!). The side
effect is that we can handle problems with known discontinuities in the solution, e.g. the
flux. Second, we looked at the manner in which the boundary conditions could be applied
more easily. While the problems essential boundary conditions were satisfied exactly, the
higher-order boundary conditions were satisfied only in the limit.

The derivation of the element equations is perhaps the most important step when
implementing the FE method. We saw how we could generate a family of trial solutions
and the associated elements. There are two important issues that we will deal with in the
next module that will make it easier to generate the element equations a more rational
way to generate the shape functions, and the isoparametric formulation.

It should be pointed out that seemingly different engineering problems are all governed by
the same differential equation. However, one still needs to contend with the physics
behind the parameters and the equations.

Better solutions are obtained by (a) increasing the number of elements in the mesh or, (b)
by increasing the order of the element (trial solution), or (c) both. With the easy
availability of FE computer programs, it is now possible to obtain accurate solutions in a
relatively short period of time. The 1DBVP

program illustrates (in its own manner) the

aspects associated with pre-processing, solution and post-processing.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 120
Review Exercises
In all the problems below, generate all the steps by hand except solve the systems
equations using a programmable calculator or a computer program. Finally check
your answers using the 1DBVP

program. Explain the differences in answers, if any.

Problem T4L2-1
Consider the " 4 bar shown below that is loaded by an axial surface traction given by the
function
2
( ) f x x lb in = .

X


The bar properties are as follows - psi E ) 10 ( 30
6
= and
2
2 in A = . Determine the stresses in
the bar using
(a) One linear element.
(b) Two linear elements.
(c) Four linear elements.

Comment on the results.

Problem T4L2-2
Resolve Problem T4L2-1 but use quadratic elements instead.
Problem T4L2-3
The structure shown in the figure below is subjected to an increase in temperature
C T

80 = A . In addition the loads are given as follows: kN P 60
1
= , kN P 75
2
= . Determine
the displacements, stresses, and support reactions using linear elements. The material
properties are as follows.




F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 121
Material
) (
2
mm A
) (GPa E
) ( C mm mm

o
Bronze 2400 83 18.9e-6
Aluminum 1200 70 23e-6
Steel 600 200 11.7e-6


P
800 mm
Bronze
Aluminum
Steel
600 mm
400 mm
P
1 2


Problem T4L3-1
Consider a brick wall of thickness cm L 30 = , C m W k

= 7 . 0 . The inner surface is at
C

28 and the outer surface is exposed to cold air at C T

15 =

. The heat transfer

coefficient associated with the outside surface is C m W h

=
2
40 . Determine the steady-
state temperature distribution within the wall and also the heat flux through the wall.
Assume a one-dimensional heat flow. Use a two linear-element model.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 122
30 cm
0
28 C
h, T
00



Problem T4L3-2
Figure shows a thin, cylindrical rod, 1m long, composed of two different materials the
two end sections each 40 cm long, are made of steel
|
.
|

\
|


=
C cm s
cm cal
k
2
12 . 0 , and the center
section is made of copper
|
.
|

\
|


=
C cm s
cm cal
k
2
92 . 0 is 20 cm long. The cross section is
circular, with a radius of 2 cm. Heat is flowing into the left at a steady rate of
2
1 . 0 cm s cal . The temperature of the right end is maintained at a constant C

0 . The rod
is in contact with air at an ambient temperature of C

20 so there is free convection from

the lateral surface. The convective coefficient is given as
C cm s
cal

2
4
10 5 . 1 .
Determine the temperature and flux distribution in the rod using (a) only linear elements
and (b) only quadratic elements.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 123

Problem T4L3-3
A variable area rectangular cross-section fin transmits heat away from a mass as indicated in the
figure. The thickness of the fin in the direction perpendicular to the paper is ten times that
shown in the plane of the paper. There is convection on the entire lateral surface.
h
L
L
M
a
s
s

a
t
t
e
m
p
e
r
a
t
u
r
e
,

T
3h T=
h , T
h , T
T
0
0
0
c
c
0
0
0
0
1
2
0
0

With cm h 5
0
= , cm L L 10
2 1
= = , C T

400
0
= , C T

100 =

, C mm W h
c

=
2 3
10 and
C mm W k

=
2
30 . 0 , find the temperature and flux distribution.
2
1 . 0
cm s
cal

= t
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 124
Problem T4L3-4
Consider a slab of thickness 0.1 m with a thermal conductivity
( )
40 k W m C =

in which
energy is generated at a constant rate of
6 3
10 W m . The boundary surface at 0 x = is insulated,
and the one at 0.1 x m = is subjected to convection with a heat transfer coefficient of
( )
2
200W m C

into an ambient at a temperature of 150 C

. Find the temperature and flux

distribution.
Problem T4L3-5
When radial heat flow takes place in cylindrical bodies, the differential equation is different
than what we saw in 1D-BVP problem. Consider steady-state radial heat conduction in a long
solid cylinder of radius r b = , in which energy is generated at a rate of
3
( ) Q r W m . The
temperature distribution ( ) T r in the cylinder is governed by the following heat conduction
equation

1 ( )
( )
d dT r
kr Q r
r dr dr
(
=
(

a r b < <
with the appropriate boundary conditions at r a = and r b = . Develop the appropriate
element equations for solving this problem using the Galerkins Method for the (a) 1D-C
0
linear
element and (b) 1D-C
0
quadratic element.
Problem T4L3-5
A 10-cm diameter solid chrome-nickel rod with thermal conductivity of
( )
20 k W m C =

is
heated electrically by the passage of an electric current which generates energy within the rod at
a uniform rate of
7 3
10 W m . The surface of the rod is subjected to convection with a heat
transfer coefficient
( )
200 h W m C =

into an ambient at 30 C

. Find the temperature

distribution in the rod.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 125
Using the OneDBVP

Program
The OneDBVP program is designed to solve the problem described by equations (T4L1-56)
and (T4L1-57). It can be used on computer systems running any version of the Windows
OS (98, NT, 2000 and XP). Install the program as you would any other Microsoft
Windows program. A few things to note before we get started with the program.

(a) The first step is to get the problem data ready. Note that the program solves the
problem as described by Eqns. (T4L1-56) and (T4L1-57). In other words, it does not
specifically solve a solid mechanics problem or a heat transfer problem etc. You must
relate the parameters from the problem to the general one-dimensional BVP as we have
done in Lessons 2-4.
(b) The program does not assume any units. The problem units must be consistent.
(c) Finally we will discuss the program terminology. The FE mesh is described in terms of
segments. A segment contains one or more elements that have exactly the same
properties. Only the nodes that make up the elements have different locations and
numbers. The problem domain is divided into one or more segments. The program
numbers the nodes and elements (starting at 1) starting with the leftmost segment. Note
that you create the segments; 1DBVP program creates the elements, nodes and loads.

Once you launch the program, you should see the following window.



Using the program is as easy as 1-2-3. Press F1 to invoke the on-line help. Go through the
on-line help on understand how to use the program. The ) (x o , ) (x | and ) (x f terms are
assumed to be described by a quadratic polynomial of the form c bx ax + +
2
and the input
to the program are the c b a , , values. The type of boundary condition is recognized by the
program using an integer code - 0 denotes an EBC, 1 denotes an NBC and 2 denotes a
mixed BC. The first two need an additional real input. The mixed BC needs two values -
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 126
the h g / and c values. Finally, an interior concentrated flux if one exists, is input
separately .

Problem Example T2L2-1 Example T4L3-1 Example T4L3-2
Units N, m, C, s N, m, C, s lb, ft, F, s
# of Segments 3 3 1
Segment 1
Element Order 1 1 1
# of elements 1 1 2
Left End Coor. 0 0 0
Right End Coor. 0.15 0.3 0.416
) (x o 0, 0, 5e7 0, 0, 20 0, 0, 24.8
) (x | 0, 0, 0 0, 0, 0 0, 0, 921.517
) (x f 0, 0, 0 0, 0, 0 0, 0, 73721.367
Segment 2
Element Order 1 1
# of elements 1 1
Left End Coor. 0.15 0.3
Right End Coor. 0.30 0.45
) (x o 0, 0, 5e7 0, 0, 20
) (x | 0, 0, 0 0, 0, 0
) (x f 0, 0, 0 0, 0, 0
Segment 3
Element Order 1 1
# of elements 1 1
Left End Coor. 0.30 0.45
Right End Coor. 0.60 0.6
) (x o 0, 0, 8e7 0, 0, 50
) (x | 0, 0, 0 0, 0, 0
) (x f 0, 0, 0 0, 0, 0
Left End BC EBC Mixed EBC
Left End BC
Value(s)
0.0 -25, 20000 150
Right End BC EBC EBC NBC
Right End BC
Value(s)
0.0 20 0
Concentrated
Flux Data

Location, Value 0.15, 300e3
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 127
Program Limitations
There are no limitations except those imposed by the operating system.

Troubleshooting
The program checks the input for a variety of errors. Most of these deal with the validity of the
input data. The error messages are shown on the screen. Apart from the input errors the most
common error message is Ill-posed problem. This usually is an indication that the boundary
conditions have not been imposed correctly.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 128
Solving Linear Algebraic
Equat ions
A class of problems in finite element analysis involves the solution of the algebraic
equations
K D F
n n n m n m
= (1)
where K
n n
is the system matrix, D
n m
is the matrix of primary unknowns and F
n m
is the
matrix of right-hand side (RHS) vectors. The number of unknowns is n and the number
of RHS vectors is m. Typically, n m >> . In the context of finite element analysis, K
n n
is
mostly symmetric and positive definite. We will assume that K
n n
is symmetric and
positive definite

Equation solvers can be categorized as being either direct or iterative. Direct solvers are
those that solve Eqns. (1) in a non-iterative fashion. The procedure involves one pass or
two passes through the n equations. On the other hand, iterative solvers transform the
problem into an equivalent problem and the solution is obtained iteratively.

There are at least four major issues in solving Eqn. (1).
(a) How much of storage space will be used?
(b) How can numerically accurate solutions be obtained?
(c) How much time will be taken to obtain the solution?
(d) How much of additional effort is needed if an additional solution (to a new right-hand
side vector) vector is to be generated?

There are other issues such as parallelizing and vectorizing the solution procedure on
specialized hardware and software systems but they are beyond the scope of the discussions
here.

Storage Schemes
The structural stiffness matrix is usually sparse. Typically (especially for larger problems), the
nonzero entries make up a few percent (1-10%) of the entire matrix. Researchers have devised
several storage schemes to minimize the amount of storage space needed by recognizing that K
is sparse and that the locations of the nonzero entries are known once the structural model is
defined.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 129
Banded and Skyline: The K matrix is symmetric and banded. A banded matrix can be though
of as a special case of a sparse matrix. Consider the planar truss shown below (Fig. 1).
4
5 3
2 1

Fig. 1 A planar truss shows the node numbers and the global degrees-of-freedom at the nodes
There are 10 degrees-of-freedom in the system. If we go through the process of constructing the
element equations and forming the system equations (symbolically not numerically), we can
generate the map for the K matrix. Fig. 2 shows only the nonzero entries in the upper
triangular portion of the matrix.
11 12 15 16
22 25 26
33 34 35 36
44 45 46
55 56 57 58 59 5,10
66 67 68 69 6,10
77 78
88
99 9,10
10,10
K K K K
K K K
K K K K
K K K
K K K K K K
K K K K K
K K
K
K K
K
(
(
(
(
(
(
(
(
(
(
(
(
(
(
(


Fig. 2 Upper triangular portion of the system stiffness matrix
As can be deduced from Fig. 2, the HBW of the given matrix is 6. We can store this matrix as a
rectangular matrix as follows.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 130
11 12 15 16
22 25 26
33 34 35 36
44 45 46
55 56
10 6
66 67 68 69 6,10
77 78
88
99 9,10
10,10
0 0
0 0 0
0 0
0 0 0
0 0 0 0
0
0 0 0 0
0 0 0 0 0
0 0 0 0
0 0 0 0 0
banded
K K K K
K K K
K K K K
K K K
K K
K K K K K
K K
K
K K
K

(
(
(
(
(
(
(
=
(
(
(
(
(
(
(
(

K
Fig. 3 Stiffness matrix stored as a banded matrix
The relationship (mapping) between the elements in the original K in Fig. 2 and the banded
form K
banded
in Fig. 3 can be derived simply as follows.
K
i j ,
= 0 if j i HBW + > 1 b g (2)
K
i j ,
= 0 if j i < b g (3)
else K K
i j i j i
banded
, ,

+1
(4)
As we can see in Fig. 4, a good number of the elements are zero.

Fig. 4 The skyline profile
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 131
We can improve the storage scheme by storing only the elements that lie within the skyline of
the matrix as shown in Fig. 4. The stiffness matrix is stored in a vector as follows.
K
35 1 11 22 12 33 44 34 10 10 9 10 8 10 7 10 6 10 5 10
=
skyline
K K K K K K K K K K K K , , , , , ,..., , , , , ,
, , , , , ,
m r
(5)
Note that each column is stored starting with the diagonal element of that column followed by
all the other elements in that column until the last nonzero entry in that column. To facilitate
mapping the original elements of the stiffness matrix, an additional indexing vector, D
n
loc
+1
is
created that has ( ) n +1 elements. These elements store the location of the diagonal element (of
each column) with the last element storing the (last element+1) in the stiffness matrix. In other
words, the last element contains one more than the total number of entries in the skyline
profile. Going back to the current example, we have the following D
n
loc
+1
vector.

D
11
1 2 4 5 71218 21 25 30 36
loc
= , , , , , , , , , , l q (6)
The relationship (mapping) between the elements in the original K in Fig. 1 and the skyline
form K
skyline
in Fig. 4 can be derived as follows.
K
i j ,
= 0 if j i D D
j
loc
j
loc
>
+
b g d i
1
(7)
K
i j ,
= 0 if j i < b g (8)
K l D j i K
i j j
loc
l
skyline
,
= + (9)
For example, to locate K
68
we note that (a) i j = = 6 8 , , (b) D
loc
8
21 = , (c) l = + = 21 8 6 23.
Hence, K K
skyline
68 23
.
Finally, let us compare the storage requirements of the three schemes full, banded, and
skyline, assuming that the stiffness matrix is stored in the double precision format, and that two
integer words make up a single double precision word ( , ) q hbw m
n
loc
= =
+
D
1
1 .
Storage
Scheme
What is to be
stored?
Equivalent integer words
Full K
n n

2
2
n
Banded K
n hbw
banded

2nq
Skyline K
m
skyline
, D
n
loc
+1

2 1 m n + +
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 132

Using our example, we have the three values in the last column as 200 120 , and 81 integer
words significant savings with increasing sophistication of the storage scheme.
Sparse: It is not evident with our simple example that the stiffness matrix is sparse. The %
sparsity of the stiffness matrix is
35
100
100 35% = . As the size of the problem increases, in most
instances, the sparsity of the stiffness matrix increases (or the number of nonzero terms
decreases). Consider the stiffness matrix shown in Fig. 4. In one of the sparse storage schemes,
the stiffness matrix is stored in a vector rowwise with only the non-zero entries being stored
11
in
K
m
sparse
1
. Using our previous example, we have the following.
K
31 1 11 12 15 16 22 25 26 88 99 9 10 10 10
=
sparse
K K K K K K K K K K K , , , , , , ,..., , , ,
, ,
m r
(10)
To access the entries in the matrix, two additional (indexing) vectors are needed. The first, C
m1

is used to store the column numbers of the nonzero entries. Again, we have with our example
C
31 1
1 2 5 6 2 5 6 3 4 5 6 8 91010

= , , , , , , , , , , ,..., , , , l q (11)
The second, R
( ) n+ 1 1
, stores the starting location of each row. Again, we have with our example
R
11 1
15 81215 21 26 28 29 3132

= , , , , , , , , , , l q (12)
Discussion of sparse equation solvers is outside the scope of this document.

Offline: There are special situations when the solution procedure operates on matrices that
are written on and retrieved from computer hard disk. This is typically done when the size
of the stiffness matrix and other associated vectors/matrices in any storage format is several
times the size of the computers random access memory (RAM). Discussion of offline
equation solvers is outside the scope of this document.

DIRECT SOLVERS
2.1 Gaussian Elimination
There are two phases to the solution forward elimination and backward substitution. In
the forward elimination phase, the basic idea is to take the set of equations from the
original form

11
It should be noted that zero entries within the skyline profile may become nonzero during the solution phase.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 133

11 12 13 1 1 1 1
21 22 23 2 2 2 2
31 32 33 3 3 3 3
1 2 3
1 2 3
i n
i n
i n
i i i ii in i i
n n n ni nn n n
K K K K K D F
K K K K K D F
K K K K K D F
K K K K K D F
K K K K K D F
(
(
(
(
(
=
` `
(

(

(

(

(
) )
(13)

to an upper triangular matrix of the form

11 12 13 1 1 1 1
(1) (1) (1) (1) (1)
22 23 2 2 2 2
(2) (2) (2) (2)
33 3 3 3 3
( 1) ( 1) ( 1)
( 1) ( )
0
0 0
0 0 0
0 0 0 0
i n
i n
i n
i i i
ii in i i
n n
nn n n
K K K K K D F
K K K K D F
K K K D F
K K D F
K D F

(
(
(
(
(
=
` `
(

(

(

(

(
) )
(14)

With each step through the equations one unknown is eliminated until only D
n
remains as
the unknown in the equation. In step k k n = 1 2 1 , ,..., b g

K K
K
K
K i j k n
ij
k
ij
k ik
k
kk
k kj
k ( ) ( )
( )
( )
( )
, ,..., = = +

1
1
1
1
1 (15)
F F
K
K
F i k n
i
k
i
k ik
k
kk
k k
k ( ) ( )
( )
( )
( )
,..., = = +

1
1
1
1
1 (16)
If K
kk
k ( )
s
1
c , where c is a small positive constant, then the system of equations is linearly
dependent.

In the backward substitution phase (dropping the superscript)

D
F
K
n
n
nn
= (17)
and D
F K D
K
i n n
i
i ij j
j i
n
ii
=

=
= +

1
1 2 1 , ,..., (18)
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 134

Algorithm
1: Forward Elimination. Loop through rows, k n = 1 1 ,..., .
2: Check if K
kk
< c . If yes, stop.
3: Loop through columns, i k n = +1,..., .
4: Compute constant, c
K
K
ik
kk
= .
5: Loop through j k n = +1,..., .
6: Set K K cK
ij ij kj
= .
7: End loop j .
8: Set F F cF
i i k
= .
9: End loop i .
10: End loop k .
11: Backward substitution. Set D F K
n n nn
= .
12: Loop through all rows, i n = 1 1 ,..., .
13: Compute sum K D
ij j
j i
n
=
= +

1
.
14: Compute D
F sum
K
i
i
ii
=

.
15: End loop i .







F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 135
Example 1
Solve the following set of equations using Gaussian Elimination method.

1
2
3
10 5 2 6
3 20 5 58
2 7 15 57
D
D
D
(

(
=
` `
(

(
) )

Solution
Forward Substitution
Noting that n = 3, the successive snapshots as we go through the algorithm are as follows.
k = 1
1
2
3
10 5 2 6
21.5 4.4 56.2
6 15.4 58.2
D
D
D
(

(
=
` `
(

(
) )


k = 2
1
2
3
10 5 2 6
21.5 4.4 56.2
14.172 42.5163
D
D
D
(

(
=
` `
(

(
) )


Backward Substitution
D
3
425163
14172
3 = =
.
.

i = 2 D
2
56 2 4 4 3
215
2 =

=
. .
.
b g

i = 1 D
1
6 4
10
1 =

=
( )

Hence the solution is
3 1
1
2
3

=
`

)
D .

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 136
2.2 Cholesky Decomposition or LDL
T

Factorization (or, Decomposition)

The matrix A in = Ax b can be factored as A LU = where L is a lower triangular matrix with
1s on the diagonals and U is an upper triangular matrix. The proof can be found in a book on
linear algebra or numerical analysis. When A is symmetric and positive definite as is K in Eqn.
(1), the matrix can be factored as K LDL
T
=

where L is a lower triangular matrix with 1s on
the diagonals and

D is a diagonal matrix with positive entries. The LDL

T

decomposition is a
variation of Cholesky Decomposition, and provides a very effective solution to the system
equilibrium equations.

The solution proceeds as follows. Starting with Eqns. (1) we first factor K LDL
T
=

by finding
the L and

D matrices.
1
11 12 1 21 1
12 22 2 21 2 2
1 2 1 2

0 . 0
. 1 0 . 0 1 .

. 1 . 0 0 . 0 0 1 .
. . . . . . . . . . . .
. . . .
. . 1 0 0 . 1

0 0 .
n n
n n
n n nn n n
n
D
K K K L L
K K K L D L
K K K L L
D
(
( ( (
(
( ( (
(
( ( (
=
(
( ( (
(
( ( (
(




1 1 21 1 31 1 1
2
2 21 1 21 31 2 32 1 21 1 2 2
2 2
1 31 2 32 3 1 31 1 2 32 2 3 3
2 2
1 1 2 2

.

.

.
. .

...
n
n n
n n n
n n n
D D L D L D L
D L D L L D L D L L D L
D L D L D D L L D L L D L
sym D L D L D
+
+
(
(
+
(
(
= + + +
(
(
(
+ + +
(

(19)
In other words, we have the following.
LDL D F
T

= (20)
Let LQ F = (21)
Then

DL D Q
T
= (22)
Note that

DL
T
is of the upper triangular form.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 137

1 1 12 1 1
2 2 2

.

.

. .

n
n
n
D D L D L
D D L
D
(
(
(
=
(
(
(

T
DL
0
(23)
We can solve Eqns. (21) for Q through a process known as forward substitution starting with
the first equation that has one unknown, Q
1
, the second equation that then has one unknown,
Q
2
and so on until the last equation that is used to find Q
n
. Once Q has been computed, we can
solve Eqns. (22) through a process known as backward substitution by starting with the last
equation that has one unknown, D
n
, then the second last equation that then has one unknown,
D
n1
, and so on until the first equation that is used to compute D
1
.

By comparing the RHS of Eqn. (19) with the LHS we have the mechanism to compute the L
and

D matrices given the K matrix.

Algorithm
Step 1: Cholesky Factorization. Loop through rows, i n =1,..., .
Step 2: Set

D K L D
i ii ij j
j
i
=
=

2
1
1
. If

D
i
< c , stop. The matrix is not positive definite.
Step 3: For j i n = +1,..., , set L
K L D L
D
ji
ji jk k ik
k
i
i
=

=

1
1
.
Step 4: End loop i. This ends the factorization phase.
Step 5: Forward Substitution. Set Q F
1 1
= .
Step 6: For i n = 2,..., , set Q F L Q
i i ij j
j
i
=
=

1
1
. This ends the Forward Substitution phase.
Step 7: Backward Substitution. Set D
Q
D
n
n
n
=

.
Step 8: For i n = 1 1 ,..., , set D
Q
D
L D
i
i
i
ij j
j i
n
=
= +

1
. This ends the Backward Substitution phase.
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 138
A careful examination of the steps will show that no extra storage is required. The storage
locations in K can be used to store both

D and L. Similarly, the storage locations in D can be

used to store the elements of Q.
Example 2
Solve the following set of equations using Cholesky Decomposition method.

1
2
3
4
5
3.5120 0.7679 0 0 0 0
0.7679 3.1520 0 2 0 0
0 0 3.5120 0.7679 0.7679 0
0 2 0.7679 3.1520 1.1520 0.04
0 0 0.7679 1.1520 3.1520 0
D
D
D
D
D
(
(

(

= (
` `
(


(

(
) )

Solution
Factorization n = 5 b g
i = 1

. D K
1 11
35120 = = .
2 j = L
K
D
21
21
1
0 7679
35120
= =

.
.
= 0.21865. Also, L L L
31 41 51
0 = = = .
i = 2

. . . . D K L D
2 22 21
2
1
2
31520 0 21865 35120 29841 = = = b g b g .
j = 3 L
K L D L
D
32
32 31 1 21
2
0 =

=

.
j = 4 L
K L D L
D
42
42 41 1 21
2
2 0
29841
=

=

=

.
0.670219. Also, L
52
0 = .
i = 3

. . D K L D L D
3 33 31
2
1 32
2
2
35120 0 0 35120 = = = .
j = 4 L
K L D L L D L
D
43
43 41 1 31 42 2 32
3
0 7679 0 0
35120
=

=

=

.
.
0.21865.
j = 5 L
K L D L L D L
D
53
53 51 1 31 52 2 32
3
0 7679 0 0
35120
=

=

=

.
.
0.21865.
i = 4

D K L D L D L D
4 44 41
2
1 42
2
2 43
2
3
=
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 139
= = 31520 0 0 670219 29841 0 21865 35120
2
2
. ( . ) . . . b g b g b g 1.64366
j = 5 L
K L D L L D L L D L
D
54
54 51 1 41 52 2 42 53 3 43
4
=

=

=
11520 0 0 0 21865 35120 0 21865
164366
. . . .
.
b gb gb g
-0.598724
i = 5

D K L D L D L D L D
5 55 51
2
1 52
2
2 53
2
3 54
2
4
=
= = 31520 0 0 0 21865 35120 0598724 164366 23949
2
2
. ( . ) . . . . b g b g b g

Forward Substitution LQ F = b g
1
2
3
4
5
1 0 0 0 0 0
0.21865 1 0 0 0 0
0 0 1 0 0 0
0 0.670219 0.21865 1 0 0.04
0 0 0.21865 0.598724 1 0
Q
Q
Q
Q
Q
(
(
(

= (
` `
(


(

(
) )


i = 1 Q F
1 1
0 = =
i = 2 Q F L Q
2 2 21 1
0 = =
i = 3 Q F L Q L Q
3 3 31 1 32 2
0 = =
i = 4 Q F L Q L Q L Q
4 4 41 1 42 2 43 3
0 04 = = .
i = 5 Q F L Q L Q L Q L Q
5 5 51 1 52 2 53 3 54 4
0 0598724 0 04 = = = . . b gb g 0.023949




F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 140
Backward Substitution

DL D Q
T
=
d i

1
2
3
4
5
3.5120 0.21865 0 0 0 0
0 2.9841 0 0.670219 0 0
0 0 3.5120 0.21865 0.21865 0
0 0 0 1.64366 0.598724 0.04
0 0 0 0 2.3949 0.023949
D
D
D
D
D
(
(

(

= (
` `
(


(

(
) )

i = 5 D
5
0 023949
23949
0 01 =

=
.
.
.
i = 4 D
Q
D
L D
4
4
4
45 5
0 04
164366
0598724 0 01 = =

=

.
.
. . b gb g 0.0303232
i = 3 D
Q
D
L D L D
3
3
3
34 4 35 5
=

= =
0
35120
0 21865 0 0303232 0 21865 0 01
.
. . . . b gb g b gb g 0.00444367
i = 2 D
Q
D
L D L D L D
2
2
2
23 3 24 4 25 5
=

= =
0
29841
0 0 670219 0 0303232 0
.
. . b gb g 0.0203232
1 i = D
Q
D
L D L D L D L D
1
1
1
12 2 13 3 14 4 15 5
=

= =
0
35120
0 21865 0 0203232 0 0 0
.
. . b gb g 0.00444367

It should also be noted that the Cholesky Decomposition requires no special effort to solve
additional RHS vectors since the factorization and the forward/backward substitutions steps
are separate. Once the factorization step is completed (once), the forward/backward
substitutions steps can be repeated as many times as required.



F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 141
ITERATIVE SOLVERS
Consider a quadratic form ( ) F d given by
( ) 2
T T
F = d d Kd d f (24)
where K is a positive definite matrix. Minimizing the quadratic form yields
= Kd f (25)

The idea of minimizing the quadratic form obviously leads to unconstrained optimization
techniques. Consider the Steepest Descent Method. Let s be a vector such that
( ) ( )
k k
F F = + d d s (26)
Using Eqn. (24), we have
( ) ( ) ( ) ( ) 2
T
T
k k k k
F + = + + + d s d s K d s d s f (27)
If s is considered to be small, the second order terms can be neglected and the above
equation can be simplified as
( ) ( ) 2
T T T
k k k k
F F + = + d s d s Kd d Ks s f (28)
However, due to Eqn. (26) we can rewrite the above equation as
( ) ( ) 0
T
T
k k
+ = s Kd f Kd f s (29)
The residual vector

k k
= g Kd f (30)
is thus orthogonal to any small s whose addition does not alter F . The steepest descent
direction is the direction of
k
g . To minimize F , we can think of moving along the
k
g
direction. In other words, we are looking for the point
k k k
o + d g by trying to find the
value of
k
o that minimizes ( )
k k k
F o + d g .

1
2
( ) ( )
( ) 2
k k k k
T T T T
k k k k k k k k k k k k
F F
F
o
o o o o
+
= +
= + + + +
d d g
d g Kd d Kg g Kg g f
(31)
Minimizing
( ) 0
k k k
k
d
F d g
d
o
o
+ = (32)
yields
0
T T T
k k k k k k
o + = g Kd g Kg g f (33)
or

T
k k
k T
k k
o =
g g
g Kg
(34)

3.1 Conjugate Gradient Method
The Steepest Descent Method while providing a good start, is not effective. It takes far too
many iterations. Error corrected in one iteration is likely to be introduced in the next iteration
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 142
since successive directions are not mutually orthogonal. The Conjugate Gradient Method
avoids this difficulty. Successive solutions are found by

1 k k k k
o
+
= + d d s (35)
where the direction vectors
k
s are computed so that successive residuals are orthogonal to each
other as

1
0
T
k k+
= g g (36)
The correction step is computed as follows

1
2
( ) ( )
( ) 2
k k k k
T T T T
k k k k k k k k k k k k
F F
F
o
o o o o
+
= +
= + + + +
d d s
d s Kd d Ks s Ks s f
(37)
As before, minimizing F leads to
0
T T T
k k k k k k
o + = s Kd s Kg s f (38)
or

T
k k
k T
k k
o =
s s
s Ks
(39)

This expression is valid for any direction
k
s . Multiplying Eqn. (35) by K and subtracting
f on both sides yields


1 k k k k
o
+
= + g g Ks (40)
Using the orthogonality relationship in Eqn. (36) and premultiplying both sides of Eqn.
(40) by
k
g yields an alternate form of Eqn. (39) (one that ensures that successive directions
are orthogonal) as

T
k k
k T
k k
o =
g g
g Kg
(41)

Rewriting Eqn. (38) we have

( ) 0
T T
k k k k k
o + = s Kd f s Kg (42)
Using Eqn. (40) we have

1
0
T
k k+
= s g (43)
To create the search direction
1 k+
s , we start by taking the residual
1 k+
g as we discussed in
the steepest descent method. Adding a component orthogonal to
1 k+
g just large enough to
satisfy Eqn. (35) we have

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 143

1 k k k k
|
+
= + s g s (44)
where the scaling constant
k
| is yet to be determined. Premultiplying both sides by
T
k
s K,
Eqn. (35) becomes

1 1
T T T
k k k k k k k
|
+ +
= + s Kd s Kg s Ks (45)
and the multiplier
k
| is given by

1 1
T T
k k k k
k T
k k
|
+ +

=
s Kg s Ks
s Ks
(46)
Orthogonality of residuals (Eqn. (35)) is the same as requiring Eqns. (39) and (41) to yield
the same value for
k
o . Premultiplying Eqn. (44) by
1
T
k+
g


1 1 1 1
T T T
k k k k k k k
|
+ + + +
= + g s g g g s (47)

Using Eqn. (43), we have
1 1
T T
k k k k + +
= g g g s (48)
or

T T
k k k k
= g g g s (49)
Using this result in Eqn. (41), we have


T
k k
k T
k k
o =
s g
g Ks
(50)
Substitute Eqn. (44) in (55)

1 1
T
k k
k T T
k k k k k
o
|

=

s g
s Ks s Ks
(51)
If
k
s and
1 k+
s are made to be conjugate to each other with respect to K as in

1
0
T
k k +
= s Ks (52)
then the second term in the denominator of Eqn. (51) vanishes and Eqns. (39) and (41) are
identical. This requirement is met by setting

1
T
k k
k T
k k
|
+
=
s Kg
s Ks
(53)
Eqn. (52) is the one that gives the method its name.

AL GORI THM ( FOR A S I NGL E RHS VECTOR)
Start with an initial guess
0
d . Set 0. k = Set c as convergence tolerance and
max
n as the
maximum number of iterations.
1. Construct
0 0
g = Kd f and
0 0
= s g .
F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 144
2. Now form
T
k k
k T
k k
o =
g g
s Ks
.
3. Update the solution vector as
1 k k k k
o
+
= + d d s .
4. Now update
1 k k k k
o
+
= + g g Ks .
5. Check if
T
k k
c s g g and if so terminate iterations. Otherwise form
1 1
T
k k
k T
k k
|
+ +
=
g g
g g
.
6. Update
1 1 k k k k
|
+ +
= + s g s .
7. Increment k . Check if
max
k n < . If so, go to Step 2. Otherwise terminate iterations.


3.2 Preconditioned Conjugate Gradient Method
The basic idea in preconditioning is to devise a scheme to accelerate the convergence process.
We would like to generate a preconditioning matrix, B such that

1
~ B K I
The area of finding the best pre-conditioner is a very active research area. One of the simplest
scheme is called Jacobi preconditioning and is discussed next.
Jacobi Preconditioning: Assume that
1
1
ii
diag
K

| |
=
|
\ .
B .
Start with an initial guess
0
0 = d . Set 0. k = Set c as convergence tolerance and
max
n as the
maximum number of iterations.
1. Construct
0
r = f ,
1
0 0

= z B r and
0 0
= s z .
2. Now form
T
k k
k T
k k
o =
r z
s Ks
.
3. Update the solution vector as
1 k k k k
o
+
= + d d s .
4. Now update
1 k k k k
o
+
= r r Ks .
5. Check if
T
k k
c s r r and if so terminate iterations. Otherwise form
1 1
T
k k
k T
k k
|
+ +
=
r z
r z
.
6. Form
1
1 1 i i

+ +
= z B r
7. Update
1 1 k k k k
|
+ +
= + s z s .
8. Increment k . Check if
max
k n s . If so, go to Step 2. Otherwise terminate iterations.

F I N I T E E L E M E N T S F O R E N G I N E E R S : M O D U L E 1
1998-2010, S. D. Raj an 145

Index
1DBVP
Input file, 115, 116
Output file, 2, 118
Program Limitations, 118
Troubleshooting, 118
1DBVP

, i, 2, 3, 107, 110, 111, 115

Input file, 115, 116
Output file, 2
Program Limitations, 118
Troubleshooting, 118
1DBVP

Output file, 118

Algebraic equations, 12, 23, 54, 63, 68
Assembly, 30, 31, 41
Boundary conditions, 21, 27, 30, 31, 43, 44, 46,
47, 52, 53, 59, 63, 70, 72, 74, 76, 79, 81, 83,
86, 87, 89, 90, 92, 93, 94, 95, 97, 104, 110,
118
EBC, 43, 44, 46, 47, 69, 72, 74, 76, 79, 83,
87, 89, 92, 93, 94, 95, 96, 97, 104, 116
Essential, 43, 110
Mixed, 21, 79, 89, 91, 92, 93, 116
Natural, 69, 72, 74, 79, 83, 89, 91, 93, 95, 96,
116
Boundary-value problem, 52, 65
Centroid, 80
Classical Technique, 13, 55, 56, 59, 60, 63, 66,
69, 99
Collocation Method, 54
Compatibility, 34, 41, 42
Conduction, 9
Convection, 9, 16, 90, 113
Convective Coefficient, 94, 113
Degrees of freedom, 40, 53, 57, 60, 76, 77, 107
Derived variable, 33, 48
Differential equations, 1, 21, 51, 53
Direct Method
Cholesky Decomposition, 44, 126, 128, 131
Direct stiffness method, 31, 52
Discrete elements, 32
Discretization, 29, 70
Displacements, 14, 26, 33, 34, 41, 45, 49, 84,
86, 87, 111
Electrical Network, 32, 39
Electromagnetics, 1, 24
Element equations, 31, 34, 36, 38, 39, 40, 42,
47, 72, 73, 77, 79, 80, 81, 85, 92, 93, 95, 101,
103, 110
Elements, 1, 4, 5, 6, 22, 24, 25, 26, 27, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
45, 46, 47, 48, 49, 50, 55, 70, 71, 72, 73, 74,
75, 76, 77, 78, 79, 80, 81, 84, 85, 86, 87, 92,
93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
104, 107, 108, 110, 111, 112, 113, 116, 118,
119, 120, 121, 122, 128
Discrete, 32
Linear, 111
Quadratic, 108, 111, 113
Elimination Approach, 44, 93
Equations
Algebraic, 12, 23, 54, 63, 68, 119
Differential, 1, 21, 51, 53, 54, 63, 64, 78, 79,
83, 89, 110
Element, 31, 34, 36, 38, 39, 40, 42, 47, 72,
73, 77, 79, 80, 81, 85, 92, 93, 95, 101, 103,
110, 120
Integral, 23
Equilibrium, 33, 34, 41, 42, 45, 46, 48, 52, 87,
127
Errors, 30, 48, 62, 81, 118
Finite element mesh, 25, 27, 28, 29, 40, 47, 70,
73, 79, 97, 106, 107, 110, 116
Fluid mechanics, 1, 10
Flux, 35, 36, 58, 61, 62, 63, 64, 66, 68, 72, 74,
76, 77, 78, 79, 81, 86, 89, 92, 94, 96, 97, 102,
104, 105, 107, 108, 110, 112, 113, 114, 116
Function
Quartic, 62
Residual, 63
h-convergence, 107
Heat transfer, 1, 9, 16, 23, 24, 32, 88, 89, 91,
112, 115
Integral equations, 23
Interpolation, 12, 70, 80, 98, 99, 102
Linear algebra, 12, 54, 126
Loading Term, 66
Mass balance, 38
Matrix
Banded, 46, 119, 121, 122
Positive Definite, 19, 44, 46, 119, 126, 128,
131
Singular, 43, 46
Symmetric, 44, 46, 67, 69, 119, 126
Mesh
refinement, 106
Method

1998-2010, S. D. Raj an 134
Direct stiffness, 31, 52
Least Squares, 55
Rayleigh-Ritz, 31
Subdomain, 54
Variational, 53
Weighted residuals, i, 6, 23, 51, 53, 54, 63
Nodes, 26, 27, 31, 33, 35, 40, 41, 45, 46, 47,
48, 50, 70, 73, 80, 85, 94, 99, 100, 101, 102,
107, 116, 120
Notation, 41, 67, 74, 79, 107
Optimizing criterion, 53, 54
p-convergence, 107
p-finite elements, 24
Polynomial, 20, 60, 63, 64, 67, 78, 99, 100,
107, 116
Post-processing, 110
Pre-processing, 110
Pressure, 37
Primary unknown, 31, 33, 47, 48, 77, 119
Residual Equations, 66, 77, 79
Review Exercises, i, 6, 7, 14, 49, 64, 78, 111
Shape Function, 71, 76, 78, 99, 100, 101, 110
Solid mechanics, 1, 23, 24, 32, 82, 115
stiffness matrix, 34, 41, 46, 48, 67, 68, 77, 84,
107
Stiffness matrix
Element, 34, 46
System, 46
Stiffness stiffness
System, 120
Stress, 8, 14, 24
Support reactions, 43, 45, 84, 87, 111
System stiffness matrix, 46
Temperature, 5, 16, 35, 36, 89, 90, 91, 92, 94,
96, 108, 111, 112, 113, 114
Theoretical Technique, 56, 59, 63, 66
Thermal Conductivity, 16, 35, 89, 94
Trial solution, 51, 52, 53, 54, 56, 57, 59, 60, 62,
63, 64, 66, 67, 69, 70, 73, 76, 77, 78, 79, 100,
107, 110
Variational approach, 23
Voltage, 39