Lectures in FEA

(Fundamentals)

Somenath Mukherjee Gangan Prathap

Scientist E2, Director,

Structural Technologies Division, National Institute of Science and
National Aerospace Laboratories (NAL), Information Resources (NISCAIR),
Bangalore. New Delhi.

1
 What everybody knows…
• The Finite Element Method is a Numerical tool to solve
differential equations (approximately) through
discretization.
• Steps of FEM
(1) Discretisation of the continuum into elements.
(2) Element formulation using approximate solution for the differential
equation.
(3) Element assembly.
(4) Solution of the system of equations, incorporating boundary
conditions.

Domain boundary

element
• Converge to more accurate results by using more
elements (i.e. through finer discretization).
2
 Element formulation
Using approximate functions and variational principles for
equilibrium, we establish {Fnodal
e
}= [K e ]{δ e }
{F }= {F
e
nodal
e
applied }+ {R e
Connectivity }
N
Fi _ node = k ij .δ j N = Total element DOF
j =1

kij=Force required at d.o.f. i for unit displacement for d.o.f. j,

with all other d.o.f. locked.

{F }= {F }
N
Assembly of elements: Global Force : G e

e =1

[K ] = [K ]
N
G e
Global Stiffness :
e =1

Solve for displacements: {F }= [K ]{δ }

G G G

Boundary conditions satisfied

Estimate element Strains/Stresses from nodal displacements

3
A machine system discretized into many finite elements

4
FEM - a tool for analysis through discretization

5
The Greeks had conceived the idea of realization through
discretization(!) . How π was discovered.

b D

Perimeter nb
Limn→∞ = Limn→∞ = 3.14159.... = π
Diameter D
How ?
nb
Limn→∞ = Limn→∞ [n. sin(π / n)] = π
D 6
 Lecture 1
Finite Element Analysis using
simple Elements

Somenath Mukherjee Gangan Prathap

Scientist E2, Director,

Structural Technologies Division, National Institute of Science and
National Aerospace Laboratories (NAL), Information Resources (NISCAIR),
Bangalore. New Delhi.

7
 Lecture 1
Finite Element Analysis using
simple Elements
Chapters
1. Analysis of a system of assembly of springs.
2. Simple bar element.
3. Simple Classical Euler beam element.
4. Transformation of co-ordinates

8
 Lecture 1

Chapter 1
Analysis of a system of assembly of
springs.

9
 1.1 A single isolated spring element
k
1 2
u1,F1 /\/\/\/\/\/\ u2,F2
Convention: Right hand direction is positive for forces/displacements.
Equilibrium :
F1 = ku1 − ku2
F2 = −ku1 + ku2
i.e.
F1 k −k u1
= or {F e } = [ K e ]{u e } (1.1)
F2 −k k u2

k −k
= [ K e ] = spring element stiffness matrix
−k k
10
 1.2 A connected spring element
k
Adjacent element 1 2 Adjacent element
or support /\/\/\/\/\/\ or support
F1( applied ) F2 ( applied )

F1( applied )
k
F2 ( applied )
/\/\/\/\/\/\
R1 R2
− R1 − R2

Convention: Right hand direction is positive for forces/displacements.

F1 = F1( applied ) + R1
F2 = F2 ( applied ) + R2 R1, R2 are end reactions
i.e.
F1( applied ) + R1 k −k u1
= or {F e } + {R e } = [ K e ]{u e } (1.2)
F2 ( applied ) + R2 −k k u2 11
 1.3 An assembly of spring elements
Fixed k1 k2 k3 4 R4
end; R 1 2 3
1 /\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\
u1=0 u4=0
F1( applied ) F2 ( applied ) F3( applied ) F4 ( applied )
Assembly by adding corresponding parameters of same degree of freedom
Assembly :
F1( applied ) + R1 = ? k1 − k1 0 0 u1 = 0
F2 ( applied ) − k1 k1 + k 2 − k2 0 u2 = ?
=
F3( applied ) 0 − k2 k 2 + k3 − k3 u3 = ?
F4 ( applied ) + R4 = ? 0 0 k3 k3 u4 = 0

