Fem class notes

Finite Element Method/Analysis Lecture Notes

1. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Finite Element Method
in Civil Engineering

Lecture Notes

Dr. A. S. Sayyad
Professor

Department of Civil Engineering
SRES’s Sanjivani College of Engineering,
Savitribai Phule Pune University,
Kopargaon-423603
Email: [email protected]
Ph. No.: ([phone]

Year-2017

2. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Unit- I
Theory of Elasticity

1.1 Fundamentals of theory of elasticity
Assumptions made in theory of elasticity:
1) Material of elastic body is continuous i.e. no sudden discontinuity such as
cracks, flaws, deep notches etc.
2) Material is homogenous, isotropic and elastic
3) strains and displacements are small
4) stress strain relationship is linear
5) higher order differential terms are neglected
6) small angle assumptions are valid such that 1sin ; cos ; tan      


1.2 Basic terms:
External forces: 1) Surface force 2) body force

1) Surface force: The force distributed over the surface of the body is called as
surface force. It is denoted by x,y & z and resolved into three components.

2) Body force: The forces distributed over the volume of the body are called as
body forces and are caused by gravity, magnetism and acceleration. Body forces
can be resolved into three components (X, Y & Z).

Internal forces: Internal forces are the stress resultants existing on the cut faces of
the isolated part of the body. Internal force distributed over an internal face may be
resolved into two components.
1) Normal component perpendicular to face (Normal stress)
2) Shear component tangential to face (Shear stress)

Normal stress: The normal force per unit area is called as normal stress. It is
denoted by ij
 where i represent the plane in which it acts and j represents the
direction of the stress.
Example:
plane and direction
plane and direction
plane and direction
xx
yy
zz
' yz' ' x'
' xz' ' y'
' xy' ' z'







3. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Sign conventions: Tensile Normal stress (Positive)
Compressive Normal stress (Negative)

Shear stress: The shear force components per unit area is called as shear stress
and denoted by ij
 . plane and direction
xy
' yz' ' y'


Displacements: As a result of a deformation of elastic body subjected to external
forces, point of a body are displaced. The displacements of a point are represented
by u, v, w in x, y, z directions respectively.

Strains: Measure of deformation of a body.
1) Linear strain: Change in dimension
2) Lateral strain: Change in lateral dimensions
3) Shear strain: Change in angle between originally perpendicular elements.
Considered as positive, if original angle is reduced and considered as negative, if
original angle is increased.

1.3 State of stress at a point:
In one plane there are three stress components. One normal stress ( ) and two
shear stresses ( ). At every point there are three mutually perpendicular planes.
i.e. orthogonal planes. Therefore, at every point there are nine stress components
(three normal stresses and six shear stresses). But since shear stresses are
complimentary (xy yx xz zx yz zy
,,        ), nine stress components are reduced to
six.  
3 3 3 3
xx xy xz xx xy xz
yx yy yz xy yy yz
zx zy zz xz yz zz
     
       
     

   
   
  
   
   
   


1.4 State of strain at a point:
In one plane there are three strain components. One normal strain ( ) and two
shear strains ( ). At every point there are three mutually perpendicular planes. i.e.
orthogonal planes. Therefore, at every point there are nine strain components (three
normal strains and six shear strains). But since shear strains are complimentary (xy yx xz zx yz zy
,,       
), nine strain components are reduced to six.

4. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 
3 3 3 3
xx xy xz xx xy xz
yx yy yz xy yy yz
zx zy zz xz yz zz
     
       
     

   
   
  
   
   
   


1.5 Equilibrium equations for 3D elasticity problem
Let consider as infinitesimal element of sides dx, dy and dz as shown in figure.
The stresses acting on the element due to external forces and body forces. These
stresses can be represented by nine components as given below.  
T
xx yy zz xy yx xz zx yz zy
, , , , , , , ,         


where, xx yy zz
,,   are the normal stresses and xy yx xz zx yz zy
, , , , ,      are the shear
stresses.


Plane Stresses on Positive faces Stresses on Negative faces Area of Plane
yz ' xx
xx xx
dx
x





xy'
xy xy
dx
x





' xz
xz xz
dx
x





xx


xy


xz


dy dz

5. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

xz yy'
yy yy
dy
y





yx'
yx yx
dy
y





yz'
yz yz
dy
y





yy


yx


yz


dx dz
xy ' zz
zz zz
dz
z





' zx
zx zx
dz
z





zy'
zy zy
dz
z





zz


zx


zy


dx dy

X = Component of body force in x-direction,
Y= Component of body force in y-direction and
Z= Component of body force in z-direction

Appling the conditions for forces and moments of static equilibrium, 0
x
F
0
yxxx
xx xx yx
zx
yx zx zx
dx dydz dydz dy dxdz
xy
dxdz dz dxdy dxdy X dxdydz
z

  

  

   
 
 

     


(1)
Simplifying the equation (1) we get
0
yxxx zx
dxdydz dxdydz dxdydz X dxdydz
x y z
 
   
  
(2)

Dividing by dxdydz to the equation (2), we will get
0
yx zxxx
X
x y z
 
   
   (First governing equation)
Similarly from 0
y
F
0
xy yy zy
Y
x y z
    
   
   (Second governing equation)

6. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

and from 0
z
F
0
yzxz zz
Z
x y z

   
   (Third governing equation)

Appling conditions of moment equilibrium to prove shear stresses are
complimentary

Now, taking moment about x-axis through the centroid of the element. The
coordinates of centroid of an element are2 2 2
dx dy dz
,,


 .
Notes:
1) Moment of 0
xx yy zz
     because they are passing through centroid.
2) Moment of 0
xy yx xz zx
      
3) Only moment due to yz zy
 will present about x-axis.
4) Anticlockwise moment positive and clockwise moment negative
22
0
22
yz
x yz yz
zy
zy zy
dy dy
M dy dxdz dxdz
y
dz dz
dz dxdy dxdy
z





  



   



(3)
Neglecting higher power of differential coefficients 
22
dy , dz
 in equation (3),
we get 20
22
yz zy
dy dz
dxdz dxdy 
yz zy

(Shear stress is complimentary)

Similarly, from0
y xz zx
M   
and from0
z xy yx
M   

7. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

1.6 Strain-Displacement Relationship
The six strain components  x y z xy xz yz
, , , , ,      are related to three displacement
component (u, v, w) at a point. Let consider elements AB and AC initially
perpendicular to each other in xy plane as shown in figure. ''AB and ''AC are the
deformed lengths of AB and AC respectively.

AB = Original line element in x-direction ''AB
= Deformed line element in x-direction
AC = Original line element in y-direction ''AC
= Deformed line element in y-direction
u = Displacement of point A in x-direction u
u dx
x



= Displacement of point B in x-direction
v = Displacement of point A in y-direction v
v dy
y



= Displacement of point C in y-direction
Change in length of the line element AB = u
u dx u
x




= u
dx
x


Linear/Normal strain of the element AB is x
 Change in length
Original length
x
u
dx
u
x
dx x




  


Similarly

8. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
y
v
y




and z
w
z




The angular displacement of line element AB =  v
dx
v
x
tan
dx x




  


Similarly the angular displacement of line element AC =  u
dy
uy
tan
dy y



  


Total shear strain xy
   =vu
xy



Similarly xz
wu
xz




and yz
wv
yz




Therefore, strain displacement relationship for the 3D elasticity problem is x
y
z
xy
xz
yz
u / x
v / y
w / z
u / y v / x
u / z w / x
v / z w / y






 
 

 
 
   
    
   
       
   
     

Example 1: An elastic body under the action of external forces has the
displacement field given by,    
2 2 2
2 5 3D x y i z y j x y k     
Evaluate components of strain at point (3, 1, 2)
Solution: 3 1 2
4 12
11
00
22
33
5 2 7
x
y
z
xy
xz
yz
, , mm/ mm
u / x x
v / y
w / z
u / y v / x y
u / z w / x
v / z w / y y







      
      
   
      
      
         
    
       
           
       
           

9. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: The general displacement field in Cartesian coordinates for stressed
body is, 22
2
0 015 0 03 0 005 0 03
0 03 0 001 0 005
u . x y . , v . y . xy
w . x . yz .
   
  

Where displacement coordinates u, v, w and x, y, z are in constant units. Find
Cartesian strain components at point (1, 0, 2)
Solution: 
2
1 0 2
0 03 0
0 01 0 03 0 03
0 001 0
0 015 0 03 0 015
0 006 0 006
0 001 0 002
x
y
z
xy
xz
yz
,,
u / x . xy
v / y . y . x .
w / z . y
u / y v / x . x . y .
u / z w / x . x .
v / z w / y . z .






     
    
  
    
    
        
     
     
         
     
         mm/ mm











Example 3: In an strained elastic body under the action of external forces has the
displacement field given by    
3 2 2
3 2 4 2 4D x y i z y j y z k     

Evaluate components of strain at point (2, 4, 1)

Example 4: In an strained elastic body under the action of external forces has the
displacement field given by    
2 2 2
23D x y i z y j x y k     

Evaluate components of strain at point (3, 1, 2)

1.7 Strain-Compatibility Conditions or Saint Venant’s Strain
compatibility conditions
In general, the elastic problem is a statically indeterminate and hence the
solution of elastic problem needs the concept of compatibility and stress-strain
relationship in addition to equilibrium equations.
The displacements are taken to be continuous single valued function and hence
the strains must be such that they do not cause any dislocations cracks or overlaps.
This means that continuity of structure must be maintained. This is the physical
interpretation of compatibility is that six components of strains must satisfy six

10. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

compatibility equations. These equations are derived from strain-displacement
relations.
23
22
x
x
uu
x y x y


  
  
   
(Differentiating x
 twice w.r.t y) (1) 2 3
22
y
y
vv
y x x y



  
   
(Differentiating y
 twice w.r.t x) (2) 2 33
22
xy
xy
v u v u
x y x y y x x y


   
    
       
(Differentiatingxy
 w.r.t x & y) (3)

Therefore, from equations (1)-(3) one can write 222
22
y xyx
y x x y


   
(I)
Similarly 222
22
z xzx
z x x z


    (II) 22 2
22
y yz z
z y y z


   
(III)
22
xy
xy
u v u v
y x z y z x z


   
    
      
(Differentiating xy
 w.r.t z) (4) 22
xz
xz
u w u w
z x y y z x y


    
    
      
(Differentiating xz
 w.r.t y) (5) 22
yz
yz
v w v w
z y x x z x y


   
    
      
(Differentiating yz
 w.r.t x) (6)

Solving equations (4) – (6) 2
2
xy yz xz
u
z y x y z

  
    
(7)
Differentiating equation (7) w.r.t x 3
2
xy yz xz
u
x z y x x y z
  
  

      
2
2
xy yz xz x
x z y x y z
   
   

     

(IV)
Similarly

11. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
2
2
xy yz y xz
y z x y x z
     
  

     

(V) 2
2
yz xyxz x
z y x z x y
  
  

     

(VI)

These are the six strain-compatibility conditions for the 3D elasticity problems.


Example 1: The general displacement field in Cartesian coordinates for stressed
body is 2
2
2
0 015 0 03
0 005 0 03
0 03 0 001 0 005
u . x y .
v . y . xy
w . x . yz .


  

Check whether the strain field given by above displacement field is compatible.

Solution: 222 2 2 2
2 2 2 2
22 2
22
0 0 0 0 0 0
0 0 0
y xyx x z xz
y yz z
, , , , , ,
y x x y z x x z
,,
z y y z
   

   
     
       
 
  
   


Above displacement field satisfy all strain displacement relationship. Therefore, it
gives possible/compatible state of strain.

Example 2: Check whether the following system of strain is possible. 23
2 2 2
23
3 3 2
4 6 4
3 34 12 9 2
x
y
xy
xy y x
x y x y
x y xy x y



    
   
    

12. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

1.8 Stress-strain relationship

3D Hooke’s law: 


1 0 0 0
1 0 0 0
1 0 0 01
0 0 0 2 1
0 0 0 2 1
0 0 0 2 1
xx
yy
zz
xy xy
xz xz
yz yz
E
 
 
 
 
 
 
   
   

   
      
   

   
    
   
      


It is also written in the form  
 
 
 
1 0 0 0
1 0 0 0
1 0 0 0
12
000
2
1 1 2
12
000
2
12
000
2
xx
yy
zz
xy xy
xz xz
yz yz
E
  
  

  



 







    
     
    

       
    

   

   
   
       




where E is Young’s modulus and  is Poisson’s ratio.

2D Hooke’s law:
Plane stress problem 2
10
10
1
0 0 1 2
xx
yy
xy xy
E
/
  
  

  
    
    
   


   
 
   

Plane strain problem  
10
10
1 1 2
0 0 1 2
xx
yy
xy xy
E
/
   
   

  
    
    
   


   
 
   

1D Hooke’s law: E

13. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example: The state of strain at a point is given by 
4
18 0 36
0 54 5 4 10
36 5 4 0
.
.




  



mm/mm
Determine the stress matrix if E = 210 GPa and 03.
Solution: 3
4
0 7 0 3 0 3 0 0 0 18
0 3 0 7 0 3 0 0 0 54
0 3 0 3 0 7 0 0 0 0210 10
10
0 0 0 0 2 0 0 01 3 0 4
0 0 0 0 0 2 0 36
0 0 0 0 0 0 2 5 4
x
y
z
xy
xz
yz
. . .
. . .
. . .
...
.
..







   
   

   
   
    

    
     
    
   
145 38
1308 04
436 6
0
290 71
43 6
x
y
z
xy
xz
yz
.
.
.
MPa
.
.






 
 

 
 
   
   
   
   
 
or 
145 38 0 290 71
0 1308 04 43 6
290 71 43 6 436 6
..
..
. . .








1.9 2D Elasticity Problems
1) Plane stress problem
2) Plane strain problem
3) Axisymmetric problem

Plane stress problem: Two dimensional elasticity problems under the following
conditions are considered as plane stress problem.
1) One dimension is very small as compared to other two dimensions
e.g. Rectangular plate (Thickness is very small as compared to length and
width)
2) The loads acting on the body are in the plane perpendicular to the thickness
of the body i.e. z-axis. The loads are uniformly distributed over the
thickness.
3) The stresses in the small direction (normally z) must be zero. 0
zz xz yz
    
but, 0
z


14. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Equations of Equilibrium 0
xyxx
X
xy

  

and 0
xy yy
Y
xy

  

State of stress at a point 
xx
yy
xy









State of strain at a point 
xx
yy
xy








and 0
z

Strain-Displacement relation x
u
x




, xy
vu
xy




Strain compatibility condition 222
22
y xyx
y x x y


   

Stress-strain relation 
10
1
10
0 0 2 1
xx
yy
xy xy
E
  
  
  
    
   
   

   
 
   
or2
10
10
1
0 0 1 2
xx
yy
xy xy
E
/
  
  

  
    
    
   


   
 
   

15. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Compatibility conditions in-terms of stresses for plane stress
problem:
As discussed in the Saint Venant’s strain compatibility conditions, substituting
plane stress conditions in the three dimensional Saint venant’s compatibility
conditions. Non-zero strain components in plane stress problems are x y z xy
, , &   
which functions of x and y. Substitute these components of strain
in compatibility equation of plane stress problem. 222
22
y xyx
y x x y


   

Other strain compatibility conditions are not imposed in plane stress problem. The
stress compatibility can be obtained either by substituting plane stress condition in
the Beltrami-Michell compatibility or directly substituting strains in-terms of
stress.    

2 2 2
2 2 2 2
21
11
x y y x xy
EE
y x x y E

    

     
   
   
         
(1)

2 2 222
2 2 2 2
21
y y xyxx
y x x y x y
  

     
       
     
   
(2)
Now, from equilibrium equations of plane stress problem neglecting body forces 22
2
0
xy xyxx
x y x x y

    
    
(3) 22
2
0
xy y y xy
x y y x y
      
    
    
(4)
Adding equations (3) and (4)
222
22
2
y xyx
x y x y

  
    (5)
Put equation (5) into the equation (2) 
2 2 22 2 2
2 2 2 2 2 2
1
y y yx x x
y x x y x y
    

         
           
     
     
2222
2 2 2 2
0
yyxx
x x y y

   
   
 
22
22
0
xy
yx


  


(6)

16. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 
2
0
xy
  
where 22
2
22
0
yx

  


This is called as stress compatibility conditions in-terms of stress. It is also written
in the form (putting Eq. (5) into Eq. (2) one can wirte)

2 2 22
22
2 2 1
y xy xyx
y x x y x y
  

   
   
     


222
22
2
y xyx
y x x y


   
(7)

This is also called as compatibility conditions in-terms of stresses. This is the plane
stress idealization of Beltrami-Michell compatibility conditions. The stress
compatibility equation is valid only for an isotropic body with constant body force
under equilibrium.

Example 1: Show that the following state of stress is in equilibrium 22
22
2
2
3 4 8 4
23
62
2
x
y
xy
x xy y
x xy y
x
xy y



   
  

   



Solution:
Differential equation of equilibrium of 2D elasticity problem is given by 0
xyx
xy



and 0
xy y
xy


 (1) 6 4 6 4
66
xyx
xy y
x y, x y
xy
x y, x y
xy



    


    


Putting these differential quantities in equations of equilibrium (1) and satisfying
those. This proves that the given state of stress is in equilibrium.

17. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: A 2D stress distribution at a point in xy coordinate system is given
below. Find constants A, B & C if system of stress is in equilibrium. The body
force is zero. 23
2
32
10
15
x
y
xy
xy Ax
. Bxy
By Cx y



  

  

Solution: 2 2 2 2
10 3 3
xyx
y Ax , By Cx
xy
 
     

23
xy y
Cxy, Bxy
xy

   


Putting these derivatives in following two governing equations   
22
3 10 3 0
xyx
A C x B y
xy

      


and 2 3 0
xy y
Cxy Bxy
xy

    


32C B /   Putting this in above equation
L.H.S. = 0 only when 3 0 and 10+3B=0AC
10 3 and A=5/3B/ 


Plane strain problem: Two dimensional elasticity problems under the following
conditions are considered as plane strain problem.
1) One dimension is infinitely long as compared to other two dimensions
e.g. Retaining wall, Dam, Bridge (Length is infinitely long as compared to
depth and width)
2) External force is perpendicular to the z-axis.
3) The strains in the long direction (normally z) must be zero. 0
zz xz yz
    
but 0
z


18. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Equations of Equilibrium 0
xyxx
X
xy

  

and 0
xy yy
Y
xy

  

State of stress at a point 
xx
yy
xy







 0
z


State of strain at a point 
xx
yy
xy









Strain-Displacement relation x
u
x




, xy
vu
xy




Strain compatibility condition 222
22
y xyx
y x x y


   


Stress-strain relation 10
1
10
0 0 2
xx
yy
xy xy
E
   

   

    
   
     

   

   
or  
10
10
1 1 2
0 0 1 2
xx
yy
xy xy
E
/
   
   

  
    
    
   


   
 
   

19. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Compatibility conditions in-terms of stresses for plane strain problem: 222
22
y xyx
y x x y


   
(1) 
 

 

22
22
2
11
11
21
x y y x
xy
y E x E
x y E

     


   
    
   

   





(2)
Now, from equilibrium equations of plane stress problem neglecting body forces 22
2
0
xy xyxx
x y x x y

    
    
(3) 22
2
0
xy y y xy
x y y x y
      
    
    
(4)
Adding equations (3) and (4) 222
22
2
y xyx
x y x y

  
   
(5)
Put equation (5) into equation (2)  
22
22
0
xy
yx


  


 
2
0
xy
  

20. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Summary of elasticity problems
Elasticity problem Equations No. of Equations No. of Unknowns
3D Equilibrium equations 0
yx zxxx
X
x y z
 
   
   0
xy yy zy
Y
x y z
    
   
   0
yzxz zz
Z
x y z

   
  

03 06 stresses
Stress-strain relations  
 
 
x x y z
y y x z
z z x y
xy xy
xz xz
yz yz
/E
/E
/E
/G
/G
/G
   
   
   



  
  
  




06 06 strains
Strain-displacement relations x
y
z
xy
xz
yz
u / x
v / y
w / z
u / y v / x
u / z w / x
v / z w / y






  
  
  
     
     
     

06 03 displacements
Compatibility conditions 222
22
y xyx
y x x y


   
222
22
z xzx
z x x z


   
22 2
22
y yz z
z y y z


   

03 ---
2D
Plane stress
Equilibrium equations 0
yxxx
X
xy

  
 0
xy yy
Y
xy

  


02 03 stresses

21. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stress-strain relations  
 
 
1
x x y
y y x
z x y
xy xy
/E
/E
/E
G
  
  
  



  


03 03 strains
Strain-displacement relations x
y
xy
u / x
v / y
u / y v / x



  
  
     

03 02 displacements
Compatibility conditions 222
22
y xyx
y x x y


   

01 ---
2D
Plane strain
Equilibrium equations 0
yxxx
X
xy

  
 0
xy yy
Y
xy

  


02 03 stresses
Stress-strain relations  
 
 
1
x x y z
y y x z
z x y
xy xy
/E
/E
G
   
   
  

  
  
  


03 03 strains
Strain-displacement relations x
y
xy
u / x
v / y
u / y v / x



  
  
     

03 02 displacements
Compatibility conditions 222
22
y xyx
y x x y


   

01 ---

22. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Unit-II
Finite Element Analysis of Spring Assembly

Springs are 1D structures subjected to axial force only. The degree of freedom at
each node is one i.e. axial displacement. Stiffness matrix for spring element having
stiffness constant k is given below which can be obtained by giving unit
displacement one by one at each node.
Let consider a two noded spring element with ui and uj displacements at each
nodes.

