Chapter 5: Finite Element Formulation (Energy Methods)

5.2. FEM formulation for Beams

Beams are slender members that are used for supporting transverse (perpendicular to axis)
loading. Long horizontal members used in buildings and bridges, and shafts supported in
bearings are some examples of beams. Beams with cross sections that are symmetrical with
respect to plane of loading are considered here. Figure below shows a general horizontal beam
with the possible types of loading. Another figure shows the cross section and the bending
stress distribution.
Pm
p
Mk
k
m

v
R dv
x
dx

Fig. 3.1 Loading and deflection of a horizontal beam

At point R, we have 2 dofs, deflection ‘v’ and slope dv
dx

‘p’ is the distributed load /unit length.

‘Pm’ is the point load acting at point ‘m’.
‘Mk’ is the moment of couple acting at point ‘k’.

y y
dA

y
Neutral
z x
Axis

Fig. 3.2 Beam section and stress distribution

1
Chapter 5: Finite Element Formulation (Energy Methods)

The Potential Energy expression for a beam:

From beam theory we have (for small deflections)
M 
    y     (3.1)
 I 
 M 
      y     (3.2)
E  EI 
d 2v M
and      (3.3)
dx 2 EI

The potential energy for a structural problem is given as:

  U  WP     (3.4)
where, U is Strain Energy, &
WP is the Work Potential

Expression for Strain Energy ‘U’:

The Strain Energy is due to the stresses and strains experienced by the beam.

The strain energy dU of the beam element of length ‘dx’ is,

1
dU    dAdx     (3.5)
2A
Substituting (3.1) & (3.2) in (3.5),
1  My   My  1 M2 
 dU      dA dx   2  y dA dx
2

2 A  I   EI  2  EI A 
But, I   y 2 dA is the area moment of inertia of beam cross section about the neutral axis.
A

1M2 
 dU    dx
2  EI 
Multiply and divide right hand side by EI
2
1  d 2v 
2
1 M 
i.e., dU  EI   dx  EI  2  dx     (3.6) [ from (3.3)]
2  EI  2  dx 
Then the total strain energy in the beam is given by integrating (3.6):
2
1
L
 d 2v 
U   EI  2  dx     (3.7)
2 0  dx 

Expression for Work Potential ‘WP’:

The Work Potential is due to the different types of loads (distributed loads, point loads
and point moments) acting on the beam.

L
The Work Potential due to Distributed Load ‘p’: WPD    pv dx
0

The Work Potential due to Point Load ‘Pm’: WPP   Pm vm

m

2
Chapter 5: Finite Element Formulation (Energy Methods)

 dv 
The Work Potential due to Point Moment ‘Mk’: WPM   M k  
k  dx  k
Putting all the Work Potential terms together,
 dv 
L
WP    pv dx   Pmvm   M k       (3.8)
0 m k  dx k
Substituting (3.7) & (3.8) in (3.4), for a beam, we can write the potential energy expression as
2
1
L
 d 2v  L
 dv 
   EI  2  dx   pv dx   Pm vm   M k       (3.9)
2 0  dx  0 m k  dx k
where ‘v’ represents the vertical deflection in general,
‘vm’ is vertical deflection at point ‘m’.
 dv 
  is the slope of the beam at point ‘k’
 dx  k

Finite element formulation of beam element:

Let us consider a beam divided into 3 elements as shown below

Q1 Q3 Q5 Q7

1 2 Q4 3 Q6 4 Q
Q2 1 2 3 8

Fig. 3.3a Global displacements

q1 q3
The global displacement vector is:
1 2
Q1  q2 e q4
Q 
 2 Fig. 3.3b Element displacements
Q3 
 
Q 
{Q}   4  The element displacement vector is:
Q5 
Q6   q1   v1 
 
Q7 
q   dv 
 2   dx 1 
{q}     
 
Q 
 8 v2 
q3   

q4   dv 
 dx 2 

