Finite Element Analysis

of a Rotating Elastic Rod

Comparison of Linear and Quadratic Elements

and Stress Smoothing for Higher Accuracy
150
Stress (MPa)

100

50 Exact
FEM
FEM (Barlow)
0
0 100 200 300 400 500
r (mm)

Linear Elements Piecewise Constant Stress approximation.

200

150
Stress (MPa)

100
Exact
FEM
50 FEM (end nodes)
FEM (interior nodes)
FEM (Gauss points)
0
0 100 200 300 400 500
r (mm)

Quadratic Elements Piecewise Discontinuous Linear Stress ap-

proximation.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 1 of 49

Centrifugal Load on Rotating Elastic Bar

aθ = ω̇r = 0

an = −ar = ω 2r
ω

Linear varying Centrifugal force per unit length fr (r) = ρAar =

ρA(ω 2r) for constant angular speed ω.

Let’s study an interesting Finite Element Analysis (FEA)

of a rotating elastic rod using both linear and quadratic
elements, perform stress smoothing for more accu-
rate stress between connected elements, and verify
the approximate solution with an exact solution for
displacement and stress.
The uniform elastic steel rod has length L = 50 cm,
section area A = 4 mm2, rotating about a fixed pivot
with constant angular speed ω in rad/sec correspond-
ing to 120 rpm.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 2 of 49

Exact Solution
Since the centrifugal force at an element of material
in the rod fr (r) = ρAar = ρAω 2r, increases linearly
with radius r, we cannot assume the internal force and
stress is constant everywhere in the rod.
Instead, since the inertial force is linear, the internal
axial force and stress will vary quadratically with r,
and for a linear elastic material, the strain also varies
quadratically.
To find displacement, integrate the quadratic strain
resulting in a cubic displacement.

From Newton’s law of motion, the strong differential

boundary value problem for the rotating elastic rod
can be summarized as: Given the centrifugal body
load fr (r) = ρAω 2r; Find u(r) such that
d du
(EA ) = −fr
dr dr
subject to the boundary conditions,
du
u(0) = 0, N (L) = EA ∣ = 0
dr r=L
LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 3 of 49
 This equation can be solved by integrating twice and
enforcing the boundary conditions.
The axial (normal) stress is a quadratic polynomial in
position r:

N (r) ρω 2 2 2
σ(r) = = (L − r )
A 2

This is the tensile stress due to centrifugal mass iner-

tia, increasing from 0 at the free end r = L and to a
maximum at the fixed end r = 0.
N (r) ρ(ωL)2
σmax = σ(0) = =
A 2
The corresponding axial displacement is the cubic poly-
nomial:

ρω 2r
u(r) = (3L2 − r2)
6E

Since the centrifugal force was linear, the internal force,

stress, and strain were quadratic. Thus, the displace-
ment varies as a cubic polynomial in the nondimen-
sional radius r/L.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 4 of 49

Finite Element Approximation
Weak form of the Boundary Value Problem
and the Principle of Virtual Displacements
The finite element approximation satisfies the corre-
sponding weak form of the boundary value problem
equivalent to the principle of virtual work for stress
analysis.
Find u(r) such that u(0) = 0 for all virtual displace-
ments δu with the constraint δu(0) = 0 (to remove un-
known reaction forces from the problem statement),
such that
L dδu du L
∫0 EA dr = ∫ δu fr dr
dr dr 0

After finite element discretization and C 0 approxima-

tion with Lagrange interpolation (shape) function for u
and δu within elements, and after enforcing displace-
ment compatibility for assembly, we solve the alge-
braic system of discrete nodal equilibrium equations:

[K]{u} = {F }

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 5 of 49

 The left side of the integral equation gives [K]{u},
where K is the assembled global symmetric stiffness
matrix and u is a vector array of nodal displacements
(degrees-of-freedom).
The right side of the integral variational equation re-
sults in the assembled load vector {F } due to work
equivalent element nodal load vector contributions
from the centrifugal load fr .