Let unit displacement at node i

Let unit displacement at node j


11
11
ij
i
j
uu
u
Kk
u






Procedure for the solution of numerical examples
1) Divide the spring assembly into number of members
2) Calculate total degrees of freedom
3) Determine stiffness matrix of each spring element
4) Assemble the global stiffness matrix
5) Impose the boundary conditions
6) Determine reduced stiffness matrix
7) Apply governing equation to determine unknown joint displacements. Kf

where f = Nodal load vector

Example 1: Determine elongations at each node and hence the forces in springs.

Solution:
Step 1: Discretization
Element k (N/m) Nodes Displacements (m) Boundary conditions
1 500 1-2 u1-u2 u1 = 0
2 100 2-3 u2-u3 ---

23. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 2: Element stiffness matrices 
12
1
11
2
1 1 1 1
500
1 1 1 1
uu
u
Kk
u
   

   
   

23
2
22
3
1 1 1 1
100
1 1 1 1
uu
u
Kk
u
   

   
   

Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix   
1 2 3
1
2
3
500 500 0
500 500 100 100
0 100 100
u u u
u
Ku
u


   

 



Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is 
23
2
3
600 100
100 100
uu
u
K
u





Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
3
600 100 0
100 100 5
u
u
    
   
    

 
23
0 01 and 0 06u . m u . m    
Step 6: Calculation of spring force
Spring 1: 
111
Kf
11
22
11
500
11
uf
uf
   
   
    

1
2
1 1 0
500
1 1 0 01
f
f.
   
   
    (12
0 and 0 01u u . )  
12
5 T and 5 Tf N f N

Spring 2:  
34
5 T and 5 Tf N f N (23
0 01 and 0 06u . u . )

24. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: Determine elongations at each node and hence the forces in springs.
Take F3 = 5000 N.

Solution:
Step 1: Discretization

Element k (N/m) Nodes Displacements (m) Boundary conditions
1 1000 1-2 u1-u2 u1 = 0
2 2000 2-3 u2-u3 ---
3 3000 3-4 u3-u4 u4 = 0

Step 2: Element stiffness matrices 
12
1
11
2
1 1 1 1
1000
1 1 1 1
uu
u
Kk
u
   

   
   

23
2
22
3
1 1 1 1
2000
1 1 1 1
uu
u
Kk
u
   

   
   

34
3
33
4
1 1 1 1
3000
1 1 1 1
uu
u
Kk
u
   

   
   

Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix 
 
 
1 2 3 4
1
2
3
4
1000 1000 0 0
1000 1000 2000 2000 0
0 2000 2000 3000 3000
0 0 3000 3000
u u u u
u
u
K
u
u


  

   




Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 and u4 = 0
Eliminate first row, first column and fourth row and fourth column.
Therefore reduced stiffness matrix is

25. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

23
2
3
3000 2000
2000 5000
uu
u
K
u






Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
3
3000 2000 0
2000 5000 5000
u
u
    
   
    

Ans.  
23
0 909 and 1 363u . m u . m    
Step 6: Calculation of spring force
Spring 1: 
111
Kf
11
22
11
1000
11
uf
uf
   
   
    

1
2
1 1 0
1000
1 1 0 909
f
f.
   
   
    (12
0 and 0 909u u . )
 
12
909 T and 909 Tf N f N

Spring 2:  
34
909 T and 909 Tf N f N (23
0 909 and 1 363u . u . )
Spring 3:  
56
4091 C and 4091 Cf N f N (34
1 363 and 0 0u . u . )


Example 3: Determine elongations at nodes 2 and 4 if u3 = 0.02 m and hence the
force at node 3.

Solution:
Step 1: Discretization
Element k (N/m) Nodes Displacements (m) Boundary conditions
1 500 1-2 u1-u2 u1 = 0
2 100 2-3 u2-u3 u3 = 0.02
3 200 3-4 u3-u4 ---

Step 2: Element stiffness matrices

26. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

12
1
11
2
1 1 1 1
500
1 1 1 1
uu
u
Kk
u
   

   
   

23
2
22
3
1 1 1 1
100
1 1 1 1
uu
u
Kk
u
   

   
   

34
3
33
4
1 1 1 1
200
1 1 1 1
uu
u
Kk
u
   

   
   

Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix 
 
 
1 2 3 4
1
2
3
4
500 500 0 0
500 500 100 100 0
0 100 100 200 200
0 0 200 200
u u u u
u
u
K
u
u


  

   




Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is 
2 3 4
2
3
4
600 100 0
100 300 200
0 200 200
u u u
u
Ku
u


  

 



Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
2
4
600 100 0 0
100 300 200 0 02
0 200 200 100
u
.f
u
    
   
      

   
 
    

  
2 4 2
0 00333 0 52 and 98 33u . m , u . m f . N      

27. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 4: Determine elongation at node 2 and pulling force (F) at node 3 for the
spring assembly given below. Take pull at node 3, 0.06m.

Solution
Step 1: Discretization
Element k (N/m) Nodes Displacements (m) Boundary conditions
1 500 1-2 u1-u2 u1 = 0
2 100 2-3 u2-u3 u1 = 0.06m

Step 2: Element stiffness matrices 
12
1
11
2
1 1 1 1
500
1 1 1 1
uu
u
Kk
u
   

   
   

23
2
22
3
1 1 1 1
100
1 1 1 1
uu
u
Kk
u
   

   
   

Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix   
1 2 3
1
2
3
500 500 0
500 500 100 100
0 100 100
u u u
u
Ku
u


   

 



Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is 
23
2
3
600 100
100 100
uu
u
K
u





Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
600 100 0
100 100 0 06
u
.F
    
   
    

 
2
0 01 and 5u . m F N    

28. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 5: Determine spring elongations and force at node 5 for the spring
assembly as shown in figure. Take stiffness of all spring elements 200 kN/m.

Solution:
Step 1: Discretization
Element k (N/m) Nodes Displacements (m) Boundary conditions
1 200 1-2 u1-u2 u1 = 0
2 200 2-3 u2-u3 ---
3 200 3-4 u3-u4 ---
4 200 4-5 u4-u5 u5 = 0.02 m

Step 2: Element stiffness matrices 
12
1
11
2
1 1 1 1
200
1 1 1 1
uu
u
Kk
u
   

   
   

23
2
22
3
1 1 1 1
200
1 1 1 1
uu
u
Kk
u
   

   
   

34
3
33
4
1 1 1 1
200
1 1 1 1
uu
u
Kk
u
   

   
   

45
4
44
5
1 1 1 1
200
1 1 1 1
uu
u
Kk
u
   

   
   


Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix 
1 2 3 4 5
1
2
3
4
5
200 200 0 0 0
200 400 200 0 0
0 200 400 200 0
0 0 200 400 200
0 0 0 200 200
u u u u u
u
u
Ku
u
u




 



 


29. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is 
2 3 4 5
2
3
4
5
400 200 0 0
200 400 200 0
0 200 400 200
0 0 200 200
u u u u
u
u
K
u
u




 



Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
3
4
400 200 0 0 0
200 400 200 0 0
0 200 400 200 0
0 0 200 200 0 02
u
u
u
.F
    
    
    
    
 
   

       

Ans.    
2 3 4
0 005 0 01 0 015 and 1u . m ,u . m ,u . m F kN       


Example 6: Figure shows three springs connected parallel. Using finite element
method determines the deflections of individual springs.


Solution: Step 1: Discretization
Element k (N/mm) Nodes Displacements (mm) Boundary conditions
1 10 1-2 u1-u2 u1 = 0, u2 =?
2 20 3-4 u3-u4 u3 = 0, u2 = u4,
3 40 5-6 u5-u6 u5 = 0, u2 = u6,
Step 2: Element stiffness matrices

30. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

12
1
11
2
1 1 1 1
10
1 1 1 1
uu
u
Kk
u
     

   
   


34
3
22
4
1 1 1 1
20
1 1 1 1
uu
u
Kk
u
     

   
   


56
5
33
6
1 1 1 1
40
1 1 1 1
uu
u
Kk
u
     

   
   


Step 3: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0, u3 = 0, u5 = 0, u2 = u4, u2 = u6.
Therefore reduced stiffness matrix is 10 20 40 70K   

Step 4: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
70 700u


246
10u u u mm   


Example 7: Figure shows cluster of four springs. One end of the spring assembly
is fixed and a force of 1000 N is applied at the other end. Using the finite element
method, determine deflection of each spring.

31. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Solution:
Step 1: Discretization
Element k (N/mm) Nodes Displacements (mm) Boundary conditions
1 4 1-2 u1-u2 u1 = 0, u2 =?
2 8 3-4 u3-u4 u3 = 0, u4 = u2
3 20 5-6 u5-u6 u5 = 0, u6 =?
4 10 7-8 u7-u8 u7 = u2, u8 = u6

Step 2: Element stiffness matrices 
12
1
11
2
1 1 1 1
4
1 1 1 1
uu
u
Kk
u
     

   
   


34
3
22
4
1 1 1 1
8
1 1 1 1
uu
u
Kk
u
     

   
   


56
5
33
6
1 1 1 1
20
1 1 1 1
uu
u
Kk
u
     

   
   


78
7
33
8
1 1 1 1
10
1 1 1 1
uu
u
Kk
u
   

   
   

Step 3: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0, u3 = 0, u5 = 0, u2 = u4, u2 = u7, u8 = u6
Therefore reduced stiffness matrix is 
26
2
6
22 10
10 30
uu
u
K
u





Step 4: Determine unknown joint displacements
Applying Equation of Equilibrium Kf
2
6
22 10 0
10 30 1000
u
u
   
   
   
 
2 4 7 8 6
17 857 39 286u u u . mm u u . mm      

32. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 8: A three spring system shown in figure has stiffnesses k1 = 40 N/mm,
k2 = 50 N/mm and k3 = 80 N/mm. The loads applied are F1 = 100 N and F2 = 50 N.
Calculate displacements at nodal points.

Example 9: The system of springs, subjected to a load of 20 kN is shown in figure.
Find the deflection of each spring.

Example 10: The figure shows cluster of five springs. One end of the assembly is
fixed while a force of 1 kN is applied at the other end. Using finite element method
determines the deflection of each spring. (Ans. 24. 39mm, 24.39mm, 34.146mm,
14.634mm)

33. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Finite Element Analysis of trusses

Stiffness matrix of a truss element
The truss may be statically determinate or indeterminate. All members are
subjected to only direct stresses (tensile or compressive). Joint displacements are
selected as unknown variables. Since there is no bending of the members we have
to ensure only displacement continuity (C
0
continuity) and there is no need to
worry about slope continuity (C
1
continuity).
Here we select two noded bar element for the formulation of stiffness matrix of
truss element. Since the members are subjected to only axial forces, the
displacements are only in the axial directions of the members. Therefore, the nodal
displacement vector for the bar element is 
1
2
e
u'
x'
u'




where, 12
and
''
uu are the displacements in axial direction of the element. The
stiffness matrix of a bar element is 12
1
2
11
11
'
u' u'
u'AE
K
u'L





Transformation matrix for the truss:
''xy
= Local coordinate
systems
x, y = global coordinate
system
12
u' ,u'
= Displacements in
local coordinate system
1 1 2 2
u ,v ,u ,v
= Displacements
in global coordinate system

=Angle measured in
anticlockwise sense w.r.t.
positive x-axis.
Since axial directions of all members of truss are not same, hence in global
coordinate system (x-y) there are two displacement components at every node.
Hence the nodal displacement vector for typical truss element is

34. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

1
1
2
2
e
u
v
x
u
v








Refereeing above figure,
At Node 1, At Node 2,
1 1 1
cos sin
'
u u v
2 2 2
cos sin
'
u u v

Therefore, in matrix form above relation are 1
11
22
2
cos sin 0 0
0 0 cos sin
'
'
u
vu
uu
v




 
   
 



'
e
x L x

where, 
'
e
x
= vector of local unknowns x
= vector of global unknowns L
= Transformation matrix = 
00
00
lm
L
lm




where, 2 1 2 1
cos or sin or
x x y y
l l m m
LL


   
Stiffness matrix of truss element in global coordinate system '
T
K L K L

0
0 1 1 0 0
0 1 1 0 0
0
l
m l mAE
K
l l mL
m


  

  
   



35. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

0
0
0
0
l
m l m l mAE
K
l l m l mL
m





 



1 1 2 2
22
1
22
1
22
2
22
2
u v u v
ul lm l lm
vlm m lm mAE
K
uLl lm l lm
vlm m lm m
 









Example 1: Analyze the truss as shown in figure. Cross-sectional area of members
are AB=1000 mm
2
, BC=800 mm
2
, CA= 800 mm
2
. Take E = 2 × 10
5
MPa



Solution: Step 1: Degrees of freedom: 06 (A A B B c c
u ,v ,u ,v ,u ,v )
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AB uA, vA, uB, vB uA = 0
A
v
2 BC uB, vB, uC, vC 0
B
v
3 CA uC, vC, uA, vA ---

Assume x-axis horizontal through point c and vertical through point A. The
coordinate of node A(0, 1.5), B(4, 1.5) and C (2, 0). Take E in GPa

Member x2-x1 y2-y1 L l m AE/L (kN/mm)
AB 4 0 4 1 0 50
BC -2 -1.5 2.5 -0.8 -0.6 64
CA -2 1.5 2.5 -0.8 0.6 64

36. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 2: Element stiffness matrices
Stiffness matrix of element AB: Stiffness matrix of element BC: 
1 0 1 0
0 0 0 0
50
1 0 1 0
0 0 0 0
A A B B
A
A
AB
B
B
u v u v
u
v
K
u
v







0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
64
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
B B c c
B
B
BC
c
c
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .







Stiffness matrix of element CA: 
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
64
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
c c A A
c
c
CA
A
A
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .








Step 3: Global stiffness matrix (Total DOF are 06, size of stiffness matrix 6×6) 
90 96 30 72 50 0 40 96 30 72
30 72 23 04 0 0 30 72 23 04
50 0 30 72
0 0 30 72 23 0
90 96 40 96 30 72
4
4 30 72 23 04
40 96 30 72 300 96 81 92 072
30 72 23 04 23 0430 72 0 46 08
A A B B c c
A
A
B
B
c
u v u v u v
u. . . .
v. . . .
u.
K
v..
. . .
..
.
..
u..
. .
.
..
  













c
v



  

Step 4: Reduced stiffness matrix (Since uA=0
A
v , 0
B
v eliminate corresponding
rows and columns from global stiffness matrix) 
90 96 40 96 30 72
40 96 81 9 0
30 72 0 46 08
B c c
B
c
c
. . .
..
..
u u v
u
Ku
v












Step 5: Equation of equilibrium Kf
90 96 40 96 30 72 0
40 96 81 9 0 0
30 72 0 46 08 120
B
c
c
. . . u
. . u
. . v
    
   
    

   

    
1 6 0 8 3 67
B c c
u . mm, u . mm, v . mm     

37. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: Figure shows a plane truss with three members. Cross-sectional area
of all members 800 mm
2
Young modulus is 200 KN/mm
2
. Determine
deflection at loaded joint.

Solution:
Step 1: Degrees of freedom: 08 (A A B B c c D D
u ,v ,u ,v ,u ,v ,u ,v )
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AD uA, vA, uD, vD 0
AA
uv
2 BD uB, vB, uD, vD 0
BB
uv
3 CD uC, vC, uD, vD 0
cc
uv

Assume origin support A (0, 0). The coordinates of other nodes B (1000, 0),
C(2000, 0) and D(1500, 1000)
Member x2-x1 y2-y1 L l m AE/L (kN/mm)
AD 1500 1000 1802.8 0.832 0.555 88.75
BD 500 1000 1118 0.447 0.894 143.112
CD -500 1000 1118 -0.447 0.894 143.112

Step 2: Element stiffness matrices
Stiffness matrix of element AD: 
61 43 40 98 61 43 40 98
40 98 27 34 40 98 27 34
6 61 43 40 981 43 4
40
0 98
40 98 27 34 98 27 34
A A D D
A
A
AD
D
D
u v u v
u. . . .
v. . . .
K
u. ..
.
.
.. .v
  

  





(0
AA
uv )

38. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix of element BD: 
28 59 57 19 28 59 57 19
57 19 114 38 57 19 114 38
28 59 57 19 28 59 57 19
57 19 114 38 57 19 114 38
B B D D
B
B
BD
D
D
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .
  

  





(0
BB
uv )
Stiffness matrix of element CD: 
28 59 57 19 28 59 57 19
57 19 114 38 57 19 114 38
28 59 57 19 28 59 57 19
57 19 114 38 57 19 114 38
C C D D
C
C
CD
D
D
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .
  

  





(0
cc
uv )
Step 3: Reduced stiffness matrix 
118 61 40 98
40 98 256 10
DD
D
D
uv
u..
K
v..





Step 4: Equation of equilibrium Kf
118 61 40 98 200
40 98 256 10 0
D
D
u..
v..
   
   
    
1 785 0 286
DD
u . mm, v . mm  

39. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 3: for the truss as shown in figure using finite element method,
determines deflections at loaded joints. The joint B is subjected to 50 kN
horizontal force towards left and 80 kN force vertically downward. Take cross-
sectional area of all members 1000 mm
2
Young modulus is 200 GPa.

Solution: Step 1: Degrees of freedom: 06 (A A B B c c D D
u ,v ,u ,v ,u ,v ,u ,v ).
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AB uA, vA, uB, vB 0
AA
uv
2 DB uD, vD, uB, vB 0
DD
uv
3 CB uC, vC, uB, vB 0
cc
uv

Assume origin point B. The coordinates of points areA (-4, 3), B (0,0), C (4,-3), D
(-4, -3)
Member x2-x1 y2-y1 L l m AE/L (kN/mm)
AB 4000 -3000 5000 0.8 -0.6 40
DB 4000 3000 5000 0.8 0.6 40
CB -4000 3000 5000 -0.8 0.6 40

Step 2: Stiffness matrix of element AB: 
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
40
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
A A B B
A
A
AB
B
B
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .
  

  





(0
AA
uv )

40. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix of element DB: 
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
40
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
D D B B
D
D
DB
B
B
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .
  

  





(0
DD
uv )
Stiffness matrix of element CB: 
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
40
0 64 0 48 0 64 0 48
0 48 0 36 0 48 0 36
C C B B
C
C
CB
B
B
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .
  

  





(0
cc
uv )
Step 3: Reduced stiffness matrix 
1 92 0 48
40
0 48 1 08
BB
B
B
uv
u..
K
v..





Step 4: Equation of equilibrium Kf
1 92 0 48 50
40
0 48 1 08 80
B
B
u..
v..
   
   
    
1 25 2 4
BB
u . mm, v . mm   


Example 3: Determine the deflections at loaded joint in two bar truss supported by
spring as shown in figure. Bar one has length of 5m and bar two a length of 10m.
The stiffness of spring is 2000 kN/m. Take A = 5×10
-4
m
2
and E = 200 GPa.

41. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Solution:
Step 1: Degrees of freedom: 06 (1 1 2 2 3 3
u ,v ,u ,v ,u ,v )
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 1-2 u1, v1, u2, v2 22
0uv
2 1-3 u1, v1, u3, v3 33
0uv
3 1-4 v1, v4 v4 = 0

Take origin node 1. The coordinates of nodes are
1(0, 0), 2(-3.535, 3.535), 3(-10, 0)
Member x2-x1 y2-y1 L l m AE/L (kN/mm)
1-2 -3.535 3.535 5 -0.707 0.707 200×10
2

1-3 -10 0 10 -1 0 100×10
2

1-4 --- --- --- --- --- ---

Step 2: Stiffness matrix of element 1-2: 
1 1 2 2
1
12
12
2
2
0 5 0 5 0 5 0 5
0 5 0 5 0 5 0 5
200 10
0 5 0 5 0 5 0 5
0 5 0 5 0 5 0 5
u v u v
u. . . .
v. . . .
K
u. . . .
v. . . .





  

  

(22
0uv )
Stiffness matrix of element 1-3: 
1 1 3 3
1
12
13
3
3
1 0 1 0
0 0 0 0
100 10
1 0 1 0
0 0 0 0
u v u v
u
v
K
u
v








(33
0uv )
Stiffness matrix of Spring element 
14
1
3
4
11
2000
11
vv
v
K
v






Step 3: Reduced stiffness matrix

42. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

11
1
1
20000 10000
10000 12000
uv
u
K
v





Step 4: Equation of equilibrium: Kf 1
1
20000 10000 0
10000 12000 40
u
v
    
   
    
11
2 857 5 714u . mm, v . mm   



Example: For the plane truss shown in figure, determine the x and y components
of displacements at node 1. Take E = 70 GPa and A = 500 mm
2
for all elements.
Length of member 1-3 is 2500mm.

Example: For the plane truss composed of three elements shown in figure
subjected to a downward force of 50 kN applied at node 1, determine the x and y
components of displacements at node 1. Take E = 200 GPa and A = 1000 mm
2
for
all elements.

Example: Figure shows a plane truss with two members. Both the members are of
cross-sectional area 70.71 mm
2
. Young’s modulus is 200 kN/mm
2
. Determine
deflections of loaded joint and hence the member forces.

43. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603



Example: A steel truss as shown in figure. The modulus of elasticity is 210 GPa.
The cross sectional area of member AB is 300 mm
2
, BC is 400 mm
2
and AC is 500
mm
2
. Calculate the horizontal and vertical displacements at point ‘A’ using finite
element method.

Example: Figure shows a plane truss with three members. All members are of
length 1000 mm and cross-sectional area 600 mm
2
. Young’s modulus is 150
kN/mm
2
. Determine unknown joint displacements of the truss.

