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: attu_sayyad@[Link]
Ph. No.: (+91) 9763567881

Year-2017
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
sin   ; cos   1; 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:
 xx  plane ' yz' and direction ' x'
 yy  plane ' xz' and direction ' y'
 zz  plane ' xy' and direction ' z'
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 .
 xy  plane ' yz' and direction ' 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.
 xx  xy  xz   xx  xy  xz 
 33   yx  yy  yz    33   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.
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 xx  xy  xz   xx  xy  xz 
 33   yx   
 yy  yz    33   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.
    xx ,  yy , zz , xy , yx , xz , zx , yz , zy 
T

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  xx dy dz
x
  xy
 'xy   xy  xy dx
x
  xz
 'xz   xz  xz dx
x
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

xz  yy  yy dx dz
 'yy   yy  dy
y
 yx  yx
 'yx   yx  dy
y
 yz  yz
 'yz   yz  dy
y
xy 
 'zz   zz  zz dz  zz dx dy
z
  zx
 'zx   zx  zx dz
z
  zy
 'zy   zy  zy dz
z

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,

 Fx  0
  xx    yx 

 xx  dx  dy dz   dy dz  
 yx  dy  dx dz
x y
xx
    (1)
  
 yx dx dz   zx  zx dz  dx dy   zx dx dy  X dx dy dz  0
 z 
Simplifying the equation (1) we get

 xx  
dxdy dz  yx dx dy dz  zx dx dy dz  X dx dy dz  0 (2)
x y z

Dividing by dxdydz to the equation (2), we will get

 xx  yx  zx
   X0 (First governing equation)
x y z
Similarly from  Fy  0
 xy  yy  zy
  Y  0 (Second governing equation)
x y z
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

and from  Fz  0
 xz  yz  zz
   Z0 (Third governing equation)
x y z

Appling conditions of moment equilibrium to prove shear stresses are

complimentary

Now, taking moment about x-axis through the centroid of the element. The
 dx dy dz 
coordinates of centroid of an element are  , ,  .
 2 2 2
Notes:
1) Moment of  xx   yy   zz  0 because they are passing through centroid.
2) Moment of  xy   yx   xz   zx  0
3) Only moment due to  yz   zy will present about x-axis.
4) Anticlockwise moment positive and clockwise moment negative

  yz  dy dy
 x  yz y dy  dx dz 2   yz dx dz 2
M   
 
(3)
  zy  dz dz
  zy  dz  dx dy   zy dx dy  0
  z  2 2
Neglecting higher power of differential coefficients  dy  , dz   in equation (3),
2 2

we get
dy dz
2 yz dx dz   zy dx dy  0
2 2
 yz   zy (Shear stress is complimentary)

Similarly, from  M y  0   xz   zx
and from  M z  0   xy   yx
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. A ' B ' and A ' C ' are the
deformed lengths of AB and AC respectively.

AB = Original line element in x-direction

A ' B ' = Deformed line element in x-direction
AC = Original line element in y-direction
A ' C ' = Deformed line element in y-direction
u = Displacement of point A in x-direction
u
u  dx = Displacement of point B in x-direction
x
v = Displacement of point A in y-direction
v
v  dy = Displacement of point C in y-direction
y
 u 
Change in length of the line element AB =  u  dx   u
 x 
u
= dx
x
Linear/Normal strain of the element AB is  x
u
dx
Change in length x u
x   
Original length dx x
Similarly
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

v w
y  and  z 
y z
The angular displacement of line element AB = 
v
dx
 v
 tan  x 
dx x
Similarly the angular displacement of line element AC = 
u
dy
y u
 tan   
dy y
v u
Total shear strain  xy     = 
x y
Similarly
w u w v
 xz   and  yz  
x z y z
Therefore, strain displacement relationship for the 3D elasticity problem is
x   u / x 
   
 y  v / y

  z   w / z 
    
 xy   u / y  v / x 
 xz  u / z  w / x 
   
 yz   v / z  w / y 
Example 1: An elastic body under the action of external forces has the
  
displacement field given by, D  2 x 2  y 2 i   5z  y  j  3x  y 2 k 
Evaluate components of strain at point (3, 1, 2)
Solution:
x   u / x   4x  12 
      1
   v / y  1
y
    
  z   w / z   0  0
      
 xy   u / y  v / x   2 y  2
 xz  u / z  w / x   3  3
       
 yz   v / z  w / y  5  2 y 3 ,1,2  7  mm / mm
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,
u  0.015 x 2 y  0.03, v  0.005 y 2  0.03xy
w  0.03x 2  0.001yz  0.005
Where displacement coordinates u, v, w and x, y, z are in constant units. Find
Cartesian strain components at point (1, 0, 2)
Solution:
x   u / x   0.03xy   0 
       0.03 
    v / y 0 . 01 y  0. 03 x
y
    
  z   w / z   0.001y   0 
      

  
xy u /  y   v / x   0. 015 x 2
 0 . 03 y   0. 015 
 xz  u / z  w / x   0.006 x   0.006 
       
 yz   v / z  w / y   0.001z 1,0 ,2 0.002  mm / mm

Example 3: In an strained elastic body under the action of external forces has the
displacement field given by
  
D  3x 3  2 y 2 i   4 z  2 y  j  4 y  z 2 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
  
D  x 2  y 2 i   2 z  y  j  3x  y 2 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
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.
u  2 x  3u
x    (Differentiating  x twice w.r.t y) (1)
x y 2 xy 2
v  2 y  3v
y    (Differentiating  y twice w.r.t x) (2)
y x 2 x 2y
v u  2 xy  3v  3u
 xy      (Differentiating  xy w.r.t x & y) (3)
x y xy yx 2 xy 2
Therefore, from equations (1)-(3) one can write
 2 x   y   xy
2 2

  (I)
y 2 x 2 xy
 2 x  2 z  2 xz
Similarly   (II)
z 2 x 2 xz
  y  2 z  2 yz
2

 2  (III)
z 2 y yz

u v  xy  2u  2v
 xy      (Differentiating  xy w.r.t z) (4)
y x z yz xz
u w  xz  2u 2w
 xz      (Differentiating  xz w.r.t y) (5)
z x y yz xy
v w  yz  2v 2w
 yz      (Differentiating  yz w.r.t x) (6)
z y x xz xy

Solving equations (4) – (6)

 xy  xz  yz  2u
  2 (7)
z y x yz
Differentiating equation (7) w.r.t x
   xy  xz  yz   3u
  
x  z x 
2
y xyz
   xy  xz  yz   2 x
   2
x  z x 
(IV)
y yz
Similarly
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

   xy  yz  xz   2 y
   2 (V)
y  z x y  xz
   xz  yz  xy   2 x
  2
z  y z 
(VI)
x xy

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
u  0.015 x 2 y  0.03
v  0.005 y 2  0.03xy
w  0.03x 2  0.001yz  0.005
Check whether the strain field given by above displacement field is compatible.

Solution:
 2 x  2 y  2 xy  2 x  2 z  2 xz
 0 ,  0 ,  0 ,  0 ,  0 ,  0,
y 2 x 2 xy z 2 x 2 xz
 2 y  2 z  2 yz
 0, 2  0, 0
z 2 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.

 x  3xy 2  3 y  x3  2
 y  4x  y 2  6x2 y 2  4
 xy  3x 2 y  34 xy 3  12 x  9 y  2
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:
x   1   0 0 0   x 
   1  0 0 0   
 y   y
  z  1      1 0 0 0   z 
    
 xy  E  0 0 0 2 1      xy 
 xz   0 0 0 2 1      xz 
    
 yz   0 0 0 2 1      yz 

It is also written in the form

1     0 0 0 
  1   0 0 0  
 x    x
     1  0 0 0   
 y   y 
 z  E  0 0 0
1  2     z 
    

  1   1  2   
2
  xy 
1  2 
xy
 xz   0 0 0   xz 
   2  
 yz  

 0 0 0
1  2    yz 
 2 
where E is Young’s modulus and  is Poisson’s ratio.

2D Hooke’s law:
Plane stress problem
 x  1  0 x 
  E    
 y   2 
 1 0
 y 
  1  
 xy   0 0 1   / 2   
 xy 
Plane strain problem
 x  1    0 x 
  E    
 y   
 1  0
 y 
  1   1  2   1  / 2  
 xy   0 0  xy 
1D Hooke’s law:   E
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

 18 0 36 
    0 54 5.4   104 mm/mm

 36 5.4 0 
Determine the stress matrix if E = 210 GPa and   0.3
Solution:
 x  0.7 0.3 0.3 0 0 0   18 
   0.3 0.7 0.3 0
 y 0 0  54 
  
 z  210  103  0.3 0.3 0.7 0 0 0  0  4
       10
 xy  1.3  0.4  0 0 0 0.2 0 0  0 
 xz   0 0 0 0 0.2 0  36 
    
 yz   0 0 0 0 0 0 .2  5 .4 
 x   145.38 
   
   1308.04
y
  145.38 0 290.71
 z   436.6  
   MPa or     0 1308.04 43.6 

 xy   0   290.71 43.6 436.6 
 xz   290.71 
   
 yz   43.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.
 zz   xz   yz  0 but,  z  0
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Equations of Equilibrium
 xx  xy   yy
  X  0 and xy  Y  0
x y x y
State of stress at a point
 xx 
 
    yy 
 
 xy 
State of strain at a point
 xx 
 
   yy  and  z  0
 
 
 xy 
Strain-Displacement relation
u v u
 x  ,  xy  
x x y
Strain compatibility condition
 2 x   y   xy
2 2

 2 
y 2 x xy
Stress-strain relation
x   1  0   x   x  1  0 x 
  1     E    
  y    1 0   y  or  y   2 
 1 0
 y 
  E 1  
 xy   0 0 2 1        
 xy   xy   0 0 1   / 2   
 xy 
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.
 2 x   y   xy
2 2

 2 
y 2 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  E  2  E   2  2 1    
y 2 1   2
 x   y   x2 1   2  y   x   xy  E  xy  (1)
    

  2 x  2 y    2 x  2 y   2 xy
 2  2  2  2 
 2 1    (2)
  y  x    x y   x y
Now, from equilibrium equations of plane stress problem neglecting body forces
 x  xy  2 x  2 xy
 0   (3)
x y x 2 xy
 xy  y  2 y  2 xy
 0   (4)
x y y 2 xy
Adding equations (3) and (4)
 2 x   y  2 xy
2

 2  2 (5)
x 2 y xy
Put equation (5) into the equation (2)
  2 x  2 y    2 x  2 y    2 x  2 y 
 2  2  2  2 
  1     2  
  y  x    x y   x y 2 
 2 x   y  2 x   y
2 2

 2  2  0
x 2 x y y 2
 2 2 
 2  2   x   y   0 (6)
  y x 
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 2 2 
  x   y   0 where   2  2   0
2 2

 y x 
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 x  2 y   2 xy  2 xy
 2  2 
 2  2 1   
  y  x  x  y xy

 2 x   y  2 xy
2

 2 (7)
y 2 x 2 xy

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

 x  3x 2  4 xy  8 y 2  4
 y  2 x 2  xy  3 y 2
 x2 
 xy     6 xy  2 y 2 
 2 
Solution:
Differential equation of equilibrium of 2D elasticity problem is given by
 x  xy  xy  y
  0 and  0 (1)
x y x y
 x  xy
 6 x  4 y,  6 x  4 y
x y
 xy  y
  x  6 y,  x  6y
x y
Putting these differential quantities in equations of equilibrium (1) and satisfying
those. This proves that the given state of stress is in equilibrium.
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.
 x  10 xy 2  Ax3
 y  1.5B xy 2
 xy   By 3  Cx 2 y
Solution:
 x  xy
 10 y 2  3 Ax 2 ,  3By 2  Cx 2
x y
 xy  y
 2Cxy,  3Bxy
x y
Putting these derivatives in following two governing equations
 
 x  xy   3 A  C  x 2  10  3B  y 2  0
x y
and
 xy  y
  2Cxy  3Bxy  0
x y
 C  3B / 2 Putting this in above equation
L.H.S. = 0 only when
3 A  C  0 and 10+3B=0
 B  10 / 3 and A=5/3

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.
 zz   xz   yz  0 but  z  0
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Equations of Equilibrium
 xx  xy   yy
  X  0 and xy  Y  0
x y x y
State of stress at a point
 xx 
 
    yy   z  0
 
 xy 
State of strain at a point
 xx 
 
    yy 
 
 xy 
Strain-Displacement relation
u v u
 x  ,  xy  
x x y
Strain compatibility condition
 2 x   y   xy
2 2

 2 
y 2 x xy

Stress-strain relation
x  1     0   x   x  1    0  x 
  1          
0    y 
E
 y    1   0   y  or  y    1 
  E  
  1   1  2   0
 xy   0 0 2   
 xy   xy   0 1   / 2   
 xy 
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:

 2 x   y   xy
2 2

 2  (1)
y 2 x xy
 2  1      2  1    
2 
y  E
 1    x   y   2   1    y   x 
 x  E 
(2)
 2  2 1    
   xy 
xy  E 
Now, from equilibrium equations of plane stress problem neglecting body forces
 x  xy  2 x  2 xy
 0   (3)
