Chapter 7 - Finite Element Method For Solid Mechanics (2016!10!10)

MECH3361/9361 Semester 2, 2016
FINITE ELEMENT METHOD FOR SOLID

MECHANICS
7.1. Introduction
7.1.1. Concept
The finite element (FE) method is based on the idea of building or dividing a complicated
object from or into small and manageable pieces, or elements. This approach to problem
solving is ubiquitous in daily life. Some examples include using bricks to construct buildings,
and developing children’s spatial awareness with Lego. A more mathematical example is the
approximation of the area of a circle by using triangles, as shown in Fig. 7.1.
R
θi
Fig 7.1 Approximation of a circular area using triangles
1
For a single triangle, the area Si can be calculated as 𝑆𝑖 = 2 𝑅2 sin𝜃. Thus, the approximate
area of the circle can be calculated as:
𝑁
1 2𝜋
𝑆𝑁 = ∑ 𝑆𝑖 = 𝑁 𝑅2 sin → 𝜋𝑅2 (when 𝑁 → ∞)
2 𝑁
𝑖=1
where N is the number of triangles. Thus, as we increase the number of triangles (N), the
approximate solution approaches the theoretical solution (𝜋𝑅2). In this case, we use
triangular elements to estimate the area. The type of element used in the FE method depends
on the geometry and physics being represented. Geometrically, elements can be 1D, 2D, or
3D as shown in Fig. 7.2.
1
node
element node
(a) 1D (b) 2D (c) 3D
Fig 7.2 Different types of elements
Once an element type has been chosen, the FE method follows the basic procedure outlined
below.
1. Divide structure into pieces (elements with nodes)

2. Describe the behaviours of the physical quantities on each element
3. Connect (assemble) the elements at the nodes to form an approximate system of
equations for the whole structure
4. Solve the system of equations involving unknown quantities at the nodes (e.g.
displacements)
5. Calculate desired quantities (e.g. strains and stresses) at selected elements
When implemented in computer software, these procedures fall into three main stages:
The pre-processor builds the FE model, material properties, loads and boundary
conditions;
The solver assembles and solves the system of equations; and
The post-processor sorts and displays the results.
7.1.2. Rationale
Analytical solutions to the equilibrium equations in 2D can be found in the form of Airy
stress functions, as shown in chapter 5. For complicated geometries and loading scenarios,
several stress functions can be superimposed to provide theoretical stress distributions, but
may still be subject to assumptions about the structure. Alternatively, experiments may be
performed to quantify the mechanical responses to specific loading cases. A properly
designed experiment may give the most realistic estimation of stress and displacement,
however when performed in isolation, experiments provide relatively limited data for the cost
and time expended.
The FE method is the most widely applied numerical method applied in mechanical
engineering. Being a numerical method, the FE method can be easily implemented in
computer software. This ties in well with the advent of computer aided design (CAD), and as
a result, it is relatively easy to perform mechanical analyses within the modern design
environment. Due to the nature of the FE method, more realistic loading and geometries can
be modelled when compared to analytical methods, however FE solutions will still be subject
to some abstraction of loading, and the inherent errors associated with discretisation.
Considering the aforementioned advantages and disadvantages of each approach, a full

mechanical analysis should be performed using a combination of theoretical analytical
solutions, experimental data, and numerical simulation. This will result in a greater
understanding of the mechanical characteristics of the problem, and enable improved design.
2
7.1.3. Brief history of FEA

Below is a short timeline outlining the development of FE analysis (FEA).
1943: Courant (Variational methods)

1956: Turner, Clough, Martin and Topp (Stiffness)
1960: Clough (“Finite Element”, plane problems)
1970s: Applications on mainframe computers
1980s: Microcomputers, pre- and postprocessors
1990s: Analysis of large structural systems
Onwards from the late 1960s, the development of FEA was closely tied to the development of
computers and specialised code bases for the FE method. Many of these codes are still
available commercially today as part of large FE software suites. Many newer software
packages have been developed to fulfil specific stages (see section 7.1.1) of the FE method.
1960s: NASTRAN = NASA STRuctural ANalysis

1970s: ANSYS, ABAQUS, MARC, DYNA3D
1980s: COMSOL, Altair, STRAND6-7
2000s: COSMOS/Solidworks, open-source
7.2. 1D spring elements

To introduce the FE method, we will consider a spring, which is the most basic of elements
modelling the linear relationship between force and displacement (Hooke’s law), discussed at
the beginning of chapter 3.
7.2.1. Single spring system
Fig 7.3 Single spring element
The single spring system illustrated in Fig. 7.3 has the following features:
Two nodes: i and j

Nodal displacements: ui and uj (displacement in the x direction in m)
Nodal forces: fi and fj (force in the x direction in N)
Spring constant (stiffness): k (N/m)
For solid mechanics problems, the displacement is known as the degree of freedom (dof).
For a spring element, there is only one displacement at each node, and hence one dof per
node. The total degrees of freedom is the total number of individual displacement
components in the system.
3
As discussed, this spring element models Hooke’s law in 1D, which is expressed by the
following spring force-displacement relationship.
𝑓𝑠𝑝𝑟𝑖𝑛𝑔 = 𝑘𝛥
where the 𝛥 = (𝑢𝑗 − 𝑢𝑖 ) is the elongation. The spring element does not model the non-linear
behaviour shown in Fig. 7.4.
fspring
Fig 7.4 The spring force-displacement relationship
Let’s formulate the equation set that defines the behaviour of this spring element. Consider
the equilibrium of the forces at nodes i and j.
At node i:
𝑓𝑖 + 𝑓𝑠𝑝𝑟𝑖𝑛𝑔 = 0
𝑓𝑖 = −𝑓𝑠𝑝𝑟𝑖𝑛𝑔 = −𝑘Δ = −𝑘(𝑢𝑗 − 𝑢𝑖 ) = 𝑘𝑢𝑖 − 𝑘𝑢𝑗
At node j:
𝑓𝑗 − 𝑓𝑠𝑝𝑟𝑖𝑛𝑔 = 0
𝑓𝑗 = 𝑓𝑠𝑝𝑟𝑖𝑛𝑔 = 𝑘Δ = 𝑘(𝑢𝑗 − 𝑢𝑖 ) = −𝑘𝑢𝑖 + 𝑘𝑢𝑗
This can also be expressed in matrix form:
𝑘 −𝑘 𝑢𝑖 𝑓𝑖
[ ] {𝑢 } = { }
−𝑘 𝑘 𝑗 𝑓𝑗
or
𝐤𝐮 = 𝐟
where:
k is the elemental stiffness matrix

u is the (element’s nodal) displacement vector
f is the (element’s nodal) force vector
Note that k is symmetric. In this case, k is singular, meaning its determinant is zero.
4
7.2.2. Multiple spring system

Now let’s consider a system of two spring connected together at a single node as shown in
Fig. 7.5. Globally, the nodes in this system are numbered 1, 2, and 3, with their
corresponding displacements (u) and forces (f).
Fig 7.5 Multiple spring system
For element 1:
(1)
𝑘 −𝑘1 𝑢1 𝑓
[ 1 ] {𝑢 } = { 1(1) }
−𝑘1 𝑘1 2 𝑓 2
For element 2:
(2)
[ 2 ] {𝑢 } = { 1(2) }
−𝑘2 𝑘2 3 𝑓 2
(𝑚)
where 𝑓i is the internal force acting on the local node i of element m (i = 1, 2).
Consider equilibrium of the spring system at nodes: (F1, F2, F3 are total forces in the nodes).
(1)
At node 1: 𝐹1 = 𝑓1
(1) (2)
At node 2: 𝐹2 = 𝑓2 + 𝑓1
(2)
At node 3: 𝐹3 = 𝑓2
The following equations can then be derived for the global equilibrium.
𝐹1 = 𝑘1 𝑢1 − 𝑘1 𝑢2
𝐹2 = (−𝑘1 𝑢1 + 𝑘1 𝑢2 ) + (𝑘2 𝑢2 − 𝑘2 𝑢3 ) = −𝑘1 𝑢1 + (𝑘1 + 𝑘2 )𝑢2 − 𝑘2 𝑢3
𝐹3 = −𝑘2 𝑢2 + 𝑘2 𝑢3
Or in matrix form:
𝑘1 −𝑘1 0 𝑢1 𝐹1
[−𝑘1 𝑘1 + 𝑘2 −𝑘2 ] {𝑢2 } = {𝐹2 }
0 −𝑘2 𝑘2 𝑢3 𝐹3
or
5
𝐊𝐔 = 𝐅
where K is the global stiffness matrix (structure matrix) for the spring system.
An alternative way of assembling the global stiffness matrix involves “enlarging” the
matrices for elements 1 and 2, and superimposing these matrices on their shared nodes.
For element 1:
(1)
[ 1 ] {𝑢 } = { 1(1) }
−𝑘1 𝑘1 2 𝑓 2
(1)
𝑘1 −𝑘1 0 𝑢1 𝑓1 → Node 1
[−𝑘1 𝑘1 0 ] { 𝑢 2 } = { (1) } → Node 2
𝑓2
0 0 0 𝑢3 0 → Node 3
For element 2:
(2)
[ 2 ] {𝑢 } = { 1(2) }
−𝑘2 𝑘2 3 𝑓 2
0 0 0 𝑢1 0 → Node 1
(2)
[0 𝑘2 −𝑘2 ] {𝑢2 } = {𝑓1 } → Node 2
0 −𝑘2 𝑘2 𝑢3 𝑓
(2)
→ Node 3
2
Adding these two expanded equations results in our global equilibrium equation.
𝑘1 −𝑘1 0 0 0 0 𝑢1 𝑓1
(1) 0
(2)
([−𝑘1 𝑘1 0] + [0 𝑘2 −𝑘2 ]) {𝑢2 } = {𝑓 (1) } + {𝑓1 }
2
0 0 0 0 −𝑘2 𝑘2 𝑢3 𝑓2
(2)
0
(1)
𝑘1 −𝑘1 0 𝑢1 𝑓1 𝐹1
[−𝑘1 𝑘1 + 𝑘2 −𝑘2 ] {𝑢2 } = {𝑓2(1) + 𝑓1(2) } = {𝐹2 }
0 −𝑘2 𝑘2 𝑢3 (2) 𝐹3
𝑓2
7.2.3. Boundary and load conditions

Now that we have a mathematical description this spring system, loads and boundary
conditions can be applied to the nodes to replicate a mechanical system.
Reaction
force F
u1=0 u2 , P u3, P
Fig 7.6 Boundary conditions and loads applied to a multiple spring system
6
Consider the problem illustrated in Fig. 7.6. At the left end, node 1 has been fixed (𝑢1 = 0)
and the system is loaded by external forces at nodes 2 and 3 (𝐹2 = 𝐹3 = 𝑃). Using the global
equilibrium equation derived for this multiple spring system, the displacements at nodes 2
and 3 (𝑢2 , 𝑢3) and the reaction force at node 1 (𝐹1 ).
The global equilibrium equation can be written as:
𝑘1 −𝑘1 0 0 𝐹1
[−𝑘1 𝑘1 + 𝑘2 −𝑘2 ] {𝑢2 } = {𝑃 }
0 −𝑘2 𝑘2 𝑢3 𝑃
Since 𝑢1 = 0, the equilibrium equation can be reduced into two separate equations:
𝑘 + 𝑘2 −𝑘2 𝑢2 𝑃
[ 1 ] {𝑢 } = { }
−𝑘2 𝑘2 3 𝑃
and
−𝑘1 𝑢2 = 𝐹1
𝑢2
The unknowns are 𝐮 = {𝑢 } and reaction force 𝐹1 . Solving the first equation set, we obtain
3
the displacements.
𝑢2 2𝑃/𝑘1
{𝑢 } = { }
3 2𝑃/𝑘1 + 𝑃/𝑘2
The reaction force can be determined from the second equation.
𝐹 = 𝐹1 = −𝑘1 𝑢2 = −2𝑃
Some remarks can be made from this fundamental example:
Both the elemental stiffness and global stiffness matrices are symmetric and positive
definite;
The calculated force satisfies the equilibrium of the external forces; and
Spring elements are easily adapted for other stiffness analyses, such as torsion.
Example 7.1
For a three-spring system shown, 𝑘1 = 100 N/mm, 𝑘2 = 200 N/mm, 𝑘3 = 100 N/mm,
𝑃 = 500 N and 𝑢1 = 𝑢4 = 0. Determine: (a) global stiffness matrix, (b) 𝑢2 and 𝑢3, (c)
reaction forces at nodes 1 and 4, (d) force in spring #2.
7
Solution
Step 1: Elemental stiffness matrices
100 −100 200 −200 100 −100

k1 = [ ], k 2 = [ ], k 3 = [ ]
−100 100 −200 200 −100 100
Step 2: Expand the elemental stiffness matrices according to the node numbers
100 −100 0 0
−100 100 0 0
k̂1 = [ ]
0 0 0 0
0 0 0 0
0 0 0 0
0 200 −200 0
k̂ 2 = [ ]
0 −200 200 0
0 0 0 0
0 0 0 0
0 0 0 0
k̂ 3 = [ ]
0 0 100 −100
0 0 −100 100
Step 3: Superposition (add) all elemental stiffness matrices for global stiffness matrix
100 100 0 0
100 100 200 200 0
K kˆ1 kˆ2 kˆ3
0 200 200 100 100
0 0 100 100
Thus the global stiffness matrix:
100 −100 0 0
−100 300 −200 0
𝐾=[ ]
0 −200 300 −100
0 0 −100 100
which is symmetric and banded.
Step 4: Global equilibrium equation
8
(1)
𝑢1 𝑓1
100 −100 0 0 𝐹1
(1) (2)
−100 300 −200 0 𝑢2 𝑓2 + 𝑓1 0
[ ] {𝑢 } = ={ }
0 −200 300 −100 3 (2)
𝑓2 + 𝑓1
(3) 𝑃
0 0 −100 100 𝑢 4 (3) 𝐹4
{ 𝑓2 }
Step 5: Apply boundary condition, 𝑢1 = 𝑢4 = 0
100 −100 0 0 0 𝐹1
−100 300 −200 0 𝑢 0
[ ] { 2} = { }
0 −200 300 −100 𝑢3 𝑃
0 0 −100 100 0 𝐹4
Delete (or separate) the 1 st and 4th rows and columns, we have non-vanishing equations
300 −200 𝑢2 0
[ ]{ } = { }
−200 300 𝑢3 𝑃
or
300𝑢2 − 200𝑢3 = 0
−200𝑢2 + 300𝑢3 = 𝑃
Step 6: Solve the linear equation set to obtain unknown displacements
𝑢2 𝑃/250 2
{𝑢 } = { } = { } mm
3 3𝑃/500 3
Step 7: Determine reaction forces from the “dropped” equations
−100𝑢2 = 𝐹1
−100𝑢3 = 𝐹4
𝐹1 = −100𝑢2 = −100×2 = −200𝑁
𝐹4 = −100𝑢3 = −100×3 = −300𝑁
Example 7.2
For the spring system shown with arbitrarily numbered nodes and elements, find the global
stiffness matrix.
9
Solution
Step 1: Construct an element connectivity table (4 elements and 5 nodes)
Element number Local node i (1) Local node j (2)
1 4 2
2 2 3
3 3 5
4 2 1
This table associates the global node numbers with the local node numbers for each
element.
𝑘 −𝑘1 𝑘 −𝑘2 𝑘 −𝑘3 𝑘 −𝑘4

k1 = [ 1 ], k 2 = [ 2 ], k 3 = [ 3 ], k 4 = [ 4 ]
⏟−𝑘1 𝑘1 ⏟−𝑘2 𝑘2 ⏟−𝑘3 𝑘3 ⏟−𝑘4 𝑘4
𝑢4 𝑢2 𝑢2 𝑢3 𝑢3 𝑢5 𝑢2 𝑢1
Step 3: Expand the elemental stiffness matrices. For the un-ordered stiffness matrices, such
as:
𝑘 −𝑘1
k1 = [ 1 ]
⏟−𝑘1 𝑘1
𝑢4 𝑢2
(1)
𝑘1 𝑢4 − 𝑘1 𝑢2 = 𝑓1
(1)
−𝑘1 𝑢4 + 𝑘1 𝑢2 = 𝑓2
If we re-order these equations to have 𝑢2 before 𝑢4 (in order of global node number)
(1)
𝑘1 𝑢2 − 𝑘1 𝑢4 = 𝑓2
(1)
𝑘1 𝑢4 − 𝑘1 𝑢2 = 𝑓1
𝑘 −𝑘1
Therefore the re-ordered stiffness matrix is k1 = [ 1 ].
⏟−𝑘1 𝑘1
𝑢2 𝑢4
10
𝑘 −𝑘4 𝑘 −𝑘4
Similarly for k 4 = [ 4 ], the re-ordered stiffness matrix is k 4 = [ 4 ].
⏟−𝑘4 𝑘4 ⏟−𝑘4 𝑘4
𝑢2 𝑢1 𝑢1 𝑢2
In expanded form, the elemental stiffness matrices are as follows.
0 0 0 0 0 0 0 0 0 0
0 𝑘1 0 −𝑘1 0 0 𝑘2 −𝑘2 0 0
k1 = 0 0 0 0 0 k 2 = 0 −𝑘2 𝑘2 0 0
0 −𝑘1 0 𝑘1 0 0 0 0 0 0
[0 0 0 0 0] [0 0 0 0 0]
𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5
0 0 0 0 0 𝑘4 −𝑘4 0 0 0
0 0 0 0 0 −𝑘4 𝑘4 0 0 0
0 0 k3 0 −k 3 0 0 0 0 0
k3 = k4 =
0 0 0 0 0 0 0 0 0 0
[0 0 −k 3 0 k3 ] [ 0 0 0 0 0]
𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5
Step 4: Superposition to form the global stiffness matrix
𝑘4 −𝑘4 0 0 0
−𝑘4 𝑘1 + 𝑘2 + 𝑘4 −𝑘2 −𝑘1 0
K= 0 −𝑘2 𝑘2 + 𝑘3 0 −𝑘3
0 −𝑘1 0 𝑘1 0
[ 0 0 −𝑘3 0 𝑘3 ]
𝑢1 𝑢2 𝑢3 𝑢4 𝑢5
7.3. 1D bar elements

Spring elements are purely 1D, and are independent of geometrical characteristics, such as
length and area. To introduce these factors into the FE elements, we will consider bar
elements.
7.3.1. Direct stiffness formulation
Fig 7.7 Single bar element
Consider the single bar element shown in Fig. 7.7. It has the following features:
11
Two nodes: i and j

Nodal displacements: ui and uj (displacement in the x direction in m)
Nodal forces: fi and fj (force in the x direction in N)
Cross-sectional area: A
Length: L
Young’s modulus: E
Strain: 𝜀 = 𝜀(𝑥)
Stress: 𝜎 = 𝜎(𝑥)
Like the spring element, the bar element models Hooke’s law, but in terms of stress and
strain:
𝜎 = 𝐸𝜀
where strain can be derived from the displacement using the strain-displacement relationship.
𝑑𝑢
𝜀=
𝑑𝑥
With these relationships in mind, let’s derive the stiffness matrix for a single bar element.
Step 1: Assuming the displacement u varies linearly along the length of the bar:
𝑥 𝑥
𝑢(𝑥) = (1 − ) 𝑢𝑖 + 𝑢𝑗
𝐿 𝐿
When 𝑥 = 0, 𝑢(𝑥 = 0) = 𝑢𝑖 ; When 𝑥 = 𝐿, 𝑢(𝑥 = 𝐿) = 𝑢𝑗
Step 2: Strain is then derived using the strain-displacement relationship, where 𝑑𝑢 is the
elongation 𝛥, defined as the relative displacement between node i and j.
𝛥 𝑢𝑗 − 𝑢𝑖
𝜀= =
𝐿 𝐿
Stress can be found from Hooke’s law.
𝐸𝛥
𝜎 = 𝐸𝜀 =
𝐿
𝐹
For a bar, 𝜎 = , therefore:
𝐴
𝐸Δ 𝐸𝐴
𝐹 = σ𝐴 = ( ) 𝐴 = ( ) Δ = 𝑘Δ
𝐿 𝐿
𝐸𝐴
where 𝑘 = is the stiffness of the bar.
𝐿
Step 3: In this case, the bar acts like the spring element discussed in section 7.2, and as a
result, the elemental stiffness matrix takes on the same form.
12
𝐸𝐴 𝐸𝐴
−
𝐤=[
𝑘 −𝑘
]=[ 𝐿 𝐿 ] = 𝐸𝐴 [ 1 −1
]
−𝑘 𝑘 𝐸𝐴 𝐸𝐴 𝐿 −1 1
−
𝐿 𝐿
Step 4: The equilibrium equation also follows the same form.
𝐸𝐴 1 −1 𝑢𝑖 𝑓𝑖
[ ] {𝑢 } = { }
𝐿 −1 1 𝑗 𝑓𝑗
From this formulation, some remarks can be made.
Like the spring element, the degree of freedom (dof) is the nodal displacement
component (u). For a 1D bar element, there is one dof per node.
The coefficients in k can be interpreted with a physical meaning. The jth column of k
(here j = 1 or 2) represents the forces applied to the bar to maintain a deformed shape
with unit displacement at node j and zero displacement at the other node.
7.3.2. Energy formulation

The direct stiffness formulation provides a simple way to define the equilibrium equation for
simple “spring-like” 1D elements such as the spring or bar elements. For more complex
elements, a formal approach is needed. This section will derive the stiffness matrix using
energy principles.
Step 1: The displacement within the bar element can be described with two linear shape
functions. These are interpolation functions 𝑁𝑖 and 𝑁𝑗 .
𝑁𝑖 (𝜉) = 1 − 𝜉
𝑁𝑗 (𝜉) = 𝜉
𝑥
𝜉 is the natural (local) co-ordinate within the element. 𝜉 = 𝐿 , 𝑤ℎ𝑒𝑟𝑒 0 ≤ 𝜉 ≤ 1. The use of a
natural co-ordinate allows us to change the co-ordinate system from [𝑥1, 𝑥2 ] to [0,1], which is
independent of the physical length of the element.
Displacement can be expressed in terms of these shape functions.
𝑢(𝑥) = 𝑢(𝜉) = 𝑁𝑖 (𝜉)𝑢𝑖 + 𝑁𝑗 (𝜉)𝑢𝑗
When 𝜉 = 0, 𝑢(𝜉 = 0) = 𝑁𝑖 (0)𝑢𝑖 + 𝑁𝑗 (0) = 𝑢𝑖