Example: For the two bar truss shown in figure determine the displacements at the
loaded joint using stiffness matrix method. Take A = 200 mm
2
and E = 70
GPa.

44. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example: Find the vertical and horizontal deflection at point C for the two
member truss as shown in figure. Area of inclined member is 2000 mm
2

whereas horizontal member is 1600 mm
2
. Take E = 200 GPa


Example: Figure shows plane truss with three members. All members are of
length 1000mm and c/s area 600mm2. E=150 KN/mm2. Determine forces in
members of truss using finite element method.


Example: Analyze the two member truss shown in figure using finite element
method. Take c/s area of each member 1000 mm
2
and E = 200 GPa. The length
of each member is 5m.


Example: For the plane truss structure shown in figure, determine the
displacements at the loaded joint using finite element method. Assume A =
2000 mm
2
and E = 200 GPa.

45. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Finite Element Analysis of Continuous Beams

A beam is a structural member which is subjected to bending deformation. There
are several methods available in the literature for the analysis of continuous beams
such as slope deflection method, moment distribution method, flexibility matrix
method, stiffness matrix method, three moment theorem etc. However all these
methods have limitations if either geometry, loading material properties or
boundary conditions. Finite element method can well handle such problems easily.

Element nodal load vector/ Equivalent load vector
In finite element method, the external forces are necessary to act at the joints
corresponding to joint displacements, which do not happen always. Beams are
often subjected to member forces, therefore these member forces we have to
convert into nodal forces. Vector of these forces is called as element nodal load
vector and apposite vector is called as equivalent load vector.

Degree of Kinematic Indeterminacy/Degrees of Freedom
Beam has two degrees of freedom at each point i.e. vertical translation and
rotation. Whereas frame has three degrees of freedom at each point i.e. two
displacements and one rotation.

Type of Support Kinematic
Unknowns for Beam
Kinematic Unknowns
for Frame
Hinge

1 ( ) 1 ( )
Roller

1 ( ) 2 (, )
Fixed

0 0
Spring

2 (, ) 2 (, )
Guided/Slider

1 ( ) 1 ( )
Internal Hinge 3 (12
,, ) 3 (12
,, )

46. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Steps for the solution of continuous (Indeterminate) beams using finite
element method:
1. Divide the beam into number of elements (Take one member as one element)
2. Identify total degrees of freedom (Two D.O.F. at each node, translation and
rotation)
3. Determine stiffness matrices of all elements ([K]1, [K]2………)
4. Assemble the global stiffness matrix [K]
5. Impose the boundary conditions and determine reduced stiffness matrix
6. Determine element nodal load vector [q] (Restrained structure)
7. Determine equivalent load vector [f]
8. Apply equation of equilibrium [K]{Δ}={f} and determine unknown joint
displacements.
9. Apply equation [K]{Δ}+[q] ={f} to determine reactions and moments

Stiffness matrix of beam
1 = Translation at node A
2= Rotation at node A
3 = Translation at node B
4= Rotation at node B


EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
K
EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
3 2 3 2
22
3 2 3 2
22
1 2 3 4
1 Reaction12 6 12 6
2 Moment6 4 6 2
3 Reaction12 6 12 6
4 Moment6 2 6 4
Reaction Moment Reaction Moment
 




   


   


Note: 1) Action corresponding to translation is reaction
2) Action corresponding to rotation is moment

47. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 1: Analyse the beam as shown in figure using finite element method.
Take EI = constant.

Solution:
Step 1: Degrees of freedom: 06 (02 DOF at each node, translation and rotation)
No. of elements: 02 (AB, BC)

Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2=3= zero (Fixed support)
2 2-3 3,4,5,6 3=5=zero (simple supports)

Step 2: Element stiffness matrices: Using standard stiffness matrix of beam
element, obtain local stiffness matrix of each element separately. (Note that the
moment of inertia of AB is 2I and BC is I).
The local stiffness matrix of element AB is: 
1 2 3 4
0 111 0 333 0 111 0 333 1
0 333 1 333 0 333 0 667 2
0 111 0 333 0 111 0 333 3
0 333 0 667 0 333 1 33
K = EI
34
AB
. . . .
. . . .
. . . .
. . . .




  



Similarly the local stiffness matrix of element BC is: 
3 4 5 6
0 1875 0 375 0 1875 0 375 3
0 375 1 0 0 375 0 5 4
0 1875 0 375 0 1875 0 375 5
0 375 0 5 0 375 1 0 6
K = EI
BC
. . . .
. . . .
. . . .
. . . .




  



Step 3: Assemble global stiffness matrix

48. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Size of global stiffness matrix will be 6×6, because total DOF are 6. Joint B is
common in both the elements; therefore elements corresponding to unknown at
joint B (3 and 4) will be added together. 
1 2 3 4 5 6
0 111 0 333 0 111 0 333 0 0 1
0 333 1 333 0 333 0 667 0 0 2
0 111 0 333 0 2985 0 042 0 1875 0 375 3
0 333 0 667 0 042 0 375 4
0 0 0 1875 0 375 0 1875 0 375 5
0 0 0 37
K = E
5 0 375 6
I
. . . .
. . . .
. . . . . .
. . . .
. . . .
..




   



    


   
2.333 0.5
0.5 1.0

Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = 5 = zero (simple supports)

Step 5: Reduced stiffness matrix
The nonzero joint displacements are 4 and 6. Therefore collect the elements
corresponding to 4 and 6 from global stiffness matrix. 
46
2 333 0 5 4
0 5 1 0 6
K = EI
. . ,
. . ,





B
C


Step 6: Element Nodal Load Vector:
The element nodal load vector is obtained by restraining the beam at all supports.
Determine fixed end moments, reactions due to external load and reactions due to
moments. Write down element nodal load vector for both the elements and then
determine reduced element nodal load vector.

49. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

AB
75 1
75 2
75 3
75 4






q

AB
50 3
50 4
50 5
50 6






q 
75 1
75 2
125 3
25 4
50 5
50 6




 
 


 


q
Step 7: Equivalent load vector
Equivalent load vector is opposite to element nodal load vector. Joint forces Fq

25 0 25 4
50 0 50 6
     
       
     
F

Step 8: Equation of equilibrium: 
B
C
KF
2 333 0 5 25
EI
0 5 1 0 50
..
..



   
   
    
BC
50
0 0 and
EI
.

Step 9: Reactions and Moments: 
A
A
B
B
C
C
f K q
R 0 111 0 333 0 111 0 333 0 0
M 0 333 1 333 0 333 0 667 0 0
R 0 111 0 333 0 2985 0 042 0 1875 0 375
M 0 333 0 667 0 042 2 333 0 375 0 5
R 0 0 0 1875 0 375 0 1875 0 375
M 0 0 0 375 0 5 0 375 1 0
. . . .
. . . .
. . . . . .
EI
. . . . . .
. . . .
. . . .
  
 



   



   


0 75
0 75
0 1251
0 25EI
0 50
50 50
    
     
     
     
   

   
   
   
     
A
A
B
B
C
C
R 0 75 75
M 0 75 75
R 18 75 125 106 25
M 25 25 0
R 18 75 50 31 25
M 50 50 0
kN
kN.m
. . kN
kN.m
. . kN
kN.m
      
      
      
      
         

       
       
       
     

50. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: Analyse the continuous beam as shown in figure using finite element
method. Take EI constant.

Solution:

Step 1: Degrees of Freedom: 06
No. of elements: 02 (AB, BC)
For simplicity, convert overhang portion into moment.
Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2=3= zero (Fixed support)
2 2-3 3,4,5,6 3=5=zero (simple supports)

Step 2: Element stiffness matrix
Using standard stiffness matrix of beam element, obtain stiffness matrix of each
element separately.
Stiffness matrix element AB 
AB
1 2 3 4
0 096 0 24 0 096 0 24 1
0 24 0 8 0 24 0 4 2
0 096 0 24 0 096 0 24 3
0 24 0 4 0 24 0 8
E
4
K = I
. . . .
. . . .
. . . .
. . . .




  



Stiffness matrix of element BC 
BC
1 2 3 4
0 1875 0 375 0 1875 0 375 1
0 375 1 0 0 375 0 5 2
0 1875 0 375 0 1875 0 375 3
0 375 0 5 0 375 1 0
K = E
4
I
. . . .
. . . .
. . . .
. . . .




  



51. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 3: Global Stiffness matrix:
Size of global stiffness matrix will be 6×6, because total DOF are 6. Joint B is
common in both the elements; therefore elements corresponding to unknown at
joint B (3 and 4) will be added together. 
1 2 3 4 5 6
0 096 0 24 0 096 0 24 0 0 1
0 24 0 8 0 24 0 4 0 0 2
0 096 0 24 0 2835 0 135 0 1875 0 375 3
0 24 0 4 0 135 0 375 4
0 0 0 1875 0 375 0 1875 0 375 5
0 0 0 375
K = EI
0 375 6
. . . .
. . . .
. . . . . .
. . . .
. . . .
..




   



    


   
1.8 0.5
0.5 1.0

Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = 5 = zero (simple supports)
Step 5: Reduced stiffness matrix
The nonzero joint displacements are 4 and 6. Therefore collect the elements
corresponding to 4 and 6 from global stiffness matrix. 
. . ,
. . ,
46
1 8 0 5 4
0 5 1
K=
0
EI
6
B
C







Step 6: Element nodal load vector
  
AB BC
17 6 1
17 6 1 60 3 24 2
24 2 40 4 92 4 3
32 4 3 60 5 4 0 4
36 4 40 6 60 5
4
qq
0
q
6
.
.
.
&
..


   

   
    
     
     
          


  


Step 7: Equivalent load vector

52. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Equivalent load vector is opposite to element nodal load vector. For simplicity
convert overhang portion into moment. (20×1.5=30kN.m clockwise) acting at joint
c. This joint moment will be considered in equivalent load vector directly. Joint forces Fq

4 0 4 0 4
40 30 10 6
.     
     
     
f

Step 8: Equation of Equilibrium: 
B
C
KF
1 8 0 5 4 0
EI
0 5 1 0 10
. . .
..



   
   
   
BC
5 806 12 903
and
EI EI
..




Step 9: Moments and Reaction Calculation f K q  
A
A
B
B
C
C
R 0 096 0 24 0 096 0 24 0 0
M 0 24 0 8 0 24 0 4 0 0
R 0 096 0 24 0 2835 0 135 0 1875 0 375 1
M 0 24 0 4 0 135 1 8 0 375 0 5 EI
R 0 0 0 1875 0 375 0 1875 0 375
M 0 0 0 375 0 5 0 375 1 0
. . . .
. . . .
. . . . . .
. . . . . .
. . . .
. . . .
 
 

 
   
 

 
    
 
 
EI
0 17 6
0 24
0 92 4
5 806 4 0
0 60
12 903 40
.
.
..
.
   
   
   
   
   

   
   
   
   
A
A
B
B
C
C
R 16 207
M 21 68
R 96 46
M 0
R 57 337
M 30
. kN
. kN.m
. kN
kN.m
. kN
kN.m
 
 
 
 
   
   
   
   



Example 3: Analyse the beam using finite element method if support B sink by
25mm. Take EI = 3800 kN.m
2


Solution: Step 1: Degrees of Freedom: 06 and No. of elements: 02 (AB, BC)

53. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603



Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2=3= zero (Fixed support)
2 2-3 3,4,5,6 3=5=zero (simple supports)

Step 2: Element stiffness matrices 
AB
1 2 3 4
0 0555 0 1667 0 0555 0 1667 1
0 1667 0 667 0 1667 0 333 2
0 0555 0 1667 0 0555 0 1667 3
0 1667 0 333 0 1667 0 66
K = EI
74
. . . .
. . . .
. . . .
. . . .




  



BC
3 4 5 6
0 0555 0 1667 0 0555 0 1667 3
0 1667 0 667 0 1667 0 333 4
0 0555 0 1667 0 0555 0 1667 5
0 1667 0 333 0 1667 0 66
K = EI
76
. . . .
. . . .
. . . .
. . . .




  



Step 3: Global Stiffness matrix 
1 2 3 4 5 6
0 0555 0 1667 0 0555 0 1667 0 0 1
0 1667 0 667 0 1667 0 333 0 0 2
0 0555 0 1667 0 111 0 0 0555 0 1667 3
0 1667 0 333 0 0 1667 4
0 0 0 0555 0 1667 0 0555 0 1667
0 0 0 1667 0 166
K
7
= EI
. . . .
. . . .
. . . . .
. . .
. . . .
..




  



   


1.33 0.333
0.333 0.667
5
6




   

Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = 5 = zero (simple supports)

54. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 5: Reduced stiffness matrix 
46
1 333 0 333 4
0 333 0 667
K = EI
6
B
C
. . ,
. . ,






Step 6: Element nodal load vector:
 
11AB BC
30 1 22 22 3
30 2 26 67 4
qq
30 3 7 78 5
30 4 13 33 6
.
.
.
.
   
   
   
   
   
      

Sinking Moments: Sinking is given in mm. Put this in m while calculating sinking
moments. Since both the element are having same length. Sinking moment of both
the elements will be same. 22
6 6 3800 0.025
Sinking moments 15.833 .
6
EI
kN m
L
  
  

 
AB BC
5 28 1 5 28 3
15 833 2 15 833 4
qq
5 28 3 5 28 5
15 833 4 15 833 6
..
..
..
..
   
   
   
   

   
       

35 28 1
45 833 2
41 66 3
q
3 33 4
13 06 5
29 163 6
.
.
.
.
.
.




 



 


Step 7: Equivalent load vector Joint forces Fq

3 33 4
F
29 163 6
.
.




55. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 8: Equation of Equilibrium: KF
B
C
1 333 0 333 3 333
EI
0 333 0 667 29 163
. . .
. . .


   
   
    
BC
9 45 48 189
EI EI
..





Step 9: Moments and Reaction Calculation 
A
A
B
B
C
C
f K q
R 0 0555 0 1667 0 0555 0 1667 0 0
M 0 1667 0 667 0 1667 0 333 0 0
R 0 0555 0 1667 0 111 0 0 0555 0 1667
M 0 1667 0 333 0 1 333 0 1667 0 333
R 0 0 0 0555 0 1667 0 0555 0 1667
M 0 0 0 1667 0 333
. . . .
. . . .
. . . . .
. . . . .
. . . .
..
  




   



   


EI
0 35 28
0 45 833
0 41 661
9 45 3 33EI
0 13 06
0 1667 0 667 48 189 29 163
.
.
.
..
.
. . . .
     
     
     
         
    

    
    
    
         
A
A
B
B
C
C
R 33 704
M 42 686
R 49 693
M 0 00
R 6 602
M 0 00
. kN
. kN.m
. kN
. kN.m
. kN
. kN.m
 
 
 
 
   
   
   
   
 


Example 4: A continuous beam ABC is loaded as shown in fig. It has constant
flexural rigidity. Fixed support at A, roller support at B and guided support at C.
Analyze the beam using finite element method.


Step 1: Degrees of Freedom: 06

Note: Guided support is having only vertical displacement. (Rotation is always
zero.) Therefore reaction at guided support is zero

56. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2=3= zero (Fixed support)
2 2-3 3,4,5,6 3=zero (simple support)
6=zero (guided support)

Step 2: Element stiffness matrices 
AB
1 2 3 4
0 0234 0 0937 0 0234 0 0937 1
0 0937 0 5 0 0937 0 25 2
0 0234 0 0937 0 0234 0 093
K =
73
0 0937 0 25 0 0937 0 5
E
4
I
. . . .
. . . .
. . . .
. . . .




  



BC
3 4 5 6
0 0234 0 0937 0 0234 0 0937 3
0 0937 0 5 0 0937 0 25 4
0 0234 0 0937 0 0234 0 093
K =
75
0 0937 0 25 0 0937 0 5
E
6
I
. . . .
. . . .
. . . .
. . . .




  



Step 3: Global Stiffness matrix 
1 2 3 4 5 6
0 0234 0 0937 0 0234 0 0937 0 0 1
0 0937 0 5 0 0937 0 25 0 0 2
0 0234 0 0937 0 0468 0 0 0234 0 0937 3
0 0937 0 25 0 1874 0 25 4
0 0 0 0234 0 0937 5
0 0 0 09
K =
37 0 25 0 0937 0 5 6
EI
. . . .
. . . .
. . . . .
. . . .
..
. . . .




  


 


1.0 -0.0937
-0.0937 0.0234


   

Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = zero (simple supports), 6 = zero (guided support)

Step 5: Reduced stiffness matrix 
45
1 0 0 0937 4
0 0937 0 0234 5
K = EI
. . ,
. . ,



B
C

57. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 6: Element nodal load vector:

AB
20 1
40 2
20 3
40 4
q








& 
BC
10 3
20 4
10 5
20 6
q







Step 7: Equivalent load vector Joint forces Fq

20 4
10 5





f

Step 8: Equation of Equilibrium: KF
B
C
20 1 0 0 0937
EI
10 0 0937 0 0234
..
..
    
   
    
BC
32 078 555 8
and
EI EI
..


  


Example 4: Analyze the continuous beam using finite element method. Take EI
constant

Solution: Step 1: Degrees of Freedom: 06
DOF at point C are 02 (rotation and translation due to spring).

58. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2=3= zero (Fixed support)
2 2-3 3,4,5,6 3=zero (simple support)
3 3-4 5, 8 8=zero (spring fixed at bottom)

Step 2: Element stiffness matrix 
AB
1 2 3 4
0 1875 0 375 0 1875 0 375 1
0 375 1 0 0 375 0 5 2
0 1875 0 375 0 1875 0 375 3
0 375 0 5 0 375 1 0
K = E
4
I
. . . .
. . . .
. . . .
. . . .




  



BC
3 4 5 6
1 5 1 5 1 5 1 5 3
1 5 2 0 1 5 1 0 4
1 5 1
K=
5 1 5 1 5 5
1 5 1 0 1 5 2 0
EI
6
. . . .
. . . .
. . . .
. . . .




  



Stiffness matrix of spring element 
11
5
5
K = EI
8
1 1 8
CD




Step 3: Global Stiffness matrix 
0 1875 0 375 0 1875 0 375 0 0
0 375 1 0 0 375 0 5 0 0
0 1875 0 375 1 6875 1 125 1 5 1 0
K = EI
0 375 0 5 1 125 3 0 1 5 1 0
0 0 1 5 1 5 2 5 1 5
0 0 1 5 1 0 1 5 2 0
1 2 3 4 5 6
1
2
3
4
5
6
. . . .
. . . .
. . . . . .
. . . . . .
. . . .
. . . .




  



   





  

Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = zero (simple supports)

59. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 5: Reduced stiffness matrix 
4 5 6
3 0 1 5 1 0 4
1 5 2 5 1 5 5
1 0 1 5

2
E
6
KI
0
=
B
C
C
. . . ,
. . . ,
. . . ,




  

 


Step 6: Element nodal load vector:

AB
20 1
20 2
20 3
20
q
4







& 
BC
03
04
05
06
q







 
20 1
20 2
20 3
20 4
05
0
q
6




 







Step 7: Equivalent load vector Joint forces Fq

Note- External moment 30 kN.m clockwise is acting at B, it is accounted in the
element corresponding to rotation at B i.e. 4

20 30 10 4
0 0 0 5
0 0 0 6
     
     
       
     
     
f

Step 8: Equation of Equilibrium: KF
B
C
C
10 3 0 1 5 1 0
0 EI 1 5 2 5 1 5
0 1 0 1 5 2 0
. . .
. . .
. . .


    
   
      

   
 
    
B C C
4 782 2 608 0 434
and
EI EI EI
. . .
,

   

60. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 5: Analyze the indeterminate beam as shown in figure using finite
element method. The beam is fixed at A, C and has internal hinge at B. Take EI
constant.

Solution: Step 1: Degrees of Freedom: 07
DOF at point B are 03 (two rotations and translation).

Discretization
Element Nodes Displacements Boundary conditions
1 1-2 1,2,3,4 1=2= zero (Fixed support)
2 2-3 3,5,6,7 6=7=zero (Fixed support)

Step 2: Element stiffness matrix 
AB
1 2 3 4
12 6 12 6 1
6 4 6 2 2
12 6 12 6 3
6
K = EI
2 6 4 4




  


 
BC
3 5 6 7
1 5 1 5 1 5 1 5 3
1 5 2 0 1 5 1 0 5
1 5 1
K =
5 1 5 1 5 6
1 5 1 0 1 5
EI
2 0 7
. . . .
. . . .
. . . .
. . . .




   




61. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 5: Reduced stiffness matrix 
13 5 6 0 1 5
K = EI 6 0 4 0 0
1 5 0 2 0
3 4 5
3
4
5
B
BA
BC
. . .
..
. ,.
,
,











Step 6: Element nodal load vector:

AB
30 1
52
30 3
54
q







& 
BC
60 3
20 5
60 6
20
q
7






 
90 3
q 5 4
20 5






Step 7: Equivalent load vector Joint forces Fq

90 3
f 5 4
20 5







Step 8: Equation of Equilibrium: KF

13 5 6 0 1 5
6 0 4 0 0
90
5
1 250 020
. . .
..
..







 
   


   
   


B
BA
BC
EI
BA BC
20 28 75 5
and
EI EI EI
B
.
,

   

62. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603



Example: For the following beam, find the vertical deflection and rotation at joint
B using finite element method. Take EI = 12×10
3
kN.m
2


Example: Determine the unknown joint displacements of the beam as shown in
figure using finite element method. Take EI constant.

Example: Analyse the beam using finite element method if support B is sink by
25mm. Take EI = 3800 kN.m
2