x y x 2 xy
 xy  y  2 y  2 xy
 0   (4)
x y y 2 xy
Adding equations (3) and (4)
 2 x   y  2 xy
2

 2  2 (5)
x 2 y xy
Put equation (5) into equation (2)
 2 2 
 2  2   x   y   0
 y x 
2  x   y   0
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 03 06 stresses
 xx  yx  zx
  X0
x y z
 xy  yy  zy
  Y  0
x y z
 xz  yz  zz
   Z0
x y z
Stress-strain relations 06 06 strains
 x    x   y   z  / E
 y    y   x   z  / E
 z    z   x   y  / E
 xy   xy / G
 xz   xz / G
 yz   yz / G
Strain-displacement relations 06 03 displacements
 x  u / x
 y   v / y
 z  w / z
 xy  u / y  v / x
 xz  u / z  w / x
 yz  v / z  w / y
Compatibility conditions 03 ---
 2 x   y   xy
2 2

 
y 2 x 2 xy
  x   z  2 xz
2 2
 
z 2 x 2 xz
  y   z  2 yz
2 2
 2 
z 2 y yz
2D Equilibrium equations 02 03 stresses
Plane stress  xx  yx
 X0
x y
 xy  yy
 Y  0
x y
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stress-strain relations 03 03 strains

 x    x   y  / E
 y    y   x  / E
 z     x   y  / E
1
 xy  
G xy
Strain-displacement relations 03 02 displacements
 x  u / x
 y  v / y
 xy  u / y  v / x
Compatibility conditions 01 ---
  x   y   xy
2 2 2

 
y 2 x 2 xy
2D Equilibrium equations 02 03 stresses
Plane strain  xx  yx
 X0
x y
 xy  yy
 Y  0
x y
Stress-strain relations 03 03 strains
 x    x   y   z  / E
 y    y   x   z  / E
 z     x   y 
1
 xy  
G xy
Strain-displacement relations 03 02 displacements
 x  u / x
 y  v / y
 xy  u / y  v / x
Compatibility conditions 01 ---
  x   y   xy
2 2 2

 
y 2 x 2 xy
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

ui u j
 1  1 u
Let unit displacement at node j  K   k  1 1  u i
  j

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.
 K    f 
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 ---
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 2: Element stiffness matrices

u1 u2
 1 1   1 1 u1
 K1   k1    500  1 1  u
 1 1    2
u2 u3
 1 1  1 1 u2
 K 2   k2    100  
 1 1   1 1  u3
Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix
u1 u2 u3
 500 500 0  u1 
 K   500  500  100  100 u2
 0 100 100  u3

Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is
u2 u3
 600 100  u2
K    
 100 100  u3
Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 600 100 u2  0
 100 100  u   5
  3  
 u2  0.01 m    and u3  0.06 m   
Step 6: Calculation of spring force
Spring 1:
 K 1 1   f 1
 1 1  u1   f1 
500      
 1 1  u2   f 2 
 1 1  0   f1 
500   0.01   f  ( u1  0 and u2  0.01 )
 1 1    2
f1  5 N  T and f 2  5 N  T 
Spring 2: f3  5N  T  and f 4  5N  T  ( u2  0.01 and u3  0.06 )
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

u1 u2
 1 1  1 1 u1
 K1   k1    1000  1 1  u
 1 1    2
u2 u3
 1 1  1 1 u2
 K 2   k2    2000  
 1 1   1 1  u3
u3 u4
 1 1  1 1 u
 K3   k3 1 1   3000 1 1  u3
    4
Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix
u1 u2 u3 u4
 1000 1000 0 0  u1 
 1000 1000  2000 0  u2
  2000
K   
 
0 2000  2000  3000  3000 u3
 
 0 0 3000 3000  u4 
 
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
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

u2 u3
 3000 2000  u2
K    
 2000 5000  u3
Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 3000 2000 u2   0 
 2000 5000  u   5000
  3  
Ans.  u2  0.909 m    and u3  1.363 m   
Step 6: Calculation of spring force
Spring 1:
 K 1 1   f 1
 1 1  u1   f1 
1000 
  u    f 
 1 1  2  2
 1 1  0   f1 
1000     ( u1  0 and u2  0.909 )
 1 1  0.909  f 2 
f1  909 N  T  and f 2  909 N  T 
Spring 2: f3  909 N  T  and f 4  909 N  T  ( u2  0.909 and u3  1.363 )
Spring 3: f5  4091N  C  and f6  4091N  C  ( u3  1.363 and u4  0.0 )

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

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

u1 u2
 1 1   1 1 u1
 K1   k1    500  1 1  u
 1 1    2
u2 u3
 1 1  1 1 u2
 K 2   k2    100  
 1 1   1 1  u3
u3 u4
 1 1  1 1 u
 K3   k3 1 1   200 1 1  u3
    4
Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix
u1 u2 u3 u4
 500 500 0 0  u1 
 500 500  100 0  u2
  100
K   
 
0 100 100  200  200 u3
 
 0 0 200 200  u4

Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is
u2 u3 u4
 600 100 0  u2
 K   100 300 200 u3

 0 200 200  u4
Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 600 100 0   u2   0 
 100 300 200  0.02    f 
    2
 0 200 200   u4  100 

 u2  0.00333 m    , u4  0.52 m    and f 2  98.33 N   

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
u1 u2
 1 1   1 1 u1
 K1   k1    500  1 1  u
 1 1    2
u2 u3
 1 1  1 1 u2
 K 2   k2    100  
 1 1   1 1  u3
Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix
u1 u2 u3
 500 500 0  u1 
 K   500  500  100  100 u2

 0 100 100  u3

Step 4: Reduced stiffness matrix
Imposing boundary conditions i.e. u1 = 0 eliminate first row and first column.
Therefore reduced stiffness matrix is
u2 u3
 600 100  u2
K    
 100 100  u3
Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 600 100  u2   0 
 100 100  0.06   F 
    
 u2  0.01 m    and F  5 N   
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
u1 u2
1 1 1 1 u1
 K1   k1  
1   1 1  u2
200
 1 
u2 u3
1 1 1 1 u2
 K 2   k2  
1   1 1  u3
200
 1 
u3 u4
1 1 1 1 u3
 K3   k3 1   200 
1  u4
 1  1
u4 u5
1 1 1 1 u4
 K 4   k4  
1   1 1  u5
200
 1 
Step 3: Global stiffness matrix
Assemble the element stiffness matrices to get the global stiffness matrix
u1 u2 u3 u4 u5
 200 200 0 0 0  u1 
 200 400 200 0 0  u2
 
 
K   0 200 400  200 0  u3
 
 0 0 200 400 200  u4
 0 0 0 200 200  u5

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
u2 u3 u4 u5
 400 200 0 0  u2
 200 400 200 0  u3
 K    
0 200 400 200  u4
 
 0 0 200 200  u5
Step 5: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 400 200 0 0   u2   0 
 200 400 200 0   u3   0 
   
 0 200 400 200   u4   0 
 
 0 0 200 200  0.02   F 
Ans.
u2  0.005 m    ,u3  0.01 m    ,u4  0.015m    and F  1 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
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

u1 u2
 1  1  1 1 u1 
 K1   k1    10  1 1  u
 1 1    2

u3 u4
 1 1  1 1 u3 
 K 2   k2    20  1 1  u
 1 1    4

u5 u6
 1 1  1 1 u 
 K3   k3 1 1   40  1 1  u5
    6

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
K  10  20  40  70
Step 4: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
70 u2  700
u2  u4  u6  10mm   

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.
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
u1 u2
 1 1  1 1 u1 
 K1   k1    4  1 1  u
 1 1    2

u3 u4
 1 1  1 1 u3 
 K 2   k2    8  1 1  u
 1 1    4

u5 u6
 1 1  1 1 u 
 K3   k3 1 1   20  1 1  u5
    6

u7 u8
 1 1  1 1 u
 K3   k3 1 1   10 1 1  u7
    8
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
u2 u6
 22 10  u2
K    
 10 30  u6
Step 4: Determine unknown joint displacements
Applying Equation of Equilibrium
 K    f 
 22 10 u2   0 
 10 30  u   1000
  6  
u2  u4  u7  17.857mm    u8  u6  39.286mm   
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)
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 (C0 continuity) and there is no need to
worry about slope continuity (C1 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
 u' 
 x'e   1 
u'2 
where, u1' and u'2 are the displacements in axial direction of the element. The
stiffness matrix of a bar element is
u'1 u'2
AE  1  1  u'1
 K '  
L  1 1  u'2
Transformation matrix for the truss:
x' y'= Local coordinate
systems
x, y = global coordinate
system

u'1 , u'2 = Displacements in

local coordinate system

u1 , v1 , u2 , v2 = 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
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 u1 
v 
 
 xe  1
u2 
 v2 
Refereeing above figure,
At Node 1, At Node 2,

u1'  u1cos  v1sin u'2  u2cos  v2sin

Therefore, in matrix form above relation are

 u1 
 u1'  cos sin 0 0   v1 
 '   
u2   0 0 cos sin  u2 
 v2 
x    Lx
'
e
where,
x  = vector of local unknowns
'
e

 x = vector of global unknowns

l m 0 0 
 L  = Transformation matrix =  L   
0 0 l m 
x2  x1 y  y1
where, l  cos or l  m  sin or m  2
L L
Stiffness matrix of truss element in global coordinate system
 K    L  K ' L
T

 l 0
m 0  1 1  l m 0 0 
 K     AE 
0 l  L  1 1  0 0 l m 
 
 0 m
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 l 0
  m l m 
AE  m 0   l
K   
L 0 l   l m l m 
 
 0 m
u1 v1 u2 v2
 l 2
lm l 2
lm  u1
 
AE  lm m 2 lm m 2  v1
K    2
L l lm l 2 lm  u2
 2 
  lm  m 2
lm m  v2

Example 1: Analyze the truss as shown in figure. Cross-sectional area of members

are AB=1000 mm2, BC=800 mm2, CA= 800 mm2. Take E = 2 × 105 MPa

Solution: Step 1: Degrees of freedom: 06 ( u A ,vA ,uB ,vB ,uc ,vc )

Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AB uA, vA, uB, vB uA = v A  0
2 BC uB, vB, uC, vC vB  0
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
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:
u A v A u B vB uB vB uc vc
 1 0 1 0  A
u  0.64 0 .48  0 .64 0 .48  uB
 0 0 0 0 v  0.48 0.36 0.48 0.36  v
 K AB  50 
  A
 K BC  64   B
1 0 1 0  uB  0.64 0.48 0.64 0.48  uc
   
 0 0 0 0  vB  0.48 0.36 0.48 0.36  vc
Stiffness matrix of element CA:
uc vc uA vA
 0.64 0.48 0.64 0.48  uc
 0.48 0.36 0.48 0.36  vc
 K CA  64   
0.64 0.48 0.64 0.48  u A
 
 0.48 0.36 0.48 0.36  v A
Step 3: Global stiffness matrix (Total DOF are 06, size of stiffness matrix 6×6)
uA vA uB vB uc vc
 90.96 30.72 50 0 40.96 30.72  u A 
 30.72 23.04 0 0 30.72 23.04  v A 
 
 50 0 90.96 30.72 40.96 30.72  uB
K    
 0 0 30 .72 23 .0 4  30.72  23 .04  vB 
 40.96 30.72 40.96 30.72 81.92 0  uc
 
 30 . 72 23 .04  30 .72  23 .04 0 46 .08  vc
  
Step 4: Reduced stiffness matrix (Since uA= vA  0 , vB  0 eliminate corresponding
rows and columns from global stiffness matrix)
uB uc vc
 90.96 40.96 30.72  uB
 K   40.96 81.9 0  uc

 30.72 0 46.08  vc
Step 5: Equation of equilibrium
 K    f 
 90.96 40.96 30.72  uB   0 
 40.96 81.9    
0   uc    0 
 
 30.72 0 46.08    
 vc  120 

uB  1.6mm, uc  0.8mm, vc  3.67mm
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 mm2 Young modulus is 200 KN/mm2. Determine
deflection at loaded joint.

Solution:
Step 1: Degrees of freedom: 08 ( u A ,vA ,uB ,vB ,uc ,vc ,uD ,vD )
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AD uA, vA, uD, vD u A  vA  0
2 BD uB, vB, uD, vD uB  vB  0
3 CD uC, vC, uD, vD uc  vc  0

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:
uA vA uD vD
 61.43 40.98 61.43 40.98  u A 
 40.98 27.34 40.98 27.34  v A 
 K AD    ( u A  vA  0 )
61.43 40.98 61.43 40.98  uD
 
 40.98 27.34 40.98 27.34  vD
 
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix of element BD:

uB vB uD vD
 28 .59 57 .19 28.59 57.19  uB 
 57.19 114.38 57.19 114.38 vB 
 K BD    ( uB  vB  0 )
28.59 57.19 28.59 57.19  uD
 
 57.19 114.38 57.19 114.38  vD
 
Stiffness matrix of element CD:
uC vC uD vD
 28 .59 57 .19 28.59 57.19  uC 
 57.19 114.38 57.19 114.38 vC 
 K CD    ( uc  vc  0 )
28.59 57.19 28.59 57.19  uD
 
 57.19 114.38 57.19 114.38  vD
 
Step 3: Reduced stiffness matrix
uD vD
118.61 40.98  uD
K    
 40.98 256.10  vD
Step 4: Equation of equilibrium
 K    f 
118.61 40.98  uD  200
 40.98 256.10   v    0 
  D   
uD  1.785 mm, v D  0.286 mm
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 mm2 Young modulus is 200 GPa.

Solution: Step 1: Degrees of freedom: 06 ( u A ,vA ,uB ,vB ,uc ,vc ,uD ,vD ).
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 AB uA, vA, uB, vB u A  vA  0
2 DB uD, vD, uB, vB uD  vD  0
3 CB uC, vC, uB, vB uc  vc  0

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:

uA vA uB vB
 0.64 0.48 0.64 0.48  u A 
 0.48 0.36 0.48 0.36  v A 
 K AB  40   ( u A  vA  0 )
0.64 0.48 0.64 0.48  uB
 
 0.48 0.36 0.48 0.36  vB
 
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Stiffness matrix of element DB:

uD vD uB vB
 0 .64 0 .48 0.64 0.48  uD 
 0.48 0.36 0.48 0.36  vD 
 K DB  40   ( uD  vD  0 )
0.64 0.48 0.64 0.48  uB
 
 0.48 0.36 0.48 0.36  vB
 
Stiffness matrix of element CB:
uC vC uB vB
 0.64 0 .48 0.64 0.48  uC 
 0.48 0.36 0.48 0.36  vC 
 K CB  40   ( uc  vc  0 )
0.64 0  .48 0.64 0.48  uB
 
 0.48 0.36 0.48 0.36  vB
 
Step 3: Reduced stiffness matrix
uB vB
 1.92 0.48 uB
 K   40  
 0.48 1.08  vB
Step 4: Equation of equilibrium
 K    f 
 1.92 0.48 uB  50
40 
   v   80
 0.48 1.08  B   
uB  1.25 mm, vB  2.4 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 m2 and E = 200 GPa.
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Solution:
Step 1: Degrees of freedom: 06 ( u1 ,v1 ,u2 ,v2 ,u3 ,v3 )
Discretization
Element Nodes Displacements (mm) Boundary conditions
1 1-2 u1, v1, u2, v2 u2  v2  0
2 1-3 u1, v1, u3, v3 u3  v3  0
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×102
1-3 -10 0 10 -1 0 100×102
1-4 --- --- --- --- --- ---

Step 2: Stiffness matrix of element 1-2:

u1 v1 u2 v2
 0.5 0.5 0.5 0.5  u1
 0.5 0.5 0.5 v1
2  0.5
 K 12  200  10   ( u2  v2  0 )
0.5 0.5 0.5 0.5 u2 
 
 0.5 0.5 0.5 0.5  v2 
 
Stiffness matrix of element 1-3:
u1 v1 u3 v3
 1 0 1 0  u1
 0 0 0  v1
2 0
 K 13  100  10   ( u3  v3  0 )
1 0 1 0  u3 
 
 0 0 0 0  v3 
 
Stiffness matrix of Spring element
v1 v4
 1 1 v
 K3   2000 1 1  v1 
  4

Step 3: Reduced stiffness matrix
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

u1 v1
 20000 10000  u1
K    
 10000 12000  v1
Step 4: Equation of equilibrium:  K    f 
 20000 10000 u1   0 
 10000 12000   v   40
  1  
u1  2.857 mm, v1  5.714 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 mm2 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 mm2 for
all elements.

Example: Figure shows a plane truss with two members. Both the members are of
cross-sectional area 70.71 mm2. Young’s modulus is 200 kN/mm2. Determine
deflections of loaded joint and hence the member forces.
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 mm2, BC is 400 mm2 and AC is 500
mm2. 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 mm2. Young’s modulus is 150
kN/mm2. 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 mm2 and E = 70
GPa.
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 mm2
whereas horizontal member is 1600 mm2. 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 mm2 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 mm2 and E = 200 GPa.
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 Kinematic Unknowns

Unknowns for Beam for Frame
Hinge 1 ( ) 1 ( )

Roller 1 ( ) 2 (  , )

Fixed 0 0

Spring 2 (  , ) 2 (  , )

Guided/Slider 1 ( ) 1 ( )

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

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

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

Note: 1) Action corresponding to translation is reaction