When 𝜉 = 1, 𝑢(𝜉 = 1) = 𝑁𝑖 (1)𝑢𝑖 + 𝑁𝑗 (1) = 𝑢𝑗
In matrix form:
𝑢𝑖
𝑢 = [𝑁𝑖 𝑁𝑗 ] {𝑢 } = 𝐍𝐮
𝑗
13
𝑢𝑖
where 𝐮 = {𝑢 } is the nodal displacement vector.
𝑗
Step 2: The strain-displacement relationship can also be expressed in terms of these shape
functions (𝐍).
𝑑𝑢 𝑑 𝑑
ε= = (𝐍𝐮) = ( 𝐍) 𝐮 = 𝐁𝐮
𝑑𝑥 𝑑𝑥 𝑑𝑥
Note that 𝐮 is independent of x, as the variation of displacement within the bar element is
𝑑
expressed within the shape functions 𝐍. 𝐁 = (𝑑𝑥 𝐍) is known as the strain-displacement
matrix, and is derived for the bar element below.
𝑑 𝑑
𝐁=( 𝐍) = [𝑁𝑖 (𝜉) 𝑁𝑗 (𝜉)]
𝑑𝑥 𝑑𝑥
𝑑 𝑑𝜉
= [𝑁𝑖 (𝜉) 𝑁𝑗 (𝜉)] ⋅
𝑑𝜉 𝑑𝑥
𝑑 𝑑 𝑥
= [1 − 𝜉 𝜉] ⋅ ( )
𝑑ξ 𝑑𝑥 𝐿
1
= [−1 1]
𝐿
1 1
= [− ]
𝐿 𝐿
Step 3: Stress can then be calculated using Hooke’s law.
σ = 𝐸ε = 𝐸𝐁𝐮
Step 4: Apply the strain energy-work principle. Here we calculate the strain energy (𝑈) and
equate it to the work done (𝑊) by the forces on the element (all work is stored as strain).
Strain energy can be calculated by:
1
𝑈= ∫ σT ε 𝑑𝑉
2
𝑉
1
= ∫(𝐸𝐁𝐮)T (𝐁𝐮) 𝑑𝑉
2
𝑉
1
= ∫(𝐮T 𝐁 T )𝐸(𝐁𝐮) 𝑑𝑉
2
𝑉
1 T
= 𝐮 [∫(𝐁 T 𝐸𝐁) 𝑑𝑉 ] 𝐮
2
𝑉
Work done for the nodal forces is:
1 1 1 1
𝑊= 𝑓𝑖 𝑢𝑖 + 𝑓𝑗 𝑢𝑗 = {𝑢𝑖 𝑢𝑗 } {𝑓𝑖 } = 𝐮T 𝐟
2 2 2 𝑓𝑗 2
14
Equating these two quantities (𝑈 = 𝑊):
1 T 1
𝐮 [∫(𝐁 T 𝐸𝐁) 𝑑𝑉] 𝐮 = 𝐮T 𝐟
2 2
𝑉
∴ [∫(𝐁 T 𝐸𝐁) 𝑑𝑉 ] 𝐮 = 𝐟
𝑉
This is in the form of 𝐤𝐮 = 𝐟, and the coefficient of 𝐮 is the elemental stiffness matrix.
𝐤 = ∫(𝐁 T 𝐸𝐁) 𝑑𝑉
𝑉
(7.1)
To show that this gives the same stiffness matrix as the direct stiffness formulation, substitute
in the 𝐁 matrix found earlier.
1 1T 1 1
𝐤 = ∫ [− ] 𝐸 [− ] (𝐴 𝑑𝑥)
𝐿 𝐿 𝐿 𝐿
𝑉
1 1
2
−
= 𝐸𝐴 𝐿 𝐿2 ∫ 𝑑𝑥
1 1
− 𝑉
[ 𝐿2 𝐿2 ]
1 1
− 2
= 𝐸𝐴 𝐿2 𝐿 𝐿
1 1
− 2
[ 𝐿 𝐿2 ]
1 1
−
𝐤 = 𝐸𝐴 [ 𝐿 𝐿]
1 1
−
𝐿 𝐿
Eq. 7.1 is a general result which can be used for the construction of other types of elements.
This expression can also be derived using other more rigorous approaches, such as the
Principle of Minimum Potential Energy, or Galerkin’s Method.
Using Eq. 7.1, strain energy in the element can be written as:
1 T
𝑈= 𝐮 𝐊𝐮
2
15
Example 7.3
The three bar structure is fully clamped at ends A and D as shown. Solve the following
problems by using finite element method. (a) Write the elemental stiffness matrices. (b) Give
the global equilibrium equation. (c) Calculate the displacement at nodes B and C. (d)
Calculate the reaction forces at points A and D.
A EA=2 B EA=1 C EA=2 D x

P 2P
1m 1m 1m
Solution
Step 1: Connectivity
Elements: AB = (1), BC = (2), CD = (3)
Nodes: A=1, B=2, C=3, D=4.
𝐸𝐴 1 −1
𝐤= [ ]
𝐿 −1 1
2 −2 1 −1 2 −2
𝐤𝟏 = [ ] 𝐤𝟐 = [ ] 𝐤𝟑 = [ ]
−2 2 −1 1 −2 2
Step 3: Expand the elemental stiffness matrices according to the global node numbers
𝑢1
𝑢2
𝐮 = {𝑢 }
3
𝑢4
2 −2 0 0 0 0 0 0
̂𝟏=[−2 2 0 0 ̂ 𝟐 = [0 1 −1 0
𝐤 ] 𝐤 ]
0 0 0 0 0 −1 1 0
0 0 0 0 0 0 0 0
0 0 0 0
̂ 𝟑 = [0
𝐤
0 0 0
]
0 0 2 −2
0 0 −2 2
Step 4: Superposition of elemental stiffness matrices to form the global stiffness matrix
16
2 2 0 0
2 2 1 1 0
K kˆ1 kˆ2 kˆ3
0 1 1 2 2
0 0 2 2
Therefore the global stiffness matrix is:
2 −2 0 0
−2 3 −1 0
𝐊=[ ]
0 −1 3 −2
0 0 −2 2
which is symmetric and banded.
Step 5: Global equilibrium equation 𝐊𝐔 = 𝐅
2 −2 0 0 𝑢1 𝐹1
−2 3 −1 0 𝑢 2 𝐹
[ ] {𝑢 } = { 2 }
0 −1 3 −2 3 𝐹3
0 0 −2 2 𝑢 4 𝐹4
Step 6: Apply the boundary conditions (𝑢1 = 𝑢4 = 0, 𝐹2 = 𝑃, 𝐹3 = 2𝑃)
2 −2 0 0 0 𝐹1
−2 3 −1 0 𝑢 𝑃
[ ] { 2} = { }
0 −1 3 −2 𝑢3 2𝑃
0 0 −2 2 0 𝐹4
Dropping the 1st and 4th rows and columns, we obtain the following linear set of equations.
3 −1 𝑢2 𝑃 3𝑢2 − 𝑢3 = 𝑃
[ ]{ } = { } or
−1 3 𝑢3 2𝑃 −𝑢2 + 3𝑢3 = 2𝑃
Step 7: Solve reduced linear equation set for displacement
𝑢2 5𝑃/8
{𝑢 } = { }m
3 7𝑃/8
Step 8: Determine reaction forces from dropped equations
5P
𝐹1 = −2𝑢2 = −
4
7𝑃
𝐹4 = −2𝑢3 = −
4
17
7.3.3. Bar elements in 2D space

So far, our treatment of bar elements has been restricted to a 1D environment, where
boundary conditions are applied only axial direction. In this section, we will consider a bar
element that can be subject to constraints and loading in two dimensions (x and y).
Step 1: Establish the generalised element parameters.
vi ’
Fig 7.8 2D bar element
Local Global
Dimensions x, y X, Y
Displacements 𝒖′𝒊 , 𝒗′𝒊 𝒖𝒊 , 𝒗𝒊

DOFs 1 DOF per node 2 DOF per node
Step 2: Transform nodal displacements and forces from the global to the local co-ordinate
system.
𝑢𝑖
𝑢𝑖′ = 𝑢𝑖 cos 𝜃 + 𝑣𝑖 cos(90° − 𝜃) = 𝑢𝑖 cos 𝜃 + 𝑣𝑖 sin 𝜃 = [𝑙 𝑚] {𝑣 }
𝑖
𝑢𝑖
𝑣𝑖′ = −𝑢𝑖 sin 𝜃 + 𝑣𝑖 cos 𝜃 = [−𝑚 𝑙 ] {𝑣 }
𝑖
𝑋𝑗 −𝑋𝑖 𝑌𝑗 −𝑌𝑖
where the direction cosines are: 𝑙 = cos 𝜃 = and 𝑚 = cos(90° − 𝜃) = sin 𝜃 =
𝐿 𝐿
In matrix form, this transformation from global to local displacements at a node can be
represented as:
𝑢′ 𝑙 𝑚 𝑢𝑖
{ 𝑖′ } = [ ]{ }
𝑣𝑖 −𝑚 𝑙 𝑣𝑖
or
̃𝐮𝐢
𝐮′𝐢 = 𝐓
̃=[ 𝑙
where the transformation matrix 𝐓
𝑚 ̃ −𝟏 = 𝐓
] is orthogonal, which means 𝐓 ̃𝑇 .
−𝑚 𝑙
18
For two nodes of a single bar element, the following transformation can be expressed.
𝑢𝑖′ 𝑙 𝑚 0 0 𝑢𝑖
𝑣𝑖′ −𝑚 𝑙 0 0 𝑣𝑖
𝑢𝑗′ = [ 0 ]{ }
0 𝑙 𝑚 𝑢𝑗
′ 0 0 −𝑚 𝑙 𝑣𝑗
{𝑣𝑗 }
or
̃
𝐮′ = 𝐓𝐮 with 𝐓 = [𝐓 0]
0 ̃
𝐓
Nodal forces can be transformed in the same way.
𝑓𝑥𝑖′ 𝑓𝑥𝑖
𝑙 𝑚 0 0
𝑓𝑦𝑖′ −𝑚 𝑙 0 0 𝑓𝑦𝑖
′ =[ ]
𝑓𝑥𝑗 0 0 𝑙 𝑚 𝑓𝑥𝑗
′ 0 0 −𝑚 𝑙
{𝑓𝑦𝑗 } {𝑓𝑦𝑗 }
or
𝐟 ′ = 𝐓𝐟 where 𝐟 ′ and 𝐟 are the local and global force vectors, respectively.
Step 3: Generate the elemental stiffness matrices. Since the lateral displacement 𝑣𝑖′ does not
contribute to the stretch of the bar within the linear theory, in the local coordinate system (x-
y) we have:
′
𝑘 −𝑘1 𝑢1 𝐸𝐴 1 −1 𝑢𝑖 𝑓𝑥𝑖′
[ 1 ] {𝑢 } = [ ] { ′ } = { ′ }
−𝑘1 𝑘1 2 𝐿 −1 1 𝑢𝑗 𝑓𝑥𝑗
We can then expand this equation to include all displacement components:
1 0 −1 0 𝑢𝑖′ 𝑓𝑥𝑖′
′
𝐸𝐴 0 0 0 0 𝑣𝑖 0
[
𝐿 −1
] 𝑢𝑗′ = {𝑓𝑥𝑗
′ } or 𝐤 ′ 𝐮′ = 𝐟 ′
0 1 0
0 0 0 0 ′
{𝑣𝑗 } 0
Using the transformations derived in step 2, we express this equation in terms of the global
co-ordinate system (X-Y).
𝐤 ′ (𝐓𝐮) = (𝐓𝐟)
As 𝐓 is orthogonal (𝐓 𝑇 𝐓 = 𝐈), multiplying both sides of this equation by 𝐓 𝑇 yields the