Example: A continuous beam has fixed support at node 1 and roller supports at
nodes 2 and 3. Analyse the beam using finite element method and draw SFD and
BMD. Take E = 200 GPa and I=4×10
6
mm
4
.

Example: Obtain rotation at B for the beam shown below using finite element
method. Consider given beam as one element. Take E = 2×10
8
kN/m
2
and I =
4×10
-6
m
4
.



Example: Analyze the continuous beam ABC as shown in Figure using finite
element method. Take EI constant.

63. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Example: Analyse the beam ABC shown in Figure 1 using finite element method.
AB = 3 m and BC = 6 m. Take EI = constant

Example: Analyse the prismatic beam ABC loaded and supported as shown in
Figure using finite element approach. Support B is sink by 25 mm. Draw SFD and
BMD. Take EI constant.

Example: Determined the prop reaction of the propped cantilever beam AB as
shown in Figure 1 using finite element method. Take EI = constant

Example: Obtain fixed end moment at support A using finite element method.
Take E = 2×10
8
kN/m
2
and I = 4×10
-6
m
4
.



Example: Determine support reactions of continuous beam ABC if support B sink
by 10 mm. Take EI = 6000 kN.m
2
. Use finite element method.

Example: Determine support reactions of continuous beam ABC as shown in
Figure 1 if support B sink by 10 mm. Take EI = 6000 kN.m
2
. Use finite element
method.

64. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Unit-III

Finite Element Analysis of Plane Frames

The plane frame is a combination of plane truss and beam. All members are
connected by rigid joints in case of frame.

Stiffness matrix of frame element in local coordinate system
Let consider a frame element of length L, flexural rigidity EI and axial rigidity AE.
A frame is having three degrees of freedom at each node i.e. displacement in x-
direction, displacement in y-direction and rotation. Therefore the size of stiffness
matrix of frame element is 66 .

D1 = Displacement in x-direction at node A
D2= Displacement in y-direction at node A
D3 = Rotation at node A
D4 = Displacement in x-direction at node B
D5= Displacement in y-direction at node B
D6 = Rotation at node B

To derive the stiffness matrix, give the unit displacements at each node one by one

Unit displacement in x-direction at node 1

Unit displacement in y-direction at node 1

Unit rotation in z-direction at node 1

Unit displacement in x-direction at node 2

65. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Unit displacement in y-direction at node 2

Unit rotation in z-direction at node 2

1 2 3 4 5 6
1
3 2 3 2
2
22
3
4
3 2 3 2
5
22
6
0 0 0 0
0 12 6 0 12 6
0 6 4 0 6 2
0 0 0 0
0 12 6 0 12 6
0 6 2 0 6 4
D D D D D D
DAE / L AE / L
DEI / L EI / L EI / L EI / L
DEI / L EI / L EI / L EI / L
K'
DAE / L AE / L
DEI / L EI / L EI / L EI / L
DEI / L EI / L EI / L EI / L




 



   




Transformation Matrix of Frame Element
In plane frame the members are oriented in different directions and hence it is
necessary to transfer stiffness matrix of individual member from local coordinate
system to global coordinate system. This is performed by using transformation
matrix.

66. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Let consider a frame element at an angle θ with respect to positive x-axis.

D1, D2 and D3 D.O.F. at each node for global co-ordinate system
??????
1
′
, ??????
2
′
and ??????
3
′
D.O.F. at each node for local co-ordinate system

Let the local DOF be expressed into global DOF

At Node 1
'
1 1 2
'
2 1 2
'
33
D Dl D m
D D m D l
DD

  

At Node 2
'
4 4 5
'
5 4 5
'
66
D D l D m
D D m D l
DD

  

In matrix form:
l = cosθ and m = sinθ are direction cosines.
1 2 3 4 5 6
'
11
'
22
'
33
'
44
'
55
'
66
0 0 0 0
0 0 0 0
0 0 1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0 0 1
D D D D D D
DlmD
DmlD
DD
DlmD
DmlD
DD
 
 

 
 
   
   
    
   
   where 
0 0 0 0
0 0 0 0
0 0 1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0 0 1
lm
ml
L
lm
ml







 

 
'
x L x

[L] = Transformation Matrix 
'
x
= Local Displacement Vector x
=Global Displacement Vector

67. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

The stiffness matrix of member is global coordinate system is obtained by
using relation (Taking θ = 90
0
, l=0, m=1)
[�]=[�]
??????
[�
′
][�]

3 2 3 2
22
3 2 3 2
22
'
0 1 0 0 0 0 / 0 0 / 0 0
1 0 0 0 0 0 0 12 / 6 / 0 12 / 6 /
0 0 1 0 0 0 0 6 / 4 / 0 6 / 2 /
0 0 0 0 1 0 / 0 0 / 0 0
0 0 0 1 0 0 0 12 / 6 / 0 12 / 6 /
0 0 0 0 0 1 0 6 / 2 / 0 6 / 4 /
T
L K L
AE L AE L
EI L EI L EI L EI L
EI L EI L EI L EI L
AE L AE L
EI L EI L EI L EI L
EI L EI L EI L EI L

  



 



   


0 1 0 0 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0 1 0
0 0 0 1 0 0
0 0 0 0 0 1
















 

















Therefore, the stiffness matrix of any member which is perpendicular
(θ = 90
0
) to reference member 
3 2 3 2
22
3 2 3 2
22
12 0 6 12 0 6
0 0 0 0
6 0 4 6 0 2
12 0 6 12 0 6
0 0 0 0
6 0 2 6 0 4
EI / L EI / L EI / L EI / L
AE / L AE / L
EI / L EI / L EI / L EI / L
K
EI / L EI / L EI / L EI / L
AE / L AE / L
EI / L EI / L EI / L EI / L
   











68. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Note:
If we neglect the axial deformation these two matrices reduced to order 4×4.
(The columns and row corresponding to axial stiffness AE/L are neglected)



















Steps for the solution of Indeterminate plane frames using finite element
method:
1. Divide the frame into number of elements (Take one member as one element)
2. Identify total degrees of freedom (Three D.O.F. at each node, two
displacements and rotation)
3. Determine stiffness matrices of all elements ([K]1, [K]2………)
4. Assemble the global stiffness matrix [K]
5. Impose the boundary conditions and determine reduced stiffness matrix
6. Determine element nodal load vector [q] (Restrained structure)
7. Determine equivalent load vector [f]
8. Apply equation of equilibrium [K]{Δ}={f} and determine unknown joint
displacements.
9. Apply equation [K]{Δ}+[q] ={f} to determine reactions and moments





Stiffness matrix for Beam Member neglecting axial deformation
(Neglect first and fourth row and columns) 
3 2 3 2
22
3 2 3 2
22
12 6 12 6
6 4 6 2
12 6 12 6
6 2 6 4
EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
K
EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
 




  



Stiffness matrix for Column Member (θ=90
0
, l = 0, m = 1) always take bottom
of column as a first node. (Neglect second and fifth row and columns) 
3 2 3 2
22
3 2 3 2
22
12 6 12 6
6 4 6 2
12 6 12 6
6 2 6 4
EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
K
EI / L EI / L EI / L EI / L
EI / L EI / L EI / L EI / L
   








69. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 1: Analyze the portal frame as shown in figure using finite element
method. Take EI constant. Neglect axial deformation.


Solution: Step 1: Total DOF = 12
(Three DOF at each node, two displacements and one rotation)
No. of elements: 03 (AB, BC, DC)
Discretization
Element Nodes Displacements Boundary conditions
1 AB 1,2,3,4,5,6 1=2=3=zero
2 BC 4,5,6,7,8,9 5=8=zero, 4=7
3 DC 10,11,12,7,8,9 10=11=12=zero

Step 2: Element Stiffness matrices
Element stiffness matrix for AB (Column member) 
   







1 3 4 6
0.048 0.24 0.048 0.24 1
0.24 1.6 0.24 0.8 3
0.04 0.048 08 0.24 4
0.2
.24
0.240 4 1.6.8 6
AB
K EI

Imposing Boundary Conditions
1=3=0
Element Stiffness Matrix for BC: (Beam member)

70. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

5 6 8 9
0.75 1.5 0.75 1.5 5
1.5 1.5 6
0.75 1.5 0.75 1.5 8
1.5
42
241.5 9
BC
K EI




   




Imposing Boundary Conditions
5=8=0
Element Stiffness Matrix for DC: (Column member) 
DC
K EI
   







10 12 7 9
0.096 0.24 0.096 0.24 10
0.24 0.8 0.24 0.4 12
0.09 0.096 0.24
0.24
6 0.24 7
0.24 0. 9 0.84

Imposing Boundary Conditions
10=12=0
Step 3: Reduced Stiffness Matrix:
Since horizontal sway at B and C are same (4=7), we can modify the above
stiffness matrix as 4 6 9
0.144 0.24 0.24 4,
[ ] 0.24 5.6 2 6,
0.24 2 4.8 9,
B
C
K EI 









Step 4: Element Nodal Load Vector

71. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

3.52 1
9.6 3
{}
4
6
65
6
{}
68
9
3.24 10
3.6 12
6.48
14.4
4
4
1.76
}
29.
{
4
7
AB
BC
DC
q
q
q


























Reduced element nodal load vector 4.72 4
{ } 10.4 6
1.6 9
q







Step 4: Equivalent Load Vector Joint forcesFq  

4.72 4
10.4 6
1.6 9
F






Step 5: Equation of Equilibrium [ ]{ } { }
0.144 0.24 0.24 4.72
0.24 5.6 2 10.4
0.24 2 4.8 1.6
B
C
KF
EI 


   
   
   

   

    

34.046 1.0419 1.803
;;
BC
m rad rad
EI EI EI
    


Step 6: Moment Calculations {f} = [K]{Δ}+{q}
Member AB

72. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
0
0.24 1.6 0.24 0.8 0 9.6 18.6041
0.24 0.8 0.24 1.6 34.046 14.4 4.562
1.0419
AB
BA
M
EI
M EI


       
         
        



Member BC 0
1.5 4 1.5 2 1.0419 4 4.5621
1.5 2 1.5 4 0 4 9.128
1.803
BC
CB
M
EI
M EI


       
         
        


Member DC 0
0.24 0.8 0.24 0.4 0 3.6 3.8491
0.24 0.4 0.24 0.8 34.046 2.4 9.128
1.803
DC
CD
M
EI
M EI


       
         
      


Example 2: Analyze the rigid frame by using finite element method. Take EI
constant. Neglect axial deformation.


Solution: Step 1: Total DOF = 12
(Three DOF at each node, two displacements and one rotation)
No. of elements: 03 (AB, BC, DC)

Discretization
Element Nodes Displacements Boundary conditions
1 AB 1,2,3,4,5,6 1=2=3=zero
2 BC 4,5,6,7,8,9 5=8=zero, 4=7
3 DC 10,11,12,7,8,9 10=11=12=zero

73. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 2: Element stiffness matrix for Column AB: 
1 3 4 6
0.1875 0.375 0.1875 0.375 1
0.375 1 0.375 0.5 3
0.1875 0.375 4
0.37
0.1875 0.375
0.375 15 0.5 6
AB
K EI
   








Element Stiffness Matrix of beam BC 
5 6 8 9
0.0469 0.1875 0.0469 0.1875 5
0.1875 0.1875 6
0.0469 0.1875 0.0469 0.1875 8
0.1875 0.187
1 0.5
0.5 1 59
BC
K EI




   




Element Stiffness Matrix for column DC 
10 12 7 9
0.1875 0.375 0.1875 0.375 10
0.375 1 0.375 0.5 12
0.1875 0.375 7
0.375 0
0.1875 0.375
0.375 1.5 9
DC
K EI
   








Imposing Boundary Conditions
1=2=3=5=8=10=11=12=0

Step 2: Reduced Stiffness Matrix
Since horizontal sway at B and C are same (4=7), we can modify the above
stiffness matrix as 4 6 9
0.375 0.375 0.375 4,
[ ] 0.375 2 0.5 6,
0.375 0.5 2 9,
B
C
K EI 









Step 3: Element Nodal Load Vector

74. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

0 1 37.5 5
0 3 6
{ } { }
4 37.5 8
0 6 9
0 10
0 12
{}
7
9
75
0
75
0
0
AB BC
DC
qq
q
   
   
   
   
   
   
   










Reduced element nodal load vector 
04
75 6
75 9
q







Step 4: Equivalent Load Vector Joint forcesFq  

0 50 50 4
75 0 75 6
75 0 75 9
F
     
     
         
     
     

Step 5: Equation of Equilibrium
0.375 0.375 0.375 50
0.375 2 0.5 75
0.375 0.5 2 75
B
C
EI 

   
   
   

   

    

??????
�=
−��.���
????????????
?????????????????? ??????
�=
��.���
????????????
?????????????????? ∆=
���.���
????????????
??????

Step 6: Moment Calculations

{f} = [K]{Δ}+{q}
Member AB 0
0.375 1 0.375 0.5 0 0 32.1431
0.375 0.5 0.375 1 190.476 0 7.143
78.571
AB
BA
M
EI
M EI


       
         
      


75. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Member BC 1 0.5 6 78.571 75 7.1431
0.5 1 9 21.428 75 92.857
BC
CB
M
EI
M EI
        
         
       

Member DC 0
0.375 1 0.375 0.5 0 0 82.1431
0.375 0.5 0.375 1 190.476 0 92.857
21.428
DC
CD
M
EI
M EI


       
         
      




Example 3: Assemble the stiffness matrix for the frame shown in Fig. using finite
element method. Take AE = 400000KN and EI = 1000 kNm
2
for both the members

Solution: Total DOF = 09 (03 at each node)
No. of elements: 02
Since c/s area of elements is given, we can not neglect axial deformation.
Therefore size of stiffness matrix is 6X6.


Element Stiffness Matrix of AB ( θ=90
0
, l =0 and m=1. Origin at A) 0 1 0 0 0 0 100 0 0 100 0 0
1 0 0 0 0 0 0 0.3 0.6 0 0.3 0.6
0 0 1 0 0 0 0 0.6 1.6 0 0.6 0.8
[ ] [ '] 1000
0 0 0 0 1 0 100 0 0 100 0 0
0 0 0 1 0 0 0 0.3 0.6 0 0.3 0.6
0 0 0 0 0 1 0 0.6 0.8 0 0.6 1.6
T
LK
   
   

   
    
   

   
      
   
   

76. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

0 0.3 0.6 0 0.3 0.6 0 1 0 0 0 0
100 0 0 100 0 0 1 0 0 0 0 0
0 0.6 1.6 0 0.6 0.8 0 0 1 0 0 0
' 1000
0 0.3 0.6 0 0.3 0.6 0 0 0 0 1 0
100 0 0 100 0 0 0 0 0 1 0 0
0 0.6 0.8 0 0.6 1.6 0 0 0 0 0 1
T
L K L
    
  

  
   
   

  
  
  
  


1 2 3 4 5 6
0.3 0 0.6 0.3 0 0.6 1
0 100 0 0 100 0 2
0.6 0 1.6 0.6 0 0.8 3
1000
0.3 0 0.6 0.3 0 0.6 4
0 100 0 0 100 0 5
0.6 0 0.8 0.6 0 1.6 6
AB
K
  




 


 



Element stiffness matrix of BC

L=5.6569 m, θ = 45
0
, l=0.7071, m=0.7071
'
0.7071 0.7071 0 0 0 0 70.71 0 0 70.71 0 0
0.7071 0.7071 0 0 0 0 0 0.106 0.30 0 0.106 0.30
0 0 1 0 0 0 0 0.30 1.1314 0 0.30 0.565
1000
0 0 0 0.7071 0.7071 0 70.71 0 0 70.71 0 0
0 0 0 0.7071 0.7071 0 0 0.106 0.30 0 0.106 0.3
0 0 0 0 0 1
T
LK 




 



   


0
0 0.30 0.565 0 0.30 1.131










[ ] [ '][ ]
50 0.075 0.2121 50 0.075 0.2121 0.7071 0.7071 0
50 0.075 0.212 50 0.075 0.2121
0 0.3 1.131 0 0.3 0.565
1000
50 0.075 0.2121 50 0.075 0.2121
50 0.075 0.2121 50 0.075 0.2121
0 0.3 0.565 0 0.3 1.131
T
L K L
   



 



   


0 0 0
0.7071 0.7071 0 0 0 0
0 0 1 0 0 0
0 0 0 0.7071 0.7071 0
0 0 0 0.7071 0.7071 0
0 0 0 0 0 1







 



77. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

4 5 6 7 8 9
35.408 35.302 0.2121 35.408 35.302 0.2121
35.302 35.408 0.2121 35.302 35.408 0.2121
0.2121 0.2121 1.131 0.2121 0.2121 0.565
1000
35.408 35.302 0.2121 35.408 35.302 0.2121
35.302 35.408 0.2121 35.302 35.408
BC
K
   




   
4
5
6
7
0.2121 8
0.2121 0.2121 0.565 0.2121 0.2121 1.131 9










Global Stiffness matrix 1 2 3 4 5 6 7 8 9
0.3 0 0.6 0.3 0 0.6 0 0 0
0 100 0 0 100 0 0 0 0
0.6 0 1.6 0.6 0 0.8 0 0 0
0.3 0 0.6 35.708 35.302 0.388 35.408 35.302 0.21 21
[ ] 10000 100 0 35.302 135.408 0.2121 35.302 35.408 0.21 21
0.6 0 0.8 0.388 0.2121 2.731 0.2121 0.2121 0.565
K
  

   
 

1
2
3
4
5
6
0 0 0 35.408 35.302 0.2121 35.408 35.302 0.2121 7
0 0 0 35.302 35.408 0.2121 35.302 35.408 0.2121 8
0 0 0 0.2121 0.2121 0.565 0.2121 0.2121 1.131 9












   




Example 4: Determine global stiffness matrix of the frame ABC shown in figure
using finite element method. Take EI constant. Neglect axial deformation.

Example 5: Analyse the frame shown in Figure using finite element method and
draw bending moment diagram. Neglect axial deformation.
(Ans.B AB BA BC CB
1.2/ EI, M 0.8kN.m,M 1.6kN.m, M 18.4kN.m, M 14.8 kN.m      )

78. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Example 6: Determine the rotation of joint B, and the horizontal displacements of
joints B and C. Take EI = 10 ×10
3
KN.m
2
. Neglect axial deformations.

Example 7: Analyze the rigid jointed portal frame shown in Figure using finite
element method. Take EI constant. Neglect axial deformation.

Example 8: Determine the unknown joint displacements of the portal frame as
shown in Figure using finite element method. Take EI constant. Neglect axial
deformation.

Example 9: Analyse the portal frame as shown in Figure using finite element
method. Neglect axial deformation.

79. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


Example 10: Determine the unknown joint displacements of the portal frame as
shown in Figure using finite element method. Take EI constant. Neglect axial
deformation.

Example 11: Derive the stiffness matrix of portal frame ABC as shown in figure
using finite element method. Neglect axial deformation.

Example 12: Analyze the rigid jointed portal frame shown in Figure 6 using finite
element method. Take EI constant. Draw BMD. Neglect axial deformation.

80. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Analysis of Grid Structures

The property of grid member is basically a combination of 2D beam with twisting
effect. The plane frame is subjected to tangential load (loaded in its own plane)
whereas grid is subjected to load perpendicular to its plane. As a result twisting
effects are included in the grid analysis. Thus grid structures withstand bending
moment, shear force and twisting moment.

Stiffness matrix for grid element
The total degrees of freedom at each node of the grid are 03 (vertical displacement,
bending rotation and twisting rotation). Therefore size of stiffness matrix is 6×6.

D1 = Displacement in y-direction at node 1
D2= Rotation in x-direction at node 1
D3 = Rotation in z-direction at node 1
D4 = Displacement in y-direction at node 2
D5= Rotation in x-direction at node 2
D6 = Rotation in z-direction at node 2

To derive the stiffness matrix, give the unit displacements at each node one by one

Unit displacement in y-direction at node 1


Unit rotation in x-direction at node 1



Unit rotation in z-direction at node 1
Moments in this figure are about
z-axis, no moment about x-axis.
Reactions are along y direction
Moments in this figure are about x-
axis i.e. twisting moment. Reactions
along y direction and moments along
z-direction are zero

81. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603



Unit displacement in y-direction at node 2


Unit rotation in x-direction at node 2




Unit rotation in z-direction at node 2


Stiffness matrix for grid element in local coordinate system from above six figures 
EI L EI L EI L EI L
GJ L GJ L
EI L EI L EI L EI L
K
EI L EI L EI L EI L
GJ L GJ L
EI L
1 2 3 4 5 6
3 2 3 2
22
3 2 3 2
2
D D D D D D
12 / 0 6 / 12 / 0 6 /
0 / 0 0 / 0
6 / 0 4 / 6 / 0 2 /
12 / 0 6 / 12 / 0 6 /
0 / 0 0 / 0
6/




  

EI L EI L EI L
2
0 2 / 6 / 0 4 /













Moments in this figure are about
z-axis, no moment about x-axis.
Reactions are along y direction
Moments in this figure are about
z-axis, no moment about x-axis.
Reactions are along y direction
Moments in this figure are about x-
axis i.e. twisting moment. Reactions
along y direction and moments along
z-direction are zero
Moments in this figure are about
z-axis, no moment about x-axis.
Reactions are along y direction

82. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Transformation Matrix of grid Element
In grid the members are orthogonally and hence it is necessary to transfer stiffness
matrix of individual member from local coordinate system to global coordinate
system. This is performed by using transformation matrix.

Let consider a grid element at an angle θ with respect to positive x-axis in
clockwise sense.

??????
1
′
, ??????
2
′
and ??????
3
′
= D.O.F. at each node for local co-ordinate system
D1, D2 and D3 = D.O.F. at each node for global co-ordinate system
Let the local DOF be expressed into global DOF

At Node 1
'
11
'
2 2 3
'
3 2 3
cos sin
sin cos
DD
D D D
D D D




  


At Node 2
'
44
'
5 5 6
'
6 5 6
cos sin
sin cos
DD
D D D
D D D




  



In matrix form: (l = cosθ and m = sinθ are direction cosines)

83. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603


'
11
'
22
'
33
'
44
'
55
'
66
1 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0
0 1 0 0 0 0 1 0 0 0
0 0 0 1 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
DD
Dl m l mD
DmmD
L
DD
Dl m l mD
Dm l m lD
    
    
    
    
        
      
      
      
     

'
x L x

[L] = Transformation Matrix 
'
x
= Local Displacement Vector x
=Global Displacement Vector
The stiffness matrix of member is global coordinate system is obtained by
using relation (θ=90
0
, l = 0, m = 1)
[�]=[�]
??????
[�
′
][�]

3 2 3 2
22
1
3 2 3 2
22
12 6 12 6
00
0 0 0 01 0 0 0 0 0
0 0 1 0 0 0
6 4 6 2
00
0 1 0 0 0 0
0 0 0 1 0 0 12 6 12 6
00
0 0 0 0 0 1
0 0 0 0 1 0
0 0 0 0
6 2 6 4
00
T
EI EI EI EI
L L L L
GJ GJ
LL
EI EI EI EI
L L L L
LK
EI EI EI EI
L L L L
GJ GJ
LL
EI EI EI EI
L L L L











 


 

  
 


 





















3 2 3 2
22
3 2 3 2
22
12 6 12 6
00
6 4 6 2
00 1 0 0 0 0 0
0 0 1 0 0 0
0 0 0 0
0 1 0 0 0 0
12 6 12 6 0 0 0 1 0 0
00
0 0 0 0 0 1
6 2 6 4
0 0 0 0 1 0
00
0 0 0 0
EI EI EI EI
L L L L
EI EI EI EI
L L L L
GJ GJ
LL

EI EI EI EI
L L L L
EI EI EI EI
L L L L
GJ GJ
LL





  



 

 




  


  













 

84. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

EI L EI L EI L EI L
EI L EI L EI L EI L
GJ L GJ L
K
EI L EI L EI L EI L
EI L EI L EI L EI L
GJ L GJ L
3 2 3 2
22
3 2 3 2
22
12 / 6 / 0 12 / 6 / 0
6 / 4 / 0 6 / 2 / 0
0 0 / 0 0 /
12 / 6 / 0 12 / 6 / 0
6 / 2 / 0 6 / 4 / 0
0 0 / 0 0 /
   



 








Example 1: Orthogonal grid is in xz plane. It consists of two prismatic members
having same EI and GJ as well as length L. Develop stiffness matrix for grid using
finite element method.

Degrees of freedom and boundary conditions
Solution: Total DOF: 09 (03 at each node)

Step 1: Discretization

Element Nodes Displacements Boundary conditions
1 AB 1,2,3,4,5,6 1=2=3=zero
2 BC 4,5,6,7,8,9 7=8=9=zero

Step 2: Element stiffness matrices
Note: First select the element to which standard stiffness matrix is applicable.
a) Direction of element must be towards right
b) Perpendicular direction must approach the observer
Stiffness matrix for member AB: (Since standard element is along x-axis,
assume twisting rotation along x and bending rotation along z)

85. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

3 2 3
22
32 32
2
1 2 3 4 5 6
12 / 0 6 / 12 / 0 6 /
0 / 0 0 / 0
6 / 0 4 / 6 / 0 2 /
12 / 0 612 / 0 6 /
0
/
0 / 0
AB
EI L EI L EI L EI L
GJ L GJ L
EI L EI L EI L
EI L E
EI
IL
L
K
EI L EI L
GL GJJ







2 2
1
2
3
4
5
66 / 0 2 /

/0
64

0/

/
L
EI L EEI L EI L IL




 









Stiffness matrix for member CB:
Rotate member CB in clockwise
direction about joint C at an angle of 90
0

w. r. t to positive x-axis. (θ=90
0
, l=0,
m=1). Therefore take ‘C’ as a first node
and ‘B’ as a second node.


3
3
2 3 2
22
32 2
7 8 9 4 5 6
12 / 6 / 0 12 / 6 / 0
6 / 4 / 0 6 / 2 / 0
0 0 / 0
1
0/
12 2 / 6 / 0/ 6 / 0
CB
EI L EI L EI L EI L
EI L EI L EI L EI L
GJ L GJ L
K
EI EI LL EI EIL L
  




2 2
7
8
9
4
56 / 2 / 0
60 0 /

6 / 4 / 0

0/

0
EI L EI L
GJ L
EI L EI L
GJ L




 






  

III) Reduced stiffness matrix    
  