2) Action corresponding to rotation is moment
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
 K AB = EI  
0.111 0.333 0.111 0.333 3
 
 0.333 0.667 0.333 1.333  4
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
 K BC = EI  
0.1875 0.375 0.1875 0.375 5
 
 0.375 0.5 0.375 1.0  6
Step 3: Assemble global stiffness matrix
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 
 K  = EI  
 0.333 0.667 0.042 2.333 0.375 0.5  4
 0 0 0.1875 0.375 0.1875 0.375 5 
 
 0 0 0.375 0.5 0.375 1.0  6
   
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.
4 6
 2.333 0.5 4, B
 K  = EI  
 0.5 1.0  6,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.
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 75  1 
 75  2 
 75  1  50  3  
 75  2  50  4 125  3 
   
q AB   q AB    q   
 75  3  50  5  25  4

754 
50 
6  50  5 
 
50  6
Step 7: Equivalent load vector
Equivalent load vector is opposite to element nodal load vector.
F   q  Joint forces
25 0 25 4
F      
50  0 50  6
Step 8: Equation of equilibrium:
 K  F
 2.333 0.5  B  25
EI      
 0.5 1.0  C  50 
50
 B  0.0 and  C 
EI
Step 9: Reactions and Moments:
f    K   q
 RA   0.111 0.333 0.111 0.333 0 0  0  75 
M   0.333 1.333 0.333 0.667 0 0  0  75 
 A       
 R B   0.111 0.333 0.2985 0.042 0.1875 0.375  1 0 125 
   EI      
M
 B  0.333 0.667 0.042 2.333 0.375 0.5  EI 0  25
 RC   0 0 0.1875 0.375 0.1875 0.375 0  50 
       
 M C   0 0 0.375 0.5 0.375 1.0  50 50
 RA   0   75   75  kN
M   0   75   75  kN .m
 A      
 RB   18.75  125  106.25 kN
       
M
 B  25   25   0  kN .m
 RC  18.75  50   31.25  kN
       
MC   50  50  0  kN .m
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
1 2 3 4
 0.096 0.24 0.096 0.24  1
 0.24 0.8 0.24 0.4  2
 K AB = EI 
 
0.096 0.24 0.096 0.24  3
 
 0.24 0.4 0.24 0.8  4
Stiffness matrix of element BC
1 2 3 4
 0.1875 0.375 0.1875 0.375  1
 0.375 1.0 0.375 0.5  2
 K BC = EI 
 
0.1875 0.375 0.1875 0.375 3
 
 0.375 0.5 0.375 1.0  4
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 
 K  = EI  
 0.24 0.4 0.135 1.8 0.375 0.5  4
 0 0 0.1875 0.375 0.1875 0.375 5 
 
 0 0 0.375 0.5 0.375 1.0  6
   
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.
4 6
1.8 0.5  4 , B
 K  = EI  
 0.5 1.0  6, C

Step 6: Element nodal load vector

17.6  1 
 24  2 
17.6  1  60  3  
 24  2  40  4 92.4  3 
q AB    & q BC 


  q    
32.4  3  60 5  4.0  4
 36  4 40 6  60  5 
 
 40  6
Step 7: Equivalent load vector
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.
F   q  Joint forces
4  0  4.0 4
f  
 40  30  10  6
Step 8: Equation of Equilibrium:
 K  F
1.8 0.5  B  4.0 
EI      
0.5 1.0   C   10 
5.806 12.903
B  and C 
EI EI
Step 9: Moments and Reaction Calculation
f  K q
 RA   0.096 0.24 0.096 0.24 0 0   0  17.6 
M   0.24 0.8 0.24 0.4 0 0   0   24 
 A     
 R B   0.096 0.24 0.2835 0.135 0.1875 0.375  1  0  92.4 
   EI      
MB   0.24 0.4 0.135 1.8 0.375 0.5  EI 5.806  4.0 
 RC   0 0 0.1875 0.375 0.1875 0.375  0   60 
       
 M C   0 0 0.375 0.5 0.375 1.0  12.903   40 
 RA  16.207  kN
M   21.68  kN .m
 A  
 RB   96.46  kN
   
M
 B  0  kN .m
 RC  57.337  kN
   
MC   30  kN .m

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

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

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

1 2 3 4
 0.0555 0.1667 0.0555 0.1667  1
 0.1667 0.667 0.1667 0.333  2
 K AB = EI 
 
0.0555 0.1667 0.0555 0.1667  3
 
 0.1667 0.333 0.1667 0.667  4
3 4 5 6
 0.0555 0.1667 0.0555 0.1667  3
 0.1667 0.667 0.1667 0.333  4
 K BC = EI 
 
0.0555 0.1667 0.0555 0.1667  5
 
 0.1667 0.333 0.1667 0.667  6
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 
 K  = EI  
 0.1667 0.333 0 1.33 0.1667 0.333 4
 0 0 0.0555 0.1667 0.0555 0.1667  5 
 
 0 0 0.1667 0.333 0.1667 0.667  6
   
Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = 5 = zero (simple supports)
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 5: Reduced stiffness matrix

4 6
1.333 0.333  4 , B
 K  = EI 0.333
0.667  6 , C

Step 6: Element nodal load vector:

 30  1  22.22  3
 30  2  26.67  4
 
q1AB    q1BC   

 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.
6 EI  6  3800  0.025
Sinking moments  2   15.833 kN .m
L 62

 35.28  1 
 45.833  2 
 5.28  1  5.28  3  
15.833 2 15.833 4  41.66  3 
qAB   
 qBC 


 q   
 5.28 3  5.28 5  3 .33 4
15.833 4 15.833 6  13.06  5 
 
 29 .163 6
Step 7: Equivalent load vector
F   q  Joint forces
 3.33  4
F   
29.163 6
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 8: Equation of Equilibrium:

 K F
1.333 0.333  B   3.333 
EI      
0.333 0.667  C  29.163
9.45 48.189
B  C 
EI EI

Step 9: Moments and Reaction Calculation

f    K   q
 RA   0.0555 0.1667 0.0555 0.1667 0 0   0   35.28 
M   0.1667 0.667 0.1667 0.333 0 0   0   45.833 
 A     
 R B   0.0555 0.1667 0.111 0 0.0555 0.1667  1  0   41.66 
   EI      
MB   0.1667 0.333 0 1.333 0.1667 0.333  EI  9.45   3.33 
 RC   0 0 0.0555 0.1667 0.0555 0.1667   0   13.06 
       
 M C   0 0 0.1667 0.333 0.1667 0.667  48.189  29.163
 RA  33.704  kN
M  42.686  kN .m
 A  
 R B   49.693 kN
   
MB   0.00  kN .m
 RC   6.602  kN
   
 M C   0.00  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
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

1 2 3 4
0.0234 0.0937 0.0234 0.0937  1
0.0937 0.5 0.0937 0.25  2
 K AB = EI  
0.0234 0.0937 0.0234 0.0937  3
 
0.0937 0.25 0.0937 0.5  4
3 4 5 6
0.0234 0.0937 0.0234 0.0937  3
0.0937 0.5 0.0937 0.25  4
 K BC = EI  
0.0234 0.0937 0.0234 0.0937  5
 
0.0937 0.25 0.0937 0.5  6
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 
 K  = EI  
0.0937 0.25 0.1874 1.0 -0.0937 0.25  4
0 0 0.0234 -0.0937 0.0234 0.0937  5
 
0 0 0.0937 0.25 0.0937 0.5  6
   
Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = zero (simple supports), 6 = zero (guided support)
Step 5: Reduced stiffness matrix
4 5
 1.0 0.0937  4, B
 K  = EI  
 0.0937 0.0234  5, C
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 6: Element nodal load vector:

20 1  10  3 
40  2  20  4 
 
q AB    & q BC   
20  3  10  5
40  4 20  6 
Step 7: Equivalent load vector
F   q  Joint forces
 20  4
 f  
10 5
Step 8: Equation of Equilibrium:
 K F
 20   1.0 0.0937   B 
   EI   
10  0.0937 0.0234  C 
32.078 555.8
B  and  C 
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).
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

1 2 3 4
0.1875 0.375 0.1875 0.375  1
0.375 1.0 0.375 0.5  2
 K AB = EI  
0.1875 0.375 0.1875 0.375 3
 
0.375 0.5 0.375 1.0  4
3 4 5 6
1.5 1.5 1.5 1.5  3
1.5 2.0 1.5 1.0  4
 K BC = EI  
1.5 1.5 1.5 1.5 5
 
1.5 1.0 1.5 2.0  6
Stiffness matrix of spring element
5 8
 1 1 5
 K CD = EI  
 1 1  8
Step 3: Global Stiffness matrix
1 2 3 4 5 6
0.1875 0.375 0.1875 0.375 0 0  1
0.375 1.0 0.375 0.5 0 0  2
 
 0.1875 0.375 1.6875 1.125 1.5 1.0  3
 K  = EI  
 0.375 0 .5 1 .125 3 .0  1 .5 1 .0  4