elemental equilibrium equation for a 2D bar element.
𝐓 𝑇 𝐤 ′ (𝐓𝐮) = 𝐟
19
(𝐓 𝑇 𝐤 ′ 𝐓)𝐮 = 𝐟
in the form of 𝐤𝐮 = 𝐅
Thus, the elemental stiffness matrix 𝐤 in the global co-ordinate system is:
𝐤 = 𝐓𝑇 𝐤′𝐓
or
𝑙2 𝑙𝑚 −𝑙 2 −𝑙𝑚
𝐸𝐴 𝑙𝑚 𝑚2 −𝑙𝑚 −𝑚2 ]
𝐤= [
𝐿 −𝑙 2 −𝑙𝑚 𝑙2 𝑙𝑚
−𝑙𝑚 −𝑚2 𝑙𝑚 𝑚2
(7.2)
The structural global stiffness matrix can then be assembled as for 1D bar elements.
Step 4: Calculate the elemental strain and stress. This is done using shape functions and the
strain-displacement matrix 𝐁.
𝑑 𝑑 1 1
𝐁=( 𝐍) = [𝑁𝑖 (𝜉) 𝑁𝑗 (𝜉)] = [− ]
𝑑𝑥 𝑑𝑥 𝐿 𝐿
𝑢𝑖′
𝜎 = 𝐸𝜀 = 𝐸 (𝐁 { ′ })
𝑢𝑗
𝑢𝑖
1 1 𝑙 𝑚 0 0 𝑣𝑖
= 𝐸 [− ][ ]{ }
𝐿 𝐿 0 0 𝑙 𝑚 𝑢𝑗
𝑣𝑗
𝑢𝑖
𝐸 𝑣𝑖
= [−𝑙 −𝑚 𝑙 𝑚] {𝑢 }
𝐿 𝑗
𝑣𝑗
Example 7.4
A simple plane truss is made of two identical bars (with E, A,
and L), and loaded as shown in the figure. Find the
displacements of node 2 and the stress in each bar.
Solution
The illustrated simple truss system is used to demonstrate the

assembly and solution process using the 1D bar element in 2D
space.
Step 1: Generate the elemental stiffness matrices in local co-ordinates.
20
For the local coordinate systems, we have:
𝐸𝐴 1 −1
𝐤 ′𝟏 = [ ] = 𝐤 ′𝟐
𝐿 −1 1
Step 2: Assemble the elemental stiffness matrix in global co-ordinates. These two elemental
stiffness matrices 𝐤 ′𝟏 and 𝐤 ′𝟐 cannot be assembled together directly, because their local co-
ordinate systems are different. Instead, both matrices need to be converted to the global co-
ordinate system (XOY).
1 1
For element 1, 𝜃 = 45°, 𝑙 = cos 45° = , 𝑚 = cos 45° = . From (7.2), we can derive the
√2 √2
elemental stiffness matrix in global co-ordinates.
𝑙2 𝑙𝑚 −𝑙 2 −𝑙𝑚 1 1 −1 −1
𝐸𝐴 𝑙𝑚 𝑚2 −𝑙𝑚 2
−𝑚 ] = 𝐸𝐴 1 1 −1 −1
𝐤𝟏 = 𝐓𝑇 𝐤′ 𝐓 = [ [ ]
𝐿 −𝑙 2 −𝑙𝑚 𝑙2 𝑙𝑚 2𝐿 −1 −1 1 1
−𝑙𝑚 −𝑚2 𝑙𝑚 𝑚2 −1 −1 1 1
In expanded form:
1 1 −1 −1 0 0
1 1 −1 −1 0 0
𝐸𝐴 −1 −1 1 1 0 0
̂𝟏=
𝐤
2𝐿 −1 −1 1 1 0 0
0 0 0 0 0 0
[0 0 0 0 0 0]
1 1
For element 2, 𝜃 = 135°, 𝑙 = cos 135° = − , 𝑚 = cos 45° =
√2 √2
𝑙2 𝑙𝑚 −𝑙 2 −𝑙𝑚 1 −1 −1 1
𝐸𝐴 𝑚2 2 𝐸𝐴 −1 1 1 −1
𝐤𝟐 = 𝐓𝑇 𝐤′ 𝐓 = [ 𝑙𝑚2 −𝑙𝑚 −𝑚 ] = [ ]
𝐿 −𝑙 −𝑙𝑚 𝑙2 𝑙𝑚 2𝐿 −1 1 1 −1
−𝑙𝑚 −𝑚2 𝑙𝑚 𝑚2 1 −1 −1 1
In expanded form:
0 0 0 0 0 0
0 0 0 0 0 0
𝐸𝐴 0 0 1 −1 −1 1
̂𝟐=
𝐤
2𝐿 0 0 −1 1 1 −1
0 0 −1 1 1 −1
[0 0 1 −1 −1 1 ]
Step 3: Assemble the global stiffness matrix and FE equilibrium equation. Superimposing the
two stiffness matrices yields the global stiffness matrix, which can then be installed into the
equilibrium equation.
21
1 1 −1 −1 0 0 𝑢1 𝐹1𝑋
1 1 −1 −1 0 0 𝑣1 𝐹1𝑌
𝐸𝐴 −1 −1 2 0 −1 1 𝑢2 𝐹
𝑣 = 2𝑋
2𝐿 −1 −1 0 2 1 −1 2 𝐹2𝑌
0 0 −1 1 1 −1 𝑢3 𝐹3𝑋
[0 0 1 −1 −1 1 ] {𝑣3 } {𝐹3𝑌 }
𝑢1 𝑣1 𝑢2 𝑣2 𝑢3 𝑣3
Step 4: Apply boundary conditions and loads.
𝑢1 = 𝑣1 = 𝑢3 = 𝑣3 = 0, 𝐹2𝑋 = 𝑃1 , 𝐹2𝑌 = 𝑃2
1 1 −1 −1 0 0 0 𝐹1𝑋
1 1 −1 −1 0 0 0 𝐹1𝑌
𝐸𝐴 −1 −1 2 0 −1 1 𝑢2 𝑃1
=
2𝐿 −1 −1 0 2 1 −1 𝑣2 𝑃2
0 0 −1 1 1 −1 0 𝐹3𝑋
[0 0 1 −1 −1 1 ] { 0 } {𝐹3𝑌 }
Step 5: Condense the FE equations by eliminating all zero-displacement rows and columns.
𝐸𝐴 2 0 𝑢2 𝑃
[ ] {𝑣 } = { 1 }
2𝐿 0 2 2 𝑃2
Step 6: Solve for the displacements at node 2.
𝑢2 𝐿 𝑃1
{𝑣 } = { }
2 𝐸𝐴 𝑃2
Step 7: Calculate the elemental stress in the two bars. From our earlier formulation of stress for
bar elements in 2D environments:
𝑢𝑖
𝐸 𝑣𝑖
𝜎 = 𝐸𝜀 = [−𝑙 −𝑚 𝑙 𝑚] {𝑢 }
𝐿 𝑗
𝑣𝑗
0
𝐸 1 𝐿 0 1
𝜎1 = [−1 −1 1 1] { }= (𝑃1 + 𝑃2 )
𝐿 √2 𝐸𝐴 𝑃1 𝐴 √2
𝑃2
𝑃1
𝐸 1 𝐿 𝑃2 1
𝜎2 = [1 −1 −1 1] { }= (𝑃1 − 𝑃2 )
𝐿 √2 𝐸𝐴 0 𝐴√2
0
22
7.4. 1D beam elements

Spring and bar elements displace only in their local axial directions. Now we will consider a
beam element, which can model transverse displacements, i.e. orthogonal to the local axial
direction.
7.4.1. Simple plane beam element
Mi EI Mj
x
(vi , i )
i j (v j , j)
L
Fi Fj
Fig 7.9 A simple plane beam element
The simple plane beam element shown in Fig. 7.9 has the following features.
Two nodes: i and j

Length: L
2nd moment of inertia of the cross-sectional area: 𝐼 = ∫𝐴 𝑦 2 𝑑𝐴
Young’s modulus: E
Deflection (transverse/lateral displacement) of the neutral axis: 𝑣 = 𝑣(𝑥)
Slope or rotation about the z-axis: 𝜃 = 𝑑𝑣/𝑑𝑥
Shear force: 𝐹 = 𝐹(𝑥)
Moment about the z-axis: 𝑀 = 𝑀(𝑥)
This beam element models the fundamental beam theory concepts as taught in Mechanics of
Solids 1, governed by the following deflection and flexure formulae.
𝑑2𝑣
𝐸𝐼 = 𝑀(𝑥)
𝑑𝑥 2
𝑀𝑦
𝜎=−
𝐼
7.4.2. Energy formulation

Due to the differences between “spring-like” and beam elements, the energy formulation is
used to derive the stiffness matrix and equilibrium equation in this case.
Step 1: Introduce the shape functions, which describe the variation in deflection in the beam
element.
3 2 2 3 2 2 1 3
𝑁1 (𝑥) = 1 − 𝑥 + 3𝑥 𝑁2 (𝑥) = 𝑥 − 𝑥 + 2𝑥
𝐿2 𝐿 𝐿 𝐿
3 2 2 3 1 2 1 3
𝑁3 (𝑥) = 2 𝑥 − 3 𝑥 𝑁4 (𝑥) = − 𝑥 + 2 𝑥
𝐿 𝐿 𝐿 𝐿
23
Step 2: Calculate the deflection. This is achieved by multiplying the shape functions (𝐍) to
the nodal displacement vector (𝐮). Due to the nature of the shape functions, the deflection is a
cubic function of x.
𝑣𝑖
𝜃𝑖
𝑣(𝑥) = 𝐍𝐮 = [𝑁1 (𝑥) 𝑁2 (𝑥) 𝑁3 (𝑥) 𝑁4 (𝑥)] {𝑣 }
𝑗
𝜃𝑗
This deflection function can be used to calculate the curvature of the beam, or slope. This is
similar to the element strain calculation performed in bar elements.
𝑑2𝑣 𝑑2 𝑑2
= (𝐍𝐮) = ( 𝐍) 𝐮 = 𝐁𝐮
𝑑𝑥 2 𝑑𝑥 2 𝑑𝑥 2
In this case, 𝐁 is the strain-displacement matrix, derived below.
𝑑2
𝐁= 𝐍 = [𝑁1′′ (𝑥) 𝑁2′′ (𝑥) 𝑁3′′ (𝑥) 𝑁4′′ (𝑥)]
𝑑𝑥 2
6 12𝑥 4 6𝑥 6 12𝑥 2 6𝑥
= [− 2 + 3 − + 2 2
− 3 − + 2]
𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿
Step 3: Calculate the strain energy stored in the beam element.
1 1 𝑀𝑦 𝑇 1 𝑀𝑦
𝑈= ∫ 𝜎 𝑇 𝜀 𝑑𝑉 = ∫ (− ) (− ) 𝑑𝐴𝑑𝑥
2 𝑉 2 𝑉 𝐼 𝐸 𝐼
1 1
= ∫𝑀𝑇 𝑀 𝑑𝑥
2 𝐿 𝐸𝐼
𝑇
1 𝑑2𝑣 𝑑2𝑣
= ∫ ( 2 ) 𝐸𝐼 ( 2 ) 𝑑𝑥
2 𝑑𝑥 𝑑𝑥
𝐿
1
= ∫(𝐁𝐮)𝑇 𝐸𝐼(𝐁𝐮) 𝑑𝑥
2
𝐿
1 𝑇 1
= 𝐮 ( ∫ 𝐁 𝑇 𝐸𝐼𝐁 𝑑𝑥) 𝐮
2 2
𝐿
This form of strain energy is similar to the form encountered in section 7.3.2. The stiffness
matrix for a simple plane beam element is therefore:
1
𝐤= ∫ 𝐁 𝑇 (𝐸𝐼)𝐁 𝑑𝑥
2
𝐿
This results in the equilibrium equation for a simple plane beam element.
24
12 6𝐿 −12 6𝐿 𝑣𝑖 𝐹𝑖
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃𝑖 𝑀𝑖
[ ] {𝑣 } = 𝐹
𝐿3 −12 −6𝐿 12 −6𝐿 𝑗 𝑗
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃𝑗 𝑀
{ 𝑗}
𝑣𝑖 𝜃𝑖 𝑣𝑗 𝜃𝑗
Combining this with the bar element, a formulation of a general 2D beam element can be
derived.
Example 7.5
The step shaft is fully clamped at the both ends and loaded in a transverse force P in the
middle as shown. Derive the global finite element equilibrium equation (Exam 2011).
P
2EI EI x
1 2 3
L L
Solution
Step 1: Discretise into 2 elements and 3 nodes, as shown in the figure.
Step 2: Write the elemental equilibrium equations for both elements.
For element 1:
(1)
𝑣1 𝑓1𝑌
12 6𝐿 −12 6𝐿 (1)
2𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃1 𝑚1
[ ]{ } =
𝐿3 −12 −6𝐿 12 −6𝐿 𝑣2 (1)
𝑓2𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃2 (1)
{𝑚2 }
For element 2:
(2)
𝑣2 𝑓2𝑌
12 6𝐿 −12 6𝐿 (2)
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃2 𝑚2
[ ]{ } =
𝐿3 −12 −6𝐿 12 −6𝐿 𝑣3 (2)
𝑓3𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃3 (2)
{𝑚3 }
Step 3: Expand the elemental equilibrium equations.
Element 1:
25
(1)
𝑓1𝑌
24 12𝐿 −24 12𝐿 0 0 𝑣1
(1)
12𝐿 8𝐿2 −12𝐿 4𝐿2 0 0 𝜃1 𝑚1
𝐸𝐼 −24 −12𝐿 24 −12𝐿 0 0 𝑣2 (1)
= 𝑓2𝑌
𝐿3 12𝐿 4𝐿2 −12𝐿 8𝐿2 0 0 𝜃2 (1)
0 0 0 0 0 0 𝑣3 𝑚2
[ 0 0 0 0 0 0] {𝜃3 } 0
{ 0 }
Element 2:
0
0 0 0 0 0 0 𝑣1
0
0 0 0 0 0 0 𝜃1 (2)
𝑓2𝑌
𝐸𝐼 0 0 12 6𝐿 −12 6𝐿 𝑣2
= 𝑚(2)
𝐿3 0 0 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃2 2
0 0 −12 −6𝐿 12 −6𝐿 𝑣3
(2)
𝑓3𝑌
[0 0 6𝐿 2𝐿2 −6𝐿 4𝐿2 ] {𝜃3 } (2)
{𝑚3 }
Step 4: Assemble the global equilibrium equation.
(1)
𝑓1𝑌
24 12𝐿 −24 12𝐿 𝑣1 (1)
𝜃1 𝑚1
12𝐿 8𝐿2 −12𝐿 4𝐿2 (1) (2)
𝐸𝐼 −24 −12𝐿 24 + 12 −12𝐿 + 6𝐿 −12 6𝐿 𝑣2 𝑓2𝑌 + 𝑓2𝑌
=
𝐿3 12𝐿 4𝐿2 −12𝐿 + 6𝐿 8𝐿2 + 4𝐿2 −6𝐿 2𝐿2 𝜃2 (1)
𝑚2 + 𝑚2
(2)
−12 −6𝐿 12 −6𝐿 𝑣3 (2)

[ 𝑓3𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 ] {𝜃3 }
(2)
{ 𝑚3 }
Step 5: Apply the boundary conditions.
24 12𝐿 −24 12𝐿 0 𝐹1

12𝐿 8𝐿2 −12𝐿 4𝐿2 0 𝑀1
𝐸𝐼 −24 −12𝐿 24 + 12 −12𝐿 + 6𝐿 −12 6𝐿 𝑣 2 = −𝑃
𝐿3 12𝐿 4𝐿2 −12𝐿 + 6𝐿 8𝐿2 + 4𝐿2 −6𝐿 2𝐿2 𝜃2 0
−12 −6𝐿 12 −6𝐿 0 𝐹3
[ 2 ]{ 0 } { 𝑀3 }
6𝐿 2𝐿2 −6𝐿 4𝐿
Step 6: Delete zero displacement equations to obtain the reduced equilibrium equation.
𝐸𝐼 36 −6𝐿 𝑣2 −𝑃
[ ] { } = { }
𝐿3 −6𝐿 12𝐿2 𝜃2 0
Step 7: Solve the reduced equation to obtain the unknown displacement vector. Reaction
forces and moments can be calculated by using the eliminated equations.
26
𝐿2
𝑣2 − 𝑃
∴ {𝜃 } = 66𝐸𝐼
2 𝐿3
{ 33𝐸𝐼 𝑃}
−
7.5. Modelling distributed loads on finite elements

In our treatment of finite elements, loading has only been accomplished on nodes. In some
cases, distributed loading along elements best describes the mechanical system. In this case,
these distributed loads can be abstracted to the nodes by the use of equivalent nodal forces.
Fig 7.10 Two statically equivalent bar elements
As an example, a uniformly distributed load q (N/mm, N/m, lb/in) can be converted to two
equivalent nodal forces of magnitude qL/2, as shown in Fig. 7.10. We verify this by
considering the work done by the load q.
𝐿
1
𝑊𝑞 = ∫ 𝑢𝑞 𝑑𝑥
0 2
1 1
= ∫ 𝑢(𝜉)𝑞(𝐿𝑑𝜉)
2 0
𝑞𝐿 1 𝑢𝑖
= ∫ [𝑁𝑖 (𝜉) 𝑁𝑗 (𝜉)] {𝑢 } 𝑑𝜉
2 0 𝑗
1 𝑢𝑖
𝑞𝐿
= (∫ [1 − 𝜉 𝜉]𝑑𝜉 ) {𝑢 }
2 0 𝑗
1 𝑞𝐿 𝑞𝐿 𝑢𝑖
= [ ]{ }
2 2 2 𝑢𝑗
1 𝑞𝐿/2
= [𝑢𝑖 𝑢𝑗 ] { }
2 𝑞𝐿/2
1 𝑞𝐿/2
That is, 𝑊𝑞 = 2 𝐮𝑇 𝐟𝑞 , where 𝐟𝑞 = { }. In general, there may be other nodal forces in
𝑞𝐿/2
addition to the distributed load q. By equating strain energy and work done in this element,
we come to the following equation.
1 𝑇 1 1
𝑈≡ 𝐮 𝐤𝐮 = 𝑊𝑞 ≡ 𝐮𝑇 𝐟 + 𝐮𝑇 𝐟𝑞
2 2 ︸ 2 ︸
other distributed
nodal load
force
This is in the form of 𝐤𝐮 = 𝐟 + 𝐟𝑞 , which provides a generalised nodal force vector below.
27
𝑓𝑖 𝑞𝐿/2 𝑓𝑖 + 𝑞𝐿/2
𝑓 + 𝑓𝑞 = { } + { }={ }
𝑓𝑗 𝑞𝐿/2 𝑓𝑗 + 𝑞𝐿/2
7.5.1. Multiple bar element system

For multiple bar elements as illustrated in Fig. 7.11, one can expand the force vector in each
element, and superimpose the expanded force vectors.
Fig 7.11 Distributed loading on two connected bar elements
(1)
𝑓1 + 𝑞𝐿/2 0
(2)
(𝑓 + 𝑓𝑞 ) + (𝑓 + 𝑓𝑞 ) = {𝑓 (1) + 𝑞𝐿/2} + {𝑓1 + 𝑞𝐿/2}
element 1 element 2 2 (1)
0 𝑓2 + 𝑞𝐿/2
𝑓1 + 𝑞𝐿/2 𝑞𝐿/2
= { 𝑓2 + 𝑞𝐿 } = { 𝑞𝐿 }
𝑓3 + 𝑞𝐿/2 𝑞𝐿/2
7.5.2. Equivalent nodal loading for beam elements
Fig 7.12 Distributed transverse loading on a beam element
In Fig. 7.12, the distributed transverse loading can be made equivalent to nodal shear forces
(qL/2) and moments (qL2/12). This can be verified by considering the work done by the
distributed load q.
Fig 7.13 Distributed loading on two beam elements
For two beam elements, the equivalent loading can be calculated as shown in Fig 7.13.
28
Example 7.6
A cantilever beam with distributed lateral load p as shown above. Find the deflection and
rotation at the right end, and the reaction force and moment at the left end.
Solution
Step 1: The work-equivalent nodal loads for the distributed load p are shown on the right end
of the figure below, where 𝑓 = 𝑝𝐿/2, 𝑚 = 𝑝𝐿2 /12
Step 2: Apply the FE equilibrium equation.
12 6𝐿 −12 6𝐿 𝑣1 𝐹1𝑌
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃1 𝑀
[ ] {𝑣 } = { 1 }
𝐿3 −12 −6𝐿 12 −6𝐿 2 𝐹2𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃2 𝑀2
Step 3: Apply boundary conditions 𝐹2𝑌 = −𝑓, 𝑀2 = 𝑚, 𝑣1 = 𝜃1 = 0, and eliminate the zero
displacement equations.
é 12 6L - 12 6 L ù ì 0 ü ì F1Y ü
ê
EI ê 6 L 4L 2
- 6L 2 L2 ú ïï 0 ïï ïï M1 ïï
úí ý = í ý
L3 ê - 12 - 6L 12 - 6 Lú ïv2 ï ï - f ï
ê ú
ë 6L 2 L2 - 6L 4 L2 û ïîq 2 ïþ ïî m ïþ
The reduced equation set is therefore:
𝐸𝐼 12 −6𝐿 𝑣2 −𝑓
[ ]{ } = { }
𝐿3 −6𝐿 4𝐿2 𝜃2 𝑚
Step 4: Solve for the displacements using the reduced equation set.
29
𝑝𝐿4
𝑣2 𝐿 −2𝐿2 𝑓 + 3𝐿𝑚 −
{𝜃 } = { }= 8𝐸𝐼
2 6𝐸𝐼 −3𝐿𝑓 + 6𝑚 𝑝𝐿3
−
{ 6𝐸𝐼 }
These nodal values are the same as from engineering beam theory. Note that the deflection
v(x) (for 0 < x< L) in the beam by the FEM is, however, different from that by the exact
solution. The exact solution by the engineering beam theory is a 4th order polynomial of x,
while the FE solution of v is only a 3rd order polynomial of x. If more elements are used,
such an approximate error can be reduced.
Step 5: Solve for the reactions forces and moments using the eliminated FE equations.
𝑝𝐿
3
𝐹 𝐿 −12 6𝐿 𝑣2 2
{ 1𝑌 } = [ 2 ] {𝜃 } = { }
𝑀1 𝐸𝐼 −6𝐿 2𝐿 2 5𝑝𝐿2
12
Example 7.7
Given P = 50 kN, k = 200 kN/m, L = 3 m, E = 210 GPa, I = 2×10–4 m4, find the deflections,
rotations and reaction forces and moments in this mixed element system.
Solution
Step 1: Connectivity analysis: this is a combined problem of beam and spring elements. The
system has four nodes, connected by two beam elements and one spring element. There are
7 DOFs, represented in the unknown displacement vector.
𝐮 = {𝑣1 , 𝜃1 , 𝑣2 , 𝜃2 , 𝑣3 , 𝜃3 , 𝑣4 }𝑇
Step 2: Elemental stiffness matrices for beam and spring elements.
For beam element 1,
12 6𝐿 −12 6𝐿 𝑣1 𝑓1𝑌
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃1 𝑚
[ ] {𝑣 } = { 1 }
𝐿3 −12 −6𝐿 12 −6𝐿 2 𝑓2𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃2 𝑚2
30
For beam element 2,
12 6𝐿 −12 6𝐿 𝑣2 𝑓2𝑌
𝐸𝐼 6𝐿 4𝐿2 −6𝐿 2𝐿2 𝜃2 𝑚
[ ] { } = { 2}
𝐿3 −12 −6𝐿 12 −6𝐿 𝑣3 𝑓3𝑌
6𝐿 2𝐿2 −6𝐿 4𝐿2 𝜃3 𝑚3
For the spring element, which is tied to displacements 𝑣3 and 𝑣4.
𝑘 −𝑘
𝐤𝑠 = [ ]
−𝑘 𝑘
𝑣3 𝑣4
Step 3: Expand the elemental equilibrium equations
For beam element 1:
12 6𝐿 −12 6𝐿 0 0 0 𝑣1 𝑓1
6𝐿 4𝐿2 −6𝐿 2𝐿2 0 0 0 𝜃1 𝑚1
𝐸𝐼 −12 −6𝐿 12 −6𝐿 0 0 0 𝑣2 𝑓2
6𝐿 2𝐿2 −6𝐿 4𝐿2 0 0 0 𝜃2 = 𝑚2
𝐿3
0 0 0 0 0 0 0 𝑣3 0
0 0 0 0 0 0 0 𝜃3 0
[ 0 0 0 0 0 0 0] {𝑣4 } { 0 }
For beam element 2:
0 0 0 0 0 0 0 𝑣1 0
0 0 0 0 0 0 0 𝜃1 0
𝑣 𝑓
𝐸𝐼 0 0 12 6𝐿 −12 6𝐿 0 2 2
0 0 6𝐿 4𝐿2 −6𝐿 2𝐿2 0 𝜃2 = 𝑚2
𝐿3
0 0 −12 −6𝐿 12 −6𝐿 0 𝑣3 𝑓3
0 0 6𝐿 2𝐿2 −6𝐿 4𝐿2 0 𝜃3 𝑚3
[0 0 0 0 0 0 0] {𝑣4 } { 0 }
For the spring element:
0 0 0 0 0 0 0 𝑣1 0
0 0 0 0 0 0 0 𝜃1 0
𝑣2 0
𝐸𝐼 0 0 0 0 0 0 0
0 0 0 0 0 0 0 𝜃2 = 0
𝐿3 𝑣3
0 0 0 0 𝑘′ 0 −𝑘 ′ 𝑓3
0 0 0 0 0 0 0 𝜃3 0
[0 0 0 0 −𝑘 ′ 0 𝑘 ′ ] {𝑣 }
4 {𝑓4 }
𝐿3
where 𝑘 ′ = 𝐸𝐼 𝑘. We do this to accommodate the elemental stiffness matrix of the spring
element into a single global stiffness matrix.
Step 4: Assemble the global FE equilibrium equation.
31
Beam element 1
Beam element 2
spring element
Step 5: Apply the boundary conditions 𝑣1 = 𝜃1 = 𝑣2 = 𝑣4 = 0, 𝑀2 = 𝑀3 = 0, 𝐹3𝑌 = −𝑃
Eliminating the zero displacement equations, a reduced set of FE equations is obtained.
𝐸𝐼 8𝐿
2
−6𝐿 2𝐿2 𝜃2 0
[−6𝐿 12 + 𝑘 ′ −6𝐿 ] {𝑣3} = {−𝑃}
𝐿3
2𝐿2 −6𝐿 4𝐿2 𝜃3 0
Step 6: Solve the reduced set of equations for the displacements (deflections and rotations).
𝜃2 3 −0.002492 rad
𝑃𝐿2
{𝑣3} = − {
′ ) 7𝐿
} = { −0.01744 m }
𝜃3 𝐸𝐼(12 + 7𝑘
9 −0.007475 rad
Step 7: Calculate the reaction moments and forces from the eliminated equations.
𝐹1𝑌 −69.78 kN
𝑀1 −69.78 kN⋅m
{ }={ }
𝐹2𝑌 116.2 kN
𝐹4𝑌 3.488 kN
A free body diagram for the beam portion of the system can be drawn as shown.
32
7.6. FE methods for vibrational analysis

Vibrational analysis is often performed on mechanical systems to evaluate their dynamic
behaviour and stability during vibrations.
7.6.1. Single degree-of-freedom systems

A system that can be described by a single co-ordinate is called a single degree-of-freedom
(DOF) system. For example, a system consisting of one mass, one spring, and one damper
vibrating in one direction is a typical single DOF system, as shown in Fig. 7.14. A single
mass vibrating in a single direction is a single DOF system.
Fig 7.14 A typical single DOF vibrating system (m-c-k system).
We are concerned with the motion of the mass m, i.e. x(t). The mass m is considered to be a
rigid body.
7.6.2. Multiple degree-of-freedom systems

A system that must be described by two or more independent co-ordinates is called a multi-
DOF system. A single mass can have a maximum of 6 DOFs, i.e. 3 DOFs in translation
(straight line motion) and 3 DOFs in rotation, as shown in Fig. 7.15.
z
θz
m y
θx θ
x y
Fig. 7.15 Translation and rotation of a mass (6 DOFs)
Examples of 2 DOF systems are shown in Fig 7.16 and 7.17. A single mass vibrating in
different directions constitutes a multi DOF system, as well as multiple masses vibrating in a
single (or multiple) direction(s).
y(t) m
x(t)
Fig 7.16 A single mass vibrating in x and y, constituting a 2 DOF system.
33
Fig 7.17 Two masses vibrating in a single direction (x), constituting a 2 DOF system.
In an n-DOF system, the number of co-ordinates required to describe the motion of the
masses is N. An N-DOF system has N number of natural frequencies (ωn). For each ωn,
there is a corresponding vibration mode (i.e. mode shape).
7.6.3. Equations of motion

In vibration analysis, the displacement of each mass (i.e. x(t)) is measured from its static
equilibrium position, as illustrated in Fig 7.16. Thus we can ignore the gravity force of each
mass and the initial deformation of the springs when we derive the vibration differential
equations using Newton’s second law.
Fig 7.17 Static equilibrium positions of 2 DOF system
Fig 7.18 A 2 DOF free vibration system and free body diagram
To derive the equations of motion, consider the free body diagram of the 2 DOF system
illustrated in Fig. 7.18. Newton’s second law (∑𝐹 = 𝑚𝑎 = 𝑚𝑥̈ ) can be applied to each mass.
For mass 1:
𝑚1 𝑥̈ 1 = −𝑘1 𝑥1 − 𝑘2 𝑥1 + 𝑘2 𝑥2
𝑚1 𝑥̈1 + (𝑘1 + 𝑘2 )𝑥1 − 𝑘2 𝑥2 = 0
34
For mass 2:
𝑚2 𝑥̈ 2 = −𝑘2 𝑥2 + 𝑘2 𝑥1
𝑚2 𝑥̈ 2 − 𝑘2 𝑥1 + 𝑘2 𝑥2 = 0
These two equations can be expressed in matrix form.
𝑚1 0 𝑥̈ 1 𝑘 + 𝑘2 −𝑘2 𝑥1 0
[ ]{ }+[ 1 ] {𝑥 } = { }
0 𝑚2 𝑥̈ 2 −𝑘2 𝑘2 2 0
or
[𝐌]{𝐱̈ } + [𝐊]{𝐱} = {0}

(7.3)
where:
𝑚1 0
[𝐌] = [ ] is the mass matrix;
0 𝑚2
𝑘 + 𝑘2 −𝑘2
[𝐊] = [ 1 ] is the stiffness matrix; and
−𝑘2 𝑘2
𝑥
{𝐱} = {𝑥1 } is the displacement vector
2
7.6.4. Characteristic free vibration (synchronous)

In the above example, each mass undergoes harmonic motion at the same frequency
(synchronous) and exhibits a certain vibration mode. Therefore, we assume the following
synchronous solutions.
𝑥1 (𝑡) = 𝐴1 sin 𝜔𝑡 𝑥̈ 1 (𝑡) = −𝜔2 𝐴1 sin 𝜔𝑡

𝑥2 (𝑡) = 𝐴2 sin 𝜔𝑡 𝑥̈ 2 (𝑡) = −𝜔2 𝐴2 sin 𝜔𝑡
(7.4)
Substituting the synchronous solutions in (7.4) into the equation of motion (Eq. 7.3), we
obtain the following equation set.
𝑚 0 𝐴1 𝑘 + 𝑘2 −𝑘2 𝐴1 0
−𝜔2 [ 1 ][ ]+ [ 1 ][ ] = [ ]
0 𝑚2 𝐴2 −𝑘2 𝑘2 𝐴2 0
𝑘1 + 𝑘2 − 𝜔2 𝑚1 −𝑘2 𝐴1 0
[ ][ ] = [ ]
−𝑘2 𝑘2 − 𝜔2 𝑚2 𝐴2 0
(7.5)
For a non-trivial solution (at least one of the A’s is non-zero) we have:
35
𝑘1 + 𝑘2 − 𝜔2 𝑚1 −𝑘2
| |=0
−𝑘2 𝑘2 − 𝜔2 𝑚2
(7.6)
which is the characteristic equation of the system.
7.6.5. Natural frequencies

The natural frequencies are the solutions to the characteristic equation. For simplicity, let’s
look at the special system where 𝑚1 = 𝑚2 = 𝑚 and 𝑘1 = 𝑘2 = 𝑘. The characteristic
equation becomes:
2
|2𝑘 − 𝜔 𝑚 −𝑘 | = 0
−𝑘 𝑘 − 𝜔2 𝑚
(2𝑘 − 𝜔2 𝑚)(𝑘 − 𝜔2 𝑚) − 𝑘 2 = 0
𝑚2 𝜔4 − 3𝑘𝑚𝜔2 + 𝑘 2 = 0
𝑘
Divide by k2, and put 𝜔0 = √𝑚 , where 𝜔0 is an assumed reference frequency.
𝜔4 𝜔2
− 3 +1=0
𝑘2 𝑘
𝑚2 𝑚
4
𝜔 𝜔2
− 3 +1=0
𝜔04 𝜔02
𝜔 2
This is a quadratic equation, in terms of (𝜔 ) . Therefore, using the quadratic equation:
0
𝜔 2 3 ± √9 − 4 3 ± √5 0.382
( ) = = ={
𝜔0 2 2 2.618
Two natural frequencies are obtained, 𝜔𝑛1 and 𝜔𝑛2 .
𝑘
1st natural circular frequency: 𝜔𝑛1 = 0.618𝜔0 = 0.618√ rad/s
𝑚
𝑘
2nd natural circular frequency: 𝜔𝑛2 = 1.618𝜔0 = 1.618√𝑚 rad/s
7.6.6. Mode shapes

For each 𝜔𝑛 there is a corresponding mode shape. For a multi DOF spring-mass system, a
mode shape is simply an amplitude ratio. Recall in Eq. 7.4, we provided two synchronous
solutions with different amplitudes 𝐴1 and 𝐴2 .
36
𝑥1 (𝑡) = 𝐴1 sin 𝜔𝑡
𝑥2 (𝑡) = 𝐴2 sin 𝜔𝑡
𝐴
Hence the mode shape is 𝐴1.
2
Expanding Eq. 7.5,
𝑘 + 𝑘2 − 𝜔2 𝑚1 −𝑘2 𝐴 0
[ 1 2 ] [ 1] = [ ]
−𝑘2 𝑘2 − 𝜔 𝑚2 𝐴2 0
(𝑘1 + 𝑘2 − 𝜔2 𝑚1 )𝐴1 − 𝑘2 𝐴2 = 0
−𝑘2 𝐴1 + (𝑘2 − 𝜔2 𝑚2 )𝐴2 = 0
Mode shape expressions can be derived from either equations.
𝐴1 𝑘2
=
𝐴2 𝑘1 + 𝑘2 − 𝜔 2 𝑚1
𝐴1 𝑘2 − 𝜔2 𝑚2
=
𝐴2 𝑘2
(7.7)
Both ratios will be equivalent. For our special case in section 7.6.5, this can be confirmed by
𝑘
substituting the natural circular frequency 𝜔𝑛1 = 0.618√𝑚 into both equations in (7.7).
𝐴1 𝑘2 𝑘 𝑘 1
( ) = = = = = 0.618
𝐴2 1 𝑘1 + 𝑘2 − 𝜔 𝑚1 2𝑘 − 0.382 𝑘 𝑚 1.618𝑘 1.618
2
𝑚
𝑘
𝐴1 𝑘2 − 𝜔2 𝑚2 𝑘 − 𝜔2 𝑚 𝑘 − 0.382 𝑚 𝑚 0.618𝑘
( ) = = = = = 0.618
𝐴2 1 𝑘2 𝑘 𝑘 𝑘
𝑘
Similarly, for 𝜔𝑛2 = 1.618√𝑚
𝐴1 𝑘2 𝑘 𝑘 1
( ) = = = =− = −1.618
𝐴2 2 𝑘1 + 𝑘2 − 𝜔 𝑚1 2𝑘 − 2.618 𝑘 𝑚 −0.618𝑘
2 0.618
𝑚
𝑘
𝐴1 𝑘2 − 𝜔2 𝑚2 𝑘 − 𝜔2 𝑚 𝑘 − 2.618 𝑚 𝑚 −1.618𝑘
( ) = = = = = −1.618
𝐴2 2 𝑘2 𝑘 𝑘 𝑘
For convenience, the amplitude A2 is taken unit (A2=1). As a result, the mode shapes can also
be expressed in terms of vector form.
37
0.618
1st mode shape: { }
1
−1.618
2nd mode shape: { }
1
These are the actually the eigenvectors corresponding to the eigenvalues (𝜔𝑛1 , 𝜔𝑛2 ) found
from the characteristic equation.
𝑘 𝐴
For 𝜔𝑛1 = 0.618√ and ( 1) = 0.618, the first order vibration mode (if A2=1, then
𝑚 𝐴2 1
A1=0.618) can be drawn as shown Fig. 7.19.
Fig 7.19 The first mode shape
𝑘 𝐴
Similarly, for 𝜔𝑛2 = 1.618√𝑚 and (𝐴1) = −1.618, the second order vibration mode (if
2 2
A2=1, then A1=-1.618) can be drawn as shown in Fig. 7.20. Notably, when vibrating at the 2 nd
order natural frequency, the two masses are out of phase, i.e. they are vibrating in different
directions.
Fig 7.20 The second mode shape
38
7.7. Two-dimensional finite elements

So far, we have only dealt with problems using one-dimensional (1D) finite elements (spring,
bar, and beam elements). In solid mechanics, two-dimensional (2D) elements are plane
elements (which can model plane stress and plane strain problems) with two displacements
(u, v) for both dimensions.
7.7.1. General formula for the stiffness matrix in 2D

Recall the general formulation of a 1D bar element in terms of shape functions from section
7.3.2. The displacement u within the bar element could be expressed in terms of the nodal
displacements 𝐝.
𝑢
𝑢 = [𝑁𝑖 𝑁𝑗 ] {𝑢𝑖 } = 𝐍𝐝
𝑗
This can be extended to two dimensions, where the displacements (u, v) in a plane element
are interpolated from nodal displacements (ui, vi) using the shape functions Ni.
𝑢 = 𝑢(𝑥, 𝑦) = 𝑁1 (𝑥, 𝑦)𝑢1 + 𝑁2 (𝑥, 𝑦)𝑢2 + ⋯
𝑣 = 𝑣(𝑥, 𝑦) = 𝑁1 (𝑥, 𝑦)𝑣1 + 𝑁2 (𝑥, 𝑦)𝑣2 + ⋯
This can be expressed in matrix form.
𝑢1
𝑣1
𝑢 𝑁 0 𝑁2 0 ⋯
𝐮={ }=[ 1 ] 𝑢2 = 𝐍𝐝
𝑣 0 𝑁1 0 𝑁2 ⋯
𝑣2
{⋮}
(7.8)
where 𝐍 is the shape function matrix, 𝐮 the displacement vector inside each element and
𝐝 is the nodal displacement vector of an element. Here we assume that 𝑢 depends on the
nodal values of 𝑢𝑖 only, and 𝑣 on nodal values of 𝑣𝑖 only.
From the strain-displacement relation, the strain vector can be derived as:
𝜀 = 𝐃𝐮 = 𝐃(𝐍𝐝) = (𝐃𝐍)𝐝
or
𝜀 = 𝐁𝐝
As in the 1D case, 𝐁 is the strain-displacement matrix, which is formed from the derivative
matrix 𝐃 and the shape function matrix 𝐍. Once again, the strain energy-work formulation
can be used to derive the expression for the elemental stiffness matrix.
39
Consider the strain energy stored in a 2D element and 𝐄 is a symmetric matrix (𝐄 𝑇 = 𝐄):
1 1
𝑈= ∫[𝜎 𝑇 𝜀] 𝑑𝑉 = ∫(𝜎𝑥𝑥 𝜀𝑥𝑥 + 𝜎𝑦𝑦 𝜀𝑦𝑦 + 𝜎𝑥𝑦 𝜀𝑥𝑦 ) 𝑑𝑉
2 2
𝑉 𝑉
1 1 1
= ∫(𝚬𝛆)𝑇 𝛆 𝑑𝑉 = ∫ 𝛆𝑇 𝚬𝛆 𝑑𝑉 = ∫(𝐁𝐝)𝑇 𝚬(𝐁𝐝) 𝑑𝑉
2 2 2
𝑉 𝑉 𝑉
Since nodal displacement vector 𝐝 is independent on the elemental coordinate 𝑑𝑉.
1 𝑇 1
𝑈= 𝐝 (∫ 𝐁 𝑇 𝚬𝚩 𝑑𝑉) 𝐝 = 𝐝𝑇 𝐤𝐝
2 2
𝑉
The general formula for the elemental stiffness matrix is:
𝐤 = ∫ 𝐁 𝑇 𝚬𝚩 𝑑𝑉
𝑉
Some remarks can be made from this formulation:
Note that unlike 1D cases, 𝐄 is a matrix that is given by the stress-strain relation
(generalised Hooke’s law);
The elemental stiffness matrix 𝐤 is symmetric as 𝐄 is symmetric;
For any given material property 𝐄, the behaviour of 𝐤 depends on the 𝐁 matrix only,
which in turn depends on the shape functions. Thus, the quality of finite elements in
representing the behaviour of a structure is entirely determined by the choice of shape
functions.
7.7.2. 2D constant strain triangle (CST)
Fig 7.21 A linear triangular element or CST
The constant strain triangle (CST), also known as T3 or linear triangular element is the
simplest 2D element, illustrated in Fig. 7.21. For this element, we have three nodes at the
vertices of the triangle, which are numbered around the element in the counter-clockwise
direction. Each node has two degrees of freedom (can move in both x and y directions). The
displacements u and v are assumed to be linear functions within the element, so that:
40
𝑢 = 𝑢(𝑥, 𝑦) = 𝑏1 + 𝑏2 𝑥 + 𝑏3 𝑦
𝑣 = 𝑣(𝑥, 𝑦) = 𝑏4 + 𝑏5 𝑥 + 𝑏6 𝑦
(7.9)
where 𝑏𝑖 (𝑖 = 1, 2, . . . , 6) are constants. In terms of these constants, the strains can be

determined from the strain-displacement relationship.
𝜕𝑢 𝜕
𝜀𝑥𝑥 = = (𝑏 + 𝑏2 𝑥 + 𝑏3 𝑦) = 𝑏2
𝜕𝑥 𝜕𝑥 1
𝜕𝑣 𝜕
𝜀𝑦𝑦 = = (𝑏 + 𝑏5 𝑥 + 𝑏6 𝑦) = 𝑏6
𝜕𝑦 𝜕𝑦 4
𝜕𝑢 𝜕𝑣 𝜕 𝜕
2𝜀𝑥𝑦 = ( + )= (𝑏4 + 𝑏5 𝑥 + 𝑏6 𝑦) + (𝑏 + 𝑏2 𝑥 + 𝑏3 𝑦) = 𝑏3 + 𝑏5
𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 1
Therefore, the strains are constant throughout the element, which gives rise to the name
constant strain triangle. The displacement functions should satisfy the following six
equations, based on the nodal displacements.
𝑢1 = 𝑢(𝑥 = 𝑥1, 𝑦 = 𝑦1 ) = 𝑏1 + 𝑏2 𝑥1 + 𝑏3 𝑦1
𝑣1 = 𝑣(𝑥 = 𝑥1, 𝑦 = 𝑦1 ) = 𝑏4 + 𝑏5 𝑥1 + 𝑏6 𝑦1
𝑢2 = 𝑢(𝑥 = 𝑥2 , 𝑦 = 𝑦2 ) = 𝑏1 + 𝑏2 𝑥2 + 𝑏3 𝑦2
𝑣2 = 𝑣(𝑥 = 𝑥2, 𝑦 = 𝑦2 ) = 𝑏4 + 𝑏5 𝑥2 + 𝑏6 𝑦2
𝑢3 = 𝑢(𝑥 = 𝑥3 , 𝑦 = 𝑦3 ) = 𝑏1 + 𝑏2 𝑥3 + 𝑏3 𝑦3
𝑣3 = 𝑣(𝑥 = 𝑥3, 𝑦 = 𝑦3 ) = 𝑏4 + 𝑏5 𝑥3 + 𝑏6 𝑦3
Solving these six equations, we can find the coefficients 𝑏1 , 𝑏2 , . . . , 𝑏6 in terms of nodal
displacements and coordinates. Substituting these six coefficients into displacement equations
and rearranging the terms, we obtain a linear distribution of displacements within the
element.
𝐮 = 𝐍𝐝
3
𝑢 = ∑ 𝑁𝑖 (𝑥, 𝑦)𝑢𝑖 = 𝑁1 𝑢1 + 𝑁2 𝑢2 + 𝑁3 𝑢3
𝑖=1
3
𝑣 = ∑ 𝑁𝑖 (𝑥, 𝑦)𝑣𝑖 = 𝑁1 𝑣1 + 𝑁2 𝑣2 + 𝑁3 𝑣3
𝑖=1
or
41
𝑢1
𝑣1
𝑢 𝑁 0 𝑁2 0 𝑁3 0 𝑢 2
{ }=[ 1 ]
𝑣 0 𝑁1 0 𝑁2 0 𝑁3 𝑣2
𝑢3
{ 𝑣3 }
The shape functions in 𝐍 are linear functions of x and y.
1
𝑁1 (𝑥, 𝑦) = [(𝑥 𝑦 − 𝑥3 𝑦2 ) + (𝑦2 − 𝑦3 )𝑥 + (𝑥3 − 𝑥2 )𝑦]
2𝐴 2 3
1
𝑁2 (𝑥, 𝑦) = [(𝑥 𝑦 − 𝑥1 𝑦3) + (𝑦3 − 𝑦1 )𝑥 + (𝑥1 − 𝑥3 )𝑦]
2𝐴 3 1
1
𝑁3 (𝑥, 𝑦) = [(𝑥 𝑦 − 𝑥2 𝑦1) + (𝑦1 − 𝑦2 )𝑥 + (𝑥2 − 𝑥1 )𝑦]
2𝐴 1 2
(7.10)
1 𝑥1 𝑦1
1
where 𝐴 = 2 det [1 𝑥2 𝑦2 ] is the area of the triangle.
1 𝑥3 𝑦3
The strain matrix can be determined using the strain-displacement 𝜀 = 𝐁𝐝, but with explicit
shape functions for this CST.
𝜕𝑢 𝜕 1
𝜀𝑥𝑥 = = (𝑁1 𝑢1 + 𝑁2 𝑢2 + 𝑁3 𝑢3 ) = [(𝑦 −𝑦 )𝑢 + (𝑦3 −𝑦1 )𝑢2 + (𝑦1 −𝑦2 )𝑢3 ]
𝜕𝑥 𝜕𝑥 2𝐴 2 3 1
𝜕𝑣 𝜕 1
𝜀𝑦𝑦 = = (𝑁1 𝑣1 + 𝑁2 𝑣2 + 𝑁3 𝑣3 ) = [(𝑥 −𝑥 )𝑣 + (𝑥1 −𝑥3)𝑣2 + (𝑥2−𝑥1 )𝑣3]
𝜕𝑦 𝜕𝑦 2𝐴 3 2 1
𝜕𝑣 𝜕𝑢 𝜕 𝜕
2𝜀𝑥𝑦 = + = (𝑁1 𝑣1 + 𝑁2 𝑣2 + 𝑁3 𝑣3 ) + (𝑁 𝑢 + 𝑁2 𝑢2 + 𝑁3 𝑢3 )
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 1 1
1
= [(𝑦 −𝑦 )𝑣 + (𝑦3 −𝑦1 )𝑣2 + (𝑦1 −𝑦2)𝑣3]
2𝐴 2 3 1
1
+ [(𝑥 −𝑥 )𝑢 + (𝑥1−𝑥3 )𝑢2 + (𝑥2 −𝑥1)𝑢3 ]
2𝐴 3 2 1
or in matrix form:
𝑢1
𝑣1
𝜀𝑥𝑥
1 𝑦23 0 𝑦31 0 𝑦12 0
𝑢2
{𝜀𝑦𝑦 } = [ 0 𝑥32 0 𝑥13 0 𝑥21 ] 𝑣
𝜀𝑥𝑦 2𝐴 𝑥 𝑦23 𝑥13 𝑦31 𝑥21 𝑦12 𝑢2
32
3
{𝑣3 }
where xij = xi – xj and yij = yi – yj (𝑖, 𝑗 = 1, 2, 3). Again, we see constant strains within the
element (because variables x and y disappeared after differentiation). By Hooke’s law, we see
that stresses obtained in the CST element should also be constant.
42
The general formula for stiffness matrix can be applied for the CST element.
𝐤 = ∫ 𝐁 𝑇 𝚬𝚩 𝑑𝑉 = 𝑡𝐴(𝐁 𝑇 𝚬𝚩)
𝑉
𝑡 is the thickness of the element (for a plane stress example). 𝐤 for a CST is a 6 by 6
symmetric matrix. The matrix multiplication can be carried out by a computer program.
Due to the constant strain nature of this CST element, it uses are often limited as per
following recommendations.
Areas where the strain gradient is small.

Mesh transition areas (fine mesh to coarse mesh) to accommodate for different mesh
topologies.
Avoid using CST in stress concentration or other crucial areas in the structure, such as
edges of holes and corners.
Recommended only for quick and preliminary FE analysis of 2D problems.
7.7.3. 2D linear strain triangle (LST)
Fig 7.22 A quadratic triangular element or LST
The linear strain triangle (LST) is also known as T6 or quadratic triangular element,
illustrated in Fig. 7.22. There are six nodes on this element: three corner nodes and three mid-
side nodes. Each node has two degrees of freedom. The displacements (u, v) are assumed to
be quadratic functions of the co-ordinates (x, y) within the element.
𝑢 = 𝑏1 + 𝑏2 𝑥 + 𝑏3 𝑦 + 𝑏4 𝑥 2 + 𝑏5 𝑥𝑦 + 𝑏6 𝑦 2
𝑣 = 𝑏7 + 𝑏8 𝑥 + 𝑏9 𝑦 + 𝑏10 𝑥 2 + 𝑏11 𝑥𝑦 + 𝑏12 𝑦 2
(7.11)
where 𝑏𝑖 (𝑖 = 1,2, … ,12) are the constants. In terms of these constants, the strains can be
determined.
𝜕𝑢 𝜕
𝜀𝑥𝑥 = = (𝑏1 + 𝑏2 𝑥 + 𝑏3 𝑦 + 𝑏4 𝑥 2 + 𝑏5 𝑥𝑦 + 𝑏6 𝑦 2 ) = 𝑏2 + 2𝑏4 𝑥 + 𝑏5 𝑦
𝜕𝑥 𝜕𝑥
43
𝜕𝑣 𝜕
𝜀𝑦𝑦 = = (𝑏 + 𝑏8 𝑥 + 𝑏9 𝑦 + 𝑏10 𝑥 2 + 𝑏11 𝑥𝑦 + 𝑏12 𝑦 2 ) = 𝑏9 + 𝑏11 𝑥 + 2𝑏12 𝑦
𝜕𝑦 𝜕𝑦 7
𝜕𝑢 𝜕𝑣
2𝜀𝑥𝑦 =( + )
𝜕𝑦 𝜕𝑥
𝜕
= (𝑏 + 𝑏2 𝑥 + 𝑏3 𝑦 + 𝑏4 𝑥 2 + 𝑏5 𝑥𝑦 + 𝑏6 𝑦 2 )
𝜕𝑦 1
𝜕
+ (𝑏 + 𝑏8 𝑥 + 𝑏9 𝑦 + 𝑏10 𝑥 2 + 𝑏11 𝑥𝑦 + 𝑏12 𝑦 2)
𝜕𝑥 7
= (𝑏3 + 𝑏8 ) + (𝑏5 + 2𝑏10 )𝑥 + (2𝑏6 + 2𝑏11 )𝑦
All these strains are the linear functions of x and y. Thus, we have the linear strain
triangle, which may provide better accuracy than the CST.
In a natural coordinate system (𝜉, 𝜂), the six shape functions for the LST element can be
defined.
𝑁1 = 𝜉(2𝜉 − 1)
𝑁2 = 𝜂(2𝜂 − 1)
(1
𝑁3 = − 𝜉 − 𝜂)[2(1 − 𝜉 − 𝜂) − 1]
𝑁4 = 4𝜉𝜂
𝑁5 = 4𝜂(1 − 𝜉 − 𝜂)
𝑁6 = 4(1 − 𝜉 − 𝜂)𝜉
(7.12)
Using these shape functions, we obtain a quadratic distribution of displacements within the
elements. The displacements are found using the following equations.
6 6
𝑢 = ∑ 𝑁𝑖 𝑢𝑖 𝑣 = ∑ 𝑁𝑖 𝑣𝑖
𝑖=1 𝑖=1
For given nodal displacements (𝑢𝑖 , 𝑣𝑖 ) of the element, one can know the displacement at any
point within the element. In other words, we use shape functions and nodal displacements to
express displacement field (function) within the entire element.
The elemental stiffness matrix is still given by 𝐤 = ∫𝑉 𝐁 𝑇 𝚬𝚩 𝑑𝑉, however given the quadratic
nature of 𝐍, 𝐁 𝑇 𝚬𝚩 is also quadratic with respect to x and y. In general, this integral has to be
computed numerically.
44
7.7.4. 2D bilinear quadrilateral (Q4)
Fig 7.23 A bi-linear quadrilateral or Q4 element
The bilinear quadrilateral element is also known as a Q4 element, as there are four nodes on
this element: one at each of the corners of the quadrilateral shape (Fig. 7.23). Each node has
two degrees of freedom. In the natural co-ordinate system (𝜉, 𝜂), four shape functions can be
defined for the Q4 element.
1
𝑁1 = (1 − 𝜉)(1 − 𝜂)
4
1
𝑁2 = (1 + 𝜉)(1 − 𝜂)
4
1
𝑁3 = (1 + 𝜉)(1 + 𝜂)
4
1
𝑁4 = (1 − 𝜉)(1 + 𝜂)
4
(7.13)
These shape functions can be used to map the Cartesian co-ordinates to the natural co-
ordinate system (Fig. 7.23) using the relationships in (7.14).
4 4
𝑥 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑥𝑖 𝑦 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑦𝑖

𝑖=1 𝑖=1
(7.14)
Displacements within the element can also be mapped using the same shape functions. These
shape functions result in a bilinear distribution of displacement within the element.
Elements using the same shape functions for both co-ordinate and displacement mapping are
known as isoparametric elements.
4 4
𝑢 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑢𝑖 𝑣 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑣𝑖

𝑖=1 𝑖=1
(7.15)
45
From the displacement functions in (7.15), the strain can be determined.
6
𝜕𝑢 𝜕
𝜀𝑥𝑥 = = (∑ 𝑁𝑖 (𝜉, 𝜂)𝑢𝑖 )
𝜕𝑥 𝜕𝑥
𝑖=1
6
𝜕𝑣 𝜕
𝜀𝑦𝑦 = = (∑ 𝑁𝑖 (𝜉, 𝜂)𝑣𝑖 )
𝜕𝑦 𝜕𝑦
𝑖=1
6 6
𝜕𝑢 𝜕𝑣 𝜕 𝜕
2𝜀𝑥𝑦 =( + )= (∑ 𝑁𝑖 (𝜉, 𝜂) 𝑣𝑖 ) + (∑ 𝑁𝑖 (𝜉, 𝜂) 𝑢𝑖 )
𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦
𝑖=1 𝑖=1
In this form, the displacements are expressed as a function of the natural co-ordinates (𝜉, 𝜂),
rather than the global Cartesian co-ordinates (x, y). To calculate the strain, the chain rule can
be applied, as shown in (7.16). Note that u and v have been expressed in terms of 𝜉 and 𝜂.
𝜕𝑢 𝜕𝑢 𝜕𝑥 𝜕𝑢 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑢
+
𝜕𝜉 𝜕𝑥 𝜕𝜉 𝜕𝑦 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝑥
= =
𝜕𝑢 𝜕𝑢 𝜕𝑥 𝜕𝑢 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑢
+
{𝜕𝜂} {𝜕𝑥 𝜕𝜂 𝜕𝑦 𝜕𝜂 } ⏟[𝜕𝜂 𝜕𝜂] {𝜕𝑦}
[𝐉]
𝜕𝑢 𝜕𝑥 𝜕𝑦 −1 𝜕𝑢 𝜕𝑢
𝜕𝑥 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝐽∗ ∗
𝐽12 𝜕𝜉
𝜕𝑢 = = [ 11
∗ ∗ ]
𝜕𝑥 𝜕𝑦 𝜕𝑢 𝐽21 𝐽22 𝜕𝑢
{𝜕𝑦} [⏟𝜕𝜂 𝜕𝜂] {𝜕𝜂} {𝜕𝜂}
[𝐉]−1
(7.16)
where [𝐽] is called the Jacobian matrix. The coefficients in [𝐉] can be obtained as:
4
𝜕𝑥 𝜕𝑁𝑖
∵ 𝑥 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑥𝑖 → = ∑( 𝑥)
𝜕𝜉 𝜕𝜉 𝑖
𝑖=1
4
𝜕𝑥 𝜕𝑁𝑖
∵ 𝑥 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑥𝑖 → = ∑( 𝑥)
𝜕𝜂 𝜕𝜂 𝑖
𝑖=1
4
𝜕𝑦 𝜕𝑁𝑖
∵ 𝑦 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑦𝑖 → = ∑( 𝑦)
𝑖=1
4
𝜕𝑦 𝜕𝑁𝑖
∵ 𝑦 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑦𝑖 → = ∑( 𝑦)
𝜕𝜂 𝜕𝜂 𝑖
𝑖=1
Thus strains can be calculated using the inverse Jacobian matrix [𝐉]−1 as shown in Eq. 7.16.
∗ ∗
For example, 𝜀𝑥 is shown below where 𝐽11 and 𝐽12 are the coefficients in the first row of
−1
[𝐽] .
46
𝜕𝑢 ∗
𝜕𝑢 ∗
𝜕𝑢
𝜀𝑥 = = 𝐽11 + 𝐽12
𝜕𝑥 𝜕𝜉 𝜕𝜂
𝜕𝑢 𝜕𝑁𝑖 𝜕𝑢 𝜕𝑁𝑖
= ∑( 𝑢 ), = ∑( 𝑢 )
𝜕𝜉 𝜕𝜉 𝑖 𝜕𝜂 𝜕𝜂 𝑖
The general formula for the stiffness matrix can be applied for the Q4 element, but in terms
of the natural co-ordinates used by the shape functions.
1 1
𝐤 = ∫ 𝐁 𝐄𝐁 𝑑𝑉 = ∫ ∫ 𝐁 𝑇 𝐄𝐁𝑡|𝐉| 𝑑𝜉𝑑𝜂
𝑇
−1 −1
𝑉
|𝐉| is the determinant of the Jacobian matrix. To evaluate this integral, numerical integration
techniques, such as Gaussian quadrature are commonly used. Using the Gaussian
quadrature, one can have:
𝑛 𝑚
1 1
∫ ∫ 𝜙(𝜉, 𝜂) 𝑑𝜉𝑑𝜂 ≈ ∑ ∑ 𝑊𝑖 𝑊𝑗 𝜙(𝜉𝑖 , 𝜂𝑗 )
−1 −1 𝑖=1 𝑗=1
where 𝜙 is the function to be integrated, 𝑊𝑖 , 𝑊𝑗 are the weights of the quadrature rule, and
𝜉𝑖 , 𝜂𝑗 are the quadrature points (or Gaussian points). For this integration, we will select a
2×2 Gaussian integration order, which is illustrated in Fig 7.24.
Fig 7.24 2×2 Gaussian quadrature points
Gaussian points and weights for various integration orders are shown in the table below.
m×n Gaussian points (η, ξ) Weights (W)
1×1 0 2
2×2 1 Include this in the cheat sheet
±1/√3 = ±0.5774
0 0.8888
3×3
±√3/5 = ±0.7746 0.5555
Stress and strain results at the Gaussian points should be the most accurate in the entire
element.
47
7.7.5. 2D quadratic quadrilateral (Q8)

The quadratic quadrilateral element is also known as a Q8 element, as there are eight nodes
on this element: one at each of the corners, and one at the midpoint of each side (Fig. 7.25).
Each node has two degrees of freedom. In the natural co-ordinate system (𝜉, 𝜂), eight shape
functions can be defined for the Q8 element.
Fig 7.25 Co-ordinate transformation for 2D quadratic quadrilateral (Q8) element
1 1
𝑁1 = − (1 − 𝜉)(1 − 𝜂)(1 + 𝜉 + 𝜂) 𝑁5 = (1 − 𝜉 2 )(1 − 𝜂)
4 2
1 1
𝑁2 = − (1 + 𝜉)(1 − 𝜂)(1 − 𝜉 + 𝜂) 𝑁6 = (1 + 𝜉)(1 − 𝜂2 )
4 2
1 1
𝑁3 = − (1 + 𝜉)(1 + 𝜂)(1 − 𝜉 − 𝜂) 𝑁7 = (1 − 𝜉 2 )(1 + 𝜂)
4 2
1 1
𝑁4 = − (1 − 𝜉)(1 + 𝜂)(1 + 𝜉 − 𝜂) 𝑁8 = (1 − 𝜉)(1 − 𝜂2 )
4 2
(7.17)
The Q8 element is also isoparametric, and hence both co-ordinate mapping and displacement
fields utilise the same shape functions, as shown below. The shape functions in (7.17) result
in a quadratic distribution of displacement within the element.
8 8

𝑖=1 𝑖=1
8 8

𝑖=1 𝑖=1
Similar to the Q4 element, the stiffness matrix is expressed in terms of the natural co-ordinate
system, and can be calculated numerically using Gaussian quadrature (2×2).
1 1
𝐤 = ∫ 𝐁 𝑇 𝐄𝐁 𝑑𝑉 = ∫ ∫ 𝐁 𝑇 𝐄𝐁𝑡|𝐉| 𝑑𝜉𝑑𝜂
−1 −1
𝑉
48
𝑛 𝑚
1 1
∫ ∫ 𝜙(𝜉, 𝜂) 𝑑𝜉𝑑𝜂 ≈ ∑ ∑ 𝑊𝑖 𝑊𝑗 𝜙(𝜉𝑖 , 𝜂𝑗 )
−1 −1 𝑖=1 𝑗=1
7.7.6. Remarks on 2D elements

In sections 7.7.2 to 7.7.5, we have described the formulation of 2D triangular and rectangular
elements. These elements have either linear (T3 and Q4) or quadratic (T6 and Q8)
displacement distributions as a result of their shape functions. Hence, T3 and Q4 elements are
classified as linear elements and T6 and Q8 elements are classified as quadratic elements.
Linear T3 and Q4 elements are usually applied together in a mesh without mid-side nodes.
Quadratic T6 and Q8 elements are preferred for stress analysis because of their higher
accuracy and their flexibility in modelling complex geometries such as curved boundaries.
Example 7.8
Sketch the displacement distribution of the following elements (CST, LST, Q4 and Q8) and
explain the difference (final exam 2010, 2011, 2012).
Solution
3-node element (T3, CST): linear distribution of displacement and constant in the strain.
4-node element (Q4): bilinear distribution of displacements and bilateral linear in the strain
3 node CST element 4 node bilinear element
6-node element (T6, LST): quadratic distribution of displacement and linear distribution in
the strain
8-node element (Q8): quadratic distribution of displacements and linear distribution in the
strain
49
For LST draw curves
6-node LST element 8-node quadratic element
Example 7.9
For the 8-node element in Example 7.8, if 2×2 Gaussian points are considered, go on to
determine the displacements of all the Gaussian points if the nodal displacement vectors are
given as below:
𝑢 1 0 −1 −1 0 1 2 1
{ }=[ ]
𝑣 2 1 −1 −1 1 0 −1 −2
Solution
The displacement functions as per shape functions Ni and nodal displacements ui and vi can
be expressed as
4 4
𝑢 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑢𝑖 , 𝑣 = ∑ 𝑁𝑖 (𝜉, 𝜂)𝑣𝑖

𝑖=1 𝑖=1
Shape functions for 8-node elements can be written in terms of natural coordinates as shown
in (7.17):
1 1
𝑁1 = − (1 − 𝜉)(1 − 𝜂)(1 + 𝜉 + 𝜂) 𝑁5 = (1 − 𝜉 2 )(1 − 𝜂)
4 2
1 1
𝑁2 = − (1 + 𝜉)(1 − 𝜂)(1 − 𝜉 + 𝜂) 𝑁6 = (1 + 𝜉)(1 − 𝜂2 )
4 2
1 1
𝑁3 = − (1 + 𝜉)(1 + 𝜂)(1 − 𝜉 − 𝜂) 𝑁7 = (1 − 𝜉 2 )(1 + 𝜂)
4 2
1 1
𝑁4 = − (1 − 𝜉)(1 + 𝜂)(1 + 𝜉 − 𝜂) 𝑁8 = (1 − 𝜉)(1 − 𝜂2 )
4 2
Gaussian points can be defined in natural coordinates as (refer to the table of Gaussian points
and weights):
Gaussian point #1: 𝜉 = 0.5774, 𝜂 = 0.5774;
Gaussian point #2: 𝜉 = −0.5774, 𝜂 = 0.5774;
Gaussian point #3: 𝜉 = −0.5774, 𝜂 = 0.5774;
50
Gaussian point #4: 𝜉 = 0.5774, 𝜂 = −0.5774;
For {𝑢} = {1 0 −1 −1 0 1 2 1}𝑇 , 𝑢 at Gaussian point #1 can be calculated:

8
𝑢(𝜉, 𝜂) = 𝑢(0.5774, 0.5774) = ∑ 𝑁𝑖 𝑢𝑖

𝑖=1
= 𝑁1 𝑢1 + 𝑁2 𝑢2 + 𝑁3 𝑢3 + 𝑁4 𝑢4 + 𝑁5 𝑢5 + 𝑁6 𝑢6 + 𝑁7 𝑢7 + 𝑁8 𝑢8
= −0.25(1 − 𝜉)(1 − 𝜂)(1 + 𝜉 + 𝜂)𝑢1 − 0.25(1 + 𝜉)(1 − 𝜂)(1 − 𝜉 + 𝜂)𝑢2
− 0.25(1 + 𝜉)(1 + 𝜂)(1 − 𝜉 − 𝜂)𝑢3 − 0.25(1 − 𝜉)(1 + 𝜂)(1 + 𝜉 − 𝜂)𝑢4
+ 0.5(1 − 𝜉 2 )(1 − 𝜂)𝑢5 + 0.5(1 + 𝜉)(1 − 𝜂2 )𝑢6 + 0.5(1 − 𝜉 2 )(1 + 𝜂)𝑢7
+ 0.5(1 − 𝜉)(1 − 𝜂2 )𝑢8
= −0.25×(1 − 0.5774)(1 − 0.5774)(1 + 0.5774 + 0.5774)×1
− 0.25×(1 + 0.5774)(1 − 0.5774)(1 − 0.5774 + 0.5774)×0
− 0.25×(1 + 0.5774)(1 + 0.5774)(1 − 0.5774 − 0.5774)×(−1)
− 0.25×(1 − 0.5774)(1 + 0.5774)(1 + 0.5774 − 0.5774)×(−1)
+ 0.5×(1 − 0.57742 )(1 − 0.5774)×(0)
+ 0.5×(1 + 0.5774)(1 − 0.57742 )×(1)
+ 0.5×(1 − 0.57742 )(1 + 0.5774)×(2)
+ 0.5×(1 − 0.5774)(1 − 0.57742 )×(1)
∴ 𝑢(0.5774, 0.5774) = 1.6923
Similarly, for
{𝑣} = {2 1 −1 −1 1 0 −1 −2}𝑇 :
8
𝑣(𝜉, 𝜂) = 𝑣(0.5774, 0.5774) = ∑ 𝑁𝑖 𝑣𝑖

𝑖=1
= 𝑁1 𝑣1 + 𝑁2 𝑣2 + 𝑁3 𝑣3 + 𝑁4 𝑣4 + 𝑁5 𝑣5 + 𝑁6 𝑣6 + 𝑁7 𝑣7 + 𝑁8 𝑣8
= −0.25×(1 − 0.5774)(1 − 0.5774)(1 + 0.5774 + 0.5774)×2
− 0.25×(1 + 0.5774)(1 − 0.5774)(1 − 0.5774 + 0.5774)×1
− 0.25×(1 + 0.5774)(1 + 0.5774)(1 − 0.5774 − 0.5774)×(−1)
− 0.25×(1 − 0.5774)(1 + 0.5774)(1 + 0.5774 − 0.5774)×(−1)
+ 0.5×(1 − 0.57742 )(1 − 0.5774)×(1)
+ 0.5×(1 + 0.5774)(1 − 0.57742 )×(0)
+ 0.5×(1 − 0.57742 )(1 + 0.5774)×(−1)
+ 0.5×(1 − 0.5774)(1 − 0.57742 )×(−2)
∴ 𝑣(0.5774, 0.5774) = −0.9553
Displacements at other Gaussian points can be found using the same procedure. MATLAB
can be used to expedite this calculation.
xi=0.5774; eta=0.5774;
u1=1; u2=0; u3=-1; u4=-1; u5=0; u6=1; u7=2; u8=1;
v1=2; v2=1; v3=-1; v4=-1; v5=1; v6=0; v7=-1; v8=-2;
N1=-(1-xi)*(1-eta)*(1+xi+eta)/4; N2=-(1+xi)*(1-eta)*(1-xi+eta)/4;
51
N3=-(1+xi)*(1+eta)*(1-xi-eta)/4; N4=-(1-xi)*(1+eta)*(1+xi-eta)/4;
N5=(1-xi*xi)*(1-eta)/2; N6=(1+xi)*(1-eta*eta)/2;
N7=(1-xi*xi)*(1+eta)/2; N8=(1-xi)*(1-eta*eta)/2;
u=N1*u1+N2*u2+N3*u3+N4*u4+N5*u5+N6*u6+N7*u7+N8*u8
v=N1*v1+N2*v2+N3*v3+N4*v4+N5*v5+N6*v6+N7*v7+N8*v8
The results are shown below.
𝑢(0.5774, 0.5774) = 1.6923 𝑣(0.5774, 0.5774) = −0.9553

𝑢(−0.5774, 0.5774) = 1.6218 𝑣(−0.5774, 0.5774) = −1.7956
𝑢(−0.5774, −0.5774) = 1.3075 𝑣(−0.5774, −0.5774) = −0.3778
𝑢(0.5774, −0.5774) = 1.0445 𝑣(0.5774, −0.5774) = 0.1290
Example 7.10
Calculate the stiffness matrix of a plane stress CST element as shown in the figure below.
y
3(0.5,0.5)
x
1(0,0) 2(1,0)
Solution
To generate the elemental stiffness matrix for a CST element, MATLAB can be used to
replicate the following matrix procedures.
Step 1: Material matrix (Generalised Hooke’s law)
For plane stress, this is:
1 𝜈 0
𝐸 𝜈 1 0
[𝐄] = 2
[ 1 − 𝜈]
1−𝜈
0 0
2
For plane strain, this is:
1−𝜈 𝜈 0
𝐸 𝜈 1−𝜈 0
[𝐄] = [ 1 − 2𝜈]
(1 + 𝜈)(1 − 2𝜈)
0 0
2
Step 2: Calculate the strain-displacement matrix, [𝐁] from
1 𝑦23 0 𝑦31 0 𝑦12 0

[𝐁𝟏 ] = [ 0 𝑥32 0 𝑥13 0 𝑥21]
2𝐴 𝑥 𝑦23 𝑥13 𝑦31 𝑥21 𝑦12
32
52
Step 3: Calculate the stiffness matrix, using 𝐤 = ∫𝑉 𝐁 𝑇 𝐄𝐁 𝑑𝑉 = 𝑡𝐴(𝐁 𝑇 𝐄𝐁)
[𝐤] = ∫[𝐁𝟏 ]𝑇 [𝐄][𝐁𝟏 ] 𝑑𝑉 = ∫ 𝑑𝑉 ([𝐁𝟏 ]𝑇 [𝐄][𝐁𝟏 ]) = 𝑡𝐴([𝐁𝟏 ]𝑇 [𝐄][𝐁𝟏 ])

𝑉 𝑉
𝑦23 0 𝑥32
0 𝑥32 𝑦23 1 𝜈 0 𝑦23 0 𝑦31 0 𝑦12 0
𝑡 𝑦31 0 𝑥13 𝐸 𝜈 1 0
=
4𝐴 0 𝑥13
(
𝑦31 1 − 𝜈 2
[ 1 − 𝜈 ]) [ 0 𝑥32 0 𝑥13 0 𝑥21 ]
0 0 𝑥32 𝑦23 𝑥13 𝑦31 𝑥21 𝑦12
𝑦12 0 𝑥21 2
[ 0 𝑥21 𝑦12 ]
Step 4: The following MATLAB code can be used to calculate the stiffness matrix, assuming
some basic values for 𝐸, 𝜈, and 𝑡, and using the co-ordinates of the nodes provided in the
diagram.
E=1;
nu=0.3;
t=1;
x1=0;
y1=0;
x2=1;
y2=0;
x3=0.5;
y3=0.5;
A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;
x21=x2-x1;
x13=x1-x3;
x32=x3-x2;
y12=y1-y2;
y31=y3-y1;
y23=y2-y3;
Bmat=[y23 0 y31 0 y12 0; 0 x32 0 x13 0 x21; x32 y23 x13 y31 x21
y12];
Emat=(E/(1-nu*nu))*[1 nu 0; nu 1 0; 0 0 (1-nu)/2];
Ke=t/(4*A)*Bmat'*Emat*Bmat
Ke =
0.3709 0.1786 -0.1786 -0.0137 -0.1923 -0.1648

0.1786 0.3709 0.0137 0.1786 -0.1923 -0.5495
-0.1786 0.0137 0.3709 -0.1786 -0.1923 0.1648
-0.0137 0.1786 -0.1786 0.3709 0.1923 -0.5495
-0.1923 -0.1923 -0.1923 0.1923 0.3846 0
-0.1648 -0.5495 0.1648 -0.5495 0 1.0989
53
Example 7.11
Sketch the distribution of v displacement and y strain in the following two connected CST
elements, given a global displacement vector {𝐝} as shown.
ìuA ü ì 1 ü v
ïv ï ï 0 ï
ï Aï ï ï
ïuB ï ï 4 ï
ï ï ï ï
ï v ï ï - 1ï A (0,0)
{d}= í B ý = í ý
ï uC ï ï 2 ï
ï vC ï ï - 2ï x D (1,0) 1
2
ï ï ï ï B (0,1)
ïu D ï ï 3 ï
ïv ï ïî 1 ïþ
î Dþ C (1,1) y
Solution
Step 1: Select v displacement at each node in the {𝐝} vector.
{𝐯} = {𝑣𝐴 𝑣𝐵 𝑣𝐶 𝑣𝐷 }𝑇 = {0 −1 −2 1}𝑇
Step 2: Draw the displacement at each corresponding node and connect them as shown. Note
that the shape functions in a CST result in a linear distribution of displacements, so straight
lines can be used to connect the displacements at each node.
A (0,0)
D (1,0) 1
x 2
B (0,1)
C (1,1)
y
Step 3: Calculate the y strains in each of the CST elements according the following equation.
𝛆 = 𝐁𝐝
𝑢1
𝑣1
𝜀𝑥𝑥
1 𝑦2 − 𝑦3 0 𝑦3 − 𝑦1 0 𝑦1 − 𝑦2 0
𝑢2
𝜀
{ 𝑦𝑦 } = [ 0 𝑥3 − 𝑥2 0 𝑥1 − 𝑥3 0 𝑥2 − 𝑥1 ] 𝑣
𝜀𝑥𝑦 2𝐴 𝑥 − 𝑥 𝑦2 − 𝑦3 𝑥1 − 𝑥3 𝑦3 − 𝑦1 𝑥2 − 𝑥1 𝑦1 − 𝑦2 𝑢2
3 2
3
{𝑣3 }
Note that the node numbering (1, 2, 3) should be in a counter-clockwise sequence, which is
consistent with the sample CST element shown in Fig. 7.21. For element 1, the local nodes
(1,2,3) can be the global nodes (A,D,C), but not nodes (A,C,D). Similarly for element 2, local
nodes (1,2,3) can be global nodes (C,B,A), but not nodes (A,B,C).
54
We can create co-ordinate vectors for the positions of the nodes for each element.
𝑥1 𝑥𝐴 0
𝑦1 𝑦𝐴 0
𝑥2 𝑥𝐷 1
For element 1: 𝐱 (1) = 𝑦 = 𝑦 =
2 𝐷 0
𝑥3 𝑥𝐶 1
{𝑦3 } {𝑦𝐶 } {1}
𝑥1 𝑥𝐶 1
𝑦1 𝑦𝐶 1
𝑥2 𝑥𝐵 0
For element 2: 𝐱 (2) = 𝑦 = 𝑦 =
2 𝐵 1
𝑥3 𝑥𝐴 0
{𝑦3 } {𝑦𝐴 } {0}
Therefore, the strain in element 1, where (Node 1 = Node A, Node 2 = Node D, Node 3 =
Node C) is:
(1) 1
𝜀𝑦𝑦 = [(𝑥 − 𝑥2 )𝑣1 + (𝑥1 − 𝑥3 )𝑣2 + (𝑥2 − 𝑥1 )𝑣3]
2𝐴 3
1
= [(𝑥 − 𝑥𝐷 )𝑣𝐴 + (𝑥𝐴 − 𝑥𝐶 )𝑣𝐷 + (𝑥𝐷 − 𝑥𝐴 )𝑣𝐶 ]
2𝐴 𝐶
1
= [(1 − 1)×(0) + (0 − 1)×(1) + (1 − 0)×(−2)] = −3
2×0.5
And the strain in element 2, where (Node 1 = Node C, Node 2 = Node B, Node 3 = Node A)
is:
(2) 1
𝜀𝑦𝑦 = [(𝑥 − 𝑥2)𝑣1 + (𝑥1 − 𝑥3 )𝑣2 + (𝑥2 − 𝑥1 )𝑣3]
2𝐴 3
1
= [(𝑥 − 𝑥𝐵 )𝑣𝐶 + (𝑥𝐶 − 𝑥𝐴 )𝑣𝐵 + (𝑥𝐵 − 𝑥𝐶 )𝑣𝐴 ]
2𝐴 𝐴
1
= [(0 − 0)×(−2) + (1 − 0)×(−1) + (0 − 1)×(0)] = −1
2×0.5
Step 4: Plot the y strain for each element. There should be no variation of strain within the
element, because they are constant strain elements!
e yy
D (1,0) A (0,0)
1
x 2
B (0,1)
C (1,1) -1
-3 y
-1
-3
55
Example 7.12
Determine the co-ordinates of the Gaussian integration points (G1, G2, G3, G4) in the
Cartesian co-ordinate system of isoparametric element Q4 as shown (the nodal coordinates
have been given in the figure). Compute the Jacobian matrix [𝐉] and det[𝐉] for the Q4
element.
Solution
Q4 and Q8 are isoparametric elements in which co-ordinates and displacement share the
same shape functions for mapping and interpolation, as follows.
Co-ordinate mapping:
4 4

𝑖=1 𝑖=1
Displacement interpolation:
4 4

𝑖=1 𝑖=1
Step 1: Node numbering: Node (1,2,3,4) = Node (A,B,C,D)
Step 2: Co-ordinates in terms of the shape functions.
The shape functions for a Q4 element are
1 1
𝑁1 = (1 − 𝜉 )(1 − 𝜂) 𝑁2 = (1 + 𝜉 )(1 − 𝜂)
4 4
1 1
𝑁3 = (1 + 𝜉)(1 + 𝜂) 𝑁4 = (1 − 𝜉)(1 + 𝜂)
4 4
Therefore, the Cartesian co-ordinates (x,y) in terms of the natural co-ordinates (ξ, η) can be
derived.
56
𝑥 = ∑ 𝑁𝑖 𝑥𝑖
𝑖=1
1 1 1
= (1 − 𝜉)(1 − 𝜂)𝑥1 + (1 + 𝜉)(1 − 𝜂)𝑥2 + (1 + 𝜉)(1 + 𝜂)𝑥3
4 4 4
1
+ (1 − 𝜉)(1 + 𝜂)𝑥4
4
1 1 1
= (1 − 𝜉)(1 − 𝜂)1 + (1 + 𝜉)(1 − 𝜂)2 + (1 + 𝜉)(1 + 𝜂)2.25
4 4 4
1
+ (1 − 𝜉)(1 + 𝜂)1.25
4
1
∴ 𝑥(𝜉, 𝜂) = (6.5 + 2𝜉 + 0.5𝜂)
4
4
𝑦 = ∑ 𝑁𝑖 𝑦𝑖
𝑖=1
1 1 1
= (1 − 𝜉)(1 − 𝜂)𝑦1 + (1 + 𝜉)(1 − 𝜂)𝑦2 + (1 + 𝜉)(1 + 𝜂)𝑦3
4 4 4
1
+ (1 − 𝜉)(1 + 𝜂)𝑦4
4
1 1 1
= (1 − 𝜉)(1 − 𝜂)0 + (1 + 𝜉)(1 − 𝜂)0 + (1 + 𝜉)(1 + 𝜂)1.5
4 4 4
1
+ (1 − 𝜉)(1 + 𝜂)1
4
1
∴ 𝑦(𝜉, 𝜂) = (2.5 + 0.5𝜉 + 2.5𝜂 + 0.5𝜉𝜂)
4
Step 3: Calculate the co-ordinates of G1.
1
𝑥(𝜉 = 0.5774, 𝜂 = 0.5774) = (6.5 + 2×0.5774 + 0.5×0.5774)
4
= 1.9859
1 2.5 + 0.5×0.5774 + 2.5×0.5774
𝑦(𝜉 = 0.5774, 𝜂 = 0.5774) = ( )
4 +0.5×0.5774×0.5774
= 1.0997
This can be repeated for the rest of the Gaussian points by simply substituting in their natural
co-ordinates.
Step 4: Calculate the Jacobian matrix.
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑢 𝜕𝑥 𝜕𝑦
𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝑥 𝜕𝜉 𝜕𝜉
= → [𝐉] =
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑢 𝜕𝑥 𝜕𝑦
{𝜕𝜂} ⏟[𝜕𝜂 𝜕𝜂] {𝜕𝑦} [𝜕𝜂 𝜕𝜂]
[𝐉]
57
𝜕 1 𝜕 1
[ (6.5 + 2𝜉 + 0.5𝜂)] [ (2.5 + 0.5𝜉 + 2.5𝜂 + 0.5𝜉𝜂)]
𝜕𝜉 4 𝜕𝜉 4
∴ [𝐉] =
𝜕 1 𝜕 1
[ (6.5 + 2𝜉 + 0.5𝜂)] [ (2.5 + 0.5𝜉 + 2.5𝜂 + 0.5𝜉𝜂)]
[𝜕𝜂 4 𝜕𝜂 4 ]
1 (2) (0.5 + 0.5𝜂)

[𝐉(𝜉, 𝜂)] = [ ]
4 (0.5) (2.5 + 0.5𝜉)
For Gaussian point G1:
1 (2) (0.5 + 0.5×0.5774) 0.5 0.1972

[𝐉(0.5774,0.5774)] = [ ]=[ ]
4 (0.5) (2.5 + 0.5×0.5774) 0.125 0.6972
This can be repeated for the other Gaussian points.
Step 5: Calculate the determinant of the Jacobian
1 (2) (0.5 + 0.5𝜂)

det[𝐉(𝜉, 𝜂)] = [ ] = 4.75 + 𝜉 − 0.25𝜂
4 (0.5) (2.5 + 0.5𝜉)
which is a linear function of the natural co-ordinates (ξ, η). At G1:
0.5 0.1972
det[𝐉(0.5774,0.5774)] = det [ ] = 0.5×0.6972 − 0.125×0.1972 = 0.3239
0.125 0.6972
Again, this can be repeated for the other Gaussian points.
7.8. Three-dimensional finite elements

For 3D problems, we introduce a z dimension to displacements, strains, and stresses. The
same methods can be applied as for 2D elements, but with larger matrices and vectors
including the necessary z components, as shown below.
3D stress vector:
𝛔 = {𝜎𝑥𝑥 𝜎𝑦𝑦 𝜎𝑧𝑧 𝜎𝑥𝑦 𝜎𝑦𝑧 𝜎𝑧𝑥 }𝑇
3D strain vector:
𝛆 = {𝜀𝑥𝑥 𝜀𝑦𝑦 𝜀𝑧𝑧 𝜀𝑥𝑦 𝜀𝑦𝑧 𝜀𝑧𝑥 }𝑇
3D Hooke’s law:
𝛔(6×1) = 𝐄(6×6) 𝛆(6×1)
58
1−𝜈 𝜈 𝜈 0 0 0
𝜈 1−𝜈 𝜈 0 0 0
𝜎𝑥𝑥 𝜀𝑥𝑥
𝜈 𝜈 1−𝜈 0 0 0
𝜎𝑦𝑦 𝜀𝑦𝑦
1 − 2𝜈
𝜎𝑧𝑧 𝐸 0 0 0 0 0 𝜀𝑧𝑧
𝜎𝑥𝑦 = 2 𝜀𝑥𝑦
(1 + 𝜈)(1 − 2𝜈) 1 − 2𝜈
𝜎𝑦𝑧 0 0 0 0 0 𝜀𝑦𝑧
{ 𝜎𝑧𝑥 } 2 { 𝜀𝑧𝑥 }
1 − 2𝜈
[ 0 0 0 0 0
2 ]
3D displacement:
𝑢(𝑥, 𝑦, 𝑧)
𝐮 = { 𝑣(𝑥, 𝑦, 𝑧) }
𝑤(𝑥, 𝑦, 𝑧)
7.8.1. 3D finite element formulation

The displacement field can be expressed using 3D shape functions.
𝑁 𝑁 𝑁
𝑢 = ∑ 𝑁𝑖 𝑢𝑖 𝑣 = ∑ 𝑁𝑖 𝑣𝑖 𝑤 = ∑ 𝑁𝑖 𝑤𝑖
𝑖=1 𝑖=1 𝑖=1
or in matrix form:
𝐮 = 𝐍𝐝
𝑢1
𝑣1
𝑢(𝑥, 𝑦, 𝑧) 𝑁1 0 0 𝑁2 0 0 ⋯ 𝑤1
{ 𝑣(𝑥, 𝑦, 𝑧) } =[0 𝑁1 0 0 𝑁2 0 ⋯] 𝑢2
𝑤(𝑥, 𝑦, 𝑧) (3×1) 0 0 𝑁1 0 0 𝑁2 ⋯ (3×3𝑁) 𝑣2
𝑤2
{ ⋮ }(3𝑁×1)
where 𝑁𝑖 are the shape functions for a 3D element and
𝐝 = {𝑢1 𝑣1 𝑢2 𝑣2 𝑢3 𝑣3 … }𝑇 is the nodal displacement vector.
Using the strain-displacement relation, we can have the following formula for strain.
𝛆(6×1) = 𝐁(6×3𝑁) 𝐝(3𝑁×1)
The stiffness matrix can once again be determined using the same formula used for 2D
elements (𝐤 = ∫𝑉 𝐁 𝑇 𝐄𝐁 𝑑𝑉).
𝐤 (3𝑁×3𝑁) = ∫ 𝐁 𝑇 (3𝑁×6) 𝐄(6×6) 𝐁(6×3𝑁) 𝑑𝑉 𝑒

𝑉𝑒
59
Numerical quadrature rules can be used to evaluate the above integration.
7.8.2. Typical 3D elements

3D elements come in a variety of shapes (tetrahedral, hexahedral, pentahedral) as shown
in Fig. 7.25. Being 3D elements, each node has 3 degrees of freedom (𝑢, 𝑣, 𝑤) in these 3D
elements.
(a) Tetrahedron
(b) Hexahedron
(c) Pentahedron
Fig 7.25 3D finite element shapes
Note that the 4 node tetrahedral element is a constant strain element. Thus we do not
recommend using 4 node tetrahedral elements for capturing high stress/strain gradients.
60
7.8.3. Linear hexahedron (8-node brick) element

Just like the 2D elements in section 7.7, 3D elements also utilise a natural co-ordinate system,
which can be mapped to the global co-ordinate system. Consider the linear hexahedron (8-
node brick) element in Fig. 7.26.
Fig 7.26 Linear hexahedron element
Co-ordinate transformation from the 3D natural co-ordinate system (ξ,η,ζ) to the global co-
ordinate system (x,y,z) can be achieved using the shape functions.
8 8 8
𝑥 = ∑ 𝑁𝑖 𝑥𝑖 𝑦 = ∑ 𝑁𝑖 𝑦𝑖 𝑧 = ∑ 𝑁𝑖 𝑧𝑖
𝑖=1 𝑖=1 𝑖=1
For a linear hexahedron element, the shape functions are:
1 1
𝑁1 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 − 𝜂)(1 − 𝜁) 𝑁2 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 − 𝜂)(1 − 𝜁)
8 8
1 1
𝑁3 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 + 𝜂)(1 − 𝜁) 𝑁4 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 + 𝜂)(1 − 𝜁)
8 8
1 1
𝑁5 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 − 𝜂)(1 + 𝜁) 𝑁6 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 − 𝜂)(1 + 𝜁)
8 8
1 1
𝑁7 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 + 𝜂)(1 + 𝜁) 𝑁8 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 + 𝜂)(1 + 𝜁)
8 8
The linear hexahedron element is isoparametric, meaning the above shape functions can be
used for the displacement field, as well as for co-ordinate mapping.
𝑁 𝑁 𝑁
𝑖=1 𝑖=1 𝑖=1
To calculate strains in the element, we need to apply the strain-displacement relationship,