3 2 2
2
2
4 5 6
24 / 6 / 6 / 4,
6 / 4 / / 0 5,
6 / 0 4 / / 6,
Bz
Bx
By
EI L EI L EI L
K EI L EI L GJ L
EI L EI L GJ L


 






86. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 2: Analyze the grid structure as shown in figure using finite element
method. Take GJ = 0.4 EI

Degrees of freedom and boundary conditions
Solution: Total DOF: 09 (03 at each node)
No. of elements: 02 (AB and CB)
Step 1: Discretization

Element Nodes Displacements Boundary conditions
1 AB 1,2,3,4,5,6 1=2=3=zero
2 CB 4,7,8,9,5,6 7=8=9=zero

Note: First select the element to which standard stiffness matrix is applicable.
a) Direction of element must be towards right
b) Perpendicular direction must approach the observer
Standard stiffness matrix for member AB: (Since standard element AB is along
x-axis, assume twisting rotation along x and bending rotation along y)

Step 2: Element stiffness matrices
Stiffness matrix of element AB (Using standard stiffness matrix) 
1 2 3 4 5 6
0.096 0 0.24 0.096 0 0.24
0 0.08 0 0 0.08 0
0.24 0 0.8 0.24 0 0.4
0.096 0 0.24 0.096 0 0.24
0 0.08 0 0 0.08 0
0
AB
K EI




  

1
2
3

4
5
.24 0 0.4 0.24 0 0.8 6





 







87. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix for member BC:
Rotate member BC in clockwise
direction about joint C at an angle of 90
0

w. r. t to positive x-axis. (θ=90
0
, l=0,
m=1). Therefore take ‘B’ as a first node
and ‘C’ as a second node.

4 5 6 7 8 9
0.444 0.667 0 0.444 0.667 0
0.667 1.333 0 0.667 0.667 0
0 0 0.133 0 0 0.133

0.444 0.667 0 0.444 0.667 0
0.667 0.6
BC
K EI
  





4
5
6

7
67 0 0.667 1.333 0 8
0 0 0.133 0 0 0.133 9








 




Reduced stiffness matrix is 
4 5 6
0.540 0.667 0.24 4,
0.667 1.413 0 5,
0.24 0 0.933 6,
Bz
Bx
By
K EI 

  






Step 3: Element nodal load vector
Since there is no load on the members, fixed end moments and corresponding
reactions are zero. Therefore element nodal load vector is null matrix.
  
0 1 0 4
0 2 0 5
04
0 3 0 6
05
0 4 0 7
06
0 5 0 8
0 6 0 9
AB BC
q and q q
   
   
   

    
        
     

   
   
   


Step 4: Equivalent load vector Joint forcesfq  

88. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

0 70 70 4
0 0 0 5
0 0 0 6
f
     
     
       
     
     

Step 5: Equation of Equilibrium
[K]{Δ} = {f}
0.540 0.667 0.24 70
0.667 1.413 0 0
0.24 0 0.933 0
Bz
Bx
By
EI 

      
   
    

   

    
Bz Bx By
EI EI EI
428.371 202.21 110.192
;;
  
   


Step 6: Moment Calculations
Element AB 
AB AB ABAB
K q f  
0
0 0.08 0 0 0.08 0 0 0
0.24 0 0.8 0.24 0 0.4 0 0 1
0 0.08 0 0 0.08 0 428.371 0
0.24 0 0.4 0.24 0 0.8 202.21 0
110.192
Ax
Ay
Bx
By
M
M
EI
MEI
M


    

   
    
      
 
     

         



Element BC
428.371
0.667 1.333 0 0.667 0.667 0 202.21 0
0 0 0.133 0 0 0.133 110.192 0 1
0.667 0.667 0 0.667 1.333 0 0 0
0 0 0.133 0 0 0.133 0 0
0
Bx
By
Cx
Cy
M
M
EI
MEI
M


      

   
    
      

     

         


Ax
Ay
Bx
By
M
M
M
M
16.176
58.732
16.176
14.655
 
 
   
   

   
   

and Bx
By
Cx
Cy
M
M
M
M
16.176
14.655
150.84
14.655
 
 
   
   
   
   


89. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 3: Analyze the balcony grid as shown in figure using finite element
method. Take EI = 1600 kN.m
2
and GJ = 800 kN.m
2

Degrees of freedom and boundary conditions
Solution:
Step 1: Total DOF = 09 (03 at each node)
No. of elements: 02 (AB and BC)

Step 1: Discretization

Element Nodes Displacements Boundary conditions
1 AB 1,2,3,4,5,6 1=2=3=zero
2 BC 4,7,8,9,5,6 7=8=9=zero

Step 2: Element stiffness matrices
Note:
1) According to standard derivation, element BC is satisfying conditions of
standard derivation (direction along the member is along right and
perpendicular direction approaching the observer). Therefore, standard
stiffness matrix is applicable to BC and 90
0
matrix is applicable to member
AB.
2) Since BC is along y-direction assume rotation in y-direction (By
 ) as second
unknown and x-direction rotation (Bx
 ) as third unknown.
3) While measuring angle for second member, measure w.r.t. positive y-axis.

90. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix of element BC 
4 5 6 7 8 9
300 0 600 300 0 600 4
0 200 0 0 200 0 5
600 0 1600 600 0 800 6
300 0 600 300 0 600 7
0 200 0 0 200 0 8
600 0 800 600 0 1600 9

BC
K




 

   

 


  


Stiffness matrix for member AB:
Rotate member AB in clockwise
direction about joint A at an angle of 90
0

w. r. t to positive y-axis. (θ=90
0
, l=0,
m=1). Therefore take ‘A’ as a first node
and ‘B’ as a second node.


1 2 3 4 5 6
300 600 0 300 600 0 1
600 1600 0 600 800 0 2
0 0 200 0 0 200 3
300 600 0 300 600 0 4
600 800 0 600 1600 0 5
0 0 200 0 0 200 6

AB
K
   



 






  

Reduced stiffness matrix is 
4 5 6
600 600 600 4,
600 1800 0 5,
600 0 1800 6,
Bz
By
Bx
K 








91. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 3: Element nodal load vector: (Fixed end moments and reactions)

 

20 4 60 1
0 5 40 2
20 6 0 3
and
20 7 60 4
0 8 40 5
20 9 0 6
80 4,
40 5,
20 6,
BC AB
Bz
By
Bx
qq
R
qM
M
   
   

   
    
   

   
   
   
   







Step 4: Equivalent load vector Joint forcesfq  

80 4
40 5
20 6
f







Step 5: Equation of Equilibrium
[K]{Δ} = {f} 600 600 600 80
600 1800 0 40
600 0 1800 20
Bz
By
Bx


   
   
   

   
 
    
0.3 ; 0.0888 ; 0.0777
Bz Bx By
m rad rad    

Step 6: Moment Calculations K q f  

Element AB 0
600 1600 0 600 800 0 0 0
0 0 200 0 0 200 0 40
600 800 0 600 1600 0 0.3 0
0 0 200 0 0 200 0.0888 40
0.0777
Ay
Ax
By
Bx
M
M
M
M


    

   
    
      

     

         



92. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
0.3
0 200 0 0 200 0 0.0888 20
600 0 1600 600 0 800 0.0777 0 1
0 200 0 0 200 0 0 20
600 0 800 600 0 1600 0 0
0
By
Bx
Cy
Cx
M
M
EI
MEI
M


    

   
    
      

     

         



Element BC Ay
Ax
By
Bx
M
M
M
M
17.76
157.84
17.76
15.68
 
 
   
   
   
   
and By
Bx
Cy
Cx
M
M
M
M
17.76
15.68
128.96
15.68
 
 
   
   

   
   


Example 4: Using finite element method, determine unknown joint displacements
of the grid as shown in figure. E=2×10
5
MPa, I = 20×10
5
mm
4
, G = 0.8×10
5
MPa, J
=50×10
5
mm
4


Degrees of freedom and boundary conditions
Solution:
Step 1: Total DOF = 09 (03 at each node)
No. of elements: 02 (BA and BC)
EI = 2×10
5
× 20×10
5
= 40×10
10
N.mm
2
= 400 kN.m
2

GJ = 0.8×10
5
× 50×10
5
= 40×10
10
N.mm
2
= 400 kN.m
2

93. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 1: Discretization

Element Nodes Displacements Boundary conditions
1 BA 4,5,6,1,2,3 1=2=3=zero
2 BC 4,7,8,9,5,6 7=8=9=zero
Note:
1) According to standard derivation, element BA is satisfying conditions of
standard derivation (direction along the member is along right and
perpendicular direction approaching the observer). Therefore, standard
stiffness matrix is applicable to BA and 90
0
matrix is applicable to member
BC.
2) Since BA is along y-direction assume rotation in y-direction (By
 ) as second
unknown and x-direction rotation (Bx
 ) as third unknown.
3) While measuring angle for second member, measure w.r.t. positive y-axis.

4 5 6 1 2 3
600 0 600 600 0 600 4
0 200 0 0 200 0 5
600 0 800 600 0 400 6
600 0 600 600 0 600 1
0 200 0 0 200 0 2
600 0 400 600 0 800 3

BA
K




 

   

 


  



Stiffness matrix for member BC:
Rotate member BC in clockwise
direction about joint B at an angle of 90
0

w. r. t to positive y-axis. (θ=90
0
, l=0,
m=1). Therefore take ‘B’ as a first node
and ‘C’ as a second node.

94. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

4 5 6 7 8 9
177.77 266.67 0 177.77 266.67 0
266.67 533.33 0 266.67 266.67 0
0 0 133.33 0 0 133.33

177.77 266.67 0 0.177 266.67 0
BC
K
  




4
5
6
7
266.6 266.67 0 266.67 533.33 0 8
0 0 133.3 0 0 133.33 9








  


  

Reduced stiffness matrix is 






Bz
ByBC
Bx
K
4 5 6
777.77 266.67 600 4
266.67 733.33 0 5
600 0 933.33 6




Step 3: Element nodal load vector: (Fixed end moments and reactions)
  
0 4 6 4
0 5 3 5
64
0 6 0 6
and 3 5
0 1 6 7
06
0 2 3 8
0 3 0 9
BA BC
q q q
   
   

   

    
         
     

   
   
   

Step 4: Equivalent load vector Joint forcesfq  

95. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

64
35
06
f






Step 5: Equation of Equilibrium
[K]{Δ} = {f}
     
 
    

   

    
Bz
By
Bx
777.77 266.67 600 6
266.67 733.33 0 3
600 0 933.33 0



0.0166 ; 0.00195 ; 0.0106
Bz By Bx
m rad rad     




Example 6: Derive the stiffness matrix for the grid elements as shown in Figure
using finite element method. Take flexural rigidity EI and torsional rigidity GJ
same for both the elements

Example 7: Analyze the grid structure ABC as shown in Figure using finite
element method. Take EI=2×10
5
kN.m
2
and GJ = 1.2×10
5
kN.m
2
.

96. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Example 8: Orthogonal grid ABC is in x-z plane. It consists of two prismatic
members having same EI and GJ as well as length ‘L’. Develop stiffness matrix for
grid using finite element method.


Example 9: Analyze the grid structure ABC as shown in Figure using finite
element method method. Take E = 210 GPa, G = 84 GPa, I = 16.6 X10
-5
m
4
, J =
4.6 X10
-5
m
4
for all elements.

97. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
1

Unit - I
Finite Element Method

1.1 Introduction to FEM
For any structural problem two types of solutions are available; (i) Analytical
solution and (ii) Numerical solutions. Analytical solutions are accurate for the
simple boundary conditions, loading conditions and linear problems. But, for the
complex geometry, irregular boundary conditions and geometric non-linearity,
analytical solutions are not effective and accurate. Therefore, various numerical
methods are developed by the researchers for solving such complex problems. But
these numerical methods give approximate solutions of the problem.
The finite element method (FEM) is a numerical technique used to find out
solutions of complex engineering problems. Originally this method was developed
for the aerospace engineering but, it is now widely used in other engineering
disciplines such as Civil Engineering, Mechanical Engineering and Electrical
Engineering. The first book on finite element method (FEM) was written by
Zienkiewicz in three volumes.

1.2 What is Finite Element Analysis (FEA)?
In finite element analysis, solution of complex problem is obtained by dividing
domain (structure) into ‘n’ number of subdomains (elements). The study of
properties of one element is called as element formulation whereas assembly of
properties of all elements (global study) to obtain solution of problem is called as
system formulation.

1.3 Principles of FEA
1) The finite element method (FEM) is a computational technique used to
obtain approximate solutions of boundary value problems in engineering.
2) Boundary value problems are also called field problems. The field is the
domain of interest and most often represents a physical structure.
3) The field variables are the dependent variables of interest governed by the
differential equation.
4) The boundary conditions are the specified values of the field variables (or
related variables such as derivatives) on the boundaries of the field.

1.4 What is Discretization?
Discretization of structure is an important task in finite element analysis and
requires some skill and knowledge. In this procedure, first, the number, shape, size
and configuration of elements have to be decided in such a manner that the real

98. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
2

structure is simulated as closely as possible. The discretization is to be in such that
the results converge to the true solution.



1. First, the domain (Rectangular slab) is presented
as a collection of finite number ‘n’ of subdomain
i. e. Rectangular element. This is called as
discretization of the domain.
2. Each domain is called as ‘Element’.
3. The collection of element is called as ‘Finite
Element Mesh’.
4. The elements are connected at points called as
‘Nodes’.
5. When elements are of same dimension it is called
as uniform mesh otherwise non-uniform mesh.

1.5 Advantages of FEM over conventional methods
1. For problems involving irregular shape, irregular boundary condition and
irregular loading conditions, conventional methods makes certain
assumptions whereas in FEM no such assumptions are made. The problems
are treated as it is.
2. For anisotropic material properties, solution by classical method is very
difficult. FEM can handle such structures without any difficulty.
3. Material non-linearity and geometric non-linear problems cannot handle by
classical methods. There is no difficulty in FEM to handle such problem.
4. FEM superior to irregular problems, for regular problems classical methods
are best solution

1.6 Disadvantages of FEM
1. In classical method exact solution is obtained whereas in FEM approximate
solution is obtained.
2. This technique depends upon skill designer in assuming element type,
number of nodes and displacement fields etc.

99. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
3

3. Still adequate for several problems such as; cracking behaviour, bond
failure, problems of composite materials.
4. It requires large amount of computer memory and computational time to
obtained results.

1.7 Finite Element Method (FEM) Vs Finite Difference Method (FDM)
F.D.M. F.E.M.
1. FDM gives values at nodes points
only. At other point interpolation is
required
1. FEM gives values at any nodes
points including node points.
2. FDM needs larger number of nodes
to get good results.
2. FEM needs fever nodes to get
good results.
3. Fairly complicated problems can be
handled by FDM.
3. All type of complicated problem
can be handled by FEM.
4. FDM makes stair type of
approximation to sloping and
curved boundaries.
4. FEM can consider sloping and
curved boundaries exactly.
5. FDM makes point wise
approximation i.e. satisfy continuity
at node points only, along the side
continuity are not ensured.
5. FEM makes piece wise
approximation i.e. Satisfy
continuity at node point as well as
along the side/edge of element.
6. It is less efficient and more
approximate.
6. It is more efficient and more
approximate.
7. Not applicable for non-linearity of
domain.
7. Applicable to non-linearity of
domain.
8. Difficult to apply FDM for unusual
boundary condition and loading
conditions.
8. FEM can be applied to any type of
boundary condition and loading
conditions.

1.8 Applications of FEM
Finite Element Method was originally developed for the aerospace engineering but,
it is now widely used in other engineering disciplines also. The applications of
FEM are as follows:
a) Civil Engineering:
1. Analysis of bars, trusses, beams, frames and grids
2. Earthquake analysis of structures
3. Analysis of bridges, dams, retaining walls and water tanks
4. Bending, buckling and vibration analysis of beams, plates and shells.
b) Mechanical Engineering:

100. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
4

1. Analysis of pressure vessels, flywheels, gears, fluid flow, wave propagation
problems
2. IC engines, turbines, blades, steam pipes, nozzles
3. Analysis of casting, forming, welding and machining process.
c) Aerospace Engineering
1. Bending, buckling and vibration analysis of aircrafts, rockets, missiles,
spacecraft
2. Transient analysis of aircraft and spacecraft
d) Medical
1. Stress analysis of bones and teeth
2. Mechanics of heart valves
e) Electrical Engineering
1. Analysis of electromechanical devices such as motors, actuators
2. Eddy current and core losses in electric machines

1.9 Co-ordinate systems used in FEM
Three different coordinate systems are used in the finite element analysis
1) Local Co-ordinate system
2) Global Co-ordinate system
3) Natural Co-ordinate system
Local Co-ordinate System
When for each element in FEM, a separate coordinate system is used for deriving
element properties it is called as Local Co-ordinate system.

Global Co-ordinate System
The coordinate system is used to define the point in the entire structure is called as
Global Co-ordinate system.

101. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
5

Natural Co-ordinate System:
Natural Co-ordinate System is a coordinate system which permits the specification
of a point within the element by a set of dimensionless numbers, whose magnitude
never exceeds unity. In this coordinate system, origin is always assumed at the
center of element.

1.10 Principle of Virtual work
When the force and displacement are unrelated to the cause and effect relation, the
work is called virtual work. Therefore, the virtual work may be caused by true
force moving through imaginary displacements or vice versa. Principle of Virtual
work stated that “A body is in equilibrium if the internal virtual work equals to the
external virtual work for every kinematically admissible displacement fields. For
the linear elasticity following equation is the principle of virtual work. ie
WW


Internal workdone
External workdone
i ij ij
dv
e
ds
W dv
W q wds
  





ij ij
dv ds
dv q wds  


1.11 Principle of minimum potential energy
The Principle of minimum potential energy stated that “for a conservative system
for all kinematically admissible displacement fields those corresponding to
equilibrium, extremes the total potential energy, if extremum condition is
minimum, the equilibrium state is stable. UV

where 
= Total potential energy U
=Strain energy (due to internal forces) V
= Work potential (due to external forces) 
1
and
2
U V qw  

102. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
6

1.12 What is Node?
Nodes are selected points at which basic unknowns are to be determined in the
finite element analysis. Nodes are those points where elements are connected.
There are two types of nodes
a) External nodes
b) Internal nodes
External Nodes
Nodes which are occur at edges or surfaces of element are called as external nodes.
They are common to two or more elements.
Primary external nodes: Nodes which are occur at corners or ends of element are
called as primary external nodes
Secondary external nodes: Nodes which are occur along the edges of element
excluding corner nodes are called as secondary external nodes
Internal Nodes:
Internal nodes are the one which occur inside an element. They are specific to the
element selected and not common to other element.



1, 2 External Nodes,
3, 4 Internal Nodes.

1.13 Effective Node Numbering
Node Numbering
Node numbering should be such that semi band width/half band width is minimum.
Semi band width
A matrix is banded when all elements are zero except those within a band on either
side of the principal diagonal. The semi-bandwidth of such matrix is the number of
terms within the band to the right of (including) diagonal. The semi-bandwidth ‘B’
is given by expression (1 )B D f

Where,
D = maximum difference in node number in an element after considering all
elements.
f = degree of freedom per node.

1 3 4 2

103. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
7

Example: Let consider a rectangular domain (Rectangular slab) discretized into 24
elements:




Solution: Divide the given rectangular slab into 24 elements. Give the node
numbering in sequence following four different ways.
For Case I and II, move horizontally
For Case III and IV, move vertically

Determine the maximum difference between two consecutive node numbers to
determine semi/half band width.
Semi/half band width (1 )B D f

For two dimensional problem, f = 02 (DOF)
Case I) D = 13, f = 02 1 13 2 28B   
Case II) D = 8, f = 02 1 8 2 18B   
Case III) D = 9, f = 02 1 9 2 20B   
Case VI) D = 6, f = 02 1 6 2 14B    
Minimum band width = 14
Note: If, B = Semi band width and order of matrix is N x N,
Need to solve elements of N x B only.

104. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
8


Examples 1: Determine half bandwidth of overall stiffness matrix of the truss as
shown in figure. Is it possible to minimize bandwidth? if yes, then suggest
alternative node numbering scheme and hence half band width.

Solution: Case I)
Half band width (1 )B D f
For plane truss, DOF at each node are ( f ) = 02
To determine ‘D’ find out difference between two consecutive node numbers. The
maximum value of the difference is the value of ‘D’.

Maximum difference between two consecutive node numbers is 6 i.e. ‘D=6’.
Therefore, half band width is (1 )
(1 6)2
14
B D f



Case II) Yes, it is possible to minimize bandwidth.

105. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
9


Half band width
Case I) D = 5, f = 02 B = (1+5) 2 = 12
Case II) D = 5, f = 02 B = (1+5) 2 = 12
Case III) D = 6, f = 02 B = (1+6) 2 = 14
Case VI) D = 6, f = 02 B = (1+6) 2 = 14
Therefore, Minimum band width is equal to 12

Examples 2: Determine minimum band width of overall stiffness matrix of truss
shown in fig.

Solution:
The node numbers are given in four different ways and then maximum difference
between two consecutive node numbers.

106. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
10



1.14 Aspect Ratio of Element
It is a ratio of largest to smallest size of element. For accuracy of solution aspect
should be as close to unity as possible.
Example: Rectangle (4 x 6) m, divide into 6 element.

Size of element = 4 x 1 Size of element = 6 x 0.667
Aspect ratio = 4/1 = 4 Aspect ratio = 6/0.667 = 9

Size of element = 2 x 2 Size of element = 3 x 1.333
Aspect ratio = 2/2 = 1 Aspect ratio = 3/1.333 = 2.25
For the accuracy of solution, size of element should be (2 x 2)m.

107. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
11


1.15 Step by Step Procedure of FEM
Various steps involved in the finite element analysis are:
1) Select suitable field variables and the elements
The basic unknowns or the field variables which are encountered in the
engineering problem are,
Displacement in solid mechanics,
Velocities in fluid mechanics,
Electric and magnetic potentials in electrical engineering,
Temperature in heat transfer problem.
2) Discretize the Continuum
3) Select Interpolation Function: The finite element procedure reduces infinite
unknowns to a finite numbers by dividing domain into sub-domains called as
elements and by expressing the unknown field variables in terms of assumed
approximating/displacement function (Interpolating function /shape function)
within each element.
4) Find Element Properties: After selecting elements and nodal unknowns next
step in finite element analysis is to assemble element properties for each
element,
Example: 
eee
KQ
Where, 
e
K = element stiffness matrix, 
e

= nodal displacement vector of element, 
e
Q
= nodal force vector of element.
Element properties can be assembled by four methods,
1. Direct approach
2. Variational approach
3. Weighted residual approach
4. Energy balance approach.
6) Assemble Global Properties: Element properties are used to assemble global
properties to get system equations,
Example: KQ
Where, K = Global stiffness matrix of structure, 
= nodal load vector of structure, Q
= force vector of structure.
7) Impose the Boundary Condition: The boundary conditions are imposed to
find the solution of system equations which gives nodal unknowns.

108. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
12

8) Make the additional calculations to get the required values: Using nodal
unknowns additional calculations are made to get the required values, e.g.
stresses, strains, moments etc.

1.16 Types of Elements

1D Elements

Applications
1D Elements: Analysis of trusses, beams


2D Elements


Applications
2D Elements: Analysis of plane stress, plane strain and plate bending. Shell
structures are discretized using shell element or curve element.
3D Elements

109. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
13


Applications
3D Elements: Three-dimensional analysis, analysis of axisymmetric solids.

C
0
Continuity Element:
Elements in which only continuity of nodal variables are to be ensured. C
0
continuity is ensured in plane stress and plane strain problem. Example: Due to
Kirchhoff assumption that plane section remains plane even after bending, we have
the relation between slopes and displacements as /
y
wx  and /
x
wy  . If
Kirchhoff’s assumption is not made, slopes are independent of deflection and
hence ,,
xy
w are nodal unknowns reduced to C
0
continuity requirement.

C
1
continuity elements:
First order continuity elements in which higher order derivative of ‘w ’ is one only.
C
1
continuity elements satisfy not only displacement continuity but also slope
continuity should be satisfied. Hence the flexural problems (beams, plates, shells)
displacement and their first derivative are selected as nodal variables i.e. ,,
ww
w
xy


.
C
2
continuity elements:
Second order continuity elements in which second derivatives of ‘w ’ are also
nodal unknowns i.e. 2 2 2
22
, , , , ,
w w w w w
w
x y x y x y
    
      .


1.17 Difference of CST and LST Elements

110. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
14

Constant Strain Triangle (CST) elements:

Displacement field in terms of generalized coordinates 1 2 3
u x y    
(Continuous across element)
Resulting strain field is 2
(Constant)
x
u
x




( NOT continuous across element)

1. The CST element has three nodes and six degrees of freedom.
2. Inside each element, all components of strain are constant: hence the name
Constant Strain Triangle.
3. This is the simplest 2D element, which is also called linear triangular element.
4. Element stresses are also constant.
5. The displacement field is continuous across element boundaries
6. The strains and stresses are NOT continuous across element boundaries
7. Avoid CST in critical areas of structures (e.g., stress concentrations, edges of
holes, corners)
8. In general CSTs are not recommended for general analysis purposes as a very
large number of these elements are required for reasonable accuracy.
Linear Strain Triangle (LST) elements:

Displacement field in terms of generalized coordinates 22
1 2 3 4 5 6
u x y x xy y          

Resulting strain field is 2 4 5
2 (Linear)
x
u
xy
x
   

   



1. The LST element has six nodes and twelve degrees of freedom.
2. This element is also called quadratic triangular element.
3. The displacement function for the triangle is quadratic.
4. Strains and Stresses are linear
5. Quadratic elements are preferred for stress analysis, because of their high
accuracy and the flexibility in modeling complex geometry, such as curved
boundaries.

For a given number of nodes, a better representation of true stress and
displacement is generally obtained using LST elements than is obtained using

111. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
15

the same number of nodes a finer subdivision of CST elements. For example, a
single LST element gives better results than four CST elements.

1.18 Convergence Requirement of Displacement Function
1) Displacement function must be continuous and compatible within elements.
Continuous: Second element should start by taking one edge common of first
element.

Compatible: When it deforms, there should not be any discontinuity between
elements, i.e. Elements must not separate or overlap and 2 There should not be
any sudden change in slope across the inter element boundaries.
2) The displacement function must be capable of representing constant strain state
within the elements. This will achieved by dividing structure into smaller and
smaller elements.
3) The displacement function must be capable of representing rigid body
displacement, i.e. when nodes are given such displacement under rigid body
motion, elements should not experience any strain and hence leads to zero nodal
force. Constant term in polynomial ensures this condition.


1.19 2D and 3D Pascal’s Triangle
Besides the convergence and compatibility requirement, one of the important
considerations in choosing proper terms in the polynomial expansion is that the
element should have no preferred direction. That is displacement shapes will not
change with change in local coordinate system. This property is known as
Geometric Isotropy or Geometric Invariance. The Geometric Invariance is
achieved by Pascal’s triangle.




2D Pascal’s triangle

112. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
16


Geometric Invariance can be achieved by selecting the corresponding order of
terms on either side of triangle.

3D Pascal’s triangle

1.20 Displacement functions
Displacement function represents the approximate variation of the displacement
within the element. On the basis of the problem to be solved, the displacement
function needs to be approximated in the form of either linear or higher-order
function. A convenient way to express it is by the use of polynomial expressions
using Pascal’s triangle.

Displacement functions for various elements is as follows
1) Two noded bar element (1D)

DOF per node: 01 (u), Total DOF: 02
Select two elements from Pascal triangle to write
displacement function. Select only x or y-coordinate
Displacement Function: 12
ux

113. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
17

2) Two noded beam element (1D) (Bending element)

DOF per node: 02 (,w ), Total DOF: 04
Select four elements from Pascal triangle to write
displacement function. Select only x coordinate
Displacement Function: 23
1 2 3 4
w x x x      


3) Three noded Triangular (CST) element (2D)

DOF per node: 02 (u, v), Total DOF: 06
03 DOF in x-direction (u1 u2 u3); 03 DOF in y-direction (v1 v2
v3). Select three elements from Pascal triangle to write
displacement function. Since it is 2D element, select x and y
coordinates both
Displacement Function: 1 2 3
4 5 6
u x y
v x y
  
  
  
  


4) Six noded Triangular (LST) element (2D)

DOF per node: 02 (u, v), Total DOF: 12
06 DOF in x-direction (u1 u2 u3 u4 u5 u6)
06 DOF in y-direction (v1 v2 v3 v4 v5 v6)
Select six elements from Pascal triangle to write
displacement function. Since it is 2D element, select x
and y coordinates both
Displacement Function: 22
1 2 3 4 5 6
u x y x xy y          
22
7 8 9 10 11 12
v x y x xy y          


5) Four noded Rectangular element (2D)

DOF per node: 02 (u, v), Total DOF: 08
04 DOF in x-direction (u1 u2 u3 u4 )
04 DOF in y-direction (v1 v2 v3 v4)
Select four elements from Pascal triangle to write
displacement function. Since it is 2D element, select x
and y coordinates both
Displacement Function: 1 2 3 4
5 6 7 8
u x y xy
v x y xy
   
   
   
   

114. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
18





6) Four noded Rectangular plate bending element (2D)

Displacement Function: 22
1 2 3 4 5 6
3 2 2 3 3 3
7 8 9 10 11 12
w x y x xy y
x x y xy y x y xy
     
     
     
     

DOF per node: 03 ,,
xy
w , Total
DOF: 12
Select 12 elements from Pascal
triangle to write displacement
function. Since it is bending
element, displacement function will
be for deflection only. select x and
y coordinates both


7) Four noded tetrahedron element (3D)

DOF per node: 03 (u, v, w), Total DOF: 12
04 DOF in x-direction (u1 u2 u3); 04 DOF in
y-direction (v1 v2 v3); 04 DOF in z-direction
(w1 w2 w3). Select four elements from Pascal
triangle to write displacement function.
Since it is 3D element, select x, y and z
coordinates
Displacement Function: 1 2 3 4
5 6 7 8
9 10 11 12
u x y z
v x y z
w x y z
   
   
   
   
   
   


1.21 Natural coordinates
Natural coordinate system is basically a local coordinate system which allows the
specification of a point within the element by a set of dimensionless numbers
whose magnitude never exceeds unity. This coordinate system is found to be very
effective in formulating the element properties in finite element formulation. This
system is defined in such that the magnitude at nodal points will have unity or zero
or a convenient set of fractions. It also facilitates the integration to calculate
element stiffness.
I) Natural coordinates of 1D bar element (x-y coordinate system)
Let us consider a two noded bar/line element (1D). Node 1 and 2 have Cartesian
coordinates 1
x and 2
x respectively. Cartesian coordinate of any point ‘P’ is ‘x ’.
Natural coordinates of any point ‘P’ are (L1, L2).

115. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
19


From the definition of line element 12
1LL
1 1 2 2
L x L x x

In matrix form, 1
2 1 2
11 1L
L x x x
   
   
  
1
1 2 2
2 1 2 1 1 21
1 1 1 ( )11 11
1 ( )
L x x x
L x x x x xxxx x L

          
            
          
2
1
xx
L
L


and 1
2
xx
L
L


At node 11
xx , 1
1L , 2
0L
At node 2 2
xx , 1
0L , 2
1L
Therefore natural coordinates at node 1 is (1, 0) and natural coordinates at node
2 is (0, 1).


Natural coordinates of 1D bar element (, coordinate system)


=It is ratio of distance of any point P from origin (OP) to its maximum distance
from origin (O2)

116. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
20
1 2 2 1 1
1
/2
222
22
22
2
OP
L
x x x x x
xx
LL
Lx
x
L




         
     
      
      
  
  



At Node 1, 1
xx 1 1 1
1
2 2 2 2 2
22
L x x L x
x
LL

       
    
 
   
1  

At Node 2, 2
xx 1 2 1
2
2 2 2 2 2
22
22
2
L x x L x
x
LL
LL
L


       
     
 
   




1

Finally Natural Coordinates of 1D bar element in , coordinate system are:



II) Natural coordinates of rectangular element (, -coordinate system)

Four Noded Eight Noded



III) Natural Coordinates for 2D Triangular (CST) Element

117. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
21

(Area coordinates of CST Element)


Let us consider three noded triangular element as shown in figure. (x1, y1), (x2, y2),
(x3, y3) are the Cartesian coordinates of nodes 1, 2, 3 respectively. (u1, v1), (u2, v2),
(u3, v3) are the displacements of nodes 1, 2, 3 respectively.
Let ‘P’ be the any point on the element having Cartesian coordinates (x, y) and
natural coordinates (L1, L2, L3)
According to definition of natural coordinates
1 2 3
1 1 2 2 3 3
1 1 2 2 3 3
1L L L
L x L x L x x
L y L y L y y
  
  
   (1)
In matrix form, 1
11
2 1 2 3 2 1 2 3
3 1 2 3 3 1 2 3
1 1 1 1 1 1 1 1LL
L x x x x L x x x x
L y y y y L y y y y

         
          
         
   
       
   
         

1 1 1 1
2 2 2 2
3 3 3 3
1
1
L a b c
L a b c x
D
L a b c y
    
   
   

   

    (2)
where
   
2 3 3 2 3 1 1 3 1 2 2 1
1 2 3 3 2 1 3 1 1 3 1 1 2 2 1
2 3 2 2 3 1 2 1 2
3 3 1 3 1 3 3 2 1
D x y x y x y x y x y x y
a x y x y b x y x y c x y x y
a y y b y y c y y
a x x b x x c x x
     
     
     
      (3)
Area of triangle
Note: It is well known that the two times area of triangle is equals to its
determinant (2DA ). The following is the proof of this.

118. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
22



Area of  = Area of 1354+ Area of 3265 – Area of 1264
Area of trapezoidal =  
distance between two parallel sides
Sum of parallel sides
2
       
 
1 3 3 1 3 2 2 3 1 2 2 1
1 3 3 1 3 2 2 3 1 2 2 1
1 1 1
2 2 2
1
2
2
A y y x x y y x x y y x x
A y x y x y x y x y x y x
AD
        
     


Putting this in equation (2), we get
1 1 1 1
2 2 2 2
3 3 3 3
1
1
2
L a b c
L a b c x
A
L a b c y
    
   
   

   

   

3 3 31 1 1 2 2 2
1 2 3
; and
2 2 2
a b x c ya b x c y a b x c y
L L L
A A A
   
   (4)

Let divide the total area of triangle into three parts (A1, A2, A3).
Now, let us consider any point P shifted to node 1: 1
23
0
AA
AA



119. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
23
   
23
23
2 3 3 2 3 3 2 2 1
1 1 1 1
1 1 1
2
2
D x x x
y y y
D x y x y xy x y xy x y A
a b x c y A

      
   

Similarly, if any point P is shifted to node 2 and 3, we get 2 2 2 2 3 3 3 3
2 and 2a b x c y A a b x c y A     

Therefore, from equation (4) natural coordinates are written as 331 1 2 2
1 2 3
222
and
2 2 2
AAA A A A
L L L
A A A A A A
       

The variation of natural coordinates at each node is represented as


1.22 Introduction to 3D elements
A three dimensional elements can be considered in the problems where field
variables are dependent of x, y, & z. An example of a 3D Solid structure under
loading is as shown in figure.

Fig: 3D Solid under loading
 3-D elements can actually be used to model all kinds of structural components
including trusses, beams, plates, shells and so on.
 Typically 3-D solid elements can be tetrahedron or hexahedron in shape with
either flat or curves surfaces.

120. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
24

Applications
3D Elements: Three-dimensional analysis, analysis of axisymmetric solids.

3D solid elements
There are two basic families of three-dimensional elements similar to two-
dimensional case. Extension of triangular elements will produce tetrahedrons in
three dimensions. Similarly, rectangular parallelepipeds are generated on the
extension of rectangular elements. Following are few commonly used 3D solid
elements for finite element analysis.

Tetrahedron parallelepiped


4 Noded 10 Noded 8 Noded 20 Noded

3D Tetrahedron element:
The simplest element of the tetrahedral family is 4 noded tetrahedron.

DOF per node: 03 (u, v, w), Total DOF: 12
04 DOF in x-direction (u1 u2 u3); 04 DOF in y-direction (v1 v2 v3); 04 DOF in z-
direction (w1 w2 w3).
Displacement Function: 1 2 3 4
5 6 7 8
9 10 11 12
u x y z
v x y z
w x y z
   
   
   
   
   
   

121. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
25


3D Brick element (Hexahedron):

Figure shows, 8 noded
brick/hexahedron element in
natural coordinate system.
Displacement Function:
DOF per node: 03 (u, v, w),
Total DOF: 24
08 DOF in x-direction (u1 u2 u3);
08 DOF in y-direction (v1 v2 v3);
08 DOF in z-direction (w1 w2 w3).
Displacement Function:
(Natural coordinates, Pascal
Triangle) 1 2 3 4 5
6 7 8 5
u       
     
    
   

122. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
26

Unit-V
Shape functions/Interpolation Function

4.1 Shape functions
In FEM analysis, the displacement model we assume the variation of
displacements within the element since the true variation of displacements are not
known. But in higher engineering mathematics, analytical solution of some
problems is either not known or difficult to find out. In such cases we replaced that
function by another function which is easy to solve mathematically. That function
is called as “Shape function” or “Interpolation Function”.
Notes:
1) The magnitude of shape function at each node is Unity.
2) Number of shape functions are equal to number of nodes
3) The sum of shape functions is always unity

Methods for deriving shape functions:
1) Shape functions using polynomials in Cartesian coordinate system
2) Shape functions using polynomials in natural coordinates
3) Shape functions using Lagrange interpolation function in natural coordinates

4.2 Shape functions using polynomials in Cartesian coordinate system (x, y)

1) Shape functions for two noded bar element in (x, y) coordinates
Let consider two noded bar element of length L. x1 and x2 are the Cartesian
coordinates of nodes 1 and 2 respectively. u1 and u2 are the displacements of nodes
1 and 2. u is the displacement of any point in x-direction.

Total DOF = 02 (one at each node)

I) Displacement function 12
ux

In matrix form 
1
2
1ux



 


uP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements

123. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
27

Express displacement function in terms of nodal displacements using the
coordinates of nodes x1 and x2. 1 1 1
2 2 2
1
1
ux
ux


    
   
    


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