Shape functions for beam element:
The linear and quadratic ‘1D’ element shape function represents both displacement ‘u’
and geometry ‘x’ along axial direction. But, for beam element the shape functions will be
different because, now the shape functions should represent displacement (deflection ‘v’
perpendicular to axis of beam) and they must also represent the slope ‘dv/dx’. The Hermite
Cubic polynomial functions are the best-suited shape functions to represent the displacement
(deflection ‘v’ & slope ‘dv/dx’) for beam element. They are represented by the letter ‘H’ (not
‘N’). However, the geometry of the beam ‘x’ can be represented by the shape functions [N]
used for 2 noded 1D element. Hence, the beam element formulation is called ‘Sub-Parametric
Formulation’.
3
Chapter 5: Finite Element Formulation (Energy Methods)

The Hermite Cubic polynomial is given by:

Hi = ai + bi + ci2 +di3 for i = 1, 2, 3, 4 --------(3.10)

H1 & H2 are the shape functions at node 1 for representing deflection and slope
respectively
H3 & H4 are the shape functions at node 2 for representing deflection and slope
respectively

To evaluate the coefficients ai, bi, ci, & di the following conditions must be met.

dH1 dH 2 dH 3 dH 4
H1 H2 H3 H4
d d d d
=1 1 0 0 1 0 0 0 0
=+1 0 0 0 0 1 0 0 1
a1 = 1 2 a2 = 1 4 a3 = 1 2 a4 =  4
1

b1=  3 4 , b2=  1 4 , b3 = 3 4 , b4=  1 4 ,

c1 = 0, c2 =  1 4 , c3 = 0, c4 = 1 4 ,
d1 = ¼ d2 = ¼ d3 = ¼ d4 = ¼

Substituting values of ai ,bi, ci & di in (3.10) we get,

H1 = ¼(2  3 + 3) or ¼(1  )2(2 + )
H2 = ¼(1     +  )
2 3
or ¼(1  )2(1 + )
---------(3.11)
H3 = ¼(2 + 3  3) or ¼(1 + )2(2  )
H4 = ¼( 1   + 2 + 3) or ¼(1 + )2(  )

The plot of variation of H1, H2, H3 & H4 with respect to ‘’ is shown below:
Slope = 0
Slope = 0
H1
H3

1 Slope = 0 Slope = 0 1

=1 =0 =+1 =1 =0 =+1

H2 H4
Slope = 0 Slope = 0 Slope = 1
Slope = 1

=1 =0 =+1 =1 =0 =+1

Fig. 3.4 Variation of shape functions of horizontal beam across element domain

So the deflection ‘v’ within the element can be represented by :

 dv   dv 
v  H1v1  H 2    H 3v2  H 4       (3.12)
 d 1  d 2
4
Chapter 5: Finite Element Formulation (Energy Methods)

The ‘x’ co-ordinate transform by the relation

 1    1  
x  x1    x2
 2   2 
x x  x x 
  1 2    2 1      (3.13)
 2   2 

dx  x2  x1  le
  
d  2  2
l 
 dx   e  d     (3.14)
2
where, le=x2x1is element length.

From chain rule of differentiation, we have,

dv  dv   dx 
        (3.15)
d  dx   d 
Substituting (14) in (15),
dv  le  dv 
        (3.16)
d  2   dx 
We have,
 dv 
At node 1, v1  q1 &    q2
 dx 1
-----------(3.17)
 dv 
At node 2, v1  q3 &    q4
 dx  2
Substituting (3.16) & (3.17) in (3.12)
l l
v  H1q1  e H 2 q2  H 3q3  e H 4 q4
2 2
Writing in matrix form,
 q1 
 
 le le  q2 
v   H1 H 2 H3 H4   
 2 2   q3 
q4 

v  {H }{q}     (3.18)

 le le 
where {H }   H1 H2 H3 H 4  is the Hermite Shape Function Matrix
 2 2 