{F G } + {R G } = [ K G ]{u G } (1.3)
Note: For every degree of freedom, you know either the displacement, or the
force. For every kinematic (displacement) boundary condition, there is an
associated unknown force (reaction). 12
 k1 k2 k3 4
R1 1 2 3 R4
/\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\ u4=a4
u1=a1
F1( applied ) F2 ( applied ) F3( applied ) F4 ( applied )

How do we incorporate arbitrary boundary conditions ?

u1=a1, u4=a4
F1( applied ) + R1 = ? k1 − k1 0 0 u1 = a1
F2( applied ) − k1 k1 + k 2 − k2 0 u2 = ?
=
F3( applied ) 0 − k2 k 2 + k3 − k3 u3 = ?
F4 ( applied ) + R4 = ? 0 0 k3 k3 u 4 = a4

Solve
F2 ( applied ) − (−k1 )u1 k1 + k 2 − k2 u2
=
F3( applied ) − (− k3)u4 − k2 k 2 + k3 u3
or
{Fnew} = [ K new ]{uunknown } u1 = a1 , u 4 = a4 (1.4)
13
 k1 k2 k3 4
R1 1 2 3 R4
/\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\ u4=a4
u1=a1
F1( applied ) F2 ( applied ) F3( applied ) F4 ( applied )

Solving for displacements Boundary Reactions

u1 = a1, u 2 = a4
{R G } = [ K G ]{u G } − {F G }
Method 1 :
R1 u1 F1( applied )
F2 ( applied ) − (−k1 )u1 k +k − k2 u2
= 1 2 0 u2 F2 ( applied )
F3( applied ) − (−k3 )u4 − k2 k 2 + k3 u3 = [K G ] −
0 u3 F3 ( applied )
R4 u4 F4 ( applied )
Method 2:
a1 1 0 0 0 u1 (1.5)
F2 ( applied ) − k1 k1 + k 2 − k2 0 u2
=
F3( applied ) 0 − k2 k 2 + k3 − k3 u3 R1 = k1u1 − k1u2 − F1( applied )
a4 0 0 0 1 u4
R4 = −k3u3 + k3u4 − F4 ( applied )
14
 k1 k2 k3 4
R1 1 2 3 R4
/\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\ u4=0
u1=0
F2 ( applied ) = 10kN F3( applied ) = 20kN
Example 1:
Solve for joint displacements and fixed support reactions for the spring
assembly under the given loading. Given data:
k1=1200 kN/m, k2=1800 kN/m, k3=1500 kN/m
F2=10 kN, F3= 20 kN
10 3000 − 1800 u2
=
Use equation (1.4) and solve: 20 − 1800 3300 u3
Displacements: u2=0.0104 m, u3= 0.0117 m
R1 = 0 − 1200u2 − 0
Use equation (1.5) : R4 = −1500u3 + 0 − 0

Support Reactions: R1= -12.43 kN ( ) , R4= -17.57 kN ( )

Note: Applied forces and displacements are rightwards (positive). Support
reactions are leftwards (negative); opposite to the direction shown in the figure.
Equilibrium: F2+F3+R1+R4=0 15
 k1 k2 k3 4
R1 1 2 3 R4
/\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\
u1=0
u4=0.05
Example 2:
Solve for joint displacements and support reactions for the spring assembly if
the left end is fixed and the right end is pulled out by u4=0.05 m. No other
forces are applied.
k1=1200 kN/m, k2=1800 kN/m, k3=1500 kN/m
0 3000 − 1800 u2
Use equation (1.4) and solve: =
1500.u4 − 1800 3300 u3

Displacements: u2=0.0203 m, u3= 0.0338 m

R1 = −1200u2
Use equation (1.5) :
R4 = −1500u3 + 1500u4
Support Reactions: R1= -24.32 kN ( ) , R4=24.32 kN ( )
Note: One should apply a force of 24.32 kN rightward at the right end to
achieve an end displacement of u4=0.05 m rightward. 16
 Lecture 1

Chapter 2
Simple bar element

17
 2.1 Classical formulation of a bar under axial distributed