Element stiffness matrix

The element stiffness matrix has the form,
r2
k=∫ B T EAB dr
r1

where B is the array relating axial strain to nodal dis-

placements – the elements of this matrix are deriva-
tives of the element shape functions.
For 2-node elements, these shape functions are lin-
ear polynomials.
For 3-node elements, these shape functions are
quadratic polynomials.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 6 of 49

Element nodal load vector
The element nodal load vector takes the form,
r2
f =∫ N T fr (r) dr
r 1

where fr (r) = ρAω 2r is the centrifugal load per length,

and N is an array of nodal shape functions.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 7 of 49

High-Order Finite Element Shape Functions

Lagrange interpolation functions of order one (linear), two

(quadratic), and three (cubic) used for approximation of displace-
ment within an element.

The element nodal shape functions N have the inter-

polation property of unit value at the defining node
and zero at other nodes, ensuring continuity between
connected elements.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 8 of 49

 Derivatives are discontinuous with a jump in values at
connecting nodes in a mesh of finite elements. As a
result, finite element approximations for strain ε = du
dr
and stress σ = E dr are discontinuous between con-
du

nected elements.
For 2-node linear elements, the strain ε = u′(r) is
constant, and if E is constant, the stress σ = Eε is
also constant within elements.
For 3-node quadratic elements, the strain and stress
(if uniform materials) vary linearly within elements.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 9 of 49

Element Stiffness Matrices for 2-node linear
and 3-node quadratic elements
The element stiffness matrices are obtained from in-
tegration over the element.
The uniform rod element stiffness matrix for the 2-
node linear element is:
1 −1
[k] = ke [ ]
−1 1

The uniform rod element stiffness matrix for the 3-

node quadratic element is:
⎡ 7 −8 1 ⎤
ke ⎢⎢ ⎥
⎥
[k] = ⎢−8 16 −8⎥
3 ⎢⎢ 1 −8 7 ⎥⎥
⎣ ⎦

where
EA
ke =
`e
is the axial ‘spring’ stiffness for an elastic bar member
of length `e defined in terms of the two end nodes of
the element.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 10 of 49

 In the case of the 2-node linear element `e = r2 − r1.
In the case of the 3-node quadratic element with
midpoint node 2 and end nodes 1 and 3: `e = r3 −r1.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 11 of 49

Nodal force vector equivalent to linear load
Load vector for 2-node linear element
For a linear body force distribution, fr (r) like our case
of Centrifugal load per unit length, after integration
over the element, the work equivalent nodal loads for
the 2-node linear element are:
`e 2fr (r1) + fr (r2)
{f } = { }
6 fr (r1) + 2fr (r2)

where fr (r1) and fr (r2) are the body force function

evaluations at the first and last node of the element
connected to other elements in the mesh.
In our case of Centrifugal load,

f (r1) = ρAω 2r1, f (r2) = ρAω 2r2

meω 2 2r1 + r2
{f } = { }
6 r1 + 2r2
where me = ρA`e is the mass for the uniform rod ele-
ment.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 12 of 49