5
Chapter 5: Finite Element Formulation (Energy Methods)

Derivation of Element Stiffness Matrix:

The element stiffness matrix is derived from the strain energy term of the potential energy
expression. The strain energy is given by:
2
1
L
 d 2v 
U   EI  2  dx (from eq. 3.7)
2 0  dx 
We have, from (3.16),
dv  2   dv 
   
dx  le   d  
d 2v  2   d 2v   d  4  d 2v   d 2 
            (3.19)   
dx 2  le   d  2   dx  le 2  d  2   dx le 
Substituting (3.18) in (3.19),
d 2v 4  d 2v  4 d 2  l l 
 2 2  2 H q  e H 2 q2  H 3 q3  e H 4 q4 
2  1 1
dx 2
le  d   le d   2 2 
4  d 2 H1 le d 2 H 2 d 2 H3 le d 2 H 4 

 q  q  q  q4      (3.20)
le 2  d  2 2 d 2 d 2 2 d 2
1 2 3

Eq. (3.20) can be written in two forms. Writing the first form of matrix notation,
 H1 
l 
 e
H 
d 2v 4 d2 2 2
 q1 q2 q3 q4 (14)  2  
dx 2 le 2 d  H 3 
 le 
 H4 
2 (41)
2 T
4 T d {H }
 {q}    (3.21a )
le 2 d 2
Rewriting (3.20) in the second form of matrix notation,
 q1 
q 
d v 4 d    2
2 2
le le
  2  1
H H H H 4  
le d  
2 3
(14)  q3 
2 2
dx 2 2
q4 
(41)

4 d 2 {H }
 2 {q}    (3.21b)
le d  2

2
 d 2v   d 2v   d 2v 
Now,  2    2    2 
 dx   dx   dx 
 4 d 2 {H }T   4 d 2 {H } 
  2 {q}T  2 {q}  eq.(3.21a)  eq.(3.21b)
 le d 2   le d 
2

16 T  d 2 {H }T d 2 {H } 
 {q}   {q}     (3.22)
le 4  d d 2 
2

 le le 
We have, {H }   H1 H2 H3 H4 
 2 2 

6
Chapter 5: Finite Element Formulation (Energy Methods)

d {H }  dH1 le dH 2 dH 3 le dH 4 
  
d  d 2 d d 2 d 
3 le 3 le 
  (1   2 ) (1  2  3 2 ) (1   2 ) ( 1  2  3 2 ) 
4 8 4 8 
d {H }  3
2
le 3 le 
  (1  3 )  (1  3 )      (3.23)
d 2
 2 4 2 4 
Substituting (3.22) & (3.23) in (3.7)
1
L
16  d 2{H }T d 2{H }  
U   EI  4 {q}T  {q} dx
2 0  le  d
2
d 2  
Changing the integral variable from ‘x’ to ‘’, limits from x = (0 to L) to  = 1 to +1 and
l 
putting dx   e d
2
1  8EI 1  d 2{H }T d 2{H }  
U  {q}T  3    d  {q}
 le 1  d d 2  
2
2

1
U  {q}T [ke ]{q}     (3.24)
2
where, [ke] is the element stiffness matrix of the beam element.
Thus, the element stiffness matrix is:
 12 6le  12 6le 
 6l 4le
2
 6le 2le 
2
EI
[ke ]  3  e -------(3.25)
le  12  6le 12  6le 
 2 
 6le 4le 
2
 6le 2le
Derivation of load vector:
Three types of load are acting on the beam:
(1) The distributed load ‘p’ per unit length
(2) Point load ‘Pm’ at point ‘m’ and
(3) Moment Mk at point ‘k’
Consideration of distributed load ‘p’:
(a) If ‘p’ is constant, i.e., Uniformly distributed load (UDL)
We have to convert the distributed load ‘p’ into equivalent point loads and moments at the
nodes. This can be done as follows.
In the equation for ‘’, the term corresponding to distributed load is:
1
 