0 0 1.5 1.5 2.5 1.5 5
 
 0 0 1 .5 1. 0  1 .5 2 . 0  6
  
Step 4: Impose the boundary conditions
1 = 2 = zero (Fixed support), 3 = zero (simple supports)
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 ,  B
 K  = EI  1.5 2.5 1.5 5 ,  C
 1.0 1.5 2.0  6 ,  C
Step 6: Element nodal load vector:

20 1 
20  2 
20 1 0  3  
20  2 0  4 20  3 
   
q AB   & q BC   q    
20  3 0  5 20  4
20  4 0  6 0  5
 
0  6
Step 7: Equivalent load vector
F   q  Joint forces
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

 f    0    0    0  5
 0   0   0 6
     
Step 8: Equation of Equilibrium:
 K F
10  3.0 1.5 1.0   B 
    
 0   EI  1.5 2.5 1.5  C 
 0   
   1.0 1.5 2.0  C 
4.782 2.608 0.434
B  , C  and C 
EI EI EI
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

1 2 3 4
12 6  12 6 1
6 4 6 2 2
 K AB = EI  
12 6 12 6  3
 
6 2 6 4 4
 
3 5 6 7
1.5 1.5 1.5 1.5  3
1.5 2.0 1.5 1.0  5
 K BC = EI  
1.5 1.5 1.5 1.5 6
 
1.5 1.0 1.5 2.0  7
 
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Step 5: Reduced stiffness matrix

3 4 5
 13.5  6.0 1.5  3 ,  B
 K  = EI  6.0 4.0 0  4 , BA
 1.5 0 2.0  5 , BC

Step 6: Element nodal load vector:

30 1 60  3
5  2 20  5 90  3
q AB    & q BC 


 q   5 4
30  3 60  6  20  5
5 4 20  7  

Step 7: Equivalent load vector

F   q  Joint forces
90  3
f    5  4
20  5
 
Step 8: Equation of Equilibrium:
 K F

13.5  6.0 1.5    B  90 

EI  6.0      5 
  BA   
4.0 0

 1.5 0 2.0   BC  20
20 28.75 5
B  ,  BA  and  BC 
EI EI EI
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×103 kN.m2

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.m2

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×106 mm4.

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

Example: Analyze the continuous beam ABC as shown in Figure using finite
element method. Take EI constant.
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×108 kN/m2and I = 4×10-6 m4.

Example: Determine support reactions of continuous beam ABC if support B sink

by 10 mm. Take EI = 6000 kN.m2. 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.m2. Use finite element
method.
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 6  6 .

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

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

D1 D2 D3 D4 D5 D6
 AE / L 0 0  AE / L 0 0  D1
 0 12 EI / L3 6 EI / L 2
0 12 EI / L 3
6 EI / L  D2
2

 0 6 EI / L2 4 EI / L 0 6 EI / L2 2 EI / L  D3
 K'     AE / L 
 0 0 AE / L 0 0  D4
 0 12 EI / L3 6 EI / L 2
0 12 EI / L3 2
6 EI / L D5
 
 0 6 EI / L2 2 EI / L 0 6 EI / L2 4 EI / L  D6

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.
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

D1'  D1l  D2 m
D2'   D1m  D2l
D3'  D3
At Node 2

D4'  D4l  D5m

D5'   D4 m  D5l
D6'  D6
In matrix form:
l = cosθ and m = sinθ are direction cosines.
D1 D2 D3 D4 D5 D6
 l m 0 0 0 0
 D1'   l m 0 0 0 0   D1   m l 0 0 0 0 
 '      
 D2   m l 0 0 0 0   D2   0 0 1 0 0 0 
 D3   0 0 1 0 0 0   D3  where  L    
'

 '    0 0 0 l m 0 
 D4   0 0 0 l m 0   D4   0 0 0 m l 0 
 D5'   0 0 0 m l 0   D5   
 '     0 0 0 0 0 1
 D6   0 0 0 0 0 1   D6 
x'    Lx
[L] = Transformation Matrix
x'  = Local Displacement Vector
 x =Global Displacement Vector
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 θ = 900, l=0, m=1)
[𝐾] = [𝐿]𝑇 [𝐾 ′ ][𝐿]

 L   K ' L  
T

0 1 0 0 0 0   AE / L 0 0  AE / L 0 0 
1 0 0 0 0  
0  0 12 EI / L3 6 EI / L2
0 12 EI / L 3 2
6 EI / L  

0 0 1 0 0 0  0 6 EI / L2 4 EI / L 0 6 EI / L2 2 EI / L 
  
0 0 0 0 1 0    AE / L 0 0 AE / L 0 0 
0 0 0 1 0 0  0 12 EI / L3 6 EI / L 2
0 12 EI / L3 6 EI / L2 
   
0 0 0 0 0 1  0 6 EI / L2 2 EI / L 0 6 EI / L2 4 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

(θ = 900) to reference member

 12 EI / L3 0 6 EI / L2 12 EI / L3 0 6 EI / L2 
 
 0 AE / L 0 0  AE / L 0 
 6 EI / L2 0 4 EI / L 6 EI / L2 0 2 EI / L 
K   
 12 EI / L
3
0 6 EI / L2 12 EI / L3 0 6 EI / L2 
 0  AE / L 0 0 AE / L 0 
 
 6 EI / L 4 EI / L 
2
0 2 EI / L 6 EI / L2 0
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)

Stiffness matrix for Beam Member neglecting axial deformation

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

 12 EI / L3 6 EI / L2 12 EI / L3 6 EI / L2 
 
 6 EI / L2 4 EI / L 6 EI / L2 2 EI / L 
K   
12 EI / L3 6 EI / L2 12 EI / L3 6 EI / L2 
 
 6 EI / L
2
2 EI / L 6 EI / L2 4 EI / L 

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
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
 K  AB  EI  0.048 0.24 0.048 0.24  4
 
 0.24 0.8 0.24 1.6  6
 
Imposing Boundary Conditions
1=3=0
Element Stiffness Matrix for BC: (Beam member)
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 4 1.5 2 6
 K BC  EI  
0.75 1.5 0.75 1.5  8 
 
 1.5 2 1.5 4 9
 
Imposing Boundary Conditions
5=8=0
Element Stiffness Matrix for DC: (Column member)
10 12 7 9
 0.096 0.24 0.096 0.24  10 
 0.24 0.8 0.24 0.4  12 
 K DC  EI  0.096 0.24 0.096 0.24  7

 
 0.24 0.4 0.24 0.8  9
 
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, 
[ K ]  EI  0.24 5.6 2  6, B
 0.24 2 4.8  9, C

Step 4: Element Nodal Load Vector

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

 3.52  1
 9.6  3
 
{q AB }   
  6.48 4
 14.4  6
 6 5
 4 6
 
{qBC }   
 6 8
 4  9
 3.24 10
3.6 12
 
{qDC }   
 1.76 7
 2.4  9
Reduced element nodal load vector
4.72  4
 
{q}  10.4  6
 1.6  9
 
Step 4: Equivalent Load Vector
F  q  Joint forces
4.72  4
F   10.4  6
 1.6  9
 
Step 5: Equation of Equilibrium
[ K ]{}  {F }
0.144 0.24 0.24     4.72 
   
EI  0.24 5.6 2   B   10.4 
 
 0.24 2 4.8  C   1.6 

34.046 1.0419 1.803

 m;  B  rad ; C   rad
EI EI EI

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

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

 0 
 
 M AB   0.24 1.6 0.24 0.8 1  0   9.6  18.604 
   EI      
 M BA   0.24 0.8 0.24 1.6  EI 34.046 14.4 4.562
1.0419 
Member BC
 0 
 
 M BC  1.5 4 1.5 2  1 1.0419   4   4.562 
   EI      
 M CB  1.5 2 1.5 4  EI  0  4  9.128

1.803

Member DC
 0 
 
 M DC   0.24 0.8 0.24 0.4 1  0  3.6 3.849
   EI   EI 34.046    2.4   9.128
 CD 
M  0.24 0.4 0.24 0.8       

1.803
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
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 
 K AB  EI  
0.1875 0.375 0.1875 0.375  4
 
 0.375 0.5 0.375 1 6
 
Element Stiffness Matrix of beam BC
5 6 8 9
 0.0469 0.1875 0.0469 0.1875  5 
 0.1875 1 0.1875 0.5  6
 K BC  EI 
 
0.0469 0.1875 0.0469 0.1875 8 
 
 0.1875 0.5 0.1875 1 9
 
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 
 K DC  EI  
0.1875 0.375 0.1875 0.375  7
 
 0.375 0.5 0.375 1 9
 
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, 
[ K ]  EI 0.375 2 0.5  6,  B
 
0.375 0.5 2  9,C

Step 3: Element Nodal Load Vector

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

0  1 37.5 5
0  3  75  6
   
{q AB }    {qBC }   
0
  4  37.5 8
0  6  75  9
0 10
0 12
 
{qDC }   
0  7
0  9
Reduced element nodal load vector
 0 4
q   75  6
75 9
 
Step 4: Equivalent Load Vector
F  q  Joint forces
 0  50  50  4
F   75   0   75 6
 75   0   75  9
     
Step 5: Equation of Equilibrium

0.375 0.375 0.375     50 

   
EI 0.375 2 0.5   B   75
 
0.375 0.5 2  C   75 
−𝟕𝟖.𝟓𝟕𝟏 𝟐𝟏.𝟒𝟐𝟖 𝟏𝟗𝟎.𝟒𝟕𝟔
𝜽𝑩 = 𝒓𝒂𝒅 𝜽𝑪 = 𝒓𝒂𝒅 ∆= 𝒎
𝑬𝑰 𝑬𝑰 𝑬𝑰

Step 6: Moment Calculations

{f} = [K]{Δ}+{q}
Member AB
 0 
 
 M AB   0.375 1 0.375 0.5 1  0  0  32.143 
   EI      
 M BA   0.375 0.5 0.375 1  EI 190.476  0 7.143
78.571
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Member BC
 M BC   1 0.5 6 1 78.571  75   7.143 
   EI   9 EI  21.428   75  92.857 
 CB 
M  0.5 1       
Member DC
 0 
 
 M DC   0.375 1 0.375 0.5 1  0  0 82.143
   EI   EI 190.476  0  92.857 
 CD 
M  0.375 0.5 0.375 1       
 21.428 