Load vector for 3-node quadratic element
For a linear body force distribution, fr (r) after inte-
gration over the element, the work equivalent nodal
loads for the 3-node quadratic element are:
⎧
⎪ fr (r1) ⎫
⎪
⎪
`e ⎪ ⎪
⎪
{f } = ⎨2fr (r1) + 2f (r3)⎬
6⎪⎪
⎪ f (r )
⎪
⎪
⎪
⎩ r 3 ⎭
where fr (r1) and fr (r3) are the body force function
evaluations at the first and last (3rd) nodes of the 3-
node element connected to other elements in the mesh.
The midpoint node 2 is internal to the quadratic el-
ement and not directly connected to other elements
in the mesh (The shape function associated with this
midpoint node is zero at the end connecting nodes).
In our case of Centrifugal load,
f (r1) = ρAω 2r1, f (r3) = ρAω 2r3
giving the nodal loads,
⎧
⎪ r ⎫
⎪
me ω2 ⎪
⎪ 1 ⎪
⎪
{f } = ⎨2(r1 + r3)⎬
6 ⎪
⎪
⎪ ⎪
⎪
⎪
⎩ r3 ⎭

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 13 of 49

Finite Element Solution using a uniform mesh
of two, 2-node linear elements
Let’s find the finite element solution using just two ele-
ments by hand and compare it to the exact solution for
displacement and stress.

The assembled load vector for the mesh of two con-

nected elements:
⎧
⎪ [1] ⎫
⎪ ⎧ ⎫
⎪
⎪ f ⎪
⎪ ⎪
⎪ 2r + r ⎪
⎪
⎪ [1] 1
[2] ⎪ m ω 2 ⎪ 1 2
⎪
⎨r1 + 4r2 + r3⎬
e
⎨ f 2 + f1 ⎬ =
⎪
⎪
⎪ ⎪
⎪
⎪ 6 ⎪ ⎪
⎪ ⎪
⎪
⎪
⎪ [2]
⎪ ⎩ r + 2r ⎭
⎩ f 2 ⎭ 2 3

⎡ ke −ke 0 ⎤ ⎧ >⎫0 ⎧ 2r1 + r2 ⎫ ⎧ F̄1⎫

⎢ ⎥⎪⎪ u
 ⎪ ⎪ ⎪
⎪ ⎪
⎪ ⎪ ⎪
⎪ ⎪
⎪ ⎪
⎥ ⎪ ⎪ meω ⎪ ⎪
1

⎢ 2


⎢ −ke 2ke −ke ⎥ ⎨u2⎬ = ⎨r1 + 4r2 + r3⎬ + ⎨ 0 ⎬

⎢ ⎥⎪ ⎪
⎢ 0 −ke ke ⎥ ⎪ ⎪ ⎪
⎪ 6 ⎪⎪
⎪ r + 2r
⎪
⎪
⎪ ⎪
⎪
⎭ ⎩ ⎪
⎪ ⎪
⎪
⎣ ⎦⎩ ⎭u3 ⎩ 2 3 0 ⎭

where F̄1 is a reaction force at the fixed pivot point at

mesh node 1.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 14 of 49

 Enforcing the essential boundary condition u(0) = u1 =
0, to remove rigid body motions and using the virtual
displacement constraint δu(0) = 0, to eliminate the
first-row equation with the unknown external reaction
force F̄1, we have the reduced, symmetric positive-
definite and invertible reduced partitioned system:

2 −1 u2 meω 2 r1 + 4r2 + r3
ke [ ]{ } = { }
−1 1 u3 6 r2 + 2r3

For the two element mesh, the coordinates are r1 = 0,

r2 = ` = L/2, and r3 = 2` = L.

2 −1 u2 meω 2` 6
ke [ ]{ } = { }
−1 1 u3 6 5

Solving by taking the inverse,

u2 1 1 6 meω 2` 11 meω 2`
{ }=[ ]{ } ={ }
u3 1 2 5 6ke 16 6ke