L
l l
0 pv dx  1 p{H }{q} 2e d  by us ing v  {H }{q} & dx  e d 
 2 
L
 ple 1 
0 pv dx   2 1{H } d {q}


 q1 
 pl  
ple  q2 
2 2
ple ple
  e      --------(3.26)
 2 12 2 12  q3 
q4 
7

 { f e }T {q}
 Chapter 5: Finite Element Formulation (Energy Methods)

(b) If ‘p’ varies linearly, i.e., Uniformly varying load (UVL)

p2
p1

le
 7 p1le 3 p2le   3 p1le 7 p2le 
 20  20   20  20 
   

 p1le 2 p2le 2   pl 2 p l 2 
    1 e  2 e 
 20 30   30 20 
Fig 3.6 Conversion of UVL into nodal loads
 q1 
 7 p1le 3 p2le   p1le 2 p2le 2   3 p1le 7 p2le   
L
 p1le 2 p2le 2   q2 
0 pv dx   20  20   20  30   20  20    30  20   q3     (3.26)
   
q4 

 { f e }T {q}
Note: (1) In case of triangular load, p1 = 0, if the load intensity varies from zero at left end to
p2 at the right end; and p2 = 0, if the load intensity varies from zero at right end to p1 at the left
end.
(2) In case of UDL, p1 = p2. It may be noted that we get the same equations what we derived
earlier for UDL.

2. Consideration of point loads Pm and Mk:

The point loads Pm and Mk are readily taken care of by introducing nodes at their points of
application.
Form the potential energy approach we get
 = (½){Q}T[K]{Q}{Q}T{F}
where, {Q} is the global displacement vector,
[K] is the global stiffness matrix &
{F} is the global force vector.
8
Chapter 5: Finite Element Formulation (Energy Methods)

Applying principle of minimum potential energy,


0 gives {F }  [ K ]{Q}
Q

Table below shows equivalent nodal forces for various loading configurations over beam

Numerical 1: 10 kN/m 20 kN

Find the deflection and slope at the point

where the concentrated load of 20 kN is
acting on the beam shown in figure. Assume 1.5 m 1m
Young’s Modulus for beam material to be 210
GPa and the area moment of inertia of the
cross section of the beam as 2 × 10-6 m4. 20 kN
+
7.5 kN 7.5 kN
1 2 2
1 3

1.875 kN m 1.875 kN m 9
Chapter 5: Finite Element Formulation (Energy Methods)

Solution:
For element (1),

The element stiffness matrix is given as:

 12 6l1  12 6l1 
 6l 4l1
2
 6l1 2l1 
2
EI
[k ](1)  13 1  1
l1  12  6l1 12  6l1 
 2 
 6l1 4l1 
2
 6l1 2l1

E1I1 (210  109 )  (2  10 6 )

  124 .44  103 N
3
l1 1.5 3 m
1 2 3 4
 1493 .33 1119 .9  1493 .33 1119 .9  1
 1119 .9 1119 .9  1119 .9 559 .9  2
3
[k ]  10
(1)
 1493 .33  1119 .9 1493 .33  1119 .9 3
 
 1119 .9 559 .9  1119 .9 1119 .9  4
For element (2),
The element stiffness matrix is given as:
 12 6l2  12 6l2 
 6l 2
 6l2 2l2 
2
E2 I 2  2 4l2
[k ]  3
( 2)

l2  12  6l2 12  6l2 
 2 
 6l2 4l2 
2
 6l2 2l2
E2 I 2 (210  109 )  (2  10 6 )
  420  103 N
3
l2 13 m
3 4 5 6
 5040 2520  5040 2520  3
 2520 1680  2520 840  4
3
[k ]( 2)  10
 5040  2520 5040  2520  5
 
 2520 840  2520 1680  6