Example 3: Assemble the stiffness matrix for the frame shown in Fig. using finite
element method. Take AE = 400000KN and EI = 1000 kNm2 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 ( θ=900, 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 
[ L]T [ K ']    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 
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
 L  K ' L   1000 
T
 
 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
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
 K AB  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

Element stiffness matrix of BC

L=5.6569 m, θ = 450, l=0.7071, m=0.7071

 L  K ' 
T

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.30 
   
 0 0 0 0 0 1  0 0.30 0.565 0 0.30 1.131 

[ L]T [ K '][ L] 
 50 0.075 0.2121 50 0.075 0.2121  0.7071 0.7071 0 0 0 0
 50 0.075  
0.212 50 0.075 0.2121 0.7071 0.7071 0 0 0 0
  
 0 0.3 1.131 0 0.3 0.565   0 0 1 0 0 0
1000   
 50 0.075 0.2121 50 0.075 0.2121   0 0 0 0.7071 0.7071 0
 50 0.075 0.2121 50 0.075 0.2121  0 0 0 0.7071 0.7071 0
  
 0 0.3 0.565 0 0.3 1.131   0 0 0 0 0 1
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 4
 35.302 35.408 0.2121 35.302 35.408 0.2121  5

 0.2121 0.2121 1.131 0.2121 0.2121 0.565  6
 K BC  1000  
 35.408 35.302 0.2121 35.408 35.302 0.2121  7
 35.302 35.408 0.2121 35.302 35.408 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 1
 0 100 0 0 100 0 0 0 0  2

 0.6 0 1.6 0.6 0 0.8 0 0 0 3
 
 0.3 0 0.6 35.708 35.302 0.388 35.408 35.302 0.2121 4
[ K ]  1000  0 100 0 35.302 135.408 0.2121 35.302 35.408 0.2121  5
 
 0.6 0 0.8 0.388 0.2121 2.731 0.2121 0.2121 0.565  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.2121 0.2121 0.2121 1.131  9
 0 0 0.565 0.2121

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  1.2 / EI, MAB  0.8kN.m, MBA  1.6kN.m, MBC  18.4kN.m, MCB  14.8kN.m )
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 ×103 KN.m2. 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.
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.
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

Moments in this figure are about

z-axis, no moment about x-axis.
Reactions are along y direction

Unit rotation in x-direction at node 1 Moments in this figure are about x-

axis i.e. twisting moment. Reactions
along y direction and moments along
z-direction are zero
Unit rotation in z-direction at node 1
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

Moments in this figure are about

z-axis, no moment about x-axis.
Reactions are along y direction

Unit displacement in y-direction at node 2

Moments in this figure are about

z-axis, no moment about x-axis.
Reactions are along y direction

Unit rotation in x-direction at node 2

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

Unit rotation in z-direction at node 2

Moments in this figure are about

z-axis, no moment about x-axis.
Reactions are along y direction

Stiffness matrix for grid element in local coordinate system from above six figures
D1 D2 D3 D4 D5 D6
 12 EI / L3 0 6 EI / L2 12 EI / L3 0 6 EI / L2 
 
 0 GJ / L 0 0 GJ / L 0 
 6 EI / L2 0 4 EI / L 6 EI / L2 0 2 EI / L 
K   
 12 EI / L 6 EI / L2 6 EI / L2 
3
0 12 EI / L3 0
 0 GJ / L 0 0 GJ / L 0 
 
 6 EI / L
2
0 2 EI / L 6 EI / L2 0 4 EI / L 
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
D1'  D1
D2'  D2 cos  D3 sin 
D3'   D2 sin   D3 cos

At Node 2
D4'  D4
D5'  D5 cos  D6 sin 
D6'   D5 sin   D6 cos

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

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

 D1'  1 0 0 0 0 0   D1  1 0 0 0 0 0 
 '    0 l m 0 0 0 
 D2  0 l m 0 0 0   D2   
 D3'  0 m 1 0 0 0   D3  0 m 1 0 0 0 
 '     L   
D
 4  0 0 0 1 0 0   D4  0 0 0 1 0 0 
 D5'  0 0 0 0 l m   D5  0 0 0 0 l m 
 '     
 D6  0 0 0 0 m l   D6  0 0 0 0 m l 

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 (θ=900, l = 0, m = 1)
[𝐾] = [𝐿]𝑇 [𝐾 ′ ][𝐿]

 12 EI 6 EI 12 EI 6 EI 
 L3 0  0
L2 L3 L2 
 
 GJ

GJ
0 
1 0 0 0 0 0  0
L
0 0
L 
0 0 1 0 0 0   6 EI 
  0
4 EI

6 EI
0
2 EI 
0 1 0 0 0 0   L2 L2 L 
 L   K    0 
L
T 1

 0 0 1 0 0 12 EI 6 EI 12 EI 6 EI 
 3 0  2 0  2 
0 0 0 0 0 
1  L L L3 L 
   GJ GJ 
 0 0 0 0 1 0   0  0 0 0 
 L L 
 6 EI 2 EI 6 EI 4 EI 
 L2 
L 
0 0
L L2
 12 EI 6 EI 12 EI 6 EI 
 L3 0  3 0
L2 L L2 
 
  6 EI 
4 EI 6 EI

2 EI 
L  
 L2
0 0 1 0 0 0 0 0
L L2
  0 0 1 0 0 0 
0  
 0 GJ GJ
0 0 
 L L  0 1 0 0 0 0
    
12 EI 6 EI 12 EI 6 EI 0 0 0 1 0 0
 3 0  0  2  
 L L2 L 3
L  0 0 0 0 0 1
 6 EI 2 EI 6 EI 4 EI   0 
0 
 2 0  0    0 0 0 1
 L L L2 L 
 GJ GJ 
 0  0 0 0 
L L 
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

 12 EI / L3 6 EI / L2 0 12 EI / L3 6 EI / L2 0 
 
 6 EI / L
2
4 EI / L 0 6 EI / L2 2 EI / L 0 
 0 0 GJ / L 0 0 GJ / L 
K   
 12 EI / L
3
6 EI / L2 0 12 EI / L3 6 EI / L2 0 
 6 EI / L2 2 EI / L 0 6 EI / L2 4 EI / L 0 
 
 0 0 GJ / L 0 0 GJ / L 

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)
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

1 2 3 4 5 6
 12 EI / L3 0 6 EI / L2 12 EI / L3 0 6 EI / L2 1
 
 0 GJ / L 0 0 GJ / L 0  2
 6 EI / L2 0 4 EI / L 6 EI / L2 0 2 EI / L  3
 K AB  
 12 EI / L 6 EI / L2 6 EI / L2
3
0 12 EI / L3 0 4
 0 GJ / L 0 0 GJ / L 0 5
 
 6 EI / L
2
0 2 EI / L 6 EI / L2 0 4 EI / L  6
  

Stiffness matrix for member CB:

Rotate member CB in clockwise
direction about joint C at an angle of 900
w. r. t to positive x-axis. (θ=900, l=0,
m=1). Therefore take ‘C’ as a first node
and ‘B’ as a second node.

7 8 9 4 5 6
 12 EI / L3 6 EI / L2 0 12 EI / L3 6 EI / L2 0 7
 
 6 EI / L 8
2
4 EI / L 0 6 EI / L2 2 EI / L 0
 0 0 GJ / L 0 0 GJ / L 9
 K CB  
 12 EI / L
3
6 EI / L2 0 12 EI / L3 6 EI / L2 0 4
 6 EI / L2 2 EI / L 0 6 EI / L2 4 EI / L 0 5
 
 0 0 GJ / L 0 0 GJ / L  6
  

III) Reduced stiffness matrix

4 5 6
 24 EI / L3 6 EI / L2 6 EI / L2  4,  Bz
 K    6 EI / L2  4 EI / L    GJ / L  0

 5,  Bx

 6 EI / L2 0  4 EI / L    GJ / L   6,  By
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  1 
 0 0.08 0 0 0.08 0  2 
 
 0.24 0 0.8 0.24 0 0.4  3 
 K AB  EI  
 0.096 0 0.24 0.096 0 0.24  4
 0 0.08 0 0 0.08 0  5
 
 0.24 0 0.4 0.24 0 0.8  6
  
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 900
w. r. t to positive x-axis. (θ=900, 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  4
 0.667 1.333 0 0.667 0.667 0  5
 
 0 0 0.133 0 0 0.133 6
 K BC  EI  
 0.444 0.667 0 0.444 0.667 0  7 
 0.667 0.667 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,  Bz
 K   EI  0.667 1.413 0  5,  Bx
 0.24 0 0.933  6,  By

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
    0  4
0  3 0  6
q AB     and qBC      q  0 5
0  4 0  7 0  6
0  5 0  8  
   
0  6 0  9

Step 4: Equivalent load vector

 f   q  Joint forces
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

0 70  70  4

 f   0   0    0  5
0   0   0  6
     
Step 5: Equation of Equilibrium
[K]{Δ} = {f}

 0.540 0.667 0.24   Bz  70 

   
EI  0.667 1.413 0   Bx    0 
 
 0.24 0 0.933   By   0 
428.371 202.21 110.192
 Bz  ;  Bx  ;  By 
EI EI EI

Step 6: Moment Calculations

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

 428.371
 0.667 1.333 0 0.667 0.667 0   202.21  0   M Bx 
 0  
 0 0.133 0 0 0.133 110.192  1 0   M By 
EI     
 0.667 0.667 0 0.667 1.333 0  0  EI 0   M Cx 
 
 0 0 0.133 0 0 0.133   0  0   M Cy 
 
 0 
 M Ax   16.176   M Bx   16.176 
M    M   
 Ay   58.732   By   14.655 
   and   
M
 Bx    16.176   Cx   150.84 
M
 M By   14.655   M Cy   14.655 
   
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.m2 and GJ = 800 kN.m2

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 900 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.
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
 K BC  
 300 0  600 300 0  600 7
 0 200 0 0 200 0 8
 
 600 0 800 600 0 1600  9 
  

Stiffness matrix for member AB:

Rotate member AB in clockwise
direction about joint A at an angle of 900
w. r. t to positive y-axis. (θ=900, 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 
 K AB   
 300 600 0 300 600 0 4
 600 800 0 600 1600 0 5
 
 0 0 200 0 0 200  6
  
Reduced stiffness matrix is
4 5 6
 600 600 600  4,  Bz
 K   600 1800 0  5,  By
 600 0 1800 6,  Bx
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
qBC     and q AB    
 20  7   60  4
 0 8   40  5
   
20  9   0 6

80  4, RBz
q  40 5, M By
20  6, M
  Bx

Step 4: Equivalent load vector

 f   q  Joint forces
80  4
 f   40 5
20  6
 
Step 5: Equation of Equilibrium
[K]{Δ} = {f}
600 600 600   Bz  80
600 1800    
0   By   40
 
600 0 1800  Bx  20
 Bz  0.3 m;  Bx  0.0888 rad ;  By  0.0777 rad
Step 6: Moment Calculations
 K   q   f 
Element AB
 0 
 600 1600 0 600 800 0   0   0   M Ay 
 
 0 0 200 0 0 200   0   40   M Ax 
    
 600 800 0 600 1600 0   0.3   0   M By 
 
 0 0 200 0 0 200  0.0888  40  M Bx 
 
0.0777 
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   M By 
 
600 0 1600 600 0 800  0.0777  1  0   M Bx 
EI      
 0 200 0 0 200 0   0  EI 20   M Cy 
 
600 0 800 600 0 1600  0   0   M Cx 
 
 0 
Element BC
 M Ay   17.76   M By   17.76 
       
 M Ax   157.84   M Bx   15.68 
   and   
M
 By   17.76   Cy   128.96 
M
 M   15.68   M   15.68 
 Bx   Cx 

Example 4: Using finite element method, determine unknown joint displacements

of the grid as shown in figure. E=2×105 MPa, I = 20×105 mm4, G = 0.8×105 MPa, J
=50×105 mm4

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×105 × 20×105 = 40×1010 N.mm2 = 400 kN.m2
GJ = 0.8×105 × 50×105 = 40×1010 N.mm2 = 400 kN.m2
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 900 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
 K BA   
 600 0 600 600 0 600  1 
 0 200 0 0 200 0  2
 
 600 0 400 600 0 800  3 
  

Stiffness matrix for member BC:

Rotate member BC in clockwise
direction about joint B at an angle of 900
w. r. t to positive y-axis. (θ=900, l=0,
m=1). Therefore take ‘B’ as a first node
and ‘C’ as a second node.
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 4
 266.67 533.33 0 266.67 266.67 0 5
 
 0 0 133.33 0 0 133.33 6
 K BC   
 177.77 266.67 0 0.177 266.67 0 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
4 5 6
 777.77 266.67 600  4  Bz
 K BC   266.67 733.33 0  5  By
 600 0 933.33  6  Bx

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

0  4  6 4
0  5 3 5
     6 4
0  6  0 6
qBA    and qBC      q  3 5
0  1  6 7  0 6
0  2  3 8  
   
0  3  0 9
Step 4: Equivalent load vector
 f   q  Joint forces
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

6  4
 f    3  5
 0 6
 
Step 5: Equation of Equilibrium
[K]{Δ} = {f}

 777.77 266.67 600    Bz  6 

 266.67 733.33    
 0    By    3 
 600 0 933.33    Bx   0 

 Bz  0.0166 m;  By  0.00195 rad ;  Bx  0.0106 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×105 kN.m2 and GJ = 1.2×105 kN.m2.
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 m4, J =
4.6 X10-5 m4 for all elements.
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
1
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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.

2
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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 1. FEM gives values at any nodes
only. At other point interpolation is points including node points.
required
2. FDM needs larger number of nodes 2. FEM needs fever nodes to get
to get good results. good results.
3. Fairly complicated problems can be 3. All type of complicated problem
handled by FDM. can be handled by FEM.
4. FDM makes stair type of 4. FEM can consider sloping and
approximation to sloping and curved boundaries exactly.
curved boundaries.
5. FDM makes point wise 5. FEM makes piece wise
approximation i.e. satisfy continuity approximation i.e. Satisfy
at node points only, along the side continuity at node point as well as
continuity are not ensured. along the side/edge of element.
6. It is less efficient and more 6. It is more efficient and more
approximate. approximate.
7. Not applicable for non-linearity of 7. Applicable to non-linearity of
domain. domain.
8. Difficult to apply FDM for unusual 8. FEM can be applied to any type of
boundary condition and loading boundary condition and loading
conditions. 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:

3
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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.

4
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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.
Wi  We
Internal workdone  Wi     ij  ij dv
dv

External workdone  We    q  wds

ds

 dv
 ij  ij dv   q  wds
ds

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.
 U V
where
 = Total potential energy
U =Strain energy (due to internal forces)
V = Work potential (due to external forces)
1
U     and V  q w
2

5
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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 3 4 2

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
B  (1  D) f
Where,
D = maximum difference in node number in an element after considering all
elements.
f = degree of freedom per node.

6
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
B  (1  D) f
For two dimensional problem, f = 02 (DOF)
Case I) D = 13, f = 02  B  1  13 2  28
Case II) D = 8, f = 02  B  1  8 2  18
Case III) D = 9, f = 02  B  1  9 2  20
Case VI) D = 6, f = 02  B  1  6 2  14
 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.

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

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  B   (1  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
B  (1 D) f
 (1  6)2
 14
Case II) Yes, it is possible to minimize bandwidth.

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

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.

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

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.
10
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

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:  K e  e  Qe
Where,  K e = element stiffness matrix,
 e = nodal displacement vector of element,
Qe = 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:  K    Q
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.

11
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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

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

Applications
3D Elements: Three-dimensional analysis, analysis of axisymmetric solids.
C0 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  w / x and  x  w / y . If
Kirchhoff’s assumption is not made, slopes are independent of deflection and
hence w,  x ,  y are nodal unknowns reduced to C0 continuity requirement.

C1 continuity elements:
First order continuity elements in which higher order derivative of ‘ w ’ is one only.
C1 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.
w w
w, , .
x y
C2 continuity elements:
Second order continuity elements in which second derivatives of ‘ w ’ are also
w w  2 w  2 w  2 w
nodal unknowns i.e. w, , , 2 , 2 , .
x y x y xy

1.17 Difference of CST and LST Elements

13
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Constant Strain Triangle (CST) elements:

Displacement field in terms of generalized coordinates

u  1   2 x  3 y (Continuous across element)
Resulting strain field is
u
x    2 (Constant) ( NOT continuous across element)
x

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

u  1   2 x  3 y   4 x 2  5 xy  6 y 2
Resulting strain field is
u
x    2  2 4 x   5 y (Linear)
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

14
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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

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

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: u  1   2 x

16
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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: w  1   2 x  3 x   4 x
2 3

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
u  1   2 x   3 y
Displacement Function:
v   4  5 x  6 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:
u  1   2 x  3 y   4 x 2  5 xy  6 y 2
v  7  8 x  9 y  10 x2  11xy  12 y 2

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
u  1   2 x   3 y   4 xy
Displacement Function:
v   5   6 x   7 y  8 xy

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

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

DOF per node: 03 w, x , y , 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
Displacement Function: y coordinates both
w  1   2 x   3 y   4 x 2   5 xy   6 y 2
 7 x3  8 x 2 y   9 xy 2  10 y 3  11x3 y  12 xy 3

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:
u  1   2 x   3 y   4 z
v  5   6 x   7 y  8 z
w   9   10 x  11 y  12 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 x1 and x2 respectively. Cartesian coordinate of any point ‘P’ is ‘ x ’.
Natural coordinates of any point ‘P’ are (L1, L2).
18
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

From the definition of line element

L1  L2  1
L1x1  L2 x2  x
In matrix form,
 L1  1 1  1 
   
 L2  x1 x2   x 
1
 L1   1 1  1  1  x2 1 1  1 ( x2  x) 
       
 L2   x1 x2   x  x2  x1   x1 1   x  L  ( x  x1 ) 

x x x  x1
 L1  2 and L2 
L L
At node 1 x  x1 ,  L1  1, L2  0
At node 2 x  x2 ,  L1  0 , L2  1
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)

19
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
OP
 
L/2
2   x  x  2   x  x  2 x1  
   x   1 2      x   2 1 
L   2  L  2 
2   L  2 x1  
   x   
L   2  
At Node 1, x  x1
2  L  2 x1   2  2 x1  L  2 x1 
   x1       
L  2  L 2 
  1
At Node 2, x  x2
2  L  2 x1   2  2 x2  L  2 x1 
  x       
L   L  
2
 2  2
2  2L  L 
 
L  2 
  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

20
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
(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
L1  L2  L3  1
L1 x1  L2 x2  L3 x3  x (1)
L1 y1  L2 y2  L3 y3  y
In matrix form,
1
 L1  1 1 1  1   L1   1 1 1 1 
       
 L2   x1 x2 x3   x    L2    x1
 x2 x3  x
  
L   y y2 y3   y  L   y y3   y
 3  1  3  1 y2  
 L1   a1 b1 c1   1 
  1   
 L2   a2 b2 c2   x  (2)
 
 L  D a b c   y 
 3  3 3 3 
where
D   x2 y3  x3 y2    x3 y1  x1 y3    x1 y2  x2 y1 
a1  x2 y3  x3 y2 b1  x3 y1  x1 y3 c1  x1 y2  x2 y1
(3)
a2  y3  y2 b2  y3  y1 c2  y1  y2
a3  x3  x1 b3  x1  x3 c3  x2  x1
Area of triangle
Note: It is well known that the two times area of triangle is equals to its
determinant ( D  2 A ). The following is the proof of this.

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

Area of  = Area of 1354+ Area of 3265 – Area of 1264

distance between two parallel sides
Area of trapezoidal = Sum of parallel sides  
2
1 1 1
A   y1  y3  x3  x1    y3  y2  x2  x3    y1  y2  x2  x1 
2 2 2
1
A   y1 x3  y3 x1  y3 x2  y2 x3  y1 x2  y2 x1 
2
2A  D
Putting this in equation (2), we get

 L1   a1 b1 c1   1 
  1   
 L2   a2 b2 c2   x 
 
L  2 A a b c   y 
 3  3 3 3 
a  b x  c1 y a  b2 x  c2 y a  b3 x  c3 y
L1  1 1 ; L2  2 and L3  3 (4)
2A 2A 2A

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

22
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
1 1 1
D x x2 x3
y y2 y3
D   x2 y3  x3 y2    xy3  x3 y    xy2  x2 y   2 A1
 a1  b1 x  c1 y  2 A1
Similarly, if any point P is shifted to node 2 and 3, we get
a2  b2 x  c2 y  2 A2 and a3  b3 x  c3 y  2 A3
Therefore, from equation (4) natural coordinates are written as
2 A1 A1 2 A2 A2 2 A3 A3
 L1    L2   and L3  
2A A 2A A 2A 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.
23
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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:
u  1   2 x   3 y   4 z
v  5   6 x   7 y  8 z
w   9   10 x  11 y  12 z

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

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)
u  1   2    3   4   5
 6   7  8   5

25
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
u  1   2 x
In matrix form
 
u  1 x  1 
 2 
u   P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements

26
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Express displacement function in terms of nodal displacements using the
coordinates of nodes x1 and x2.
u1  1 x1  1 
   
u2  1 x2   2 
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe
where [N] = Shape functions
 N    P A
1

III) Shape functions

1
1 x 
 N    P A  1 x   1 
1

1 x2 
Inverse is obtained by using method of adjoin
1  x  x1 
 N    P A  1 x  2
1

2  1 1 
 N1  1  x2  x 
   
 N 2  2  x  x1 
x x x  x1
N1  2 and N 2 
L L
Sum of shape functions is always unity
x  x x  x1 x2  x1 L
N1  N 2  2    1
L L L L

At node 1, x=x1
x x L x x
N1  2 1   1 and N 2  1 1  0
L L L
At node 2, x=x2
x  x2 x x L
N1  2  0 and N 2  2 1   1
L L L

27
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
w  1   2 x  3 x 2   4 x3 (From Pascal Triangle)
dw
   2  2 3 x  3 4 x 2
dx
In matrix form
1 
 
 w 1 x x 2
x   2 
3

  2 
  0 1 2 x 3 x   3 

 4 

    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.
 w1  1 0 0 0  1 
   0 1 0 0  
 1    2 
   
   3 
2 3
w
  2 1 L L L
  2
 2   0 1 2 L 3L    4 

xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
w   P A xe
1

w   N xe
where [N] = Shape functions
28
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

III) Shape functions

 1 0 0 0 
 0 1 0 0 
 N    P  A  1 x x x   3  
1 2

3 / L 2 / L 3 / L
2 2
1 / L 
 
 2/ L
3
1 / L2 2 / L3 1 / L2 
Inverse of [A] is obtained using elementary operations of matrix algebra.
3x 2 2 x3 2 x 2 x3
N1  1  2  3 , N2  x   2,
L L L L
2 3 2 3
3x 2x x x
N3  2  3 , N4    2
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
u  1   2 x  3 y
In matrix form
1 
 
u  1 x y   2 
 
 3
u   P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements

29
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Express displacement function in terms of nodal displacements using the
coordinates of nodes (x1, y1), (x2, y2), (x3, y3).
 u1  1 x1 y1  1 
    
u2   1 x2 y2   2 
u  1 x y   
 3  3 3 3

xe   A  (2)

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

u   N xe
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
1 x1 y1 
 N    P  A  1 x y  1 x2 y2 
1

1 x3 y 
 a1 b1 c1 
 N   1 x y  a2 b2 c2 
1 
2A
 a3 b3 c3 
a  a2 x  a 3 y b  b x  b3 y c  c x  c3 y
N1  1 N2  1 2 N3  1 2
2A 2A 2A
where
a1  x2 y3  x3 y2 a2  y2  y3 a3  x3  x2
b1  x3 y1  x1 y3 b2  y3  y1 b3  x1  x3
c1  x1 y2  x2 y1 c2  y1  y2 c3  x2  x1

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
 N    P A
1

where
30
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
1 x1 y1 
 P   1 x y  A  1 x2 y2 

1 x3 y3 
1
1 0 0 
 N   1 x y  1 4 0 
1 2 2 
Inverse is obtained using method of adjoin
 8 0 0
 N   1 x y  2 2 0
1
8
 2 2 4 
8  2x  2 y 2x  2 y 4y
N1  N2  N3 
8 8 8

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
 N    P A
1

1
1 1 2 
 N    P  A  1 x y  1 5 3 
1
 
1 4 6 

18 2 7 
 N   1 x y  3 4 1
1
13
 1 3 4 
18  3x  y 2  4x  3 y 7  x  4 y
N1  N2  N3 
13 13 13
Now, at x = 3.3 N1 = 0.3
18  3  3.3  y
0.3   y  4.2
13
Using x = 3.3 and y = 4.2
2  4  3.3  3  4.2 7  3.3  4  4.2
N2   0.2 N3   0.5
13 13

31
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
u  1   2 x  3 y   4 xy
In matrix form
1 
 
 
u  1 x y xy   2 
 3 

 4 

u   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), (a, 0), (a, b), (0, b).
 u1  1 0 0 0  1 
u  1 a 0 0   
 2   2
   
u3  1 a b ab   3 
  
u4  1 0 b 0    4 

xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
32
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603

u   P A xe

1

u   N xe
where [N] = Shape functions

III) Shape functions

1
1 0 0 0 
1 a 0 0 
 N    P  A  1 x y xy   
1

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

1 / b 0 0 1/ b 
 
1 / ab 1 / ab 1 / ab 1 / ab 
x xy y x xy
N1  1    ; N2   ;
a ab b a ab
xy y xy
N3  ; N2  
ab b ab

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

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

33
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 1 0 0 0 
 1 / 6 1 / 6 0 0 
 N    P  A  1 x y xy   
1

1 / 4 0 0 1/ 4 
 
1 / 24 1 / 24 1 / 24 1 / 24 
x xy y x xy
N1  1    ; N2   ;
6 24 4 6 24
xy y xy
N3  ; N2  
24 4 24

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 coordinate   1 and node 2
has coordinate   1

I) Displacement function in-terms of natural coordinate

(use Pascal triangle in natural coordinates)
u  1   2
 
u  1    1 
 2 
u   P (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.
u1  1  1 1 
   
u2  1 1   2 
xe   A 
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe
34
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
where [N] = Shape functions
III) Shape functions
1
1 1
 N    P A  1   
1

1 1 
 N1  1 1   
   
 N 2  2 1   
1 1
N1  and N 2 
2 2
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 coordinate   1 , node 2 has
coordinate   0 , and node 3 has coordinate   1

I) Displacement function in-terms of natural coordinate

(Use Pascal triangle in natural coordinates)
u  1   2  3 2
1 
u  1   2   2 
 
 3
u   P (1)
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of node.
u1  1 1 1  1 
    
u2   1 0 0   2 
u  1 1 1   
 3   3
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe
35
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
where [N] = Shape functions
III) Shape functions
1
1 1 1 
 N    P  A  1   2  1 0 0
1

1 1 1 
 0 2 0
 N   1   2  1 0 1 
1
2
 1 2 1 

N1 
   1  ; N 2   1   2  and N3 
1    
2 2

2. Four noded rectangular element

I) Displacement function
u  1   2  3   4
1 
 
u  1      2 
 3 

 4 

u   P (1)
II) Displacement function in-terms of nodal displacements
Express displacement function in terms of nodal displacements using the
coordinates of node.
 u1  1 1 1 1  1 
u  1 1 1 1  
 2   2
   
u3  1 1 1 1   3 
  
u4  1 1 1 1   4 

36
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe
where [N] = Shape functions
III) Shape functions
1 1 1 1
 1 1
1  1 1
 N    P  A  1      
1

4 1 1 1 1
 
 1 1 1 1

N1 
 1   1    N2 
1   1   
4 4

N3 
 1    1    N2 
1   1   
4 4

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
n
  m x  xm
Nk   or
m1,m k  m   k xm  xk
1. Two noded bar element

  2   1 1     1   1 1  
N1    N2   
1  2 1  1 2 2  1 1  1 2

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

   2   3   0   1    1
N1     
1   2 1  3 1  0 1  1 2
  1   3   1   1
N2      1   1   
 2  1  2  3 0  1 0  1
   2   1   0   1    1
N3     
3   2 3  1 1  0 1  1 2

3. Four noded rectangular element

   2   4   1   1 1   1   
N1     
1   2 1  4 1  1 1  1 4
  1   3   1   1 1   1   
N2     
 2  1 2  3 1  1 1  1 4
   4   2   1   1 1   1   
N3     
 3   4 3   2 1  1 1  1 4
  3   1   1   1 1   1   
N4     
 4  3 4  1 1  1 1  1 4
Note: Generally shape functions for four noded rectangular element are given by
1  k 1  k 
Nk  
4
38
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
4. Nine noded rectangular element

Shape functions for corner nodes (1, 2, 3, 4):

   5    2    4   8   0   1   0   1
N1        
1  5 1   2 1  4 1  8 1  0 1  1 1  0 1  1

 N1 
1   1   
4
  5   1   6   3 1     1
N2     
2  5 2  1 2  6 2  3 4
  7   4   6   2 1   1   
N3     
3  7 3  4 3  6 3  2 4
  7   3   8   1   11   
N4     
4  7 4  3 4  8 4  1 4
In general shape functions for corner nodes are given by
Nk  
1  k 1  k kk
4
Shape functions for middle nodes (5, 6, 7, 8):
  1    2   9   7 1      1
2

N5     
5  1 5   2 5  9 5  7 2

1    1   2  1      1 ;
2
  1 1   2 
N6  ; N7  N8 
2 2 2
In general shape functions for middle nodes are given by

Nk  
1   2  1  k k
 at nodes   0 axis 
2
1  k  1   2 k
Nk    at nodes   0 axis 
2
39
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Shape function for central node:
   6   8    7   5   1   1   1   1
N9        
 9   6  9   8  9   7 9  5 0  1 0  1 0  1 0  1
 N9  1   2 1   2 
5. Eight noded hexahedron element

Shape functions using Lagrange interpolation function

  5   2    4   1   1   1 1   1   1   
N1       
1  5 1  2  1   4 1  1 1  1 1  1 8
  6   1    3 1   1   1   
N2    
 2  6 2  1  2   3 8
  7   4    2 1   1   1   
N3    
 3   7 3   4  3   2 8
  8   3    1 1   1   1   
N4    
 4  8  4  3  4   1 8
  1   6    8 1   1   1   
N5    
5  1 5  6  5   8 8
  2   5    7 1   1   1   
N6    
6  2 6  5  6   7 8
  3   8    6 1   1   1   
N7    
7  3 7  8  7   6 8

40
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
  4   7    5 1   1   1   
N8    
8   4 8  7  8   5 8
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

4 8 12

Shape functions for serendipity elements

1. Four noded rectangular element

Shape function for the node 1:

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

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

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

Shape functions for the corner nodes:

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

 N1  
1   1   1     
4
Similarly
N2  
1   1   1      N   1   1   1     
3
4 4

N4  
1   1   1     
4
Shape functions for the middle nodes:
N5  C 1   1   1   
N5  1 at   0 &   1
 C 1 / 2
 1    1    1     1     1   
2

 N5  
2 2
Similarly
1     1   2   1    1   
2
1     1   2 
N6  N7  N8 
2 2 2

42
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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 Mapped element in

Natural coordinate system Global coordinate system

43
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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.
x  N1 x1  N 2 x2  ..............  N8 x8
(Geometry)
y  N1 y1  N 2 y2  ..............  N8 y8
u  N1u1  N 2u2  ..............  N8u8
(Displacements)
v  N1v1  N 2v2  ..............  N8v8

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.
x  N1 x1  N 2 x2  N3 x3  N 4 x4
(Geometry)
y  N1 y1  N 2 y2  N3 y3  N 4 y4
u  N1u1  N 2u2  ..............  N8u8
(Displacements)
v  N1v1  N 2v2  ..............  N8v8

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.
x  N1 x1  N 2 x2  .................  N8 x8
(Geometry)
y  N1 y1  N 2 y2  .................  N8 y8
u  N1u1  N 2u2  N 3u3  N 4u4
(Displacements)
v  N1v1  N 2v2  N3v3  N 4v4

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

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
x  N1 x1  N 2 x2  .................  N n xn
y  N1 y1  N 2 y2  .................  N n yn
where [N] = Shape function
 x = Coordinates of any point
 xe = Coordinates of nodal points
45
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
N1 
1   1    ; N  1   1   
2
4 4

N3 
1   1    ; N  1   1   
4
4 4
Values of these shape functions at given coordinates (   0.5,  0.6 )

N1 
1  0.51  0.6   0.05 N2 
1  0.5 1  0.6   0.15
4 4

N3 
1  0.51  0.6   0.6 N4 
1  0.5 1  0.6   0.2
4 4
Given Nodal coordinates
x1  2, y1  1, x2  8, y2  3,
x3  7 , y3  7 , x4  3, y4  5
Therefore, coordinates of point ‘P’ in Cartesian system are P (x, y) where
x  N 1 x1  N 2 x2  N 3 x3  N 4 x4  6.1
y  N 1 y1  N 2 y2  N 3 y3  N 4 y4  5.7

46
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
N1 
1   1    ; N  1   1   
2
4 4

N3 
1   1    ; N  1   1   
4
4 4
Cartesian coordinates of point ‘P’ are
x  N 1 x1  N 2 x2  N 3 x3  N 4 x4
y  N 1 y1  N 2 y2  N 3 y3  N 4 y4
Given P ( x, y ) = P (7, 4)

7
1   1    3  1   1    6  1   1    8  1   1    2
4 4 4 4
 9  9    3 (1)

4
1   1    1  1   1    1  1   1    6  1   1    5
4 4 4 4
 3    9   (2)
3  9
From equation (1),  (3)
1
Put into the equation (2)
47
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 3  9   3  9 
9  9     3 
 1   1 
26 2  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

x  N 1 x1  N 2 x2  N 3 x3  N 4 x4
y  N 1 y1  N 2 y2  N 3 y3  N 4 y4
Displacements in-terms of shape functions is represented as
u  N 1u1  N 2u2  N 3u3  N 4u4
v  N 1v1  N 2v2  N 3v3  N 4v4
Shape functions for the four noded rectangular element in natural coordinate
system are

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

N1 
1   1    ; N2 
1   1   
4 4

N3 
1   1    ; N4 
1   1   
4 4
Strain displacement relationship
 u1 
 
 u   N N 2 N 3 N 4  u2 
   1 0 0 0 0  u 
 x   x x x x  3
 v   N1 N 2 N 3 N 4  u4 
     0  
y   v1 
0 0 0
 y   y y y

 u v   N1 N 2 N 3 N 4 N1 N 2 N 3 N 4   v2 
     
 y x   y y y y x x x x   v3 
v 
 4
    Bxe
where [B] = Strain displacement matrix
N i N i  N i 
 
x  x  x
(i = 1, 2, 3, 4)
N i N i  N i 
 
y  y  y
In matrix form
 N i       N i 
 x   x x    
   
 i   
 N     N i 
 y   y y    
where [J] = Jacobian Matrix

     x x 
 x x     
J 
1
 and J 
     y y 
 y y     


Elements of Jacobian matrix are

49
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
x N1 N 2 N 3 N 4 4
N i
 x1  x2  x3  x4   xi
     i 1 
x N1 N 2 N 3 N 4 4
N i
 x1  x2  x3  x4   xi
     i 1 
y N1 N N N 4
N i
 y1  2 y2  3 y3  4 y4   yi
     i 1 
y N1 N N N 4
N i
 y1  2 y2  3 y3  4 y4  yi
     i 1 
Finally Jacobian matrix is obtained as
  x  x   4 N i 4
N i 
       x i   i 
x
J    i 1 i 1

 y y   4
N i 4
N i 
         y i  y
 i 
   i 1 i 1

6.6 Strain displacement matrix for four noded quadrilateral element

Strain displacement matrix is obtained from strain displacement relationship
explained in Jacobian matrix.
 N1 N 2 N3 N 4 
 0 0 0 0 
 x x x x 
 N1 N 2 N 3 N 4 
 B   0
y 
0 0 0
y y y
 
 N1 N 2 N3 N 4 N1 N 2 N 3 N 4 
 y y y y x x x x 

Elements of Strain displacement matrix are as follows
N1 N1  N1  N1 N1  N1 
   
x  x  x y  y  y
N 2 N 2  N 2  N 2 N 2  N 2 
   
x  x  x y  y  y
N3 N3  N3  N3 N3  N3 
   
x  x  x y  y  y
N 4 N 4  N 4  N 4 N 4  N 4 
   
x  x  x y  y  y

50
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
N1 
1   1    N2 
1   1   
4 4

N3 
1   1    N4 
1   1   
4 4
Elements of Jacobian matrix are as follows
x N1 N N N
 x1  2 x2  3 x3  4 x4
    
1    1    1    1     6  2 
  2   4   5  1 
44 4 4 4
x N1 N N N
 x1  2 x2  3 x3  4 x4
    
1    1    1    1    
  2   4   5  1 
4 4 4 4 2
y N1 N N N
 y1  2 y2  3 y3  4 y4
    
1    1    1    1    1
  4   5   2  1 
4 4 4 4 2
51
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
y N1 N N N
 y1  2 y2  3 y3  4 y4
    
1    1    1    1    3
  4   5   2  1 
4 4 4 4 2

Therefore, Jacobian matrix

 x x    6  2   
     
J  4 2 
 y y   1 3 
     2 
 2

Example: Obtain strain displacement matrix for the quadrilateral element as

shown in figure using isoparametric formulation.

Solution: Elements of Jacobian matrix from previous example

 x x    6  2   
     
J  4 2 
 y y   1 3 
     2 
  2
Using Jacobian matrix
 4  2   2
   2 
x  6  2  x  y y 3
Elements of Strain-displacement matrix are as follows
Shape functions for the four noded rectangular element in natural coordinate
system are
N1 
1   1    N2 
1   1   
4 4

N3 
1   1    N4 
1   1   
4 4
52
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Elements of strain-displacement matrix are
 N1 N 2 N3 N 4 
 0 0 0 0 
 x x x x 
 N1 N 2 N 3 N 4 
 B   0
y 
0 0 0
y y y
 
 N1 N 2 N3 N 4 N1 N 2 N 3 N 4 
 y y y y x x x x 


N1 N1  N1   1    4 1     2 

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

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

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

54
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
u  1   2 x
In matrix form
 
u  1 x  1 
 2 
u   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 x1 and x2.
55
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
u1  1 x1  1 
   
u2  1 x2   2 
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe (3)

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

1 x2 
Inverse is obtained by using method of adjoin
1  x  x1 
 N    P A  1 x  2
1

2  1 1 
 N1  1  x2  x 
   
 N 2  2  x  x1 
x x x  x1
N1  2 and N 2 
L L
IV) Strain-Displacement relationship
Strain-displacement relations from linear theory of elasticity are used. Using Eq.
(4) one can write
u  N1u1  N2u2
du dN1 dN
x   u1  2 u2
dx dx dx
dN dN 2   u1 
 x    1  
 dx dx  u2 
 x    Bxe (4)
where [B] = Strain-displacement matrix
Elements of strain displacement matrix are
 B     L L   dN1
1 1 1 dN 2 1 
  ,  
   dx L dx L 

V) Stress-Strain relationship
56
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Using one dimensional Hooke’s law from linear theory of elasticity
    D 
From Eq. (4)
    D Bxe
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
Q xe    
T T

Q xe    D  B  xe  B   xe

T T T

Q   K  xe
L
where  K     D B  B 
T

0
VII) Elements of Stiffness matrix
L
dN i dN j
K ij  AE  dx
0
dx dx
 1  1 
L L
dN1 dN1 AE
K11  AE  dx  AE       dx 
0 
0
dx dx L  L  L
 1  1 
L L
dN1 dN 2 AE
K12  K 21  AE  dx  AE      dx  
0 
0
dx dx L  L  L
 1  1 
L L
dN 2 dN 2 AE
K 22  AE  dx  AE     dx 
0 
0
dx dx L  L  L
Therefore, stiffness matrix of two noded bar element is
AE  1 1
K
L  1 1 

2. Stiffness matrix for two noded beam element

57
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
w  1   2 x  3 x 2   4 x3 (From Pascal Triangle)
dw
   2  2 3 x  3 4 x 2
dx
In matrix form
1 
 
 w 1 x x 2 x 3   2 
  2 
  0 1 2 x 3 x   3 

 4 

    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.
 w1  1 0 0 0  1 
   0 1 0 0  
 1    2 
   
 w2  1 L L
2
L3   3 
  2
 2 
  0 1 2 L 3L    4 

xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
w   P A xe
1

w   N xe (3)

where [N] = Shape functions

III) Shape functions

58
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 1 0 0 0 
 0 1 0 0 
 N    P  A  1 x x x   3  
1 2

3 / L 2 / L 3 / L
2 2
1 / L 
 
 2/ L
3
1 / L 2 / L 1 / L2 
2 3

Inverse of [A] is obtained using elementary operations of matrix algebra.

3x 2 2 x3 2 x 2 x3
N1  1  2  3 , N2  x   2,
L L L L
3x 2 2 x3 x 2 x3
N3  2  3 , N4    2
L L L L
IV) Strain-Displacement relationship

d 2w  d 2 N1 d 2 N2 d 2 N3 d 2 N4 
   2   2 x
dx  dx dx 2 dx 2 dx 2  e
    Bxe (4)
where [B] = Strain-displacement matrix
Elements of strain displacement matrix are
 6 12 x 2 6x 
 B   2  3    2    2  3    2 
4 6x 6 12 x
 L L  L L   L L   L L 
V) Stress-Strain relationship
    D 
From Eq. (4)
    D Bxe
where [D] = Elasticity matrix
VI) Principle of virtual work
Q xe    
T T

Q xe    D  B  xe  B   xe

T T T

Q   K  xe
L
where  K     D B  B 
T

0
VII) Elements of Stiffness matrix

59
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
L 2
d 2 Ni d N j
Kij  EI  dx
0
dx 2 dx 2
Let i = 1 and j = 1
 6 12 x  6 12 x 
L L
d 2 N1 d 2 N1
K11  EI  dx  EI    2  3   2  3  dx
0 
0
2
dx dx 2
L L  L L 
12 EI
K11  3
L
Similarly other elements can be determined.
12 / L3 6 / L2 12 / L3 6 / L2  w1
 
6 / L2 4/ L 6 / L2 2 / L  1
 K   EI 

12 / L3 6 / L2 12 / L3 6 / L2  w2
 
6/ L
2
2/ L 6 / L2 4 / L  2

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
u  1   2 x   3 y
v   4  5 x  6 y
In matrix form
1   4 
u  1 x y   2  and v  1 x y   5 
 
   
 3  6
u   P (1)
where, [P] = Parametric matrix
II) Displacement function in-terms of nodal displacements
60
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
 u1  1 x1 y1  1 
    
u2   1 x2 y2   2 
u  1 x y   
 3  3 3 3

xe   A  (2)

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

u   N xe (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
1 x1 y1 
 N    P  A  1 x y  1 x2 y2 
1

1 x3 y 
 a1 b1 c1 
 N   1 x y  a2 b2 c2 
1 
2A
 a3 b3 c3 
a  a2 x  a 3 y b  b2 x  b3 y c  c2 x  c3 y
N1  1 N2  1 N3  1
2A 2A 2A
where
a1  x2 y3  x3 y2 a2  y2  y3 a3  x3  x2
b1  x3 y1  x1 y3 b2  y3  y1 b3  x1  x3
c1  x1 y2  x2 y1 c2  y1  y2 c3  x2  x1
Substituting shape functions in the equation (3), one can write
u  N1u1  N 2u2  N3u3
v  N1v1  N 2v2  N3v3
IV) Strain-Displacement relationship
61
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
From the linear theory of elasticity, strain-displacement relations for 2D elasticity
problem are
 u   N1 N 2 N 3   u1 
   0 0 0  u 
x   x    x x x  2
   v   N1 N 2 N 3  u3 

 y  
   
y   v1 
0 0 0
   y   y y

 xy   u v   N N N N N N   v 
     
1 2 3 1 2 3 2

 y x   y y y x x x   v3 
    Bxe (4)
Elements of strain-displacement matrix are
 a2 b2 c2 
 2A 2A 2A 0 0 0 
 
 B    0 0 0 3 
a 3 b3 c
2A 2A 2A
a b3 c3 a2 b2 c2 
 3 
 2 A 2 A 2 A 2 A 2 A 2 A 
V) Stress-Strain relationship
Using stress-strain relationship of for 2D elasticity problem
 x  1  0 x 
  E    
 y   
1 0
 y 
  1  
2

 xy   0 0 1   / 2   xy 
    D 
From Eq. (4)
    D Bxe
where [D] = Elasticity matrix
VI) Principle of virtual work
Q xe    
T T

Q xe    D  B  xe  B  xe

T T T
dA
dA

Q   K  xe
where  K     D  B  B  dA
T

dA

62
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
Example : Derive strain displacement matrix [B], elasticity matrix [D] for the
three noded triangular element as shown in figure.

Solution: I) Displacement function

u  1   2 x   3 y
v   4  5 x  6 y
In matrix form
1   4 
u  1 x y   2  and u  1 x y   5 
 
   
 3  6
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
 u1  1 1 1 1 
    
u2   1 4 3  2 
u  1 2 5  
 3   3
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe
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)

63
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
1
1 1 1
 N    P  A  1 x y  1 4 3
1

1 2 5
14 3 1
 N   1 x y  2 4 2
1
10
 2 1 3 
14  2 x  2 y 3  4 x  y 1  2 x  3 y
N1  N2  N3 
10 10 10
one can write
u  N1u1  N 2u2  N3u3
v  N1v1  N 2v2  N3v3

IV) Strain-Displacement relationship

From the linear theory of elasticity, strain-displacement relations for 2D elasticity
problem are
 u   N1 N 2 N 3   u1 
   0 0 0  u 
x   x    x x x  2
   v   N1 N 2 N 3  u3 

 y  
   
y   v1 
0 0 0
   y   y y

 xy   u v   N N N N N N   v 
     
1 2 3 1 2 3 2

 y x   y y y x x x   v3 
    Bxe
Elements of strain-displacement matrix are
 1 2 1 
5 0 0
0
5 5
 
1 1 3 
 B    0 0 0
5 10 10 
 1 1 3 1 2 1 

 5 10 10 5 5 5 

64
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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
u  1   2 x   3 y   4 xy
v   5   6 x   7 y  8 xy
In matrix form
1   5 
   
u  1 x y xy    and v  1 x y xy   6 
 2
 3   7 

 4 
 
 8 

u   P  and 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), (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

65
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 u1  1 0 0 0  1 
u  1 0   2 
 2  a 0

   
u3  1 a b ab   3 
  
u4  1 0 b 0   4 
xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

u   N xe (3)

where [N] = Shape functions
III) Shape functions
1
1 0 0 0 
1 a 0 0 
 N    P  A  1 x y xy   
1

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

1 / b 0 0 1/ b 
 
1 / ab 1 / ab 1 / ab 1 / ab 
x xy y x xy
N1  1   ; N2   ;
a ab b a ab
xy y xy
N3  ; N2  
ab b ab
Substituting shape functions in the equation (3), one can write
u  N1u1  N 2u2  N3u3  N 4u4
v  N1v1  N 2v2  N3v3  N 4v4
IV) Strain-Displacement relationship

66
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 u1 
 
 u   N1 N 2 N 3 N 4  u2 
   0 0 0 0  u 
x   x   x x x x  3
   v   N1 N 2 N 3 N 4  u4 
 y     0  
y   v1 
0 0 0
   y   y y y

 xy   u v   N N 2 N 3 N 4 N1 N 2 N 3 N 4   v2 
     
1

 y x   y y y y x x x x   v3 
v 
 4
    Bxe (4)
Elements of strain-displacement matrix are
 1 y 1 y y y 
  
 a ab a ab ab  0 0 0 0 
ab
 
 B    0
x 1 x x x 1
0 0 0    
ab b ab ab ab b 
 x 1 x x x 1 1 y 1 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
 x  1  0 x 
  E    
 y   
1 0
 y 
  1  
2

 xy   0 0 1   / 2   xy 
    D 
From Eq. (4)
    D Bxe
where [D] = Elasticity matrix

VI) Principle of virtual work

Q xe    
T T

Q xe    D  B  xe  B  xe

T T T
dA
dA

Q   K  xe

67
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
where  K     D B B
T
dA
dA

Example : Derive strain displacement matrix [B], elasticity matrix [D] for the four
noded rectangular element as shown in figure.

Solution: I) Displacement function

u  1   2 x   3 y   4 xy
v   5   6 x   7 y  8 xy
In matrix form
1   5 
   
   
u  1 x y xy   2  and v  1 x y xy   6 
 3   7 

 4 
 
 8 

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).
 u1  1 0 0 0  1 
u  1 6 0 0   
 2   2
   
u3  1 6 4 24   3 
  
u4  1 0 4 0    4 

xe   A  (2)
where, [A] = Connectivity matrix
Obtained   from Eq. (2) and put into the Eq. (1), we get
u   P A xe
1

68
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
u   N xe
where [N] = Shape functions

III) Shape functions

1
1 0 0 0 
1 6 0 0 
 N    P  A  1 x y xy   
1

1 6 4 24 
 
1 0 4 0 
Inverse is obtained using elementary operation of matrix algebra.

 1 0 0 0 
 1 / 6 1 / 6 0 0 
 N    P  A  1 x y xy   
1

 1 / 4 0 0 1/ 4 
 
1 / 24 1 / 24 1 / 24 1 / 24 
x xy y x xy
N1  1    ; N2   ;
6 24 4 6 24
xy y xy
N3  ; N2  
24 4 24
One can write
u  N1u1  N 2u2  N3u3  N 4u4
v  N1v1  N 2v2  N3v3  N 4v4

IV) Strain-Displacement relationship

 u1 
 
 u   N1 N 2 N 3 N 4  u2 
   0 0 0 0  u 
  x   x   x x x x  3
   v   N1 N 2 N 3 N 4  u4 
 y     0  
y   v1 
0 0 0
   y   y y y

 xy   u v   N N 2 N 3 N 4 N1 N 2 N 3 N 4   v2 
     
1

 y x   y y y y x x x x   v3 
v 
 4
    Bxe
Elements of strain-displacement matrix are
69
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
 1 y   1 y  y y 
  6  24   6  24   0 0 0 0
    24 24

  x 1 x x  x 1 
 B   0 0 0 0        
  24 4  24 24  24 4  
  x 1 x x  x 1  1 y  1 y  y y 
              
  24 4  24 24  24 4   6 24   6 24  24 24 

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 coordinate   1 and node 2
has coordinate   1.
x1, x2 are the Cartesian coordinates of nodes 1 and 2.

I) Shape functions in natural coordinates:

1  1 
N1  and N 2  (1)
2 2
II) Geometry:
x  N1x1  N2 x2 (2)
III) Displacements:
u  N1u1  N2u2 (3)
IV) Strain-displacement relationship
du dN1 dN 2
x   u1  u2
dx dx dx
In matrix form
dN dN 2   u1 
    1  
 dx dx  u2 
    Bxe (4)
70
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
where [B] = Strain-displacement matrix
Elements of [B] matrix
dN1 dN1 d dN 2 dN 2 d 
  (5)
dx d dx dx d  dx
From Eq. (1)
dN1 1 dN 2 1
 
d 2 d 2
From Eq. (2)
dx dN1 dN 2
 x1  x2
d d d
dx x x x x L
 1  2  2 1 
d 2 2 2 2
d 2
 
dx L
From Eq. (5), elements of strain-displacement matrix are
dN1 1 2 1 dN 2 1 2 1
     
dx 2 L L dx 2 L L
 B     L L 
1 1
 
V) Stress-Strain relationship
    D 
From Eq. (4)
    D Bxe
where [D] = Elasticity matrix
VI) Principle of virtual work
Q xe    
T T

Q xe    D  B  xe  B   xe

T T T

Q   K  xe
L

where  K     D B  B 
T
dx
0

Formulation of stiffness matrix for 2D four noded quadrilateral isoparametric

element

71
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
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:

N1 
1   1    N2 
1   1   
4 4
(1)
N3 
1   1    N4 
1   1   
4 4
II) Geometry:
x  N1 x1  N 2 x2  N3 x3  N 4 x4
(2)
y  N1 y1  N 2 y2  N3 y3  N 4 y4
III) Displacements:
u  N1u1  N 2u2  N3u3  N 4u4
(3)
v  N1v1  N 2v2   N3v3  N 4v4
IV) Strain-displacement relationship
 u1 
 
 u   N N 2 N 3 N 4  u2 
   1 0 0 0 0  u 
 x   x x x x  3
 v   N1 N 2 N 3 N 4  u4 
     0  
y   v1 
0 0 0
 y   y y y

 u v   N1 N 2 N 3 N 4 N1 N 2 N 3 N 4   v2 
     
 y x   y y y y x x x x   v3 
v 
 4
In matrix form
    Bxe (4)
where [B] = Strain-displacement matrix
72
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
N i N i  N i 
 
x  x  x
(i = 1, 2, 3, 4)
N i N i  N i 
 
y  y  y
In matrix form
 N i       N i 
 x   x x    
   
 i   
 N     N i 
 y   y y    
where [J] = Jacobian Matrix

     x x 
 x x     
J 
1
 and J 
       y y 
 y y     
 
Elements of Jacobian matrix are
x N1 N 2 N 3 N 4 4
N i
 x1  x2  x3  x4   xi
     i 1 
x N1 N N N 4
N i
 x1  2 x2  3 x3  4 x4   xi
     i 1 
y N1 N N N 4
N i
 y1  2 y2  3 y3  4 y4   yi
     i 1 
y N1 N 2 N 3 N 4 4
N i
 y1  y2  y3  y4   yi
     i 1 
Finally Jacobian matrix is obtained as
  x  x   4 N i 4
N i 
       x i   x i 
J    i 1 i 1

 y y   4
N i 4
N i 
         y i   y i 
   i 1 i 1 
Elements of Strain displacement matrix [B] are as follows
N1 N1  N1  N1 N1  N1 
   
x  x  x y  y  y

73
Finite Element Method
Lecture Notes
Dr. Atteshamuddin S. Sayyad, SRES’s Sanjivani College of Engineering, Kopargaon-423603
N 2 N 2  N 2  N 2 N 2  N 2 
   
x  x  x y  y  y
N3 N3  N3  N3 N3  N3 
   
x  x  x y  y  y
N 4 N 4  N 4  N 4 N 4  N 4 
   
x  x  x y  y  y
V) Stress-Strain relationship
    D 
From Eq. (4)
    D Bxe
where [D] = Elasticity matrix
VI) Principle of virtual work
Q xe    
T T

Q xe    D  B  xe  B  xe

T T T
dA
dA

Q   K  xe
where  K     D  B  B  dA
T

dA

74