which requires the Jacobian matrix containing derivatives with respect to different co-
ordinate systems.
For 𝑢 we have:
61
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑢 𝜕𝑢
𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝑥 𝜕𝑥
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑢 𝜕𝑢
= =𝐉
𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝑦 𝜕𝑦
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑢 𝜕𝑢
{ 𝜕𝜁} [𝜕𝜁 𝜕𝜁 𝜕𝜁] { 𝜕𝑧 } { 𝜕𝑧 }
or
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 −1 𝜕𝑢 𝜕𝑢
𝜕𝑥 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑢 𝜕𝑢
= = 𝐉 −1
𝜕𝑦 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂
𝜕𝑢 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑢 𝜕𝑢
{ 𝜕𝑧 } [𝜕𝜁 𝜕𝜁 𝜕𝜁] { 𝜕𝜁} { 𝜕𝜁}
where 𝐉 is the Jacobian matrix. Similarly for 𝑣 and 𝑤:
𝜕𝑣 𝜕𝑥 𝜕𝑦 𝜕𝑧 −1 𝜕𝑣 𝜕𝑣
𝜕𝑥 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉
𝜕𝑣 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑣 𝜕𝑣
= = 𝐉 −1
𝜕𝑦 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂
𝜕𝑣 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑣 𝜕𝑣
{𝜕𝑧 } [𝜕𝜁 𝜕𝜁 𝜕𝜁] {𝜕𝜁} {𝜕𝜁}
and
𝜕𝑤 𝜕𝑥 𝜕𝑦 𝜕𝑧 −1 𝜕𝑤 𝜕𝑤
𝜕𝑥 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉
𝜕𝑤 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑤 𝜕𝑤
= = 𝐉 −1
𝜕𝑦 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂
𝜕𝑤 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑤 𝜕𝑤
{ 𝜕𝑧 } [𝜕𝜁 𝜕𝜁 𝜕𝜁] { 𝜕𝜁 } { 𝜕𝜁 }
Strain can then be calculated in terms of the natural co-ordinate system, and then transformed
into the global co-ordinate system using the equations above. For the natural co-ordinate strain
𝜕𝑢/𝜕𝜉, the calculation is as follows.
8
𝜕𝑢 𝜕𝑁𝑖
=∑ 𝑢
𝑖=1
𝜕𝑢 𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6 𝜕𝑁7 𝜕𝑁8
= 𝑢1 + 𝑢2 + 𝑢3 + 𝑢4 + 𝑢5 + 𝑢6 + 𝑢7 + 𝑢
𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 8
62
The Jacobian matrix is then calculated using the shape functions given the isoparametric
nature of the linear hexahedron element:
8 8 8
𝑥 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑥𝑖 𝑦 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑦𝑖 𝑧 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑧𝑖

𝑖=1 𝑖=1 𝑖=1
Therefore, the Jacobian matrix is:
8 8 8
𝜕𝑁𝑖 𝜕𝑁𝑖 𝜕𝑁𝑖
∑ 𝑥 ∑ 𝑦 ∑ 𝑧
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝜉 𝑖 𝜕𝜉 𝑖 𝜕𝜉 𝑖
𝑖=1 𝑖=1 𝑖=1
𝜕𝜉 𝜕𝜉 𝜕𝜉 8 8 8
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑁 𝜕𝑁𝑖 𝜕𝑁𝑖
𝐉= = ∑ 𝑖 𝑥𝑖 ∑ 𝑦 ∑ 𝑧
𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝑖 𝜕𝜂 𝑖
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝑖=1 𝑖=1 𝑖=1
8 8 8
[𝜕𝜁 𝜕𝜁 𝜕𝜁] 𝜕𝑁𝑖 𝜕𝑁𝑖 𝜕𝑁𝑖
∑ 𝑥 ∑ 𝑦 ∑ 𝑧
𝜕𝜁 𝑖 𝜕𝜁 𝑖 𝜕𝜁 𝑖
[ 𝑖=1 𝑖=1 𝑖=1 ]
𝑥1 𝑦1 𝑧1
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6 𝜕𝑁7 𝜕𝑁8 𝑥 𝑦2 𝑧2
2
𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝑥3 𝑦3 𝑧3
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6 𝜕𝑁7 𝜕𝑁8 𝑥4 𝑦4 𝑧4
=
𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝑥5 𝑦5 𝑧5
𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝜕𝑁4 𝜕𝑁5 𝜕𝑁6 𝜕𝑁7 𝜕𝑁8 𝑥6 𝑦6 𝑧6
𝑥 𝑦7 𝑧7
[ 𝜕𝜁 𝜕𝜁 𝜕𝜁 𝜕𝜁 𝜕𝜁 𝜕𝜁 𝜕𝜁 𝜕𝜁 ] 7
[𝑥 8 𝑦8 𝑧8 ]
Use the strain-displacement relationship to find the global co-ordinate strains.
𝜀𝑥𝑥 𝜕𝑢⁄𝜕𝑥
𝜀𝑦𝑦 𝜕𝑣⁄𝜕𝑦
𝜀𝑧𝑧 𝜕𝑤 ⁄𝜕𝑧
𝛆(6×1) = 𝜀 = = 𝐁(6×24) 𝐝(24×1)
𝑥𝑦 (𝜕𝑣⁄𝜕𝑥 + 𝜕𝑢⁄𝜕𝑦)/2
𝜀𝑦𝑧 (𝜕𝑤⁄𝜕𝑦 + 𝜕𝑣⁄𝜕𝑧)/2
{ 𝜀𝑧𝑥 } {(𝜕𝑢⁄𝜕𝑧 + 𝜕𝑤 ⁄𝜕𝑥)/2}
where 𝐁(6×24) = [𝐁1 𝐁2 𝐁3 𝐁4 𝐁5 𝐁6 𝐁7 𝐁8 ] is the strain-displacement matrix