With me = ρA`, and ke = EA/`, and ` = L/2,

u2 11 ρω 2L3
{ }={ }
u3 16 48E

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 15 of 49

 Simplifying, and including the boundary condition u1 =
0, the FEM nodal displacements are:
⎧
⎪ u1⎫⎪ ⎧
⎪ 0 ⎫⎪
⎪
⎪ ⎪ ⎪⎪ ⎪ ⎪
⎪ ρω 2L3
⎨u2⎬ = ⎨11/48⎬
⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ E
u
⎩ ⎭ ⎩
3 1/3 ⎭

Let’s check these FEM nodal displacements with the

exact cubic polynomial solution u(r) evaluated at the
node coordinates.
u(r1) = u(0) = 0
ρω 2(L/2) 11ρω 2L3
u(r2) = u(L/2) = (3L − (L/2) ) =
2 2
6E 48E
2
ρω L 2
ρω L 3
u(r3) = u(L) = (3L − L ) =
2 2
6E 3E

We just found that (even with just two elements) the

finite element nodal displacements are EXACT for this
problem!
⎧
⎪ u ⎫
⎪ ⎧
⎪ u(r ) ⎫
⎪ ⎧
⎪ 0 ⎫
⎪
⎪
⎪ ⎪1
⎪ ⎪
⎪ 1 ⎪
⎪ ⎪
⎪ ⎪
⎪ ρω 2L3
⎨u2⎬ = ⎨u(r2)⎬ = ⎨11/48⎬
⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ ⎪
⎪
⎪ E
u
⎩ ⎭ ⎩
3 u(r3 ) ⎭ ⎩ 1/3 ⎭

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 16 of 49

Superconvergence
While this might appear surprising since, in general,
the FEM is an approximate numerical solution to the
weak form of the boundary value problem for our lin-
ear elastic model or a rotating elastic rod with con-
stant angular speed, this result was to be expected
from a mathematical analysis of this type of problem.
In a rigorous mathematical proof that the finite ele-
ment solution is exact at the nodes for any load like
our linear varying Centrifugal force under the condi-
tion that the rod is uniform with constant EA, where
E is Young’s modulus and A is the section area.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 17 of 49

 Strang and Fix (yes, that famous Gilbert Strang from
MIT linear algebra YouTube videos and textbooks) at-
tribute this proof to Douglas and Dupoint.
As stated by T.J.R. Hughes in his book The Finite El-
ement Method: Linear Static and Dynamic Finite Ele-
ment Analysis,

”results of this kind, embodying exceptional accuracy

characteristics, are often referred to as superconver-
ence phenomena. However, the reader should ap-
preciate that, in more complicated situations, we will
not be able, in practice, to guarantee nodal exact-
ness. Nevertheless, as we shall see later on, weighted
residual procedures provide a framework within which
optimal accuracy properties of some sort may often
be guaranteed.”
I don’t know about you, but this is reassuring to know
FEM solutions often guarantee optimal accuracy in
some sense. How is it that most finite element books
never show this? Do the authors know this?
Exercise for the reader: Show that the finite element
nodal displacement is exact using just a single ele-
ment.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 18 of 49

Accurate Stress Recovery
While displacement is exact at the nodes and 2nd-
order accurate at points within the element, the deriva-
tives, and thus strain and stress, are in general only
1st-order accurate.
In general, derivatives of a function are an order less
accurate than the function.
However, using the Taylor series and the mean value
theorem of Calculus, it can be shown that at the mid-
points of the 2-node linear element, the derivates are
2nd-order and higher-order accurate than the 1st-order
at other points.
If the exact displacement solution u(r) is linear or quadratic,
then the derivative and stress are exact at midpoints;
this is the case when a distributed body force is con-
stant.
However, for our Centrifugal force fr = ρω 2r, the load
is linear, not constant, and the displacement is cubic.
As a result, the stress will not be exact at the midpoint,
but the stress at the midpoint is second-order accurate
compared to only first-order at other points.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 19 of 49

 As pointed out by Hughes, ’The midpoints of linear
elements are called the Barlow stress points, after
Barlow, who first noted that points of optimal accuracy
existed within elements.’

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 20 of 49

Verification of Higher-Order Stress Accuracy
at 1-point Gauss Quadrature (Barlow stress)
points
Let’s use our nodal displacement solutions to recover
stresses.
The displacement within an element is the linear inter-
polation of nodal values.
The strain at a point is the derivative (gradient) of dis-
placement at that point.
Since the displacement is linear, the derivative gives
constant strain, and if E is constant from Hooke’s law,
the stress is also constant.
due ue2 − ue1
εe = = Be ue = ,
dx ` e
u2 − ue1
σe = Eεe = Be ue = E ( )
`

The strains and stresses are piecewise constant and

discontinuous at connected mesh nodes.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 21 of 49

 For our two elements,
u2 − u1 11
σ1 = E ( ) = ρ(ωL)2
` 24
u3 − u2 5
σ2 = E ( ) = ρ(ωL)2
` 24

Let’s check these FEM stress approximations with the

exact stress solution given by the quadratic polyno-
mial solution σ(r) evaluated at the node coordinates.
For Element [1], at the midpoint r = `/2 = L/4,