where [N] = Shape functions

1
N P A


III) Shape functions 
1
1 1
2
1
1
1
x
N P A x
x

 




Inverse is obtained by using method of adjoin 
1 211
1
112
xx
N P A x
 



12
21
1
2
N x x
N x x
   
   

   
21
12
and
x x x x
NN
LL



Sum of shape functions is always unity 2 1 2 1
12
1
x x x x x x L
NN
L L L L
  
     


At node 1, x=x1 2 1 1 1
12
1 and 0
x x L x x
NN
L L L

    

At node 2, x=x2 2 2 2 1
12
0 and 1
x x x x L
NN
L L L

    

124. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
28

2. Shape functions for two noded beam element (Bending element)

Let consider two noded beam element of length L, flexural rigidity EI. Degrees of
freedom at each of the beam element are two (translation and rotation). Total DOF
are 04.
w = Translation 
= Rotation
I) Displacement function
23
1 2 3 4
w x x x       (From Pascal Triangle)
2
2 3 4
23
dw
xx
dx
      
In matrix form 1
23
2
2
3
4
1
0 1 2 3
w x x x
xx






 
   
 



P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes x=0 at node 1 and x = L at node 2. 11
12
23
23
2
24
1 0 0 0
0 1 0 0
1
0 1 2 3
w
w L L L
LL




   
   
   
   

   

   
   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
w P A x



e
w N x

where [N] = Shape functions

125. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
29


III) Shape functions 
1
23
22
3 2 3 2
1 0 0 0
0 1 0 0
1
3 2 3 1
2 1 2 1
N P A x x x
/ L / L / L / L
/ L / L / L / L





  



Inverse of [A] is obtained using elementary operations of matrix algebra. 2 3 2 3
12 2 3 2
2 3 2 3
34 2 3 2
3 2 2
1
32
x x x x
N , N x ,
L L L L
x x x x
N , N
L L L L
     
    

At node 1; x = 0 N1 = 1, N2 = 0, N3 = 0, N4 = 0
At node 1; x = L N1 = 0, N2 = 0, N3 = 1, N4 = 0

3. Shape functions for three noded CST element

Let us consider three noded constant strain triangular (CST) element as shown in
figure. (x1, y1), (x2, y2), (x3, y3) are the Cartesian coordinates of nodes 1, 2, 3
respectively. (u1, v1), (u2, v2), (u3, v3) are the displacements of nodes 1, 2, 3
respectively. u, v are the displacements of any point on the element having
Cartesian coordinates (x, y).

I) Displacement function 1 2 3
u x y    

In matrix form  
1
2
3
1u x y





 



uP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements

126. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
30

Express displacement function in terms of nodal displacements using the
coordinates of nodes (x1, y1), (x2, y2), (x3, y3). 1 1 1 1
2 2 2 2
3 3 3 3
1
1
1
u x y
u x y
u x y



    
   
   

   

    


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

where [N] = Shape functions

III) Shape functions
Considering displacements in x-direction only (i.e. u, same shape functions are
applicable to displacement in y-direction i.e. v)  
1
11
1
22
3
1
11
1
xy
N P A x y x y
xy








 
1 1 1
2 2 2
3 3 3
1
1
2
a b c
N x y a b c
A
a b c






1 2 3 1 2 3 1 2 3
1 2 3
2 2 2
a a x a y b b x b y c c x c y
N N N
A A A
     
  

where 1 2 3 3 2 2 2 3 3 3 2
1 3 1 1 3 2 3 1 3 1 3
1 1 2 2 1 2 1 2 3 2 1
a x y x y a y y a x x
b x y x y b y y b x x
c x y x y c y y c x x
     
     
     


Example 1: A three noded triangular element is used in plane elasticity problem.
coordinates at nodes are 1(0,0), 2(4,0), 3(2,2). If u1, u2, u3 are nodal displacements
find shape functions.
Solutions: Shape functions are obtained by using 
1
N P A



where

127. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
31
 
11
22
33
1
11
1
xy
P x y A x y
xy






 
1
1 0 0
1 1 4 0
1 2 2
N x y








Inverse is obtained using method of adjoin  
8 0 0
1
1 2 2 0
8
2 2 4
N x y






1 2 3
8 2 2 2 2 4
8 8 8
x y x y y
N N N
  
  


Example 2: Coordinates of nodes of CST element are node 1(1, 2), 2(5, 3), 3(4, 6).
At interior point P if x = 3.3 and value of N1 = 0.3. Find coordinate of point P and
values of N2 and N3.
Solutions: Shape functions are obtained by using 
1
N P A


 
1
1
1 1 2
1 1 5 3
1 4 6
N P A x y








 
18 2 7
1
1 3 4 1
13
1 3 4
N x y


  



1 2 3
18 3 2 4 3 7 4
13 13 13
x y x y x y
N N N
      
  

Now, at x = 3.3 N1 = 0.3 18 3 3 3
0 3 4 2
13
.y
. y .
  
  

Using x = 3.3 and y = 4.2 23
2 4 3 3 3 4 2 7 3 3 4 4 2
0 2 0 5
13 13
. . . .
N . N .
       
   

128. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
32

4. Shape functions for four noded rectangular element
Let consider four noded rectangular element in Cartesian coordinate system. (x1,
y1), (x2, y2), (x3, y3), (x4, y4) are the Cartesian coordinates of nodes 1, 2, 3, 4
respectively. (u1, v1), (u2, v2), (u3, v3), (u4, v4) are the displacements of nodes 1, 2,
3, 4 respectively. u, v are the displacements of any point on the element having
Cartesian coordinates (x, y).

DOF: 02 at each node
Total DOF: 08 (04 in x-direction and 04 in y-direction)

I) Displacement function 1 2 3 4
u x y xy      

In matrix form  
1
2
3
4
1u x y xy







 




uP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes (0, 0), (a, 0), (a, b), (0, b). 11
22
33
44
1 0 0 0
1 0 0
1
1 0 0
u
u a
u a b ab
u b




   
   
   
   

   

   
   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get

129. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
33

1
e
u P A x



e
u N x

where [N] = Shape functions

III) Shape functions  
1
1
1 0 0 0
1 0 0
1
1
1 0 0
a
N P A x y xy
a b ab
b









Inverse is obtained using elementary operation of matrix algebra.  
1
1 0 0 0
1 1 0 0
1
1 0 0 1
1 1 1 1
/ a / a
N P A x y xy
/ b / b
/ ab / ab / ab / ab









12
32
1
x xy y x xy
N ; N ;
a ab b a ab
xy y xy
N ; N
ab b ab
     
  



Example: Derive shape functions for four noded rectangular element as shown in
figure.

Solution:
Shape functions are obtained by using 
1
N P A



130. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
34
 
1
1 0 0 0
1 6 1 6 0 0
1
1 4 0 0 1 4
1 24 1 24 1 24 1 24
//
N P A x y xy
//
/ / / /








12
32
1
6 24 4 6 24
24 4 24
x xy y x xy
N ; N ;
xy y xy
N ; N
     
  


4.2 Shape functions using polynomials in natural coordinate system (, )

1. Two noded bar element
Let consider a two noded bar element in natural coordinate system. Center of
element is assumed as origin (0 ). The node 1 has coordinate1 and node 2
has coordinate 1

I) Displacement function in-terms of natural coordinate
(use Pascal triangle in natural coordinates)
12
u  

1
2
1u




 

uP (1)
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes 1 and 2. 11
22
11
11
u
u


   
   
   

e
xA 

where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

131. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
35

where [N] = Shape functions

III) Shape functions 
1
1 11
1
11
N P A 

 



1
2
11
12
N
N


 
   

12
11
and
22
NN



2. Three noded bar element
Let consider a three noded bar element in natural coordinate system. Center of
element is assumed as origin. The node 1 has coordinate1 , node 2 has
coordinate0 , and node 3 has coordinate 1

I) Displacement function in-terms of natural coordinate
(Use Pascal triangle in natural coordinates)
2
1 2 3
u      

1
2
2
3
1u

  



 


uP (1)
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of node. 11
22
33
1 1 1
1 0 0
1 1 1
u
u
u



   
   
   

   

   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

132. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
36

where [N] = Shape functions
III) Shape functions 
1
1
2
1 1 1
1 1 0 0
1 1 1
N P A 





 



2
0 2 0
1
1 1 0 1
2
1 2 1
N 



 



 

2
1 2 3
11
1 and
22
N ; N N
   


   


2. Four noded rectangular element

I) Displacement function 1 2 3 4
u        
 
1
2
3
4
1u


  





 




uP (1)
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of node. 11
22
33
44
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
u
u
u
u




   
   
   
   

   

      

133. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
37


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

where [N] = Shape functions

III) Shape functions  
1
1 1 1 1
1 1 1 11
1
1 1 1 14
1 1 1 1
N P A   









 
 
12
32
1 1 1 1
44
1 1 1 1
44
NN
NN
   
   
   

   



4.3 Shape functions using Lagrange interpolation function in natural
coordinate system
If only continuity of basic unknown is to be satisfied, Lagrange polynomials can
be used to derive shape functions. Lagrange polynomials in one dimension is
defined by 1
or
n
mm
k
m ,m k m k m k
xx
N
xx








1. Two noded bar element
21
12
1 2 2 1
1 1 1 1
1 1 2 1 1 2
NN
       
   
     
     
    

134. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
38

2. Three noded bar element



23
1
1 2 1 3
13
2
2 1 2 3
21
3
3 2 3 1
101
1 0 1 1 2
11
11
0 1 0 1
101
1 0 1 1 2
N
N
N
     
   
     

   
     
   
   
    
     
   
      
   
   
    
   


3. Four noded rectangular element





24
1
1 2 1 4
13
2
2 1 2 3
42
3
3 4 3 2
31
4
4 3 4 1
1111
1 1 1 1 4
1111
1 1 1 1 4
1111
1 1 1 1 4
1111
1 1 1 1 4
N
N
N
N
     
   
     
   
     
   
     
   
   
    
     
   
    
    
   
    
   
   
    
    

Note: Generally shape functions for four noded rectangular element are given by   11
4
kk
k
N
 


135. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
39

4. Nine noded rectangular element


Shape functions for corner nodes (1, 2, 3, 4): 
5 2 4 8
1
1 5 1 2 1 4 1 8
1
0 1 0 1
1 0 1 1 1 0 1 1
11
4
N
N
           
       
  
       
       
           



5 1 6 3
2
2 5 2 1 2 6 2 3
11
4
N
         
       
   
    
   

7 4 6 2
3
3 7 3 4 3 6 3 2
11
4
N
         
       
   
    
   

7 3 8 1
4
4 7 4 3 4 8 4 1
11
4
N
         
       
   
    
   

In general shape functions for corner nodes are given by   11
4
k k k k
k
N
   


Shape functions for middle nodes (5, 6, 7, 8):  
2
1 2 9 7
5
5 1 5 2 5 9 5 7
11
2
N
         
       
   
    
   

     
2 2 2
6 7 8
1 1 1 1 1 1
2 2 2
N ; N ; N
             
  

In general shape functions for middle nodes are given by   
 
  
 
2
2
11
at nodes 0 axis
2
11
at nodes 0 axis
2
kk
k
kk
k
N
N
  

  


  

  

136. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
40

Shape function for central node:   
6 8 7 5
9
9 6 9 8 9 7 9 5
22
9
1 1 1 1
0 1 0 1 0 1 0 1
11
N
N
           
       

       
       
       
   

5. Eight noded hexahedron element

Shape functions using Lagrange interpolation function 
5 2 4
1
1 5 1 2 1 4
1 1 11 1 1
1 1 1 1 1 1 8
N
          
     
       
      
       




6 1 3
2
2 6 2 1 2 3
7 4 2
3
3 7 3 4 3 2
8 3 1
4
4 8 4 3 4 1
1 6 8
5
5 1 5 6 5 8
111
8
111
8
1 1 1
8
111
8
N
N
N
N
       
     
       
     
       
     
       
     
  
   
  
  
   
  
    
   
  
  
   
  

2 5 7
6
6 2 6 5 6 7
1 1 1
8
N
       
     
    
   
  

3 8 6
7
7 3 7 8 7 6
1 1 1
8
N
       
     
    
   
  

137. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
41

4 7 5
8
8 4 8 7 8 5
111
8
N
       
     
  
   
  

Serendipity elements
Higher order Lagrange elements contains internal nodes, which do not
contribute to the interelement connectivity. Internal nodes require extra
computational work. To avoid this, internal nodes can remove. Therefore, elements
which are having nodes only along the external boundaries are called as
serendipity elements. The elimination of these internal nodes results in reduction in
size of the element matrices.
Example





Shape functions for serendipity elements

1. Four noded rectangular element

Shape function for the node 1: 
1
11NC   

Since at any node value of shape function is unity, 1
1 at 1 1N&    14C/

1
11
4
N



Similarly   
2 3 4
1 1 1 1 1 1
4 4 4
N N N
          
  

4 8 12

138. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
42

2. Eight noded rectangular element

Shape functions for the corner nodes:  
1
111NC        

Since at any node value of shape function is unity, 1
1 at 1 1N&    14C/  
 
1
111
4
N
      
  

Similarly    
 
23
4
1 1 1 1 1 1
44
1 1 1
4
NN
N
       
   
       
   
   


Shape functions for the middle nodes: 
5
111NC      
5
1 0 1N at &  
12C/
 
2
5
11111
22
N
   
  

Similarly      
2 2 2
6 7 8
1 1 1 1 1 1
2 2 2
N N N
          
  

139. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
43

Unit-VI
Isoparametric formulation

6.1 Isoparametric elements:
For the analysis of structural problems of complex shapes involving curved
boundaries or surfaces, simple triangular or rectangular elements are no longer
sufficient. This has led to the development of elements of more arbitrary shapes
known as “Isoparametric elements”.
Concept of mapping in isoparametric element

Parent element in
Natural coordinate system
Mapped element in
Global coordinate system

140. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
44

There are three sub classes in isoparametric elements
1) Isoparametric elements
2) Sub-parametric elements
3) Super-parametric elements

Isoparametric elements: If the same number of elements is used to define
geometry as well as displacements, the element is called as isoparametric elements.
e.g.
1 1 2 2 8 8
1 1 2 2 8 8
x N x N x .............. N x
y N y N y .............. N y
   
    (Geometry) 1 1 2 2 8 8
1 1 2 2 8 8
u N u N u .............. N u
v N v N v .............. N v
   
   
(Displacements)

Sub-parametric elements: The elements in which less number of nodes is used to
define geometry compared to the number of nodes used to define displacements,
the element is called as sub-parametric elements.
e.g.
1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
x N x N x N x N x
y N y N y N y N y
   
    (Geometry) 1 1 2 2 8 8
1 1 2 2 8 8
u N u N u .............. N u
v N v N v .............. N v
   
   
(Displacements)

Super-parametric elements: The elements in which less number of nodes is used
to define displacements compared to the number of nodes used to define geometry,
the element is called as sub-parametric elements. This element is used in problems
of stress analysis where boundaries are highly curved but stress gradient is not
high.
e.g.
1 1 2 2 8 8
1 1 2 2 8 8
x N x N x ................. N x
y N y N y ................. N y
   
    (Geometry) 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
u N u N u N u N u
v N v N v N v N v
   
   
(Displacements)

141. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
45



6.2 Basic theorems of isoparametric formulation
Isoparametric concept is developed based on the following three basic theorems.

Theorem I: If two adjusted elements are generated using shape functions, then
there is continuity at the common edge.
Theorem II: It states, if the shape functions used are such that continuity of
displacement is represented in the parent coordinates then the continuity
requirement will be satisfied in the isoparametric elements also.
Theorem III: The constant strain condition and constant derivative condition
satisfied by all isoparametric elements.

6.3 Advantages of isoparametric elements
1. Shape functions are used for defining the geometry as well as displacements
2. Suitable for structures having complex shapes or curved boundaries
3. When computation of [A]
-1
to obtain stiffness matrix of each elements
requires more computational time.
4. Increase accuracy of results
5. Satisfy pascal’s triangle and required less computational work
6.4 Coordinate transformation
Shape functions are used for transformation of natural coordinate system into
global coordinate system. Thus the (Cartesian coordinate) coordinates of a point in
global coordinate system may be expressed as 1 1 2 2
1 1 2 2
nn
nn
x N x N x ................. N x
y N y N y ................. N y
   
   

where [N] = Shape function x
= Coordinates of any point 
e
x
= Coordinates of nodal points

142. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
46

Example 1: Determine the Cartesian coordinate (x, y) of the any point P (0 5 0 6. , .
) as shown in figure.

Solution: Parent element for the quadrilateral in natural coordinate system is four
noded rectangle.

Given, P (, ) = P (0 5 0 6. , . )
Shape functions for the four noded rectangular element in natural coordinate
system are  
 
12
34
1 1 1 1
44
1 1 1 1
44
N ; N
N ; N
   
   
   

   


Values of these shape functions at given coordinates (0 5 0 6. , . )      
     
12
34
1 0 5 1 0 6 1 0 5 1 0 6
0 05 0 15
44
1 0 5 1 0 6 1 0 5 1 0 6
0 6 0 2
44
. . . .
N . N .
. . . .
N . N .
   
   
   
   

Given Nodal coordinates 1 1 2 2
3 3 4 4
2 1 8 3
7 7 3 5
x , y , x , y ,
x , y , x , y
   
   

Therefore, coordinates of point ‘P’ in Cartesian system are P (x, y) where 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
61
57
x N x N x N x N x .
y N y N y N y N y .
    
    

143. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
47

Example 2: For the isoparametric quadrilateral element shown in figure, determine
local coordinates of the point ‘Q’ which has Cartesian coordinates (7, 4).

Solution: Parent element for the quadrilateral in natural coordinate system is four
noded rectangle.

Shape functions for the four noded rectangular element in natural coordinate
system are  
 
12
34
1 1 1 1
44
1 1 1 1
44
N ; N
N ; N
   
   
   

   


Cartesian coordinates of point ‘P’ are 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
x N x N x N x N x
y N y N y N y N y
   
   

Given P (x,y ) = P (7, 4) 1 1 1 1 1 1 1 1
7 3 6 8 2
4 4 4 4
              
   

9 9 3      (1) 1 1 1 1 1 1 1 1
4 1 1 6 5
4 4 4 4
              
   

39      (2)
From equation (1), 39
1





 (3)
Put into the equation (2)

144. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
48
3 9 3 9
9 9 3
11



   
  
   
   
2
26 80 18 0  
0 2105.

Put this value into the equation (3) 0 91325.


6.5 Jacobian matrix for four noded quadrilateral element
Let consider a four noded quadrilateral element in Cartesian coordinate system.
The coordinates of nodes are (x1, y1), (x2, y2), (x3, y3), (x4, y4).

Parent element for the quadrilateral in natural coordinate system is four noded
rectangle.

Geometry in-terms of shape functions is represented as 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
x N x N x N x N x
y N y N y N y N y
   
   

Displacements in-terms of shape functions is represented as 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
u N u N u N u N u
v N v N v N v N v
   
   

Shape functions for the four noded rectangular element in natural coordinate
system are

145. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
49
 
 
12
34
1 1 1 1
44
1 1 1 1
44
N ; N
N ; N
   
   
   

   


Strain displacement relationship 
1
2
1 2 3 4
3
41 2 3 4
1
21 2 3 4 1 2 3 4
3
4
0 0 0 0
0 0 0 0
u
u
u N N N N
u
x x x x x
uv N N N N
vy y y y y
vN N N N N N N Nuv
y y y y x x x xyx v
v



      
 

     

      
   
    
   
          
   
         



e
Bx

where [B] = Strain displacement matrix
i i i
i i i
N N N
x x x
N N N
y y y




    

    
    

     (i = 1, 2, 3, 4)
In matrix form ii
i i
NN
x x x
N N
yyy




  
 
     
   
      
     

where [J] = Jacobian Matrix
1 xx
J
yy










and xx
J
yy










Elements of Jacobian matrix are

146. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
50
4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
i
i
i
i
i
i
i
i
i
i
i
i
x N N N N N
x x x x x
x N N N N N
x x x x x
y N N N N N
y y y y y
y N N N N N
y y y y y
     
     
     
     




     
    
     
     
    
     
     
    
     
     
    
     





Finally Jacobian matrix is obtained as 44
11
44
11
ii
ii
ii
ii
ii
ii
NNxx
xx
J
yy NN
yy
   
 


   

   

 

  




6.6 Strain displacement matrix for four noded quadrilateral element
Strain displacement matrix is obtained from strain displacement relationship
explained in Jacobian matrix. 
1 2 3 4
1 2 3 4
1 2 3 4 1 2 3 4
0 0 0 0
0 0 0 0
N N N N
x x x x
N N N N
B
y y y y
N N N N N N N N
y y y y x x x x
   

   

    


   

       

       

Elements of Strain displacement matrix are as follows 1 1 1 1 1 1
N N N N N N
x x x y y y
   
   
         
   
         
2 2 2 2 2 2
N N N N N N
x x x y y y
   
   
         
   
         
3 3 3 3 3 3
N N N N N N
x x x y y y
   
   
         
   
         
4 4 4 4 4 4
N N N N N N
x x x y y y
   
   
         
   
         

147. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
51

Example: Obtain Jacobian matrix for the quadrilateral element as shown in figure
using isoparametric formulation.

Solution: Parent element for the quadrilateral in natural coordinate system is four
noded rectangle.

Shape functions for the four noded rectangular element in natural coordinate
system are  
 
12
34
1 1 1 1
44
1 1 1 1
44
NN
NN
   
   
   

   


Elements of Jacobian matrix are as follows 







 
1 2 3 4
1 2 3 4
1 1 1 1 6 2
2 4 5 1
4 4 4 4 4
x N N N N
x x x x
    
    
    
   
    
    
     








1 2 3 4
1 2 3 4
1 1 1 1
2 4 5 1
4 4 4 4 2
x N N N N
x x x x
    
    
    
   
    
    
    








1 2 3 4
1 2 3 4
1 1 1 1 1
4 5 2 1
4 4 4 4 2
y N N N N
y y y y
    
   
    
   
    
   
     

148. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
52








1 2 3 4
1 2 3 4
1 1 1 1 3
4 5 2 1
4 4 4 4 2
y N N N N
y y y y
    
   
    
   
    
   
    


Therefore, Jacobian matrix  62
42
13
22
xx
J
yy



  
 

 
 
 
 


Example: Obtain strain displacement matrix for the quadrilateral element as
shown in figure using isoparametric formulation.


Solution: Elements of Jacobian matrix from previous example  62
42
13
22
xx
J
yy



  
 

 
 
 
 

Using Jacobian matrix  
4 2 2
2
6 2 3x x y y
   

   
     
    

Elements of Strain-displacement matrix are as follows
Shape functions for the four noded rectangular element in natural coordinate
system are  
 
12
34
1 1 1 1
44
1 1 1 1
44
NN
NN
   
   
   

   


149. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
53

Elements of strain-displacement matrix are 
1 2 3 4
1 2 3 4
1 2 3 4 1 2 3 4
0 0 0 0
0 0 0 0
N N N N
x x x x
N N N N
B
y y y y
N N N N N N N N
y y y y x x x x
   

   

    


   

       

       


 


 



1 1 1
1 1 1
11 42
4 6 2 4
11
6 2 2
11 2
2
4 4 3
11
26
N N N
x x x
N N N
y y y

   





        
   

      
  


       
   

     
  


 


 



2 2 2
2 2 2
11 42
4 6 2 4
11
6 2 2
11 2
2
4 4 3
11
26
N N N
x x x
N N N
y y y

   





      
   

      



     
   

     



 


 



3 3 3
3 3 3
11 42
4 6 2 4
11
6 2 2
11 2
2
4 4 3
11
26
N N N
x x x
N N N
y y y

   





      
   

      



     
   

     



150. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
54

 


 



4 4 4
4 4 4
11 42
4 6 2 4
11
6 2 2
11 2
2
4 4 3
11
26
N N N
x x x
N N N
y y y

   





        
   

      
  


       
   

     
  


151. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
55

Unit-VII
Stiffness Matrices of various elements

6.1 Stiffness Matrix
Element Stiffness Matrix
The stiffness matrix is an inherent property of the structure. Element stiffness is
obtained with respect to its axes and then transformed this stiffness to structure
axes. The characteristics of stiffness matrix are as follows:
1. Stiffness matrix is symmetric and square.
2. In stiffness matrix, all diagonal elements are positive.
3. Stiffness matrix is positive definite

Global Stiffness Matrix
A structural system is an assemblage of number of elements. These elements are
interconnected together to form the whole structure. Therefore, the element
stiffness of all the elements are first need to be calculated and then assembled
together in systematic manner. This matrix is called as global stiffness matrix.

6.2 Stiffness Matrices for various elements

1. Stiffness matrix for two noded bar element
Let consider two noded bar element of length L. x1 and x2 are the Cartesian
coordinates of nodes 1 and 2 respectively. u1 and u2 are the displacements of nodes
1 and 2. u is the displacement of any point in x-direction.

Total DOF = 02 (one at each node)

I) Displacement function 12
ux

In matrix form 
1
2
1ux



 


uP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes x1 and x2.

152. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
56
1 1 1
2 2 2
1
1
ux
ux


    
   
    


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x




e
u N x (3)
where [N] = Shape functions
III) Shape functions 
1
1 1
2
1
1
1
x
N P A x
x

 




Inverse is obtained by using method of adjoin 
1 211
1
112
xx
N P A x
 



12
21
1
2
N x x
N x x
   
   

   
21
12
and
x x x x
NN
LL



IV) Strain-Displacement relationship
Strain-displacement relations from linear theory of elasticity are used. Using Eq.
(4) one can write 1 1 2 2
u N u N u

12
12
112
2
x
x
du dN dN
uu
dx dx dx
udN dN
udx dx


  

 
 


x e
Bx (4)
where [B] = Strain-displacement matrix
Elements of strain displacement matrix are 
11
B
LL




12
11dN dN
,
dx L dx L

  



V) Stress-Strain relationship

153. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
57

Using one dimensional Hooke’s law from linear theory of elasticity D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix
VI) Principle of virtual work
The principle of virtual work states that the internal work done is equal to external
work done 
 

0
TT
e
L
TTT
e e e
e
Qx
Q x D B x B x
Q K x





where 
0
L
T
K D B B
VII) Elements of Stiffness matrix 0
L
ji
ij
dNdN
K AE dx
dx dx

11
11
00
12
12 21
00
22
22
00
11
11
11
LL
LL
LL
dN dN AE
K AE dx AE dx
dx dx L L L
dN dN AE
K K AE dx AE dx
dx dx L L L
dN dN AE
K AE dx AE dx
dx dx L L L
  
    
  
  
  
     
  
  
  
  
  
  




Therefore, stiffness matrix of two noded bar element is 11
11
AE
K
L






2. Stiffness matrix for two noded beam element

154. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
58

Let consider two noded beam element of length L, flexural rigidity EI. Degrees of
freedom at each of the beam element are two (translation and rotation). Total DOF
are 04.
w = Translation 
= Rotation
I) Displacement function
23
1 2 3 4
w x x x       (From Pascal Triangle)
2
2 3 4
23
dw
xx
dx
      
In matrix form 1
23
2
2
3
4
1
0 1 2 3
w x x x
xx






 
   
 



P (1)
where, [P] = Parametric matrix

II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes x=0 at node 1 and x = L at node 2. 11
12
23
23
2
24
1 0 0 0
0 1 0 0
1
0 1 2 3
w
w L L L
LL




   
   
   
   

   

   
   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
w P A x




e
w N x (3)
where [N] = Shape functions

III) Shape functions

155. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
59

1
23
22
3 2 3 2
1 0 0 0
0 1 0 0
1
3 2 3 1
2 1 2 1
N P A x x x
/ L / L / L / L
/ L / L / L / L





  



Inverse of [A] is obtained using elementary operations of matrix algebra. 2 3 2 3
12 2 3 2
2 3 2 3
34 2 3 2
3 2 2
1
32
x x x x
N , N x ,
L L L L
x x x x
N , N
L L L L
     
    

IV) Strain-Displacement relationship

2 2 2 2 2
1 2 3 4
2 2 2 2 2 e
d w d N d N d N d N
x
dx dx dx dx dx


   




e
Bx (4)
where [B] = Strain-displacement matrix
Elements of strain displacement matrix are 
2 3 2 2 3 2
6 12 4 6 6 12 2 6x x x x
B
L L L L L L L L
       
     
       
       

V) Stress-Strain relationship D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix
VI) Principle of virtual work 
 

0
TT
e
L
TTT
e e e
e
Qx
Q x D B x B x
Q K x





where 
0
L
T
K D B B
VII) Elements of Stiffness matrix

156. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
60
22
22
0
L
ji
ij
dNdN
K EI dx
dx dx


Let i = 1 and j = 1 22
11
11 2 2 2 3 2 3
00
11 3
6 12 6 12
12
LL
d N d N x x
K EI dx EI dx
dx dx L L L L
EI
K
L
  
     
  
  



Similarly other elements can be determined. 
3 2 3 2
1
22
1
3 2 3 2
2
22
2
12 6 12 6
6 4 6 2
12 6 12 6
6 2 6 4
w/ L / L / L / L
/ L / L / L / L
K EI
w/ L / L / L / L
/ L / L / L / L


 




   




3. Stiffness matrix for three noded CST element

Let us consider three noded constant strain triangular (CST) element as shown in
figure. (x1, y1), (x2, y2), (x3, y3) are the Cartesian coordinates of nodes 1, 2, 3
respectively. (u1, v1), (u2, v2), (u3, v3) are the displacements of nodes 1, 2, 3
respectively. u, v are the displacements of any point on the element having
Cartesian coordinates (x, y).
I) Displacement function 1 2 3
4 5 6
u x y
v x y
  
  
  
  

In matrix form  
1
2
3
1u x y





 


and  
4
5
6
1v x y





 


uP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements

157. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
61

Express displacement function in terms of nodal displacements using the
coordinates of nodes (x1, y1), (x2, y2), (x3, y3).
Note: Shape functions used to represents the displacements in x and y directions
are same. Therefore consider nodal displacements of x-direction to determine
shape functions 1 1 1 1
2 2 2 2
3 3 3 3
1
1
1
u x y
u x y
u x y



    
   
   

   

    


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x




e
u N x (3)
where [N] = Shape functions

III) Shape functions
Considering displacements in x-direction only (i.e. u, same shape functions are
applicable to displacement in y-direction i.e. v)  
1
11
1
22
3
1
11
1
xy
N P A x y x y
xy








 
1 1 1
2 2 2
3 3 3
1
1
2
a b c
N x y a b c
A
a b c






1 2 3 1 2 3 1 2 3
1 2 3
2 2 2
a a x a y b b x b y c c x c y
N N N
A A A
     
  

where 1 2 3 3 2 2 2 3 3 3 2
1 3 1 1 3 2 3 1 3 1 3
1 1 2 2 1 2 1 2 3 2 1
a x y x y a y y a x x
b x y x y b y y b x x
c x y x y c y y c x x
     
     
     

Substituting shape functions in the equation (3), one can write 1 1 2 2 3 3
1 1 2 2 3 3
u N u N u N u
v N v N v N v
  
  

IV) Strain-Displacement relationship

158. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
62

From the linear theory of elasticity, strain-displacement relations for 2D elasticity
problem are 1
1 2 3
2
31 2 3
1
21 2 3 1 2 3
3
0 0 0
0 0 0
x
y
xy
u
u N N N
u
x x x x
uv N N N
vy y y y
vu v N N N N N N
y x y y y x x x v



      
   
       
      
      
   
      
           
   
           


e
Bx (4)
Elements of strain-displacement matrix are 
2 2 2
3 3 3
3 3 3 2 2 2
000
222
000
222
222222
a b c
AAA
a b c
B
AAA
a b c a b c
AAAAAA










V) Stress-Strain relationship
Using stress-strain relationship of for 2D elasticity problem 2
10
10
1
0 0 1 2
xx
yy
xy xy
E
/
  
  

  
    
    
   


   
 
   
D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix
VI) Principle of virtual work 
 

TT
e
TTT
e e e
dA
e
Qx
Q x D B x B x dA
Q K x





where 
T
dA
K D B B dA

159. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
63

Example : Derive strain displacement matrix [B], elasticity matrix [D] for the
three noded triangular element as shown in figure.

Solution: I) Displacement function 1 2 3
4 5 6
u x y
v x y
  
  
  
  

In matrix form  
1
2
3
1u x y





 


and  
4
5
6
1u x y





 


u,v P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes (1, 1), (4, 3), (2, 5).
Consider only x-directional displacement 11
22
33
1 1 1
1 4 3
1 2 5
u
u
u



   
   
   

   

   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



e
u N x

where [N] = Shape functions

III) Shape functions
Considering displacements in x-direction only (i.e. u, same shape functions are
applicable to displacement in y-direction i.e. v)

160. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
64
 
1
1
1 1 1
1 1 4 3
1 2 5
N P A x y








 
14 3 1
1
1 2 4 2
10
2 1 3
N x y


  



1 2 3
14 2 2 3 4 1 2 3
10 10 10
x y x y x y
N N N
       
  

one can write 1 1 2 2 3 3
1 1 2 2 3 3
u N u N u N u
v N v N v N v
  
  


IV) Strain-Displacement relationship
From the linear theory of elasticity, strain-displacement relations for 2D elasticity
problem are 1
1 2 3
2
31 2 3
1
21 2 3 1 2 3
3
0 0 0
0 0 0
x
y
xy
u
u N N N
u
x x x x
uv N N N
vy y y y
vu v N N N N N N
y x y y y x x x v



      
   
       
      
      
   
      
           
   
           

e
Bx

Elements of strain-displacement matrix are 
1 2 1
000
555
1 1 3
000
5 10 10
1 1 3 1 2 1
5 10 10 5 5 5
B








   



161. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
65

4. Stiffness matrix for four noded rectangular element
Let consider four noded rectangular element in Cartesian coordinate system. (x1,
y1), (x2, y2), (x3, y3), (x4, y4) are the Cartesian coordinates of nodes 1, 2, 3, 4
respectively. (u1, v1), (u2, v2), (u3, v3), (u4, v4) are the displacements of nodes 1, 2,
3, 4 respectively. u, v are the displacements of any point on the element having
Cartesian coordinates (x, y).

DOF: 02 at each node
Total DOF: 08 (04 in x-direction and 04 in y-direction)

I) Displacement function 1 2 3 4
5 6 7 8
u x y xy
v x y xy
   
   
   
   

In matrix form  
1
2
3
4
1u x y xy







 



and  
5
6
7
8
1v x y xy







 



uP and vP (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes (0, 0), (a, 0), (a, b), (0, b).
Note: Shape functions used to represents the displacements in x and y directions
are same. Therefore consider nodal displacements of x-direction to determine
shape functions

162. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
66
11
22
33
44
1 0 0 0
1 0 0
1
1 0 0
u
u a
u a b ab
u b




   
   
   
   

   

   
   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x




e
u N x (3)
where [N] = Shape functions
III) Shape functions  
1
1
1 0 0 0
1 0 0
1
1
1 0 0
a
N P A x y xy
a b ab
b









Inverse is obtained using elementary operation of matrix algebra.  
1
1 0 0 0
1 1 0 0
1
1 0 0 1
1 1 1 1
/ a / a
N P A x y xy
/ b / b
/ ab / ab / ab / ab









12
32
1
x xy y x xy
N ; N ;
a ab b a ab
xy y xy
N ; N
ab b ab
     
  

Substituting shape functions in the equation (3), one can write 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
u N u N u N u N u
v N v N v N v N v
   
   

IV) Strain-Displacement relationship

163. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
67
1
2
1 2 3 4
3
41 2 3 4
1
21 2 3 4 1 2 3 4
3
4
0 0 0 0
0 0 0 0
x
y
xy
u
u
u N N N N
u
x x x x x
uv N N N N
vy y y y y
vu v N N N N N N N N
y x y y y y x x x x v
v





        
   

    
    
          
      
    
      
             
   
             






e
Bx (4)
Elements of strain-displacement matrix are 
11
0 0 0 0
11
0 0 0 0
1 1 1 1
y y y y
a ab a ab ab ab
x x x x
B
ab b ab ab ab b
x x x x y y y y
ab b ab ab ab b a ab a ab ab ab

   



    


       


IV) Stress-Strain relationship
Using stress-strain relationship of for 2D elasticity problem 2
10
10
1
0 0 1 2
xx
yy
xy xy
E
/
  
  

  
    
    
   


   
 
   
D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix

VI) Principle of virtual work 
 

TT
e
TTT
e e e
dA
e
Qx
Q x D B x B x dA
Q K x





164. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
68

where 
T
dA
K D B B dA

Example : Derive strain displacement matrix [B], elasticity matrix [D] for the four
noded rectangular element as shown in figure.

Solution: I) Displacement function 1 2 3 4
5 6 7 8
u x y xy
v x y xy
   
   
   
   

In matrix form  
1
2
3
4
1u x y xy







 



and  
5
6
7
8
1v x y xy







 



u,v P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of nodes (0, 0), (6, 0), (6, 4), (0, 4). 11
22
33
44
1 0 0 0
1 6 0 0
1 6 4 24
1 0 4 0
u
u
u
u




   
   
   
   

   

   
   


e
xA  (2)
where, [A] = Connectivity matrix
Obtained  from Eq. (2) and put into the Eq. (1), we get 
1
e
u P A x



165. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
69

e
u N x

where [N] = Shape functions

III) Shape functions  
1
1
1 0 0 0
1 6 0 0
1
1 6 4 24
1 0 4 0
N P A x y xy









Inverse is obtained using elementary operation of matrix algebra.
 
1
1 0 0 0
1 6 1 6 0 0
1
1 4 0 0 1 4
1 24 1 24 1 24 1 24
//
N P A x y xy
//
/ / / /








12
32
1
6 24 4 6 24
24 4 24
x xy y x xy
N ; N ;
xy y xy
N ; N
     
  

One can write 1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
u N u N u N u N u
v N v N v N v N v
   
   


IV) Strain-Displacement relationship 1
2
1 2 3 4
3
41 2 3 4
1
21 2 3 4 1 2 3 4
3
4
0 0 0 0
0 0 0 0
x
y
xy
u
u
u N N N N
u
x x x x x
uv N N N N
vy y y y y
vu v N N N N N N N N
y x y y y y x x x x v
v





        
   

    
    
          
      
    
      
             
   
             





e
Bx

Elements of strain-displacement matrix are

166. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
70

11
0 0 0 0
6 24 6 24 24 24
11
0 0 0 0
24 4 24 24 24 4
1 1 1 1
24 4 24 24 24 4 6 24 6 24 24 24
y y y y
x x x x
B
x x x x y y y y
   
   
   
   

    
    
    
   

       
              
       









Formulation of stiffness matrix for 1D isoparametric element

Let consider a two noded bar element in natural coordinate system. Center of
element is assumed as origin (0 ). The node 1 has coordinate1 and node 2
has coordinate 1 .
x1, x2 are the Cartesian coordinates of nodes 1 and 2.

I) Shape functions in natural coordinates:
12
11
and
22
NN

 (1)
II) Geometry:
1 1 2 2
x N x N x (2)
III) Displacements:
1 1 2 2
u N u N u (3)
IV) Strain-displacement relationship 12
12x
du dN dN
uu
dx dx dx
  

In matrix form 
112
2
udN dN
udx dx


 
 


e
Bx (4)

167. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
71

where [B] = Strain-displacement matrix
Elements of [B] matrix
1 1 2 2
dN dN d dN dN d
dx d dx dx d dx


 (5)
From Eq. (1) 12
11
22
dN dN
dd



From Eq. (2) 12
12
dx dN dN
xx
d d d  

1 2 2 1
2 2 2 2
dx x x x x L
d

    
2d
dx L



From Eq. (5), elements of strain-displacement matrix are 12
1 2 1 1 2 1
22
dN dN
dx L L dx L L
       

11
B
LL





V) Stress-Strain relationship D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix
VI) Principle of virtual work 
 

0
TT
e
L
TTT
e e e
e
Qx
Q x D B x B x
Q K x





where 
0
L
T
K D B B dx

Formulation of stiffness matrix for 2D four noded quadrilateral isoparametric
element

168. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
72

Let consider four noded quadrilateral element in Cartesian coordinate system.
(x1, y1), (x2, y2), (x3, y3), (x4, y4) are the Cartesian coordinates of nodes 1, 2, 3, 4
respectively. (u1, v1), (u2, v2), (u3, v3), (u4, v4) are the displacements of nodes 1, 2,
3, 4 respectively. u, v are the displacements of any point on the element having
Cartesian coordinates (x, y).

I) Shape functions in natural coordinates:
 
 
12
34
1 1 1 1
44
1 1 1 1
44
NN
NN
   
   
   

   
 (1)
II) Geometry:
1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
x N x N x N x N x
y N y N y N y N y
   
    (2)
III) Displacements:
1 1 2 2 3 3 4 4
1 1 2 2 3 3 4 4
u N u N u N u N u
v N v N v N v N v
   
     (3)
IV) Strain-displacement relationship 
1
2
1 2 3 4
3
41 2 3 4
1
21 2 3 4 1 2 3 4
3
4
0 0 0 0
0 0 0 0
u
u
u N N N N
u
x x x x x
uv N N N N
vy y y y y
vN N N N N N N Nuv
y y y y x x x xyx v
v



      
 

     

      
   
    
   
          
   
         



In matrix form

e
Bx (4)
where [B] = Strain-displacement matrix

169. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
73

i i i
i i i
N N N
x x x
N N N
y y y




    

    
    

     (i = 1, 2, 3, 4)
In matrix form ii
i i
NN
x x x
N N
yyy




  
 
     
   
      
     

where [J] = Jacobian Matrix
1 xx
J
yy










and xx
J
yy









Elements of Jacobian matrix are 4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
4
1 2 3 4
1 2 3 4
1
i
i
i
i
i
i
i
i
i
i
i
i
x N N N N N
x x x x x
x N N N N N
x x x x x
y N N N N N
y y y y y
y N N N N N
y y y y y
     
     
     
     




     
    
     
     
    
     
     
    
     
     
    
     





Finally Jacobian matrix is obtained as 44
11
44
11
ii
ii
ii
ii
ii
ii
NNxx
xx
J
yy NN
yy
   
 


   

   

 

  



Elements of Strain displacement matrix [B] are as follows 1 1 1 1 1 1
N N N N N N
x x x y y y
   
   
         
   
         

170. Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
74
2 2 2 2 2 2
N N N N N N
x x x y y y
   
   
         
   
         
3 3 3 3 3 3
N N N N N N
x x x y y y
   
   
         
   
         
4 4 4 4 4 4
N N N N N N
x x x y y y
   
   
         
   
         

V) Stress-Strain relationship D

From Eq. (4) 
e
D B x

where [D] = Elasticity matrix
VI) Principle of virtual work 
 

TT
e
TTT
e e e
dA
e
Qx
Q x D B x B x dA
Q K x





where 
T
dA
K D B B dA