and the displacement vector 𝐝(24×1) = {𝑢1 𝑣1 𝑤1 𝑢2 𝑣2 𝑤2 ⋯ 𝑢8 𝑣8 𝑤8 }𝑇 .
𝜕𝑁𝑖 ⁄𝜕𝑥 0 0
0 𝜕𝑁𝑖 ⁄𝜕𝑦 0
0 0 𝜕𝑁𝑖 ⁄𝜕𝑧
Each term in the strain-displacement matrix is 𝐁i =
0 𝜕𝑁𝑖 ⁄𝜕𝑧 𝜕𝑁𝑖 ⁄𝜕𝑦
𝜕𝑁𝑖 ⁄𝜕𝑧 0 𝜕𝑁𝑖 ⁄𝜕𝑥
[𝜕𝑁𝑖 ⁄𝜕𝑦 𝜕𝑁𝑖 ⁄𝜕𝑥 0 ](6×3)
63
The Jacobian transformation for chain rule of differentiation is needed to fill the strain-
displacement matrix terms.
𝜕𝑁𝑖 𝜕𝑥 𝜕𝑦 𝜕𝑧 −1 𝜕𝑁𝑖 𝜕𝑁𝑖

𝜕𝑥 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉
𝜕𝑁𝑖 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑁𝑖 𝜕𝑁𝑖
= = 𝐉 −1
𝜕𝑦 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂 𝜕𝜂
𝜕𝑁𝑖 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑁𝑖 𝜕𝑁𝑖
{ 𝜕𝑧 } [𝜕𝜁 𝜕𝜁 𝜕𝜁] { 𝜕𝜁 } { 𝜕𝜁 }
Thus, the Jacobian matrix needs to be positive definite to enable a non-trivial invertible
matrix.
The stiffness matrix is once again derived using the energy formulation equating strain
energy (𝑈) and work done
1 1 1 1
𝑈= ∫ 𝛔𝛵 𝛆 𝑑𝑉 = ∫(𝐄𝛆)𝛵 𝜺 𝑑𝑉 = ∫ 𝛆𝛵 𝐄𝛆 𝑑𝑉 = ∫(𝚩𝐝)𝛵 𝐄(𝚩𝐝) 𝑑𝑉
2 2 2 2
𝑉 𝑉 𝑉 𝑉
1
= 𝐝𝑇 [ ∫ 𝚩 𝛵 𝐄𝚩 𝑑𝑉] 𝐝
2
𝑉
1
∴ 𝐤 (24×24) = ∫ 𝚩 𝛵 (24×6) 𝐄(6×6) 𝚩(6×24) 𝑑𝑉
2
𝑉
In the natural co-ordinate system, 𝑑𝑉 = (det 𝐉)𝑑𝜉𝑑𝜂𝑑𝜁, so therefore the elemental stiffness
matrix can be expressed in natural co-ordinates, and can be evaluated numerically.
1 1 1
1
𝐤 (24×24) = ∫ ∫ ∫ 𝚩 𝛵 (24×6) 𝐄(6×6) 𝚩(6×24) (det 𝐉) 𝑑𝜉𝑑𝜂𝑑𝜁
2
−1 −1 −1
Example 7.13
The 8-node isoparametric brick element shown in Cartesian coordinate system (x,y,z) below
is mapped to the natural coordinate (ξ,η,ζ). Use shape functions to determine (a) the Cartesian
co-ordinate of centroid, and (b) the displacements 𝑢, 𝑣, 𝑤 at the centroid if the nodal
displacement is:
𝐝𝑇 = [0 0 𝐝3 𝐝4 0 0 𝐝7 𝐝8 ]1×24
The non-zero nodal displacements are:
𝐝3 = {1,2, −1}, 𝐝4 = {−2,0, −2}, 𝐝7 = {2,2,2}, 𝐝8 = {−1,1,2}
64
ζ
Solution
Nodal positions
𝐱 = {0,0,0, 2,0,0, 2,1,0, 0,3,0, 0,0,4, 3,0,3, 3,3,3, 0,4,4}𝑇
Nodal displacements
𝐮 = {0,0,0, 0,0,0, 1,2, −1, − 2,0, −2, 0,0,0, 0,0,0, 2,2,2, − 1,1,2}𝑇
Shape functions
1 1
𝑁1 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 − 𝜂)(1 − 𝜁) 𝑁2 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 − 𝜂)(1 − 𝜁)
8 8
1 1
𝑁3 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 + 𝜂)(1 − 𝜁) 𝑁4 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 + 𝜂)(1 − 𝜁)
8 8
1 1
𝑁5 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 − 𝜂)(1 + 𝜁) 𝑁6 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 − 𝜂)(1 + 𝜁)
8 8
1 1
𝑁7 (𝜉, 𝜂, 𝜁) = (1 + 𝜉)(1 + 𝜂)(1 + 𝜁) 𝑁8 (𝜉, 𝜂, 𝜁) = (1 − 𝜉)(1 + 𝜂)(1 + 𝜁)
8 8
Co-ordinate mapping to find global co-ordinates of the centroid.
8 8 8
𝑥 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑥𝑖 𝑦 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑦𝑖 𝑧 = ∑ 𝑁𝑖 (𝜉, 𝜂, 𝜁) 𝑧𝑖