ρω 2 2 15
σ(L/4) = (L − (L/4)2) = ρ(ωL)2
2 32
For Element [2], with midpoint r = 3`/2 = 3L/4,

ρω 2 2 7
σ(3L/4) = (L − (3L/4) ) = ρ(ωL)2
2
2 32
These don’t match the finite element stress values as
expected for a linear body load with cubic displace-
ment field.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 22 of 49

 Let’s check the error at the midpoints of the two ele-
ments.
At the midpoints of elements [1] and [2], the percent
error is:
σ1 1
( − 1) × 100% = − × 100% = 2.2%
σ(L/4) 45
σ2 1
( − 1) × 100% = − × 100% = 4.8%
σ(3L/4) 21

Wow, Barlow was right; this is incredible stress accu-

racy with a mesh of only two 2-node linear elements!

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 23 of 49

Nodal Stress Smoothing by Interpolation
Since the finite element stress solution is discontinu-
ous with a jump, there is no unique solution at the
connecting node.
Let’s do some stress smoothing using simple interpo-
lation; this results in a simple average of the two-element
constant solutions averaged at the connecting node.

σ1 + σ2 ρ(ωL)2
σavg = =
2 3

Let’s compare it to the exact stress evaluated at node

2.
3ρ(ωL)2
σ(r2) = σ(L/2) =
8
This gives an error,

σavg 1
( − 1) × 100% = − × 100% = 11%
σ(L/2) 9

It’s not too bad, but we need mesh refinement with

more elements for improved accuracy.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 24 of 49

 Finite Element Solution using a single 3-
node quadratic element
Now, let’s examine the accuracy using just a single 3-
node quadratic element.
⎡ 7 −8 1 ⎤ ⎧ ⎪ >⎫0
⎪ ⎧
⎪ ⎫
⎪ ⎧
⎪ ⎫
⎪
⎢
ke ⎢ ⎥ ⎪
⎪
u 1

⎪
⎪ ⎪
⎪
r1 ⎪
⎪ ⎪
⎪
F̄1 ⎪
⎪
⎥ 2


m ω
⎢ −8 16 −8 ⎥ ⎨u2⎬ = ⎨2(r1 + r3)⎬ + ⎨ 0 ⎬
e
3 ⎢⎢ ⎥⎪
⎥⎪⎪ ⎪
⎪
⎪ 6 ⎪ ⎪
⎪ ⎪
⎪
⎪ ⎪
⎪ ⎪
⎪
⎣ 1 −8 7 ⎦⎩ ⎭u 3 ⎩ r3 ⎭ ⎩ ⎪
⎪ 0 ⎭

ke 16 −8 u2 meω 2 2(r1 + r3)

[ ]{ } = { }
3 −8 7 u 3 6 r3

For the 3-node quadratic element with middle node

2, the element length is ` = r3 − r1 = L.

ke 16 −8 u2 meω 2` 2
[ ]{ } = { }
3 −8 7 u 3 6 1

Solving this, we verify that the nodal displacements

are EXACT.
u2 1 7 8 2 meω 2` 11 ρω 2L3
{ }= [ ]{ } ={ }
u3 48 8 16 1 2ke 16 48E

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 25 of 49

 We leave it as an exercise for the reader to show that
the finite element stress solution for this 3-node quadratic
element is higher-order accurate at the 2-point Gauss-
Legendre quadrature points within the element.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 26 of 49

Deciding on an optimal finite element mesh
of connected elements
0.6
Stress Gradient

0.4
= ;r! 2

0.2
d<
dr

0
0 100 200 300 400 500
r (mm)
Analytical Stress Gradient

The exact stress gradient is proportional to the linear