Assembling [k](1) and [k](2) to get the global stiffness matrix [K]
(Size of the Global Stiffness Matrix) = (Total no. of nodes)×(No. of dof per node)
= 3×2 = 6
1 2 3 4 5 6
 1493 .33 1119 .9  1493 .33 1119 .9 0 0 1
 1119 .9 1119 .9  1119 .9 559 .9 0 0  2

 1493 .33  1119 .9 1493 .33     1119 .9  
     5040 2520  3
  5040   2520  
[ K ]  103  
 1119 .9   1119 .9   1119 .9  
559 .9      2520 840  4
  2520   1680  
 0 0  5040  2520 5040  2520  5
 
 0 0 2520 840  2520 1680  6
10
Chapter 5: Finite Element Formulation (Energy Methods)

Assembling the Global force vector:

 F1   P1   pl1 2 
 F   M   pl 2 12 
 2  1  1 
 F3   P2   pl1 2  Pi is the Point Load acting at node i
     Mi is the Point Moment acting at node i
 F4  M 2   pl1 12 
2
p =  10 kN/m (negative, acting down ward)
 F5   P3   0 
     
 F6   M 3   0 
 F1  0   7. 5    7.5 
F  0  1.875   1.875 
 2      
 F3  
 
20  
  7 . 5 
 
3  27 . 5 
    103    103    10   N
 F4  0  1.875   1.875 
 F5  0  0   0 
       
 F6   0   0   0 

Now writing the {F} = [K]{Q} relation, where,

{F} is the Global Force Vector
[K] is the Global Stiffness Matrix
{Q} is the Global Displacement Vector
1 2 3 4 5 6
 F1   1493 .33 1119 .9  1493 .33 1119 .9 0 0  1 Q1 
F   1119 .9 1119 .9  1119 .9 559 .9 0 0  2 Q2 
 2 
 F3  
3  1493 .33  1119 .9 6533 .33 1400  5040 2520  3 Q3 
   10    
 F4   1119 .9 559 .9 1400 2800  2520 840  4 Q4 
 F5   0 0  5040  2520 5040  2520  5 Q5 
     
 F6   0 0 2520 840  2520 1680  6 Q6 

Applying Boundary Conditions

Nodes 1 is fixed. Hence, both deflection & slope are zero.
 Q1 = 0, Q2 = 0
Node 3 is a roller support. Hence, deflection is zero and slope is expected here.
 Q5 = 0, Q6 ≠ 0
Apply a method of handling Boundary Conditions.
Applying Elimination Approach of handling Boundary Conditions.
Delete rows and columns corresponding to dof nos. which are fixed (or specified).
i.e., delete rows and columns 1, 2, & 5.
The modified (reduced) {F} = [K]{Q} relationship becomes:
3 4 6
27.5 6533.33 1400 2520  3 Q3 
   
103 1.875   103  1400 2800 840  4 Q4 
 0  840 1680  6 Q6 
   2520
Solving the set of linear algebraic equations,
Q3 = 0.01005 m; Q4 = 0.001379 rad; Q6 = 0.0144 rad;

11
Chapter 5: Finite Element Formulation (Energy Methods)

Numerical 2:
For the beam structure shown in figure, find the slopes and displacements at the supports, if
any and at the point where the 20 kN-m moment is acting. The beam diameter is 23 cm,
uniform throughout its length and the beam material is having a Modulus of elasticity of 200
GPa.

Solution:
The given beam can be modelled with minimum two elements. We have element
stiffness matrix for beam as,
 12 6l1  12 6l1 
 2
 6l1 2l1 
2
E1I1  6l1 4l1
[k ]  3
(1)

l1  12  6l1 12  6l1 
 2 
 6l1 4l1 
2
 6l1 2l1
Area Moment of Inertia of beam cross section,

12
Chapter 5: Finite Element Formulation (Energy Methods)

Evaluating the global force vector:

13
Chapter 5: Finite Element Formulation (Energy Methods)

By solving these equations, we get,

14
Chapter 5: Finite Element Formulation (Energy Methods)

Two Dimensional (2D) or Plane Frames

We consider here plane structures with rigidly connected members. These members
will be similar to the beams except that axial loads and axial deformations are present. The
elements also have different orientations. Figure shows a frame element. We have two
displacements and a rotational deformation for each node. The nodal displacement vector is
given by

 q1 
q 
 2
 q 
q   3  ------(3.27)
q4 
 q5 
 
q6 

We also define the local or

body coordinate system x1, y1such
that x1 is oriented along 1-2, with
direction cosines C, S (where C =
cosθ, S = sinθ). These are evaluated
using relationships given for the trussFig. 3.8 Local & Global displacements of a Plane Frame Element
element. The nodal displacement
vector in the local system is
 q11 
 1
q 2 
 q 1  ----------------(3.28)
q   31 
1

q 4 
q 1 
 51 
q6 
Recognizing that q31=q3 and q61=q6, which are rotations with respect to the body. We
obtain the local-global transformation
{q1}=[L]{q} where,

C S 0 0 0 0
S C 0 0 0 0
 
0 0 1 0 0 0
 L    --------(3.29)
0 0 0 C S 0
0 0 0 S C 0
 
0 0 0 0 0 1 

It is now observed that q21, q31, q51 and q61 are like the beam degrees of freedom, while
q11 and q41 are similar to the displacements of a rod elements as discussed. Combining the two
stiffnesses and arranging in proper locations. We get the element stiffness for a frame element
as
15
Chapter 5: Finite Element Formulation (Energy Methods)

 EA EA 
 l 0 0 - 0 0 
le
 e

 12 EI 6 EI 12 EI 6 EI 
 0 le 3 le 2
0 -
le 3 le 2 
 
 6 EI 4 EI 6 EI 2 EI 
 0 0 - 
le 2 le le 2 le 
 k1e    ---------(3.30)
 EA EA 
- 0 0 0 0 
 le le 
 12 EI 6 EI 12 EI 6 EI 
 0 - 3 - 0 - 2 
 le le 2 le 3 le 
 6 EI 2 EI 6 EI 4 EI 
 0 0 - 
 le 2 le le 2 le 

As discussed in the development of a truss element, we recognize that the element

strain energy is given by
U e  q1  k1e  q1
1 T

2
1
 q  L   k1e   L q
T T

2
1
 q  k e q
T

2
where [k ] = [L]T[ke1][L]---------------------------------------------------(3.31)
e

1 2 3 4 5 6
 P Q TS P Q
TS  1
 Q R TC Q  R TC  2

E  TS TC 4 I TS TC 2 I  3
[k e ]   
le   P Q TS P Q TS  4
 Q R TC Q R TC  5
 
 TS TC 2 I TS TC 4 I  6

12 I 2  12 I 
where, P  AC 2  S , Q   A  2  CS ,
le2  le 
12 I 2 6I
R  AS 2  C , T
l2 le
e
C = cos() S = sin()
If there is distributed load on a member, as shown

16
Chapter 5: Finite Element Formulation (Energy Methods)

Fig 3.9 Distributed load on a frame element

{q1}T {f1} = {q}T[L]T {f1} ------------------------------(3.32)

where  0 
 pl 
 e 
 2 
 pl 2 
 e 
 12 
{ f 1}    ........................................(3.33)
 0 
 ple 
 
 2 
  ple2 
 
 12 
The nodal loads due to the distributed load p are govern by

{f} = [L]T{f1} ------------------------------(3.34)

The values of are added to the global load vector. On gathering stiffness and loads, we get
system of equations
[K]{Q} = {F}

17
Chapter 5: Finite Element Formulation (Energy Methods)

Numerical for Practice:

Question statement:
For the beams configurations shown in following figures, determine the
displacements and slopes at the nodes, and the reactions at the supports through
finite element method.

**************************************

18