𝑖=1 𝑖=1 𝑖=1
For the x co-ordinate:

𝑥0 (𝜉 = 0, 𝜂 = 0, 𝜁 = 0) = 𝑁1 𝑥1 + 𝑁2 𝑥2 + 𝑁3 𝑥3 + 𝑁4 𝑥4 + 𝑁5 𝑥5 + 𝑁6 𝑥6 + 𝑁7 𝑥7 + 𝑁8 𝑥8
1 1
= (1 − 0)(1 − 0)(1 − 0)(0) + (1 + 0)(1 − 0)(1 − 0)(2)
8 8
1 1
+ (1 + 0)(1 + 0)(1 − 0)(2) + (1 − 0)(1 + 0)(1 − 0)(0)
8 8
1 1
+ (1 − 0)(1 − 0)(1 + 0)(0) + (1 + 0)(1 − 0)(1 + 0)(3)
8 8
1 1
+ (1 + 0)(1 + 0)(1 + 0)(3) + (1 − 0)(1 + 0)(1 + 0)(0)
8 8
2 2 3 3 10
𝑥0 = 0 + + + 0 + 0 + + + 0 = = 1.25
8 8 8 8 8
65
Similarly, the y and z co-ordinates of the centroid can be found.
𝑥0 = 1.25, 𝑦0 = 1.375, 𝑧0 = 1.75
Determine the displacement at the centroid using the shape functions.
8 8 8
𝑖=1 𝑖=1 𝑖=1
For u:
𝑢0 (𝜉 = 0, 𝜂 = 0, 𝜁 = 0) = 𝑁1 𝑢1 + 𝑁2 𝑢2 + 𝑁3 𝑢3 + 𝑁4 𝑢4 + 𝑁5 𝑢5 + 𝑁6 𝑢6 + 𝑁7 𝑢7 + 𝑁8 𝑢8
1 1
= (1 − 0)(1 − 0)(1 − 0)(0) + (1 + 0)(1 − 0)(1 − 0)(0)
8 8
1 1
+ (1 + 0)(1 + 0)(1 − 0)(1) + (1 − 0)(1 + 0)(1 − 0)(−2)
8 8
1 1
+ (1 − 0)(1 − 0)(1 + 0)(0) + (1 + 0)(1 − 0)(1 + 0)(0)
8 8
1 1
+ (1 + 0)(1 + 0)(1 + 0)(2) + (1 − 0)(1 + 0)(1 + 0)(−1)
8 8
𝑢0 (𝜉 = 0, 𝜂 = 0, 𝜁 = 0) = 0
Similarly, the v and w displacements can be calculated at the centroid.
𝑢0 = 0, 𝑣0 = 0.625, 𝑤0 = 0.125
Alternatively, this calculation can be performed in MATLAB.
xi=0; nu=0; zt=0;

x1=0; x2=2; x3=2; x4=0; x5=0; x6=3; x7=3; x8=0;
y1=0; y2=0; y3=1; y4=3; y5=0; y6=0; y7=3; y8=4;
z1=0; z2=0; z3=0; z4=0; z5=4; z6=3; z7=3; z8=4;
u1=0; u2=0; u3=1; u4=-2; u5=0; u6=0; u7=2; u8=-1;
v1=0; v2=0; v3=2; v4=0; v5=0; v6=0; v7=2; v8=1;
w1=0; w2=0; w3=-1;w4=-2; w5=0; w6=0; w7=2; w8=2;
N1=(1-xi)*(1-nu)*(1-zt)/8; N2=(1+xi)*(1-nu)*(1-zt)/8;
N3=(1+xi)*(1+nu)*(1-zt)/8; N4=(1-xi)*(1+nu)*(1-zt)/8;
N5=(1-xi)*(1-nu)*(1+zt)/8; N6=(1+xi)*(1-nu)*(1+zt)/8;
N7=(1+xi)*(1+nu)*(1+zt)/8; N8=(1-xi)*(1+nu)*(1+zt)/8;
x=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8
y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8
u=N1*u1+N2*u2+N3*u3+N4*u4+N5*u5+N6*u6+N7*u7+N8*u8
v=N1*v1+N2*v2+N3*v3+N4*v4+N5*v5+N6*v6+N7*v7+N8*v8
w=N1*w1+N2*w2+N3*w3+N4*w4+N5*w5+N6*w6+N7*w7+N8*w8
66
7.9. Practical FE modelling techniques

In this section, we will briefly discuss some issues and techniques associated with FE models
which may be relevant in practice. Understanding these points might help you when your
model doesn’t work, or does something you don’t expect.
7.9.1. Symmetry
Symmetry was discussed in chapter 4 to reduce the complexity of the analytical solution. In
the same way, introducing symmetrical boundary conditions will reduce the size of the
problem, and hence make the stiffness matrix smaller and lower the total DOFs in the model.
Different types of symmetry are listed below, including the reflective symmetry discussed in
chapter 4.
Reflective (mirror, bilateral) symmetry;

Rotational (cyclic) symmetry;
Translational symmetry; and
Axisymmetry
Axisymmetry is a commonly used and very effective technique for simplifying the
complexity of models. Many typical engineering problems exhibit axisymmetric
characteristics. For axisymmetric problems, there are two conditions: (1) rotational
structure and (2) axisymmetric loading/constraints. Examples of rotational structures
suitable for axisymmetric modelling are shown in Fig. 7.27.
Fig 7.27 Rotational structures suitable for axisymmetric modelling
As these structures are rotational, the cylindrical co-ordinate system (r,θ,z) is used instead
of a Cartesian co-ordinate system (x,y,z) as shown in Fig. 7.28.
Fig 7.28 Cylindrical co-ordinates for axisymmetric problems
67
In axisymmetric problems, a 2D circumferential cross-section (r-z plane) is modelled, and

rotated about the z-axis to produce the full 3D solution.
The displacement fields are 𝑢(𝑟, 𝑧), 𝑤(𝑟, 𝑧). There is no v displacement (circumferential
displacement), as it vanishes in an axisymmetric assumption (similar to how the displacement
normal to the symmetry plane in planar symmetry is constrained).
The strain-displacement relationship for axisymmetric problems can be written in 2D.
𝜕𝑢(𝑟, 𝑧)
𝜀𝑟𝑟 =
𝜕𝑟
𝑢(𝑟, 𝑧)
𝜀𝜃𝜃 =
𝑟
𝜕𝑤(𝑟, 𝑧)
𝜀𝑧𝑧 =
𝜕𝑧
𝜕𝑤(𝑟, 𝑧) 𝜕𝑢(𝑟, 𝑧)
𝜀𝑟𝑧 = +
𝜕𝑟 𝜕𝑧
Note that 𝜀𝑟𝜃 = 𝜀𝜃𝑧 = 0.
The stress-strain relationships can also be derived, with a non-zero 𝜎𝜃𝜃 .
𝜎𝑟𝑟 1−𝜈 𝜈 𝜈 0 𝜀𝑟𝑟

𝜎𝜃𝜃 𝐸 𝜈 1−𝜈 𝜈 0 𝜀𝜃𝜃
{𝜎 } = 𝜈 𝜈 1−𝜈 0 {𝜀 }
𝑧𝑧 (1 + 𝜈)(1 − 2𝜈) 1 − 2𝜈 𝑧𝑧
𝜎𝑟𝑧 𝜀𝑟𝑧
[ 0 0 0
2 ]
Include this in the cheat sheet
In the context of finite elements, special 2D axisymmetric elements exist (Fig. 7.29), which will
represent a ring when swept around an axis.
Fig 7.29 2D axisymmetric finite elements
The elemental stiffness matrix for such an element is derived in cylindrical co-ordinates.
2𝜋 1 1
𝐤 = ∫ 𝐁 𝑇 𝐄𝐁 𝑟𝑑𝜃 𝑑𝑟𝑑𝑧 = ∫ ∫ ∫ 𝐁 𝑇 𝐄𝐁 𝑟𝑑𝜃 (det 𝐉)𝑑𝜉𝑑𝜂

𝑉
0 −1 −1
1 1
𝐤 = 2𝜋 ∫ ∫ 𝐁 𝑇 𝐄𝐁 𝑟(det 𝐉) 𝑑𝜉 𝑑𝜂
−1 −1
68
An example of an axisymmetric problem is the rotating flywheel in Fig. 7.30. This rotating
wheel can be classified as an axisymmetric problem because it satisfies both conditions, that
is: (1) the structure is rotational, and (2) the body force generated by the angular velocity is
axisymmetric. Thus flywheels are easily modelled in 2D.
Fig 7.30 An axisymmetric rotating flywheel
7.9.2. Error
The main categories of error associated with create FE models include:
Precision error from solving the FE equations;

Modelling error from approximating the real problem to simply solved analogues
(beam, bar, plate, loading assumptions); and
Discretisation error from the use of finite elements and piecewise functions.
7.9.3. The inherent FE assumption

In solid mechanics, the motivation behind the use of the FE method is to create mathematical
models of real structural components and their responses. This is enabled by the
approximations discussed at the beginning of this chapter. FE models have some finite
number of nodes, and therefore a limited number of DOFs. Additionally, the DOFs can be
interpolated within an element from each node using shape functions, however the accuracy
still depends on the nodal results (Fig. 7.31). In other words, the displacement field is
controlled (or constrained) by the values at a limited number of nodes. Compare this to a real
structure, which can be considered to have an infinite number of nodes (or physical points or
particles), and hence an infinite number of DOFs.
Fig 7.31 Displacement field dependence on nodal results
7.9.4. Stiffening effect

In displacement based FE methods, the FE model is always stiffer than the real structure.
This means in general, displacement results are smaller in magnitude than the exact value.
69
Hence, the exact solution provides an upper bound for the FE solution, as shown on the
mesh convergence plot in Fig. 7.32. Note that the mesh convergence plot has the total
number of DOFs on the x-axis, and displacement on the y-axis. As the number of DOFs
increases (which may be a result of mesh refinement), the displacement increases and
approaches the exact solution from below. Such a plot is useful in evaluating the
discretisation error.
Fig 7.32 A typical mesh convergence plot for a displacement based FE analysis
Assignment 3 Q3 Convergence plot
7.9.5. Finite element refinements
To reduce the effect of errors in the FE method, various refinements can be implemented.
These techniques are formally known as h-method, p-method, r-method, and hp-method.
h-method involves reducing the size of the finite elements (element size is typically referred
to as h). This is the most commonly used in commercial software such as ANSYS.
p-method involves increasing the order of the polynomials used for the shape functions in an
element without changing the element size (p refers to the higher order polynomial). This is
less common in commercial software, and is no longer available in ANSYS (now the only
option is to go from linear to quadratic elements).
r-method involves rearranging the nodal positions for a given mesh to give more accurate
solutions in areas of high stress concentration.
hp-method is a combination of both h-method and p-method, which enables the solution to
converge to the exact solution at a faster rate.
Newer algorithms have been developed, including adaptive mesh refinements which
automate the discretisation process based on measures of error in the model, but these
generally utilise one of the above methods for producing a new mesh.
7.9.6. Mesh quality

Mesh quality refers to both the type of the element used for discretisation, and the shape of
the element used. Therefore, it is imperative for the engineer to understand the behaviour of
the elements. For example, T3 and Q4 elements calculate linear displacement fields and give
constant stress and strain, whereas T6 and Q8 elements calculate quadratic displacement
70
fields and give linear stress and strain. By understanding these behaviours, an engineer will
be able to choose the most suitable element for a given problem. When in doubt, use higher
order elements or a finer mesh—this does resolve a lot of issues when everything else in
the model seems correct.
The shape of the elements may also have a large impact on the accuracy of the solution. As a
general rule, avoid elements with large aspect ratios and corner angles. Aspect ratios of
elements can be calculated as a measure of element quality.
𝐿max
Aspect ratio =
𝐿min
𝐿max and 𝐿min are the largest and smallest characteristic lengths of an element, respectively.
Other measures of element quality may be calculated, and are frequently available in finite
element software suites (especially dedicated pre-processors). Fig. 7.33 shows the difference
between elements with poor quality and good quality.
Fig 7.33 Low and high quality 2D finite elements
Singular elements should also be avoided in most problems, especially when continuous
fields are expected. Examples of singular elements are shown in Fig. 7.34. These elements
will produce singularities when numerical quadrature is applied. Such elements may have
some utility in modelling fracture mechanics, where the stress at the crack tip can be
considered singular.
Fig 7.34 Singular quadrilateral elements
7.9.7. Ill-conditioned matrices

Consider the basic 1D spring system shown in Fig. 7.35. After the application of boundary
conditions, the reduced global equilibrium equation can be expressed as follows:
𝑘 −𝑘1 𝑢1 𝑃
[ 1 ] {𝑢 } = { }
−𝑘1 𝑘1 + 𝑘2 2 0
71
Fig 7.35 A simple system of 1D spring elements
The determinant of the reduced stiffness matrix is:
𝑘 −𝑘1
| 1 | = 𝑘1 𝑘2
−𝑘1 𝑘1 + 𝑘2
which means that the matrix is positive definite (𝑘1 𝑘2 > 0) and there is a non-trivial solution
to this set of equations. The equations can be expressed as:
𝑃 𝑘1
𝑢2 = 𝑢1 − 𝑢2 = 𝑢
𝑘1 𝑘1 + 𝑘2 1
These linear equations can be plotted as in Fig. 7.36, and their intersection is the solution.
Fig 7.36 Ill-conditioned and well-conditioned scenarios for the spring system
On the left, the system is ill-conditioned when the lines have very similar gradients. This
means that for very small changes in 𝑘1 , 𝑘2 , or 𝑃, there will be a very large change in the
solutions for 𝑢1 and 𝑢2. These problems are very difficult to solve computationally because
they can be easily affected by precision (rounding) error. This specific system will become
ill-conditioned when 𝑘2 ≪ 𝑘1 , thus causing the gradients of both lines to converge (Fig. 7.36
left). Conversely, if 𝑘2 ≫ 𝑘1 , then the gradients of both lines are more distinct, and the
solution will be less susceptible to any precision error (Fig. 7.36 right).
In general practice for solid mechanics, ill-conditioned problems may arise when:
Structures are insufficiently constrained;

There is a large mismatch in the stiffness (as shown by the above example); or
Mesh quality is low (elements are overly distorted).
72
Example 7.13
For an axisymmetric problem, if a 3-node triangular element is used, derive the formula for
the strain-displacement matrix, 𝐁.
Solution
Shape functions: after replacing x and y with r and z, respectively, the shape functions for a
T3 element can be used for this axisymmetric problem.
1
𝑁1 = [(𝑥 𝑦 − 𝑥3𝑦2) + (𝑦2 − 𝑦3 )𝑥 + (𝑥3 − 𝑥2 )𝑦]
2𝐴 2 3
1
→ [(𝑟 𝑧 − 𝑟3 𝑧2 ) + (𝑧2 − 𝑧3 )𝑟 + (𝑟3 − 𝑟2 )𝑧]
2𝐴 2 3
1
𝑁2 = [(𝑥 𝑦 − 𝑥1 𝑦3 ) + (𝑦3 − 𝑦1 )𝑥 + (𝑥1 − 𝑥3 )𝑦]
2𝐴 3 1
1
→ [(𝑟 𝑧 − 𝑟1 𝑧3 ) + (𝑧3 − 𝑧1 )𝑟 + (𝑟1 − 𝑟3 )𝑧]
2𝐴 3 1
1
𝑁3 = [(𝑥 𝑦 − 𝑥2 𝑦1 ) + (𝑦1 − 𝑦2 )𝑥 + (𝑥2 − 𝑥1 )𝑦]
2𝐴 1 2
1
→ [(𝑟 𝑧 − 𝑟2 𝑧1 ) + (𝑧1 − 𝑧2 )𝑟 + (𝑟2 − 𝑟1 )𝑧]
2𝐴 1 2
Alternatively, we can write:
1 1 1
𝑁1 = [𝛼1 + 𝛽1 𝑟 + 𝛾1 𝑧], 𝑁2 = [𝛼2 + 𝛽2 𝑟 + 𝛾2 𝑧], 𝑁3 = [𝛼 + 𝛽3 𝑟 + 𝛾3 𝑧]
2𝐴 2𝐴 2𝐴 3
where:
𝛼1 = 𝑟2 𝑧3 − 𝑟3 𝑧2 𝛼2 = 𝑟3 𝑧1 − 𝑟1 𝑧3 𝛼3 = 𝑟1 𝑧2 − 𝑟2 𝑧1
𝛽1 = 𝑧2 − 𝑧3 𝛽2 = 𝑧3 − 𝑧2 𝛽3 = 𝑧1 − 𝑧2
𝛾1 = 𝑟3 − 𝑟2 𝛾2 = 𝑟1 − 𝑟3 𝛾3 = 𝑟2 − 𝑟1
and
1 𝑟1 𝑧1
2𝐴 = |1 𝑟2 𝑧2 |
1 𝑟3 𝑧3
The displacement field can be expressed using these shape functions.
73
3 3
𝑢 = ∑ 𝑁𝑖 𝑢𝑖 = 𝑁1 𝑢1 + 𝑁2 𝑢2 + 𝑁3 𝑢3 𝑤 = ∑ 𝑁𝑖 𝑤𝑖 = 𝑁1 𝑤1 + 𝑁2 𝑤2 + 𝑁3 𝑤3
𝑖=1 𝑖=1
or in matrix form:
𝑢1
𝑤1
𝑢 𝑁 0 𝑁2 0 𝑁3 0 𝑢2
{ }=[ 1 ] 𝑤
𝑤 0 𝑁1 0 𝑁2 0 𝑁3 2
𝑢3
{𝑤3 }
Applying the strain-displacement relationships for axisymmetric problems, we arrive at the

equation for strain in the form of 𝛆 = 𝐁𝐝.
𝜕𝑢(𝑟, 𝑧) 𝑢(𝑟, 𝑧) 𝜕𝑤(𝑟, 𝑧) 𝜕𝑤(𝑟, 𝑧) 𝜕𝑢(𝑟, 𝑧)

𝜀𝑟𝑟 = , 𝜀𝜃𝜃 = , 𝜀𝑧𝑧 = , 𝜀𝑟𝑧 = +
𝜕𝑟 𝑟 𝜕𝑧 𝜕𝑟 𝜕𝑧
𝜕𝑢(𝑟, 𝑧) 𝜕𝑁1 𝜕𝑁2 𝜕𝑁3

0 0 0 𝑢1
𝜕𝑟 𝜕𝑟 𝜕𝑟 𝜕𝑟
𝜀𝑟𝑟 𝑢(𝑟, 𝑧) 𝑢 𝑢 𝑢 𝑤1
1 0 0 0
𝜀𝜃𝜃 𝑟 𝑟 𝑟 𝑟 𝑢2
{𝜀 } = = 𝜕𝑁1 𝜕𝑁2 𝜕𝑁3 𝑤2
𝑧𝑧 𝜕𝑤(𝑟, 𝑧) 2𝐴
𝜀𝑟𝑧 0 0 0 𝑢3
𝜕𝑧 𝜕𝑧 𝜕𝑧 𝜕𝑧
𝜕𝑤(𝑟, 𝑧) 𝜕𝑢(𝑟, 𝑧) 𝜕𝑁1 𝜕𝑁1 𝜕𝑁2 𝜕𝑁2 𝜕𝑁3 𝜕𝑁3 { 𝑤3}
+
[ 𝜕𝑧
{ 𝜕𝑟 𝜕𝑧 } 𝜕𝑟 𝜕𝑧 𝜕𝑟 𝜕𝑧 𝜕𝑟 ]
Substituting in our shape functions and displacement fields, we have:
𝑢1
𝜀𝑟𝑟 𝛽1 0 𝛽2 0 𝛽3 0 𝑤1
𝛼 𝛾 𝛼2 𝛾2 𝛼3 𝛾3
𝜀𝜃𝜃 1 1 + 𝛽1 + 1 𝑧 0 + 𝛽2 + 𝑧 0 + 𝛽3 + 𝑧 0 𝑢2
{𝜀 } = 𝑟 𝑟 𝑟 𝑟 𝑟 𝑟 𝑤2
𝑧𝑧 2𝐴 0 𝛾1 0 𝛾2 0 𝛾3
𝜀𝑟𝑧 𝑢
[ 𝛾1 𝛽1 𝛾2 𝛽2 𝛾3 𝛽3 ] {𝑤3 }
3
Therefore, the strain-displacement matrix 𝐁 is:
𝛽1 0 𝛽2 0 𝛽3 0
𝛼1 𝛾1 𝛼2 𝛾2 𝛼3 𝛾3
1 + 𝛽1 + 𝑧 0 + 𝛽2 + 𝑧 0 + 𝛽3 + 𝑧 0
𝐁= 𝑟 𝑟 𝑟 𝑟 𝑟 𝑟
2𝐴 0 𝛾1 0 𝛾2 0 𝛾3
[ 𝛾1 𝛽1 𝛾2 𝛽2 𝛾3 𝛽3 ]
Note that 𝐁 is not a constant matrix, but is a function of the co-ordinates (r and z).
74
Furthermore, the stiffness matrix for an axisymmetric element can be derived by using Hooke’s
law and substituting it into the strain energy formulation. Since 𝐁 is not a constant matrix, the
evaluation of the final integral may require the use of numerical integration techniques.
𝜎𝑟𝑟 1−𝜈 𝜈 𝜈 0 𝜀𝑟𝑟