centrifugal load
dσ
= ρω 2r
dr
increases linearly, with the second gradient a constant.
This implies a linear gradient finite element mesh is op-
timal. However, for our studies, let’s keep things sim-
ple and only consider uniform meshes of equally spaced
nodes.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 27 of 49

Uniform Mesh of Linear Elements
Let’s do a mesh refinement study of 2, 4, and 8 linear
elements with 3, 5, and 9 mesh nodes, respectively.

Two Linear Elements

Consider a uniform mesh of a very coarse mesh of only
two 2-node linear elements.
0.3
Exact
Displacement (mm)

FEM
0.2

0.1

0
0 100 200 300 400 500
r (mm)
Piecewise Linear Displacement FEM approximation

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 28 of 49

Stress
150
Stress (MPa)

100

50 Exact
FEM
FEM (Barlow)
0
0 100 200 300 400 500
r (mm)
Piecewise Constant Stress FEM approximation.

200

150
Stress (MPa)

100

50 Exact
FEM (Nodal Avg)
0
0 100 200 300 400 500
r (mm)
FEM Nodal Stress Smoothing

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 29 of 49

Four Linear Elements
Consider a uniform mesh of a relatively coarse mesh of
four 2-node linear elements.

Displacement

0.3
Exact
Displacement (mm)

FEM
0.2

0.1

0
0 100 200 300 400 500
r (mm)
Piecewise Linear Displacement FEM approximation

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 30 of 49

Stress
Stress (MPa) 150

100

50 Exact
FEM
FEM (Barlow)
0
0 100 200 300 400 500
r (mm)
Piecewise Constant Stress FEM approximation.

200

150
Stress (MPa)

100

50 Exact
FEM (Nodal Avg)
0
0 100 200 300 400 500
r (mm)
FEM Nodal Stress Smoothing

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 31 of 49

Eight Linear Elements
Consider a uniform mesh of eight two-node linear ele-
ments.

Displacement

0.3
Exact
Displacement (mm)

FEM
0.2

0.1

0
0 100 200 300 400 500
r (mm)
Piecewise Linear Displacement FEM approximation

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 32 of 49

Stress
150
Stress (MPa)

100

50 Exact
FEM
FEM (Barlow)
0
0 100 200 300 400 500
r (mm)
Piecewise Constant Stress FEM approximation.

150
Stress (MPa)

100

50
Exact
FEM (Nodal Avg)
0
0 100 200 300 400 500
x (mm)
FEM Nodal Stress Smoothing

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 33 of 49

Results Observations
1. As expected, the FEM stress approximation is less ac-
curate than the displacement approximation.
2. The stress evaluated at the center midpoint equiva-
lent to the one-point Gauss-Legendre Points (Barlow
superconvergent points) is more accurate than other
points within the element, especially at end nodes.
3. Improved stress approximation is obtained using smooth-
ing. In this case, we used nodal interpolation of the av-
erage stress between connected elements. The nodes
at the ends are extrapolated. Another approach is to
use patch recovery for accurate nodal stress values.
4. The extrapolation at the last node tends to undershoot
or overshoot with less accuracy than the interior inter-
polated points.
5. This simple averaging is not just for plotting; It gives
a more accurate nodal stress than the constant values
from adjacent elements at connected nodes.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 34 of 49

6. A generalization of this idea of interpolating the opti-
mal Barlow stress points can be used with other stress
post-processing methods, such as global interpola-
tion and Superconvergent Patch Recovery (SPR) for ac-
curate nodal stress values.
7. Since the rod is uniform with a constant cross-section
area and material EA = constant, the FEM solution
nodal displacements are proven to be exact regard-
less of the body load, in this case, the linear varying
Centrifugal force. Stresses are shown to be higher-
order accurate compared to typical points, especially
at nodes. The exact nodal displacements are a spe-
cial case not expected in general two and three-
dimensional applications of FEA. The phenomena of
exact or higher accuracy than expected at other points
is called superconvergence.
8. The FEM approximation converges to the exact ana-
lytical solution everywhere with mesh refinement. This
can also be proved mathematically with estimates of
convergence rate.
9. Since the solution is smooth (high regularity), with a
linear stress gradient, a uniform mesh with equally spaced
nodes gave accurate solutions everywhere in the prob-
lem domain.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 35 of 49

Uniform Mesh of Quadratic Elements
Let’s do a mesh refinement study of 1, 2, and 4 quadratic
elements with the same number of 3, 5, and 9 mesh
nodes, respectively. By keeping the number of nodes
constant, we are able to make a fair comparison between
linear and quadratic elements with the same number of
degrees-of-freedom and assembled equilibrium equations
to solve.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 36 of 49

One Quadratic Element
Consider a mesh of a single 3-node quadratic element.
0.3

0.2
u(r) (mm)

0.1 Exact
FEM
FEM (end nodes)
FEM (interior nodes)
0
0 100 200 300 400 500
r (mm)

Quadratic Displacement FEM approximation

200

150
Stress (MPa)

100
Exact
FEM
50 FEM (end nodes)
FEM (interior nodes)
FEM (Gauss points)
0
0 100 200 300 400 500
r (mm)

Linear Stress FEM approximation.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 37 of 49

Two Quadratic Elements
Consider a uniform mesh of a relatively coarse mesh of
two 3-node quadratic elements.
0.3

0.2
u(r) (mm)

0.1 Exact
FEM
FEM (end nodes)
FEM (interior nodes)
0
0 100 200 300 400 500
r (mm)

Piecewise Quadratic Displacement FEM approximation

200

150
Stress (MPa)

100
Exact
FEM
50 FEM (end nodes)
FEM (interior nodes)
FEM (Gauss points)
0
0 100 200 300 400 500
r (mm)

Piecewise Linear Stress FEM approximation.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 38 of 49

Four Quadratic Elements
Consider a uniform mesh of four three-node quadratic
elements.
0.3

0.2
u(r) (mm)

0.1 Exact
FEM
FEM (end nodes)
FEM (interior nodes)
0
0 100 200 300 400 500
r (mm)

Piecewise Quadratic Displacement FEM approximation

150
Stress (MPa)

100

Exact
50 FEM
FEM (end nodes)
FEM (interior nodes)
FEM (Gauss points)
0
0 100 200 300 400 500
r (mm)

Piecewise Linear Stress FEM approximation.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 39 of 49

Observations of results for Quadratic Element
Meshes
1. Since the solution is smooth (high regularity), the quadratic
element mesh solutions are more accurate than the
linear mesh solutions for the same number of equa-
tions and unknowns.
2. The FEM piecewise constant stress approximation for
2-node linear displacement elements is less accurate
than the piecewise linear approximation of the 3-node
quadratic displacement elements.
3. The 3-node quadratic elements give the best stress
approximation at the 2-point Gauss-Legendre Quadra-
ture points as found by Barlow. The stress at these
points can be used with a smoothing technique such
as the superconvergent patch recovery (SPR) method
to recover accurate stresses at nodes.
4. All commercial software such as Ansys, Abaqus, and
Nastran use some sort of stress smoothing for nodal
stresses (something that most user software user man-
uals do not explain, likely for competitive reasons).

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 40 of 49

Comparison of Stiffness Matrix Sparsity be-
tween 2-node linear and 3-node quadratic el-
ement meshes
1. Each row in the stiffness matrix represents force equi-
librium at each node in the mesh.
2. For the linear element mesh, since the shape functions
for the nodes connect with two adjacent 2-node ele-
ments, there are three nonzero coefficients per row.
3. For the quadratic element mesh, since the shape func-
tions for connected nodes connect with two adjacent
3-node elements, there are five nonzero coefficients
per row at these connected nodes.
4. For the quadratic mesh, since the shape function for
the interior nodes is local to the element and not con-
nected with adjacent elements, there are three non-
zero coefficients per frow for these interior nodes.
5. There are more nonzero coefficients in the quadratic
element mesh, requiring additional memory storage
and compute time when solving equations.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 41 of 49

Summary and Conclusions
This detailed comparison between FEA and analytical
solutions for a rotating elastic rod illustrates the practi-
cal considerations of mesh density and element order
in achieving accurate results.
The comparison between refined finite element meshes
highlights the importance of mesh density and ele-
ment order in achieving accurate stress and displace-
ment predictions. It’s fascinating to see how nodal
interpolation (stress smoothing) can significantly en-
hance accuracy.
The insights on nodal interpolation techniques and
their impact on stress accuracy provide valuable guid-
ance for optimizing FEA models.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 42 of 49

Further observations
1. For this model problem studied, since the analytical
solution for the linear Centrifugal force has quadratic
stress and cubic displacement, a single hierarchical
basis p-version finite element or a global Rayleigh-Ritz
basis approximation for displacement with a 3rd-order
polynomial would capture the displacement and stress
exactly everywhere.
2. However, the point of this study was to examine mesh
refinement of standard linear and quadradic elements
with standard Lagrange basis shape functions used in
Commercial FEA software to give insights into the use
of these types of elements in multi-dimensional com-
plex problems.
3. The global basis Rayleigh-Ritz method can be difficult
to apply to complex geometries and does not pro-
vide a framework for general-purpose software. For
these and other compelling reasons, the Rayleigh-Ritz
method is not used in practice in industry except in
special cases.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 43 of 49

4. Analytical solutions are generally available for simple
shapes, such as one-dimensional, two-dimensional, and
three-dimensional problems with separable coordinates.
5. There are fewer nonlinear solutions to compare to,
and in general, they are not available for complex non-
linear problems. Nevertheless, a great deal can be
gained by studying how well finite element solutions
behave with mesh density relative to analytical solu-
tions where they exist.
6. Of course, in practice, we should never solve the FEA
solution with one mesh; mesh refinement and solu-
tion iterations are needed to determine mesh sensi-
tivity and convergence (grid independence) up to an
acceptable tolerance within constraints on computer
memory and processing time.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 44 of 49

 7. While this study showed quadratic elements outper-
formed linear elements, we should also note that lin-
ear elements are not always inferior to quadratic el-
ements; there are applications when linear elements
perform better, and the reasons why can be explained.
8. Contrary to simple answers, distributed loads are not
”lumped” to the nodes in an ad-hoc way; they are ju-
diciously applied to give the equivalent work of the
distributed load consistent within the finite element
approximation and virtual work principle.
9. Stress is less accurate because it is defined in terms
of the strain, which is a derivative of the displacement
function obtained after applying nodal displacement
compatibility and solving the assembled discrete nodal
equilibrium equations for nodal displacements, with
the derivatives being less accurate than the function.
10. Most of the formulations and methods used in closed
and open-source finite element codes were obtained
by the software developers reading and studying the
60 years of open literature on finite element methods
in research journal papers and research codes.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 45 of 49

11. There are valuable lessons learned that can acceler-
ate mastery of practical FEA by investing some time
in learning the fundamental principles and concepts
of FEM behind our FEA.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 46 of 49

Correction
The plots of stress are labeled as MPa;
the correct units are kPa.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 47 of 49

References
1. Thomas J.R. Hughes, The Finite Element Method: Lin-
ear Static and Dynamic Finite Element Analysis, 1987.
Republished by Dover in 2000.
2. G. Strang and G.J. Fix, An Analysis of the Finite Ele-
ment Method. 1973.
3. J. Barlow, ‘Optimal Stress Locations in Finite Element
Models’, International Journal for Numerical Methods
in Engineering, 10 (1976), 243-251.

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 48 of 49

We love the FEM
The Finite Element Method (FEM) blends elegance with
utility, rewarding mastery with the power to solve real-
world problems with FEA.

Please share this knowledge of the

FEM with others so they can also ben-
efit from its insights.
–Thanks!

LinkedIn � - Dr. Lonny Thompson, Clemson University, April 10, 2025. 49 of 49