loading intensity q(x)
x,u(x)
L
Fixed end q
d du
Governing differential equation: − EA = q (x )
dx dx (2.1)
du
Strain: ε=
dx
du
Stress: σ = Eε = E
dx
du (2.2)
Stress resultant : P = σA = EA
dx
L 2
1 du (2.3)
Strain Energy: U= EA dx
20 dx
18
2.2 A simple bar element as an equivalent spring element
Le
u1,F1 u2,F2
EAe=Elastic rigidity
1 /\/\/\/\/\/\ 2
u1,F1 u2,F2
k=EAe/Le
E=Elastic Modulus Ae=Constant section area
Equilibrium :
F1 EAe / Le − EAe / Le u1 e u1
= = [K ]
F2 e
− EA / L e e
EA / L e u2 u2 (2.4)
[ K e ] = bar element stiffness matrix

For connected element: {F e } + {R e } = [ K e ]{u e } (2.5)

19
 2.3 An assembly of bar elements as spring elements
1 2 3 4 R4
R1
L1 L3
F1( applied ) F2 ( applied ) L2 F3( applied ) F4 ( applied )
Fixed
k k2 k3 R4
end; R 1 1 2 3 4
1 /\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\
u1=0 u4=0
F1( applied ) F2 ( applied ) F3( applied ) F4 ( applied )
Equivalent spring stiffnesses k1=EA1/L1 k2=EA2/L2 k3=EA3/L3
Assembly :
F1( applied ) + R1 = ? k1 − k1 0 0 u1 = 0
F2 ( applied ) − k1 k1 + k 2 − k2 0 u2 = ?
=
F3( applied ) 0 − k2 k 2 + k3 − k3 u3 = ?
F4 ( applied ) + R4 = ? 0 0 k3 k3 u4 = 0

{F G } + {R G } = [ K G ]{u G } (2.6)
20
 R1 1 L1 2 L2 3 L3 4
R4

Fixed F2 ( applied ) = 1(kN ) F3( applied ) = 2(kN )

Fixed
end; end;
k1 k2 k3
u1=0 R 1 2 3 4 R u4=0
1 4
/\/\/\/\/\/\ /\/\/\/\/\/\ /\/\/\/\/\/\

Example 3: 1(kN ) 2(kN )

Solve for joint displacements and fixed support reactions for the bar
assembly under the given loading. F2=1 kN, F3= 2 kN.
Young’s Modulus E=200 GPa,
A1=6(10-6)m2, A2=9(10-6)m2, A3=7.5(10-6)m2, L1=L2=L3= 0.1m.

Equivalent spring stiffnesses

k1=EA1/L1= 12000 kN/m , k2=EA2/L2=18000 kN/m k3=EA3/L3 =15000 kN/m
Solution as in Example 1 (for springs)

u2=0.000104 m, u3= 0.000117 m R1= -1.243 kN ( ) , R4= -1.757 kN ( )

21
 2.4 Strain and stress resultant in bar element
Element Strain
h u2 − u1 −1 1 u1
ε = e
=
L Le Le u2
{ }
ε h = [ B] δ e (2.7)
[ B] = Strain − displacement matrix
Le
u1 u2
Element Stress Resultant
{P }= EA .[B]{δ }
e e e EAe=Elastic rigidity
In general
{σ }= [D][ B]{δ }
h e (2.8)
[ D ] = Element elastic rigidity matrix
22
 2.5 Analysis of a bar under u.d.l. q using FEM
Fixed end; L
u1=0
q
L e=L/3 L e L e
Discretization
R 1 = −q(3Le ) 1 2 3 4
= −qL
e
F1( applied ) = qL / 2 F2 ( applied ) = qLe F3( applied ) = qLe F e
4 ( applied ) = qL / 2

THE STRAIN VARIATION DIAGRAM

Analytical h
Strain ε FEA Strain ε
(linear) (constant)
qL/(EA)
qLe/(2EA)

FEA gives the best-fit to the analytical strain. Why ? 23

 h
p

h2=b2+p2

Pythagoras
24
 2.6 FEA as best fit ?
The best fit strain vector is the orthogonal projection of
the analytical strain vector, onto a given subspace.

A best fit satisfies the following

Theorem (Pythagoras)
2 2 2
ε −ε = ε −ε
apostiori analysis of FEA results shows that
FEA strain matches the best-fit strain at the element level.
ε = Best − fit
εh =ε
Thus h 2 2 h 2
ε −ε = ε −ε
(2.9)
i.e. The Energy of the Error= Error of the Energies
We need a paradigm to show -
How FEA turns out as the best fit to the analytical solution.25
 The Bilinear Symmetric Form and the inner product

• Element Strain Vector by FEA: {ε } = [B]{δ }

h e

• Analytical element strain vector : {ε }

• Bilinear form and Inner product definition :

h e
a (u , u ) = {ε } [D]{ε }dx =< ε
h T h
,ε >
e

a(u , u ) =
h h e
e
{ε } [D]{ε }dx =< ε
h T h h
, ε >= ε
h h 2

Here [D] is the element rigidity matrix

26
 Lecture 1

Chapter 3
Simple Classical Euler beam
element

27
Galileo was the first person to Euler derived the equation for
perform experiments on a the thin beam
cantilever beam 28
 3.1 Classical formulation of a the Euler beam under
transverse distributed loading intensity q(x)
R q(x)
y x
M
z,w V
M d z z M V+dV
M+dM
dx
NA Equilibrium:
dM dV d 2M
V = , q=− − 2
=q
dx dx dx
du ( R + z )dθ − Rdθ z
ε= = =
Strain: dx R.dθ R

Stress: σ = Eε = E
z
R
Governing differential equation:
d2 d 2w
NA passes through section centroid:
2
EI 2
=q
P = σ .dA =
E
RA
z.dA = 0 z=0 dx dx
A
(3.1)
Bending moment and M = σ .z.dA =
E 2
z .dA =
EI
σ =E
z M
= z I = z 2 .dA
bending stress: A
RA R R I A
29
 3.2 Basic formulation, properties and use of
the Euler beam element (using direct method)
w1,F1 w2,F2
EIe=Elastic
bending rigidity Le

1,M1
F1 k11 k12 k13 k14 w1
EIe 2,M2
M1 k 21 k 22 k 23 k 24 θ1
=
F2 k 31 k 32 k 33 k 34 w2 {F e applied } + {R e } = [ K e ]{δ e } (3.2)
M2 k 41 k 42 k 43 k 44 θ2

e e3 e e2 e e3 e e2
12( EI / L ) 6( EI / L ) − 12( EI / L ) 6( EI / L )
2 2
[K ] =
e 6( EI e / Le )
3
4( EI e / Le )
2
− 6( EI e / Le )
3
2( EI e / Le )
2 (3.3)
− 12( EI e / Le ) − 6( EI e / Le ) 12( EI e / Le ) − 6( EI e / Le )
2 2
6( EI e / Le ) 2( EI e / Le ) − 6( EI e / Le ) 4( EI e / Le )
• The stiffness matrix is symmetric. kij=kji
• C1= - C3 and R1= - R3
• C1=(C2+C4)/Le
• R1=(R2+R4)/Le
• Det [Ke]=0 WHY ? 30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 3.3 Equivalent nodal loads for the Euler beam element
under distributed loading
How do we take care of distributed loading over the element?
We need to find equivalent nodal loads that replace the distributed loading.

46
 q

Le=L/2 Le=L/2

Fixed end reactions on element ends

-q(L/2)/2 -q(L/2)/2 -q(L/2)/2 -q(L/2)/2
q q

-q(L/2)2 /12 q(L/2)2 /12 -q(L/2)2 /12 q(L/2)2 /12

Equivalent nodal loads

q(L/2)/2 q(L/2)/2 q(L/2)/2 q(L/2)/2
-q(L/2)2 /12
q(L/2)2 /12 q(L/2)2 /12 -q(L/2)2 /12

q(L/2)/2 q(L/2) q(L/2)/2

-q(L/2)2 /12
q(L/2)2 /12 47
 3.4 Strain and stress resultant in the linear Euler beam element

Element Strain z
bending strain in fibre( z ), ε=
R
2 h
1 d w
bending curvature of neutral − axis, εh = =
R dx 2
{ }
ε h = [ B] δ e (3.4)
[ B] = Strain − displacement matrix
Element Stress Resultant
{M }= EI .[B]{δ }
e e e w1,F1 w2,F2

In general {σ }= [D][B]{δ }
h e Le
(3.5)
[ D] = Element elastic rigidity matrix
1,M1 2,M2

48
 Governing differential equation d 4 wh
EI =0 (3.6)
for nodal (point) loads in FEA : dx 4

wh = ax3 + bx 2 + cx + d
w1,F1 w2,F2
2 h Le
1 d w
εh = =
R dx 2
1,M1 2,M2

Displacement (cubic polynomial of C1 continuity):

x
wh (ξ ) = N1 (ξ ) w1 + N 2 (ξ )θ1 + N 3 (ξ ) w2 + N 4 (ξ )θ 2 ξ=
Le
wh (ξ ) = [N1 (ξ ) N 2 (ξ ) N 3 (ξ ) N 4 (ξ )] δ e = [N ] δ e { } { } (3.7)
Shape functions: N1 (ξ ) = 1 − 3ξ 2 + 2ξ 3 N 2 (ξ ) = Leξ (ξ − 1) 2
N 3 (ξ ) = 3ξ 2 − 2ξ 3 N 4 (ξ ) = Leξ 2 (ξ − 1)

Bending curvature (strain):

w1
h
ε =
d 2 wh (ξ )
dx 2
=
1 d 2 wh (ξ )
e2 dξ 2
=
1
e2
[(−6 + 12ξ ) Le (6ξ − 4) (6 − 12ξ ) Le (6ξ − 2) ] wθ
1

L L 2

ε = [ B] δ e
h
{ } (3.8) θ2
49
 An example of a best fit function to a given function
Given Quadratic Curve is Linear Best-fit to y is
2 4
y = p3 = 1+ 2ξ + ξ = P1 + 2P2 + 1 P3
* * −1 ≤ ξ * ≤ 1
3 3
2
P = 1 P = ξ * P = (3ξ * −1) are Legendre Polynomials y = 4 + 2ξ * = 4 + 2P2
1 2 3 3 3
1
Orthogonality : Pi P j dξ * = 0 for i ≠ j
−1

50
 3.5 An indication of the best-fit rule in FEA …

CANTILEVER BEAM ANALYSIS USING A SINGLE EULER BEAM

ELEMENT OF LENGTH L
Uniformly distributed loading is q per unit length.
q is such that fixed end curvature (bending strain) is qL2/2EI=4 (m-1)
q

51
CANTILEVER BEAM ANALYSIS USING TWO EULER BEAM ELEMENTS
Uniformly distributed loading is q per unit length.
q is such that fixed end curvature (bending strain) is qL2/2EI=4 (m-1)
q

Le=L/2 Le=L/2

q ( Le ) 2 Bending strain
12 EI Analytical (Curvature)
(quadratic)
FEM q ( Le ) 2
(linear) 12 EI
4

Shear Force
qL
52
 Lecture 1

Chapter 4
Transformation of co-ordinates

53
4.1 What is meant by transformation of co-ordinates ?
Consider a vector A in the plane, being y A
observed from two reference frames y* a2
S and S*
A = a1i + a 2 j = a1*i * + a 2* j* a2* a1* x*

A = [i j]
a1
a2
= i [ * *
j ]aa1*
*
a1 x
2 α
A = [i [
j]{ A} = i * ]
j* { A* }
Invariance of magnitude of a
Using trigonometry:
vector upon transformation
a1* cos α sin α a1 { A* }T { A* } = { A}T { A}
=
a2* − sin α cos α a2
{ A}T [T ]T [T ]{ A} = { A}T { A}
{ A*} = [T ]{ A} (4.1) [T ]T [T ] = [ I ]

T −1
The transformation matrix is orthogonal [T ] = [T ] 54 (4.2)
 4.2 Transformation of element matrices
It is often required to transform the equilibrium equations of
an element from its own local co-ordinates to a global co-ordinate
system, common to all elements of the domain.
{F e L } = [ K e L ] δ e L{ } in local co-ordinates
[Te] is the element transformation matrix
Local (L)
{F e L } = [T e ]{F e G } {δ e L } = [T e ]{δ e G }
[T e ]{F e G } = [ K e L ][T e ]{δ e G } α
Pre-multiply both sides by [Te]T Global (G)
[T e ]T [T e ]{F e G } = [T e ]T [ K e L ][T e ]{δ e G }
e e T e e e
Transformation rule
{F G } = [T ] [K L ][T ]{δ G}
[ K eG ] = [T e ]T [ K e L ][T e ]
because [T e ]T = [T e ]−1
{F eG } = [T e ]−1{F e L } = [T e ]T {F e L }
{F e G } = [ K e G ]]{δ e G } in global co-ordinates 55 (4.3)
 4.3 Transformation of the bar (truss) element
(in the plane)
v
{F e L } = [ K e L ] δ e L{ } in local co-ordinates Local (L) uL
F e1, L 1 0 1 0 u e1, L 2
0 EAe 0 0 0 0 v e1, L 1 α
=
F e 2, L Le − 1 0 − 1 0 u e 2, L u
[Te] : element
0 0 0 0 0 v e1, L transformation matrix Global (G)
Transformation rule cos α sin α 0 0
e e T e e − sin α cos α 0 0
[K G ] = [T ] [K L ][T ] [T e ] =
0 0 cos α sin α
{F eG } = [T e ]−1{F e L } = [T e ]T {F e L } 0 0 − sin α cos α

u e1, L u e1,G
{F e G } = [ K e G ]]{δ e G } v e1, L e
e v 1,G
= [T ] e
e
u 2, L u 2,G (4.4)
in global co-ordinates e T e −1 v e 2, L v e 2,G
[T ] = [T ] 56
4.4 Transformation of the frame element (in the plane)
{F e L } = [ K e L ] δ e L { } in local co-ordinates v
Fx1, L EAe / Le 0
e e3
0
e e2
− EAe / Le 0
e e3 e
0
e2
u1, L wL Local (L) uL
Fy1, L 0 12( EI / L ) 6( EI / L ) 0 − 12( EI / L ) 6( EI / L ) w1, L
M1, L
=
0 e
6( EI / L ) e2 e
4( EI / L ) e
0 e
− 6( EI / L ) e2
2( EI e / Le ) θ1, L 2
Fx 2, L
Fy 2, L
− EAe / Le
0
0
e e3
− 12( EI / L ) − 6( EI / L )
0
e e2
EA / Le
e

0
0
e
12( EI / L ) e3
0
e
u 2, L
e2
− 6( EI / L ) w2, L
1 α
e2 θ 2, L
u
M 2, L e e e 2
0 6( EI / L ) 2( EI / L ) 0 − 6( EI e / Le ) 4( EI e / Le )

Global (G)
[Te] : element transformation matrix
Transformation rule cos α sin α 0 0 0 0
− sin α cos α
[ K eG ] = [T e ]T [ K e L ][T e ] 0 0
0
1
0
0
0
0
0
0
[T e ] =
e e −1 e e T e cos α sin α
G } = [T L } = [T
0 0 0 0
{F ] {F ] {F L}
0 0 0 − sin α cos α 0
0 0 0 0 0 1
u1, L u1,G
e e e
{F G } = [ K G ]]{δ G } w1, L v1,G
θ1, L θ1,G
= [T e ]
in global co-ordinates u 2, L u 2, G
(4.5)
e T e −1
[T ] = [T ] w2, L v2,G 57
θ 2, L θ 2, G
 4.5 Analysis of trusses and frames

Transformation matrices are obtained from element orientations

2 3
3
el=2
el=2
el=3
P el=1
1 α2 1 4
2
el=1 α3
A truss: A frame:
Assembly of axially Assembly of bending
loaded members. and axially loaded members.
58