𝜎𝜃𝜃 𝐸 𝜈 1−𝜈 𝜈 0 𝜀𝜃𝜃
𝛔 = {𝜎 } = 𝜈 𝜈 1−𝜈 0 { 𝜀 } = 𝚬𝛆
𝑧𝑧 (1 + 𝜈)(1 − 2𝜈) 1 − 2𝜈 𝑧𝑧
𝜎𝑟𝑧 𝜀𝑟𝑧
[ 0 0 0
2 ]
1 1 1
𝑈= ∫ 𝛔𝛵 𝛆 𝑑𝑉 = ∫(𝐄𝛆)𝛵 𝛆 𝑑𝑉 = ∫ 𝛆𝛵 𝐄𝛆 𝑑𝑉
2 2 2
𝑉 𝑉 𝑉
1 1
= ∫(𝚩𝐝)𝛵 𝐄(𝚩𝐝) 𝑑𝑉 = 𝐝𝑇 [ ∫ 𝚩 𝛵 𝐄𝚩 𝑑𝑉]𝐝
2 2
𝑉 𝑉
1
∴𝐤= ∫ 𝚩 𝛵 𝐄𝚩 𝑑𝑉
2
𝑉
Example 7.14
y Fy=1
Write the elemental and global equilibrium equations
D(0,1)
for the following two CST finite element model. Go Fx=1
on to determine the nodal displacements and strain in C(1,1)
2
both elements.(assume: E=1, ν=0.3, thickness=1 for a
plane stress problem) 1
Solution x
A(0,0) B(1,0)
Step 1: Connectivity of the FE model
Number of elements is 2, number of nodes is 4, total number of DOF is 8 (2 DOF per node).
Step 2: Elemental stiffness matrix calculation for a plane stress problem
[𝐤 𝑒 ] = ∫[𝐁]𝑇 [𝐄][𝐁] 𝑑𝑉 = ∫ 𝑑𝑉 ([𝐁]𝑇 [𝐄][𝐁]) = 𝑡𝐴[𝐁]𝑇 [𝐄][𝐁]

𝑉 𝑉
𝑦23 0 𝑥32
0 𝑥32 𝑦23 1 𝜈 0 𝑦23 0 𝑦31 0 𝑦12 0
𝑡 𝑦31 0 𝑥13 𝐸 𝜈 1 0
[𝐤 𝑒 ] = ( [ 1 − 𝜈 ]) [ 0 𝑥32 0 𝑥13 0 𝑥21]
4𝐴 0 𝑥13 𝑦31 1 − 𝜈 2
𝑥32 𝑦23 𝑥13 𝑦31 𝑥21 𝑦12
0 0
𝑦12 0 𝑥21 2
[ 0 𝑥21 𝑦12 ]
For element 1, the node numbering should be counter-clockwise (Node 1 = Node A, Node 2
= Node B, Node 3 = Node C)
75
E=1;
nu=0.3;
t=1;
x1=0; %Node A
y1=0; %Node A
x2=1; %Node B
y2=0; %Node B
x3=1; %Node C
y3=1; %Node C
A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;
x21=x2-x1;
x13=x1-x3;
x32=x3-x2;
y12=y1-y2;
y31=y3-y1;
y23=y2-y3;
Bmat=[y23 0 y31 0 y12 0; 0 x32 0 x13 0 x21; x32 y23 x13 y31 x21 y12];
Ke=
0.5495 0 -0.5495 0.1648 0 -0.1648
0 0.1923 0.1923 -0.1923 -0.1923 0
-0.5495 0.1923 0.7418 -0.3571 -0.1923 0.1648
0.1648 -0.1923 -0.3571 0.7418 0.1923 -0.5495
0 -0.1923 -0.1923 0.1923 0.1923 0
-0.1648 0 0.1648 -0.5495 0 0.5495
Similarly for element 2 (Node 1 = Node A, Node 2 = Node C, Node 3 = Node C)
E=1;
nu=0.3;
t=1;
x1=0; %Node A
y1=0; %Node A
x2=1; %Node B
y2=1; %Node B
x3=0; %Node C
y3=1; %Node C
A=(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2;
x21=x2-x1;
x13=x1-x3;
x32=x3-x2;
y12=y1-y2;
y31=y3-y1;
y23=y2-y3;
Bmat=[y23 0 y31 0 y12 0; 0 x32 0 x13 0 x21; x32 y23 x13 y31 x21 y12];
Ke=
0.1923 0 0 -0.1923 -0.1923 0.1923
0 0.5495 -0.1648 0 0.1648 -0.5495
0 -0.1648 0.5495 0 -0.5495 0.1648
-0.1923 0 0 0.1923 0.1923 -0.1923
-0.1923 0.1648 -0.5495 0.1923 0.7418 -0.3571
0.1923 -0.5495 0.1648 -0.1923 -0.3571 0.7418
Step 3: Write the elemental equilibrium equations
76
For element 1
(1)
𝑓𝑥𝐴
0.5495 0 −0.5495 0.1648 0 −0.1648 𝑢𝐴 𝑓𝑦𝐴
(1)
0 0.1923 0.1923 −0.1923 −0.1923 0 𝑣𝐴
(1)
−0.5495 0.1923 0.7418 −0.3571 −0.1923 0.1648 𝑢𝐵 𝑓𝑥𝐵
=
0.1648 −0.1923 −0.3571 0.7418 0.1923 −0.5495 𝑣𝐵 𝑓𝑦𝐵
(1)
0 −0.1923 −0.1923 0.1923 0.1923 0 𝑢𝐶 (1)

[−0.1648 0 0.1648 −0.5495 0 0.5495 ] { 𝑣𝐶 } 𝑓𝑥𝐶
(1)
{𝑓𝑦𝐶 }
For element 2
(2)
𝑓𝑥𝐴
0.1923 0 0 −0.1923 −0.1923 0.1923 𝑢𝐴 (2)
𝑓𝑦𝐴
0 0.5495 −0.1648 0 0.1648 −0.5495 𝑣𝐴 (2)
0 −0.1648 0.5495 0 −0.5495 0.1648 𝑢𝐶 𝑓𝑥𝐶
=
−0.1923 0 0 0.1923 0.1923 −0.1923 𝑣𝐶 𝑓𝑦𝐶
(2)
−0.1923 −0.1648 −0.5495 0.1923 0.7418 −0.3571 𝑢𝐷 (2)

[ 0.1923 −0.5495 0.1648 −0.1923 −0.3571 0.7418 ] {𝑣𝐷 } 𝑓𝑥𝐷
(2)
{𝑓𝑦𝐷 }
Step 4: Global equilibrium equations
Expand the elemental equilibrium equations according to the global displacement vector 𝐝.
Thus each elemental equilibrium equation will be in the form 𝐤𝐝 = 𝐟.
𝐝 = {𝑢𝐴 𝑣𝐴 𝑢𝐵 𝑣𝐵 𝑢𝐶 𝑣𝐶 𝑢𝐷 𝑣𝐷 }𝑇
For element 1
(1)
𝑓𝑥𝐴
0.5495 0 −0.5495 0.1648 0 −0.1648 0 0 𝑢𝐴 𝑓𝑦𝐴
(1)
0 0.1923 0.1923 −0.1923 −0.1923 0 0 0 𝑣𝐴 (1)
−0.5495 0.1923 0.7418 −0.3571 −0.1923 0.1648 0 0 𝑢𝐵 𝑓𝑥𝐵
0.1648 −0.1923 −0.3571 0.7418 0.1923 −0.5495 0 0 𝑣𝐵 𝑓𝑦𝐵
(1)
=
0 −0.1923 −0.1923 0.1923 0.1923 0 0 0 𝑢𝐶 (1)
−0.1648 0 0.1648 −0.5495 0 0.5495 0 0 𝑣𝐶 𝑓𝑥𝐶
0 0 0 0 0 0 0 0 𝑢𝐷 𝑓𝑦𝐶
(1)
[ 0 0 0 0 0 0 0 0 ] { 𝑣𝐷 } 0
{ 0 }
For element 2
(2)
𝑓𝑥𝐴
0.1923 0 0 0 0 −0.1923 −0.1923 0.1923 𝑢𝐴 (2)
𝑓𝑦𝐴
0 0.5495 0 0 −0.1648 0 0.1648 −0.5495 𝑣𝐴
𝑢𝐵 0
0 0 0 0 0 0 0 0
𝑣𝐵 0
0 0 0 0 0 0 0 0
𝑢𝐶 =
(2)
0 −0.1648 0 0 0.5495 0 −0.5495 0.1648 𝑓𝑥𝐶
−0.1923 0 0 0 0 0.1923 0.1923 −0.1923 𝑣𝐶 𝑓𝑦𝐶
(2)
−0.1923 −0.1648 0 0 −0.5495 0.1923 0.7418 −0.3571 𝑢𝐷 (2)

[ 0.1923 −0.5495 0 0 0.1648 −0.1923 −0.3571 0.7418 ] {𝑣𝐷 } 𝑓𝑥𝐷
(2)
{𝑓𝑦𝐷 }
The global equilibrium equations can be compiled by adding these expanded elemental
equations.
77
(1) (2)
𝑓𝑥𝐴 𝑓𝑥𝐴
0.7418 0 −0.5495 0.1648 0 −0.3571 −0.1923 0.1923 𝑢𝐴 (1) (2) 𝐹𝑥𝐴
𝑓𝑦𝐴 𝑓𝑦𝐴 𝐹𝑦𝐴
0.7418 0.1923 −0.1923 −0.3571 0 0.1648 −0.5495 𝑣𝐴
0.7418 −0.3571 −0.1923 0.1648 0 0 𝑢𝐵 (1)
𝑓𝑥𝐵 0 𝐹𝑥𝐵
0.7418 0.1923 −0.5495 0 0 𝑣𝐵 (1) 0 𝐹𝑦𝐵
𝑓𝑦𝐵 + 𝑓 (2)
0.7418 0 −0.5495 0.1648 𝑢𝐶 = 𝑥𝐶
=
𝐹𝑥𝐶
(1)
𝑆𝑦𝑚 0.7418 0.1923 −0.1923 𝑣𝐶 𝑓𝑥𝐶 (2)
𝑓𝑦𝐶 𝐹𝑦𝐶
0.7418 −0.3571 𝑢𝐷 (1)
𝑓𝑦𝐶 (2) 𝐹𝑥𝐷
𝑓𝑥𝐷
[ 0.7418 ] {𝑣𝐷 } 0 {𝐹𝑦𝐷 }
(2)
{ 0 } {𝑓𝑦𝐷 }
Step 5: We can go further by applying boundary conditions to solve for the displacement,
strain, and stress distributions.
0.7418 0 −0.5495 0.1648 0 −0.3571 −0.1923 0.1923 0 𝐹𝑥𝐴

0.7418 0.1923 −0.1923 −0.3571 0 0.1648 −0.5495 0 𝐹𝑦𝐴
0.7418 −0.3571 −0.1923 0.1648 0 0 0 𝐹𝑥𝐵
0.7418 0.1923 −0.5495 0 0 0 𝐹
= 𝑦𝐵
0.7418 0 −0.5495 0.1648 𝑢𝐶 1
𝑆𝑦𝑚 0.7418 0.1923 −0.1923 𝑣𝐶 1
0.7418 −0.3571 0 𝐹𝑥𝐷
[ 0.7418 ] { 0 } {𝐹𝑦𝐷 }
Step 6: Reducing the equation set to find non-zero displacements, we arrive at a two equation
set.
0.7418 0 𝑢𝐶 1 𝑢𝐶 1⁄0.7428 1.3463

[ ] {𝑣 } = { } {𝑣 } = { }={ }
0 0.7418 𝐶 1 𝐶 1⁄0.7428 1.3463
Now we have solved the global displacement vector, and the displacement distribution can be
plotted, given we have CST elements (linear displacement distributions).
𝐝 = {0 0 0 0 1.3463 1.3463 0 0}𝑇
u v
A(0,0) A(0,0)
x B(1,0) x B(1,0)
D(0,1) D(0,1)
C(1,1) y C(1,1) y
78
Step 7: The strains can be calculated given the displacement vectors.
𝑢1
𝑣1
𝜀𝑥𝑥
1 𝑦2 − 𝑦3 0 𝑦3 − 𝑦1 0 𝑦1 − 𝑦2 0
𝑢2
𝜀
{ 𝑦𝑦 } = 𝐁𝐝 = [ 0 𝑥3 − 𝑥2 0 𝑥1 − 𝑥3 0 𝑥2 − 𝑥1 ] 𝑣
𝜀𝑥𝑦 2𝐴 𝑥 − 𝑥 𝑦2 − 𝑦3 𝑥1 − 𝑥3 𝑦3 − 𝑦1 𝑥2 − 𝑥1 𝑦1 − 𝑦2 𝑢2
3 2
3
{𝑣3 }
In element 1, the node numbering should be counter-clockwise (Node 1 = Node A, Node 2 =

Node B, Node 3 = Node C)
(1) 1
𝜀𝑥𝑥 = [(𝑦2 − 𝑦3 )𝑢1 + (𝑦3 − 𝑦1 )𝑢2 + (𝑦1 − 𝑦2 )𝑢3 ]
2𝐴
1
= [(𝑦𝐵 − 𝑦𝐶 )𝑢𝐴 + (𝑦𝐶 − 𝑦𝐴 )𝑢𝐵 + (𝑦𝐴 − 𝑦𝐵 )𝑢𝐶 ]
2𝐴
1
= [(0 − 1)×(0) + (1 − 0)×(0) + (0 − 0)×(1.3463)] = 0
2×0.5
(1) 1
𝜀𝑦𝑦 = [(𝑥 − 𝑥2 )𝑣1 + (𝑥1 − 𝑥3 )𝑣2 + (𝑥2 − 𝑥1 )𝑣3 ]
2𝐴 3
1
= [(𝑥 − 𝑥𝐵 )𝑣𝐴 + (𝑥𝐴 − 𝑥𝐶 )𝑣𝐵 + (𝑥𝐵 − 𝑥𝐴 )𝑣𝐶 ]
2𝐴 𝐶
1
= [(1 − 1)×(0) + (0 − 1)×(0) + (1 − 0)×(1.3463)] = 1.3463
2×0.5
(1) 1
𝜀𝑥𝑦 = [(𝑥 − 𝑥2 )𝑢1 + (𝑦2 − 𝑦3 )𝑣1 + (𝑥1 − 𝑥3 )𝑢2 + (𝑦3 − 𝑦1 )𝑣2 + (𝑥2 − 𝑥1 )𝑢3 + (𝑦1 − 𝑦2 )𝑣3 ]
2𝐴 3
1
= [(𝑥 − 𝑥𝐵 )𝑢𝐴 + (𝑦𝐵 − 𝑦𝐶 )𝑣𝐴 + (𝑥𝐴 − 𝑥𝐶 )𝑢𝐵 + (𝑦𝐶 − 𝑦𝐴 )𝑣𝐵 + (𝑥𝐵 − 𝑥𝐴 )𝑢𝐶 + (𝑦𝐴 − 𝑦𝐵 )𝑣𝐶 ]
2𝐴 𝐶
1
= [(1 − 1)0 + (0 − 1)0 + (0 − 1)0 + (1 − 0)0 + (1 − 0)1.3463 + (0 − 0)1] = 1.3463
2×0.5
In element 2, the node numbering should be counter-clockwise (Node 1 = Node A, Node 2 =

Node C, Node 3 = Node D)
(2) 1
𝜀𝑥𝑥 = [(𝑦 − 𝑦3 )𝑢1 + (𝑦3 − 𝑦1 )𝑢2 + (𝑦1 − 𝑦2 )𝑢3 ]
2𝐴 2
1
= [(𝑦 − 𝑦𝐷 )𝑢𝐴 + (𝑦𝐷 − 𝑦𝐴 )𝑢𝐶 + (𝑦𝐴 − 𝑦𝐶 )𝑢𝐷 ]
2𝐴 𝐶
1
= [(1 − 1)×(0) + (1 − 0)×(1.3463) + (0 − 1)×(0)] = 1.3463
2×0.5
(2) 1
𝜀𝑦𝑦 = [(𝑥 − 𝑥2 )𝑣1 + (𝑥1 − 𝑥3 )𝑣2 + (𝑥2 − 𝑥1 )𝑣3 ]
2𝐴 3
1
= [(𝑥 − 𝑥𝐶 )𝑣𝐴 + (𝑥𝐴 − 𝑥𝐷 )𝑣𝐶 + (𝑥𝐶 − 𝑥𝐴 )𝑣𝐷 ]
2𝐴 𝐷
1
= [(0 − 1)×(0) + (0 − 0)×(1.3463) + (1 − 0)×(0)] = 0
2×0.5
(2) 1
𝜀𝑥𝑦 = [(𝑥 − 𝑥2 )𝑢1 + (𝑦2 − 𝑦3 )𝑣1 + (𝑥1 − 𝑥3 )𝑢2 + (𝑦3 − 𝑦1 )𝑣2 + (𝑥2 − 𝑥1 )𝑢3 + (𝑦1 − 𝑦2 )𝑣3 ]
2𝐴 3
1
= [(𝑥 − 𝑥𝐶 )𝑢𝐴 + (𝑦𝐶 − 𝑦𝐷 )𝑣𝐴 + (𝑥𝐴 − 𝑥𝐷 )𝑢𝐶 + (𝑦𝐷 − 𝑦𝐴 )𝑣𝐶 + (𝑥𝐶 − 𝑥𝐴 )𝑢𝐷 + (𝑦𝐴 − 𝑦𝐶 )𝑣𝐷 ]
2𝐴 𝐷
1
= [(0 − 1)0 + (1 − 1)0 + (0 − 0)1.3463 + (1 − 0)1.3463 + (1 − 0)0 + (0 − 1)0]
2×0.5
= 1.3463
79
ε ε
Step 8: Calculate and plot the elemental stresses using Hooke’s law.
𝜎𝑥𝑥 1 𝜈 0 𝜀𝑥𝑥
𝐸 𝜈 1 0
{𝜎𝑦𝑦 } = 𝚬𝛆 = 2[
𝜀
1 − 𝜈 ] { 𝑦𝑦 }
𝜎𝑥𝑦 1 − 𝜈 𝜀𝑥𝑦
0 0
2
1 0.3 0 𝜀𝑥𝑥 𝜀𝑥𝑥
1 1.098 0.3294 0
0.3 1 0 𝜀𝑦𝑦 } = [0.3294
= [ 1 − 0.3 ] { 1.098 0 ] {𝜀𝑦𝑦 }
1 − 0.32 𝜀𝑥𝑦
0 0 0 0 0.3843 𝜀𝑥𝑦
2
The stress in element 1

𝜎𝑥𝑥 1.098 0.3294 0 0 0.4494
{𝜎𝑦𝑦 } = [0.3294 1.098 0 ] {1.3643} = {1.4980}
𝜎𝑥𝑦 0 0 0.3843 1.3643 0.5243
The stress in element 2
𝜎𝑥𝑥 1.098 0.3294 0 1.3643 1.4980

𝜎
{ 𝑦𝑦 } = [0.3294 1.098 0 ] { 0 } = {0.4494}
𝜎𝑥𝑦 0 0 0.3843 1.3643 0.5243
σ σ
80

Chapter 7 - Finite Element Method For Solid Mechanics (2016!10!10)

Uploaded by

Copyright:

Available Formats

Chapter 7 - Finite Element Method For Solid Mechanics (2016!10!10)

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Chapter 7 - Finite Element Method For Solid Mechanics (2016!10!10)

Uploaded by

Copyright:

Available Formats

MECH3361/9361 Semester 2, 2016

FINITE ELEMENT METHOD FOR SOLID

Fig 7.1 Approximation of a circular area using triangles

Fig 7.2 Different types of elements

1. Divide structure into pieces (elements with nodes)

Considering the aforementioned advantages and disadvantages of each approach, a full

7.1.3. Brief history of FEA

1943: Courant (Variational methods)

1960s: NASTRAN = NASA STRuctural ANalysis

7.2. 1D spring elements

7.2.1. Single spring system

Fig 7.3 Single spring element

Two nodes: i and j

Fig 7.4 The spring force-displacement relationship

This can also be expressed in matrix form:

k is the elemental stiffness matrix

7.2.2. Multiple spring system

Fig 7.5 Multiple spring system

7.2.3. Boundary and load conditions

The global equilibrium equation can be written as:

The reaction force can be determined from the second equation.

Some remarks can be made from this fundamental example:

Step 1: Elemental stiffness matrices

100 −100 200 −200 100 −100

Thus the global stiffness matrix:

which is symmetric and banded.

Step 4: Global equilibrium equation

Step 5: Apply boundary condition, 𝑢1 = 𝑢4 = 0

Step 6: Solve the linear equation set to obtain unknown displacements

Step 7: Determine reaction forces from the “dropped” equations

𝐹1 = −100𝑢2 = −100×2 = −200𝑁

𝐹4 = −100𝑢3 = −100×3 = −300𝑁

Step 1: Construct an element connectivity table (4 elements and 5 nodes)

Element number Local node i (1) Local node j (2)

Step 2: Elemental stiffness matrices

𝑘 −𝑘1 𝑘 −𝑘2 𝑘 −𝑘3 𝑘 −𝑘4

In expanded form, the elemental stiffness matrices are as follows.

Step 4: Superposition to form the global stiffness matrix

7.3. 1D bar elements

7.3.1. Direct stiffness formulation

Fig 7.7 Single bar element

Two nodes: i and j

When 𝑥 = 0, 𝑢(𝑥 = 0) = 𝑢𝑖 ; When 𝑥 = 𝐿, 𝑢(𝑥 = 𝐿) = 𝑢𝑗

Stress can be found from Hooke’s law.

Step 4: The equilibrium equation also follows the same form.

From this formulation, some remarks can be made.

7.3.2. Energy formulation

Displacement can be expressed in terms of these shape functions.

𝑢(𝑥) = 𝑢(𝜉) = 𝑁𝑖 (𝜉)𝑢𝑖 + 𝑁𝑗 (𝜉)𝑢𝑗

When 𝜉 = 0, 𝑢(𝜉 = 0) = 𝑁𝑖 (0)𝑢𝑖 + 𝑁𝑗 (0) = 𝑢𝑖

Step 3: Stress can then be calculated using Hooke’s law.

Work done for the nodal forces is:

Equating these two quantities (𝑈 = 𝑊):

A EA=2 B EA=1 C EA=2 D x

Elements: AB = (1), BC = (2), CD = (3)

Nodes: A=1, B=2, C=3, D=4.

Step 2: Elemental stiffness matrices

Therefore the global stiffness matrix is: