Sa-Iii 6TH Sem Sukrit

GATEWAY EDUCATION, SONEPAT
ENOTES, Subject: Structural analysis—III, Course code –CE-308B Course: BTech ,Branch: Civil Engineering ,
Sem-6th ,
(Prepared By: Mr. SUKRIT JAIN, Assistant Professor, CE)
UNIT- 1
Matrix Method
INTRODUCTION TO MATRIX METHODS
MATRIX METHODS OF ANALYSIS :
Broadly the methods of analysis are categorised in two ways
1. Force Methods : Methods in which forces are made unknowns i.e Method
of consistent deformation and strain energy method. In both these
methods solution of number of simultaneous equations is involved.
2. Displacement Methods in which displacements are made unknowns i.e
slope deflection method, Moment distribution method and Kani’s
Method (In disguise). In slope deflection method also, the solution of
number of simultaneous equations is involved.
In both of the above methods, for the solution of simultaneous equations
matrix approach can be employed & such Method is called Matrix method
of analysis.
FORCE METHOD :
Method of consistent deformation is the base and forces are made unknown
∆b = Upward Deflection of point B on primary structure due to all causes

∆bo = Upward Deflection of point B on primary structure due to applied
load(Redundant removed i.e condition Xb = 0)
∆bb = Upward Deflection of point B on primary structure due to Xb
( i.e Redundant )
δbb = Upward Deflection of point B on primary structure due to Xb = 1
∆bb = δbb . Xb
Sem-6th ,
 ∆b = ∆bo + ∆bb Substituting for ∆bb –

∆b = ∆bo + δbb . Xb Called Super position equation
Using the compatibility condition that the net displacement at B = 0 i.e
∆b = 0 we get Xb = – ∆bo / δbb

To Conclude we can say, [ ∆ ] = [∆L ] +[ ∆ R]
DISPLACEMENT METHOD:
This method is based on slope deflection method and displacements are made
unkowns
which are computed by matrix approch instead of solving simultaneous equations
and finaly unknown forces are calculated using slpoe deflection equations.
Mab = Mab + 2EI / L ( 2 a + b + 3 / L)
Mab = Final Moment and may be considered as net force P at the joint
Mab = Fixed end moment i.e Force required for the condition of zero
displacements & is called locking force. ( i.e. P’)
The second term may be considred as the force required to produce the required
displacements at the joints. (i.e Pd )Therefore the above equation may be
written as [P] = [P’] + [Pd]
Thus, there are Two Methods in matrix methods
MATRIX METHODS
FORCE METHOD DISPLACEMENT METHOD
The force method is also called by the names 1) Flexibility Method 2)Static
Method 3)Compatibility.
Similarly the displacement method is also called by the names 1)Stiffness Method
2) Kinematic method 3) Equilibrium Method.
In both force method & displacement method there are two different approaches 1)
System Approach 2) Element Approach.
To study matrix methods there are some pre-requisites :

i) Matrix Algebra - Addition, subtraction ,Multiplication & inversion of
matrices (Adjoint Method )
ii) Methods of finding out Displacements i.e. slope & deflection at any
point in a structure, such as a)Unit load method or Strain energy method
b) Moment area method etc.
According to unit load method the displacement at any point ‘j’ is given by
∆j = ∫0 Mmjds / EI
Where
Sem-6th ,
M – B M due to applied loads & mj – B M due to unit load at j

Sem-6th ,
δij = ∫0 mi.mj.ds / EI
Sem-6th ,
When unit load is applied at i and is called flexibility coefficient.

The values of ∆j and δij can be directly read from the table depending upon the
combinations of B M diagrams & these tables are called Diagram Multipliers.
iii) Study of Indeterminacies – Static indeterminacy & kinematic indeterminacy
Sem-6th ,
INDETERMINATE STRUCTURES
1. Statically Indeterminate Structure
2. Kinematically Indeterminate
Structure INDETERMINATE
STRUCTURES
Statically Indeterminate Structure : Any structure whose reaction components
or internal stresses cannot be determined by using equations of static equilibrium
alone, (i.e.Fx = 0, Fy = 0, Mz = 0) is a statically Indeterminate Structure.
The additional equations to solve statically indeterminate structure come from
the conditions of compatibility or consistent displacements.
Roller Support : No. of reactions,r = 1
Hinged Support : No. of reactions, r = 2
Fixed Support : No. of reactions, r = 3

1. Pin Jointed Structures i.e. Trusses
Internal static indeterminacy : (Dsi) No. of members required for stability is
given by – 3 joints – 3 members – every additional joint requires two additional
members.
 m’ = 2 (j – 3)+3 j = No. of joints
 m’ = 2j – 3 Stable and statically determinate
Dsi = m – m’ Where m = No. of members in a
structure Dsi = m – (2j – 3)
External static indeterminacy (Dse)
r = No. of reaction components
Equations of static equilibrium = 3 (i.e.Fx = 0, Fy = 0, Mz = 0)
 Dse = r – 3
Total static indeterminacy Ds = Dsi + Dse  Ds = m –(2j – 3)+(r – 3)
 Ds = (m+r) – 2j
Rigid Jointed Structures : No. of reaction components over and above the no. of
equations of static equilibrium is called a degree of static indeterminacy.
Sem-6th ,
Ds = r-3
Equations of static equilibrium = 3
(i.e.Fx = 0, Fy = 0, Mz = 0)
Example 1
No. of reaction components r = 5 (as shown)

 Ds = r – 3 = 5 – 3 = 2 Ds = 2
Introduce cut in the member BC as shown. At the cut the internal stresses are
introduced i.e. shear force and bending moment as shown.
Left part : No. of unknowns = 5 Equations of equilibrium = 3
 Ds = 5 – 3 = 2
Right Part : No. of unknowns = 4 Equations of equilibrium = 2
 Ds = 4 – 2 = 2
 Ds = Static Indeterminacy = 2
Example 2
Fig. (A) Fig. (B) Fig. (C)
Sem-6th ,
Another Approach : For Every member in a rigid jointed structure there will be 3
unknowns i.e. shear force, bending moment, axial force.
Let r be the no. of reaction components and m be the no. of members
Total no. of unknowns = 3m + r
At every joint three equations of static equilibrium are available
 no. of static equations of equilibrium = 3j (where j is no. of joints)
 Ds = (3m + r) – 3j
 In the example r = 6, m = 6, j = 6
Ds = (3 x 6 + 6) – (3 x 6 ) = 6
Sem-6th ,
Kinematic Indeterminacy :
A structure is said to be kinematically indeterminate if the displacement
components of its joints cannot be determined by compatibility conditions alone.
In order to evaluate displacement components at the joints of these structures, it is
necessary to consider the equations of static equilibrium. i.e. no. of unknown joint
displacements over and above the compatibility conditions will give the degree of
kinematic indeterminacy.
Fixed beam : Kinematically determinate :
Simply supported beam Kinematically indeterminate
Any joint – Moves in three directions in a plane structure

Two displacements x in x direction, y in y direction, 
rotation about z axis as shown.
Roller Support :
r = 1, y = 0,  & x exist – DOF = 2 e=1
Hinged Support :
r = 2, x = 0, y = 0,  exists – DOF = 1 e = 2
Fixed Support :
r = 3, x = 0, y = 0,  = 0 DOF = 0 e=3
i.e. reaction components prevent the displacements  no. of restraints = no. of

reaction components.
Degree of kinematic indeterminacy :
Pin jointed structure :Every joint – two displacements components and no
rotation
Sem-6th ,
 Dk = 2j – e
where,
e = no. of equations of compatibility
j = no. of reaction components
Rigid Jointed Structure : Every joint will have three displacement components,
two displacements and one rotation.
Since, axial force is neglected in case of rigid jointed structures, it is assumed
that the members are inextensible & the conditions due to inextensibility of
members will add to the numbers of restraints. i.e to the ‘e’ value.
 Dk = 3j – e
where,
e = no. of equations of compatibility
j= no. of reaction components + constraints due to in
extensibility
Example 1 : Find the static and kinematic indeterminacies
r = 4, m = 2, j = 3
Ds = (3m + r) –
3j
= (3 x 2 + 4) – 3 x 3 = 1
Dk = 3j – e
=3x3–6=3
i.e. rotations at A,B, & C i.e. a, b & c are
the displacements.
(e = reaction components + inextensibility conditions = 4 + 2 = 6)
Example 2 :
Ds = (3m+r) – 3j
m = 3, r = 6, j = 4
 Ds = (3 x 3 + 6) – 3 x 4 = 3
Dk = 3j – e e = no. of reaction
components + conditions of
inextensibility
Sem-6th ,
= 6+3 = 9
Dk = 3 x 4 – 9 = 3 i.e. rotation b, c & sway.
Sem-6th ,
Example 3 :
Ds = (3m + r) – 3j
r = 6, m = 10, j = 9
Ds = (3 x 10 + 6) – 3 x 9 = 9
Conditions of inextensibility
: Joint : B C E F H I
1 1 2 2 2 2 Total = 10
Reaction components r=6
 e = 10 + 6 = 16
 Dk = 3j – e
= 3 x 9 – 16 = 11
Sem-6th ,
FORCE METHOD :
This method is also known as flexibility method
or compatibility method. In this method the degree of static indeterminacy of the
structure is determined and the redundants are identified. A coordinate is assigned
to each redundant. Thus,P1, P2-----------Pn are the redundants at the coordinates
1,2,---------n.If all the redundants are removed , the resulting structure known as
released structure, is statically determinate. This released structure is also known
as basic determinate structure. From the principle of super position the net
displacement at any point in statically indeterminate structure is some of the
displacements in the basic structure due to the applied loads and the redundants.
This is known as the compatibility condition and may be expressed by the
equation;
∆1 = ∆1L + ∆1R Where ∆1 - - - - ∆n = Displ. At Co-ord.at 1,2 - -n

∆2 = ∆2L + ∆2R ∆1L ---- ∆nL = Displ.At Co-ord.at 1,2---------n
| | | Due to aplied loads
| | | ∆1 R ----∆nR = Displ.At Co-ord.at 1,2--------n
∆n
∆n = L + R ∆n Due to Redudants
The above equations may be return as [∆] = [∆L] + [∆R]---------(1)

∆1 = ∆1L + δ11 P1 + δ12 P2 +----------δ1nPn
∆2 = ∆2L + δ21 P1 + δ22 P2 +---------δ2nPn
| | | | |
| | | | | - - - - - (2)
∆n = ∆nL + δn1 P1 + δn2 P2 +---------δnnPn

Sem-6th ,
∆ = [∆ L] + [δ] [P]----------------------(3)
[P]= [δ] {[∆] – [∆ L]}-----------------(4)
-1
If the net displacements at the redundants are zero then

∆1, ∆2-------∆n =0,
Then  [P] = - [δ] -1 [∆ L]--------------(5)
The redundants P1,P2,---------Pn are Thus determined
DISPLACEMENT METHOD :
This method is also known as stiffness or equilibrium. In this method the
degree of kinematic indeterminacy (D.O.F) of the structure is determined
and the coordinate is assigned to each independent displacement
component.
In general, The displacement components at the supports and joints are treated
as independent displacement components. Let 1,2,------n be the coordinates
assigned to these independent displacement components ∆1, ∆2-------∆n.
In the first instance lock all the supports and the joints to obtain the
restrained structure in which no displacement is possible at the coordinates.
Let P’1, P’2 , - - - - P’n be the forces required at the coordinates 1,2,----n in the
restrained structure in which the displacements ∆1, ∆2-------∆n are zero. Next, Let
the supports and joints be unlocked permitting displacements ∆1, ∆2-------∆n at
the coordinates. Let these displacements require forces in P1d, P2d,-------Pnd at
coordinates 1,2,-------n respectively.
If P1, P2- - - - - - Pn are the external forces at the coordinates 1,2,-------n, then the
conditions of equilibrium of the structure may be expressed as:
P1 = P’1 +P1d
P2 = P’2 +P2d-------------------------------------(1)
| | |
| | |
Pn = P’n + Pnd
or [P ] = [P’] + [ Pd]-------------------(2)
P1 = P’1 +K11 ∆1+ K12, ∆2+ K13∆3 +-----------K1n ∆n

P2 = P’2 +K21 ∆1+ K22 ∆2+ K23∆3 +------------K2n ∆n
| | | | | |
| | | | | | - - - - - - - - - (3)
Pn = P’n +Kn1 ∆1+ Kn2, ∆2+ Kn3∆3 + - - - - - - - - -Knn ∆n 
i.e [P ] = [P’ ] + [ K ] [ ∆ ]----------------------(4)
 [ ∆ ]= [ K ] –1 {[P] – [P’]}----------------------------------------(5)
Where P = External forces
Sem-6th ,
P’ = Locking forces
Pd=Forces due to displacements
If the external forces act only at the coordinates the terms P’1, P’2 ,--------P’n
vanish. i. e the Locking forces are zero,then
[ ∆ ]= [ K ] –1 [ P ]-------------------------(6)
On the other hand if there are no external forces at the coordinates then [P]=0
then [ ∆ ]= – [ K ] –1 [ P’ ]----------------------(7)
Thus the displacements can be found out. Knowing the displacements the forces
are computed using slope deflection equations:
Mab= Mab+ 2EI / L (2a+ b+3 / L)
Mba=Mba+ 2EI / L (a+ 2b+3 / L)
Where Mab& Mba are the fixed end moments for the member AB due to external
loading
FLEXIBILITY AND STIFFNESS MATRICES : SINGLE CO-ORD.
D= δ x P
D=PL3/3EI = δ x P
P=K x D
 δ=L3 /3EI
P=K x PL3 / 3EI
δ=Flexibility Coeff.
K=3EI / L3
K=Stiffness Coeff.
M=K x D =K x ML/EI
D= ML / EI
 K=EI / L
D= δ x M=ML/EI
K=Stiffness Coeff.
 δ=L / EI
δ X K= 1
δ=Flexibility Coeff.
TWO CO-ORDINATE SYSTEM

D1= δ11P1 + δ12 X P2 & D2= δ21P1+ δ22P2
 D1 δ11 δ12 P1
D2 δ21 δ22 P2
L3 / 3EI L2 / 2EI
[δ ]=
L2 / 2EI L / EI
Unit Force At Co-ord.(1)
Sem-6th ,
δ11=L3 / 3EI δ21=L2 / 2EI
Unit Force At Co-ord.(2) δ12=δ21 =L2 / 2EI δ22=L / EI

Sem-6th ,
Develop the Flexibility and stiffness matrices for frame ABCD with reference to
Coordinates shown
The Flexibility matrix can be developed by

applying unit force successively at
coordinates (1),(2) &(3) and evaluating the
displacements at all the coordinates
δij = ∫ mi mj / EI x ds δij =displacement at I due to unit load at j

Sem-6th ,
INVERSING THE FLEXIBILITY MATRIX [ δ ]

THE STIFENESS MATRIX [ K ] CAN BE OBTAINED
Sem-6th ,
Sem-6th ,
UNIT-2
U
Sem-6th ,
Sem-6th ,
77777
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
Sem-6th ,
ENOTES, Subject: Structural analysis—III, Course code –CE-308B Course: BTech ,Branch:
Civil Engineering , Sem-6th ,
Introduction 1
th
Civil Engineering , Sem-6 ,
UNIT 3
Direct stiffness method and the global stiffness matrix
Although there are several finite element methods, we analyse the Direct Stiffness
Method here, since it is a good starting point for understanding the finite element
formulation. We consider first the simplest possible element – a 1-dimensional
elastic spring which can accommodate only tensile and compressive forces. For the
spring system shown, we accept the following conditions:
 Condition of Compatibility – connected ends (nodes) of adjacent

springs have the same displacements
 Condition of Static Equilibrium – the resultant force at each node is
zero
 Constitutive Relation – that describes how the material (spring)
responds to the applied loads
Model spring system

The constitutive relation can be obtained from the governing equation
for an elastic bar loaded axially along its length:
ddu(AEΔll0)+k=0(1)
Δll0=ε(2)
ddu(AEε)+k=0(3)
ddu(Aσ)+k=0(4)
dFdu+k=0(5)
dFdu=−k(6)
dF=−kdu(7)
The spring stiffness equation relates the nodal displacements to the
applied forces via the spring (element) stiffness. The minus sign
denotes that the force is a restoring one, but from here on in we use
the scalar version of Eqn.7.
Derivation of the Stiffness Matrix for a Single Spring Element

From inspection, we can see that there are two degrees of freedom in
this model, ui and uj. We can write the force equilibrium equations:
k(e)ui−k(e)uj=F(e)i(8)
−k(e)ui+k(e)uj=F(e)j(9)
In matrix form
[ke−ke−keke]{uiuj}={F(e)iF(e)j}(10)
The order of the matrix is [2×2] because there are 2 degrees of
freedom. Note also that the matrix is symmetrical. The ‘element’
stiffness relation is:
[K(e)]{u(e)}={F(e)}(11)
Where Κ(e) is the element stiffness matrix, u(e) the nodal
displacement vector and F(e) the nodal force vector. (The element
stiffness relation is important because it can be used as a building
block for more complex systems. An example of this is provided
later.)
Derivation of a Global Stiffness Matrix
For a more complex spring system, a ‘global’ stiffness matrix is
required – i.e. one that describes the behaviour of the complete
system, and not just the individual springs.
From inspection, we can see that there are two springs (elements) and
three degrees of freedom in this model, u1, u2 and u3. As with the
single spring model above, we can write the force equilibrium
equations:
k1u1−k1u2=F1(12)
−k1u1+(k1+k2)u2−k2u3=F2(13)
k2u3−k2u2=F3(14)
In matrix form
⎡⎣⎢k1−k10−k1k1+k2−k20−k2k2⎤⎦⎥⎧⎩⎨⎪⎪u1u2u3⎫⎭⎬⎪⎪=⎧⎩⎨
⎪⎪F1F2F3⎫⎭⎬⎪⎪(15)
The ‘global’ stiffness relation is written in Eqn.16, which we

distinguish from the ‘element’ stiffness relation in Eqn.11.
[K]{u}={F}(16)
Note the shared k1 and k2 at k22 because of the compatibility
condition at u2. We return to this important feature later on.
Assembling the Global Stiffness Matrix from the Element Stiffness
Matrices
Although it isn’t apparent for the simple two-spring model above,
generating the global stiffness matrix (directly) for a complex system
of springs is impractical. A more efficient method involves the
assembly of the individual element stiffness matrices. For instance, if
you take the 2-element spring system shown,
split it into its component parts in the following way

and derive the force equilibrium equations
k1u1−k1u2=F1(17)
k1u2−k1u1=k2u2−k2u3=F2(18)
k2u3−k2u2=F3(19)
then the individual element stiffness matrices are:
[k1−k1−k1k1]{u1u2}={F1F2}and[k2−k2−k2k2]{u2u3}={F2F3}(20)
such that the global stiffness matrix is the same as that derived
directly in Eqn.15:
(Note that, to create the global stiffness matrix by assembling the

element stiffness matrices, k22 is given by the sum of the direct
stiffnesses acting on node 2 – which is the compatibility criterion.
Note also that the indirect cells kij are either zero (no load transfer
between nodes i and j), or negative to indicate a reaction force.)
For this simple case the benefits of assembling the element stiffness
matrices (as opposed to deriving the global stiffness matrix directly)
aren’t immediately obvious. We consider therefore the following
(more complex) system which contains 5 springs (elements) and 5
degrees of freedom (problems of practical interest can have tens or
hundreds of thousands of degrees of freedom (and more!)). Since
there are 5 degrees of freedom we know the matrix order is 5×5. We
also know that it’s symmetrical, so it takes the form shown below:
We want to populate the cells to generate the global stiffness matrix.

From our observation of simpler systems, e.g. the two spring system
above, the following rules emerge:
 The term in location ii consists of the sum of the direct stiffnesses of

all the elements meeting at node i
 The term in location ij consists of the sum of the indirect stiffnesses
relating to nodes i and j of all the elements joining node i to j
 Add a negative for reaction terms (–kij)

 Add a zero for node combinations that don’t interact
By following these rules, we can generate the global stiffness matrix:
This type of assembly process is handled automatically by

commercial FEM codes
Drag the springs into position and click 'Build matrix', then apply a
force to node 5. You will then see the force equilibrium equations, the
equivalent spring stiffness and the displacement at node 5.
Solving for (u)

The unknowns (degrees of freedom) in the spring systems presented
are the displacements uij. Our global system of equations takes the
following form:
To find {u} solve

u=F[K]−1(22)
Recall that [k][k]−1=I=IdentitiyMatrix=[1001].
Recall also that, in order for a matrix to have an inverse, its
determinant must be non-zero. If the determinant is zero, the matrix is
said to be singular and no unique solution for Eqn.22 exists. For
instance, consider once more the following spring system:
We know that the global stiffness matrix takes the following form
⎡⎣⎢k1−k10−k1k1+k2−k20−k2k2⎤⎦⎥⎧⎩⎨⎪⎪u1u2u3⎫⎭⎬⎪⎪=⎧⎩⎨
⎪⎪F1F2F3⎫⎭⎬⎪⎪(23)
The determinant of [K] can be found from:
det⎡⎣⎢adgbehcfi⎤⎦⎥=(aei+bfg+cdh)−(ceg+bdi+afh)(24)
Such that:
(k1(k1+k2)k2+0+0)−(0+(−k1−k1k2)+(k1−k2−k2))(25)
det[K]=(k12k2+k1k22)−(k12k2+k1k22)=0(26)
Since the determinant of [K] is zero it is not invertible, but singular.
There are no unique solutions and {u} cannot be found. If this is the
case in your own model, then you are likely to receive an error
message!
1
Instructional Objectives
After reading this chapter the student will be able to
1. Derive member stiffness matrix of a beam element.
2. Assemble member stiffness matrices to obtain the global stiffness matrix for a beam.
3. Write down global load vector for the beam problem.
4. Write the global load-displacement relation for the beam.
Introduction.
In chapter 23, a few problems were solved using stiffness method from
fundamentals. The procedure adopted therein is not suitable for computer
implementation. In fact the load displacement relation for the entire structure was
derived from fundamentals. This procedure runs into trouble when the structure is
large and complex. However this can be much simplified provided we follow the
procedure adopted for trusses. In the case of truss, the stiffness matrix of the entire
truss was obtained by assembling the member stiffness matrices of individual
members.
In a similar way, one could obtain the global stiffness matrix of a continuous beam
from assembling member stiffness matrix of individual beam elements. Towards
this end, we break the given beam into a number of beam elements. The stiffness
matrix of each individual beam element can be written very easily.
For example, consider a continuous beam ABC as shown in Fig. 27.1a. The
D
given continuous beam is divided into three beam elements as shown in Fig.
27.1b. It is noticed that, in this case, nodes are located at the supports. Thus each
span is treated as an individual beam. However sometimes it is required to
consider a node between support points. This is done whenever the cross sectional
area changes suddenly or if it is required to calculate vertical or rotational
displacements at an intermediate point. Such a division is shown in Fig. 27.1c. If
the axial deformations are neglected then each node of the beam will have two
degrees of freedom: a vertical displacement (corresponding to shear) and a
rotation (corresponding to bending moment). In Fig. 27.1b, numbers
enclosed in a circle represents beam numbers. The beam ABCD is divided into
three beam members. Hence, there are four nodes and eight degrees of freedom.
The possible displacement degrees of freedom of the beam are also shown in the
figure. Let us use lower numbers to denote unknown degrees of freedom
(unconstrained degrees of freedom) and higher numbers to denote known
(constrained) degrees of freedom. Such a method of identification is adopted in
this course for the ease of imposing boundary conditions directly on the structure
stiffness matrix. However, one could number sequentially as shown in Fig. 27.1d.
This is preferred while solving the problem on a computer.
In the above figures, single headed arrows are used to indicate translational and
double headed arrows are used to indicate rotational degrees of freedom.
Beam Stiffness Matrix.

Fig. 27.2 shows a prismatic beam of a constant cross section that is fully restrained
at ends in local orthogonal co-ordinate system x' y' z' . The beam ends
are denoted by nodes j and k . The x' axis coincides with the centroidal axis of
the member with the positive sense being defined from j to k . Let L be the length
of the member, A area of cross section of the member and I zz is the moment of
inertia about z' axis.
Two degrees of freedom (one translation and one rotation) are considered at each
end of the member. Hence, there are four possible degrees of freedom for this
member and hence the resulting stiffness matrix is of the order 4  4 . In this
method counterclockwise moments and counterclockwise rotations are taken as
positive. The positive sense of the translation and rotation are also shown in the
figure. Displacements are considered as positive in the direction of the coordinate
axis. The elements of the stiffness matrix indicate the forces exerted on
the member by the restraints at the ends of the member when unit displacements
are imposed at each end of the member. Let us calculate the forces developed in
the above beam member when unit displacement is imposed along each degree of
freedom holding all other displacements to zero. Now impose a unit
displacement along y' axis at j end of the member while holding all other
displacements to zero as shown in Fig. 27.3a. This displacement causes both shear
and moment in the beam. The restraint actions are also shown in the figure. By
definition they are elements of the member stiffness matrix. In particular they form
the first column of element stiffness matrix.
In Fig. 27.3b, the unit rotation in the positive sense is imposed at j end of the beam
while holding all other displacements to zero. The restraint actions are shown in
the figure. The restraint actions at ends are calculated referring to
tables given in lesson …
In Fig. 27.3c, unit displacement along y' axis at end k is imposed and
corresponding restraint actions are calculated. Similarly in Fig. 27.3d, unit
rotation about z' axis at end k is imposed and corresponding stiffness
coefficients are calculated. Hence the member stiffness matrix for the beam
member is
1 2 3 4
⎡ 12 EI z 6 EI z 12 EI z 6 EI z ⎤ 1
⎢ L3 L2 L3 L2 ⎥
⎢ 6 EI 4 EI 2 EIz ⎥ ⎥ 2
z z 6zEI
⎢
k  ⎢ L2
⎢  12 EI z
L
6 EI z
L2 L ⎥
6 EI z ⎥ 3
(27.1)
 12 EI z
⎢ L3 L2 L3 L2 ⎥
⎢ 6 EI L2z 2 EI z 6 EI 4 EI z ⎥⎥ 4
⎢  z
⎣ L L2 L ⎦
The stiffness matrix is symmetrical. The stiffness matrix is partitioned to separate

the actions associated with two ends of the member. For continuous beam
problem, if the supports are unyielding, then only rotational degree of freedom
shown in Fig. 27.4 is possible. In such a case the first and the third rows and
columns will be deleted. The reduced stiffness matrix will be,
⎡ 4EIz 2EIz ⎤
⎢ L L ⎥⎥
k   ⎢ 2EI 4EI (27.2)
⎢  z z⎥
⎣L L⎦
Instead of imposing unit displacement along y' at j end of the member in Fig.
a, apply a displacement u'1 along y' at j end of the member as shown in
Fig. 27.5a, holding all other displacements to zero. Let the restraining forces
developed be denoted by q11 , q21 , q31 and q41 .
The forces are equal to,
q11  k11u'1 ; q21  k21u'1 ; q31  k31u'1 ; q41  k41u'1 (27.3)
Now, give displacements u'1 , u'2 , u'3 and u'4 simultaneously along displacement
degrees of freedom 1,2,3 and 4 respectively. Let the restraining forces developed
at member ends be q1 , q2 , q3 q4 respectively as shown in Fig. 27.5b along
and
respective degrees of freedom. Then by the principle of superposition, the force
displacement relationship can be written as,
⎡q ⎤
⎡ 12EIz 6EIz  12EIz 6EIz ⎤ ⎡u' ⎤
⎢ ⎥
1 ⎢ L2 ⎥
1
⎢ L
3
L2 L3 ⎢ ⎥
⎢ ⎥ ⎢ I6E 4E 6E 2EI ⎥ ⎢
I I
⎢⎢ ⎥⎥ ⎢z L2 z
L z L2 z ⎥ ⎢u' ⎥L ⎥ ⎢ ⎥
⎢q2 ⎥  ⎢ ⎢
2 (27.4)
⎢ ⎥
⎢ ⎥ ⎢ 12EIz  6EIz 12EIz 6EIz ⎥⎥⎢ ⎥
q  u'
⎢ 3⎥ ⎢ 3⎥
⎢ ⎥ ⎢ L2 L3 L ⎥ ⎢ ⎥
2
⎢ ⎥
⎢ ⎥ L 3
2EI 6E 4EI ⎥⎥
I
⎢⎢ ⎥
6EI
⎢⎣ q 4 ⎥⎦ ⎢⎣u' 4 ⎥⎦
2z
z z
⎢
⎣
z
L 2
L L L ⎦
This may also be written in compact form as, q  k  u'

(27.5)
Beam (global) Stiffness Matrix.

The formation of structure (beam) stiffness matrix from its member stiffness
matrices is explained with help of two span continuous beam shown in Fig. 27.6a.
Note that no loading is shown on the beam. The orthogonal co-ordinate system xyz
denotes the global co-ordinate system.
For the case of continuous beam, the x - and x' - axes are collinear and other
axes ( y and y' , z and z' ) are parallel to each other. Hence it is not required to
transform member stiffness matrix from local co-ordinate system to global co
ordinate system as done in the case of trusses. For obtaining the global stiffness
matrix, first assume that all joints are restrained. The node and member numbering
for the possible degrees of freedom are shown in Fig 27.6b. The continuous beam
is divided into two beam members. For this member there are six possible degrees
of freedom. Also in the figure, each beam member with its displacement degrees
of freedom (in local co ordinate system) is also shown. Since the continuous beam
has the same moment of inertia and span, the member stiffness matrix of element 1
and 2 are the same. They are,
Global d 1 2 3 4
.o. f Local 1 2 3 4
d .o. f ⎡k k k '13 k '14 ⎤ 1 1
'11 '1 ⎥
2
⎢
k '  ⎢ k21 k k k '24 ⎥ 2 (27.6a)
' 2
'22 '23
⎢k '31 k k k '34 ⎥ 3 3
'3 '33
2
⎢k ⎥
⎣ k k k '44 ⎦ 4 4
' 41 '42 '43
Global d .o. f 3 4 5 6
Local d 1 2 3 4
.o. f
⎡k k k 214 ⎥⎤ 21 43
kk 224
211
⎢221
k 21 21
2
k22 3
k22
2 3
k 2  ⎢ ⎥ (27.6b)
k 231 k k k 234 ⎥ 3
⎢ 23 23 5
2 3
⎢2 ⎥
⎣ k41 k k k 244 ⎥⎦ 4 6
24 24
2 3
The local and the global degrees of freedom are also indicated on the top and side
of the element stiffness matrix. This will help us to place the elements of the
element stiffness matrix at the appropriate locations of the global stiffness matrix.
The continuous beam has six degrees of freedom and hence the stiffness matrix
is of the order 6  6 . Let K  denotes the continuous beam stiffness matrix of
order 6  6 . From Fig. 27.6b, K  may be written as,

Membe AB (1)
⎡ r 111 0 0 ⎤
⎢ k k121 13 k1 k141 ⎥
⎢⎢ k 1 k221 k 1
k 1 0 0 ⎥
 21 23
k 1
k 2
k 1
24
k 2 ⎥ (27.7)
K ⎢ k311 k321 33 11 34 12 13 k2 k 142 ⎥
⎢ 1 k 43  k21
1 2
k 44  k22
1 2
⎥
⎢ k41 k421 23 k2 2
k 24 ⎥
⎢ 0 0 31 k2 k322 33 k2 k 342 ⎥
⎣⎢ 0 0 k2 k2 k2 k2
⎥⎦
Member BC (2)
The 4  4 upper left hand section receives contribution from member 1 ( AB and
)
4  4 lower right hand section of global stiffness matrix receives contribution from
member 2. The element of the global stiffness matrix corresponding to global
degrees of freedom 3 and 4 [overlapping portion of equation 27.7] receives
element from both members 1 and 2.
Formation of load vector.

Consider a continuous beam ABC as shown in Fig. 27.7.
We have two types of load: member loads and joint loads. Joint loads could be
handled very easily as done in case of trusses. Note that stiffness matrix of each
member was developed for end loading only. Thus it is required to replace the
member loads by equivalent joint loads. The equivalent joint loads must be
evaluated such that the displacements produced by them in the beam should be the
same as the displacements produced by the actual loading on the beam. This is
evaluated by invoking the method of superposition.
The loading on the beam shown in Fig. 27.8(a), is equal to the sum of Fig. 27.8(b)
and Fig. 27.8(c). In Fig. 27.8(c), the joints are restrained against displacements and
fixed end forces are calculated. In Fig. 27.8(c) these fixed end actions are shown in
reverse direction on the actual beam without any load. Since the beam in Fig.
27.8(b) is restrained (fixed) against any displacement, the displacements produced
by the joint loads in Fig. 27.8(c) must be equal to the displacement produced by
the actual beam in Fig. 27.8(a). Thus the loads shown
in Fig. 27.8(c) are the equivalent joint loads .Let, p1 , p2 , p3 , p4 , p5 and p6 be the
equivalent joint loads acting on the continuous beam along displacement degrees
of freedom 1,2,3,4,5 and 6 respectively as shown in Fig. 27.8(b). Thus the global
load vector is,
⎧ Pb ⎫
⎪  ⎪
⎪ L ⎪
⎧ p1 ⎫ ⎪ ⎪
 Pab2
⎪ ⎪ ⎪⎪ L2
⎪
⎪
⎪ p2 ⎪
⎪ ⎪ ⎪ ⎛ Pa wL ⎞ ⎪
 
⎪ ⎪ ⎪⎪ ⎜⎝ L
⎪ p3 ⎪ 2 ⎠⎟ ⎪⎪
⎬ (27.8)

⎨p ⎬ ⎨ ⎛ wL Pba ⎞
2 2
⎪ 4 ⎪ ⎪ ⎜  ⎟⎪
⎪ ⎪ ⎜⎝ L2 ⎟⎠ ⎪
⎪ ⎪ 12 ⎪
⎪ p5 ⎛ wL ⎞
⎪
⎪ ⎪ ⎪
⎪
⎪ ⎜
⎪  2P ⎟
⎪⎩ p 6 ⎪ ⎝2 ⎠
⎪
⎪
⎭
⎪ wL2 ⎪
⎪ 12 ⎪
⎩ ⎭
Solution of equilibrium equations

After establishing the global stiffness matrix and load vector of the beam, the load
displacement relationship for the beam can be written as,
P  K u
(27.9)
where
Pis the global load vector, u is displacement vector and K is the
 
global stiffness matrix. This equation is solved exactly in the similar manner as
discussed in the lesson 24. In the above equation some joint displacements are
known from support conditions. The above equation may be written as
⎧⎪pk ⎫⎪ ⎡k11  k12 ⎤⎧⎪uu ⎫⎪

⎨ ⎬⎢ ⎥⎨ (27.10)
⎬
⎪⎩pu ⎪⎭ ⎢⎣k 21  k 22 ⎥⎦ ⎪⎩u k ⎪⎭
where p k  and u k  denote respectively vector of known forces and known
displacements. And pu , u u  denote respectively vector of unknown forces and
unknown displacements respectively. Now expanding equation (27.10),
{ p k }  k11 {u u }  k12 {u k (27.11a)
} (27.11b)
{ p u }  k 21 {u u }  k 22
{u k }
Since
u k is known, from equation 27.11(a), the unknown joint displacements

can be evaluated. And support reactions are evaluated from equation (27.11b),
after evaluating unknown displacement vector.
Le R1 , R3 and R be the reactions along the constrained degrees of freedom as

5
t
shown in Fig. 27.9a. Since equivalent joint loads are directly applied at the
supports, they also need to be considered while calculating the actual reactions.
Thus,
⎧R1 ⎫ ⎧ p1 ⎫
⎪R ⎪  ⎪p ⎪
K u (27.12)
⎨ 3⎬ ⎨3 21 u
⎪ ⎪ ⎬
⎪
⎪
⎪⎩ R5 ⎪⎭ ⎪⎩ p5 ⎪⎭
The reactions may be calculated as follows. The reactions of the beam shown in
Fig. 27.9a are equal to the sum of reactions shown in Fig. 27.9b, Fig. 27.9c and
Fig. 27.9d.
From the method of superposition,

Pb
R Ku
1 14 4  K16 u6 (27.13a)
L
Pa
R K
3
34  K 36 u6 (27.13b)
L u4
wL
R  2P  K
5
54  K56 u6 (27.13c)
2 u4
or
⎧R ⎫⎪ ⎧ ⎫
⎪ Pb ⎪ ⎡K K 16⎤ ⎧u ⎫
⎪ ⎪
1
⎪ ⎪ L ⎪  14 ⎪ 4⎪ (27.14a)
Pa 3 ⎬ ⎨
⎨ R ⎢K
⎬ ⎢ 34 K 36⎥ ⎥ ⎨ ⎬
L
⎪ ⎪ ⎪ wl ⎪ ⎢ ⎥⎪⎩u 6 ⎪⎭
⎢ ⎥
⎪⎩ R5 ⎪⎭ ⎪  ⎣ K 54 K56 ⎦
2P⎪
⎩2
Equation (27.14a) may be written as,
⎧ K16 ⎤
⎧R1 ⎫ ⎡ K14 ⎧u ⎫
⎫
⎪ Pb L ⎪
⎪
⎪ ⎪
R  Pa L ⎢  K ⎥⎨ ⎬
4 (27.14b)
⎪ K
 3⎬
⎨ ⎬ ⎢ 34 36 ⎥u
⎨
⎪ ⎪ ⎪ ⎪ ⎢K
wl R ⎩ 6⎭
⎩ 5⎭ ⎪  ⎣ 54 K 56⎥⎦
⎩
2
Member end actions q , q , q , q are calculated as follows. For example consider

1 2 3 4
the first element 1.
⎧ Pb
⎧q1 ⎫ ⎫ ⎧0⎫
⎪ L
⎪
⎪ ⎪ ⎪ Pab2 ⎪ ⎪ ⎪
⎪q2 ⎪ ⎪ L 2 ⎪
  K
⎪
⎪ ⎪u2 (27.16)
⎪ ⎬
⎨ ⎬ ⎨ Pa ⎪ element1 ⎨
⎪ q3 ⎪
⎪ ⎬ ⎪0⎪
⎪ ⎪ ⎪ L ⎪ ⎪ ⎪
⎪⎩q 4 ⎪ Pa 2b ⎪⎩u 4 ⎪⎭
⎪ ⎪
⎪
⎭ ⎪
In the next
lesson few
problems are
solved to
illustrate the
method so far
discussed.
ENOTES, Subject: Structural analysis—III, Course code –CE-308B Course: BTech ,Branch: Civil Engineering , Sem-6 th ,
UNIT - 4
Introduction
In engineering problems there are some basic unknowns. If they are found, the
behaviour of the entire structure can be predicted. The basic unknowns or the
Field variables which are encountered in the engineering problems are
displacements in solid mechanics, velocities in fluid mechanics, electric and
magnetic potentials in electrical engineering and temperatures in heat flow
problems.
In a continuum, these unknowns are infinite. The finite element procedure
reduces such unknowns to a finite number by dividing the solution region into
small parts called elements and by expressing the unknown field variables in
terms of assumed approximating functions (Interpolating functions/Shape
functions) within each element. The approximating functions are defined in terms
of field variables of specified points called nodes or nodal points. Thus in the
finite element analysis the unknowns are the field variables of the nodal points.
Once these are found the field variables at any point can be found by using
interpolation functions.
After selecting elements and nodal unknowns next step in finite element
analysis is to assemble element properties for each element. For example, in solid
mechanics, we have to find the force-displacement i.e. stiffness characteristics of
each individual element. Mathematically this relationship is of the form
[k]e {}e  {F}e
where [k]e is element stiffness matrix, { }e is nodal displacement vector of the

element and {F}e is nodal force vector. The element of stiffness matrix kij represent
the force in coordinate direction ‘i’ due to a unit displacement in coordinate
direction ‘j’. Four methods are available for formulating these element properties
viz. direct approach, variational approach, weighted residual approach and energy
balance approach. Any one of these methods can be used for assembling element
properties. In solid mechanics variational approach is commonly employed to
assemble stiffness matrix and nodal force vector (consistant loads).
Element properties are used to assemble global properties/structure properties to get system
equations
[k] { }  {F}. Then the boundary conditions are imposed. The solution
of these simultaneous equations give the nodal unknowns. Using these nodal
values additional calculations are made to get the required values e.g. stresses,
strains, moments, etc. in solid mechanics problems.
Thus the various steps involved in the finite element analysis are:
(i) Select suitable field variables and the elements.
(ii) Discritise the continua.
(iii) Select interpolation functions.
(iv) Find the element properties.
(v) Assemble element properties to get global properties.
(vi) Impose the boundary conditions.
(vii)Solve the system equations to get the nodal unknowns.
(viii)Make the additional calculations to get the required values.
(ix) A BRIEF EVPLANATION OF FEA FOR A STRESS ANALYSIS PROBLEM
The steps involved in finite element analysis are clarified by taking the stress
analysis of a tension strip with fillets (refer Fig.1.1). In this problem stress
concentration is to be studies in the fillet zone. Since the problem is having
symmetry about both x and y axes, only one quarter of the tension strip may be
considered as shown in Fig.1.2. About the symmetric axes, transverse
displacements of all nodes are to be made zero. The various steps involved in the
finite element analysis of this problem are discussed below:
Step 1: Four noded isoparametric element (refer Fig 1.3) is selected for the
analysis (However note that 8 noded isoparametric element is ideal for this
analysis). The four noded isoparametric element can take quadrilateral shape also
as required for elements 12, 15, 18, etc. As there is no bending of strip, only
displacement continuity is to be ensured but not the slope continuity. Hence
displacements of nodes in x and y directions are taken as basic unknowns in the
problem.
Fillet t
P
b2 C
b1
A
B D
Fig. 1.1 †ypical tension flat

Introduction 3
A 1 5 9 13 17 21 24 29 33 37 41 45
1 4 7 10 13 16 19 22 25 28 31 C
2
2 5 8 11 14
17 20 23 26 29 32 P
3
3 6 9 12 15
B 4 18
8 12 16 20 24 21 24 27 30 33 D
32 28
36 40 44 48
Fig. 1.2 Discretisation of quarter of tension flat
6 10 4 3
xP
7
(a) Element no. 5
11 1 2
(b) Typical element
Fig. 1.3
Step 2: The portion to be analysed is to be discretised. Fig. 1.2 shows discretised

portion. For this 33 elements have been used. There are 48 nodes. At each node
unknowns are x and y components of displacements. Hence in this problem total
unknowns (displacements) to be determined are 48 × 2 = 96.
Step 3: The displacement of any point inside the element is approximated by
suitable functions in terms of the nodal displacements of the element. For the
typical element (Fig. 1.3 b), displacements at P are
u  N u N u N u N u N u
i i 1 1 2 2 3 3 4 4
and v   Nivi  N1v1  N2v2  N3v3  N4v4 …(1.2)
The approximating functions Ni are called shape functions or interpolation

functions. Usually they are derived using polynomials. The methods of deriving
these functions for various elements are discussed in this text in latter chapters.
Step 4: Now the stiffness characters and consistant loads are to be found for each
element. There are four nodes and at each node degree of freedom is 2. Hence
degree of freedom in each element is 4 × 2 = 8. The relationship between the nodal
displacements and nodal forces is called element stiffness characteristics. It is of
the form
[k]e {}e  {F}e , as explained earlier.
For the element under consideration, ke is 8 × 8 matrix and  e and Fe are vectors
of 8 values. In solid mechanics element stiffness matrix is assembled using
variational approach i.e. by minimizing potential energy. If the load is acting in the
body of element or on the surface of element, its equivalent at nodal points are to
be found using variational approach, so that right hand side of the above
expression is assembled. This process is called finding consistant loads.
4 Finite Element Analysis
Step 5: The structure is having 48 × 2 = 96 displacement and load vector

components to be determined. Hence global stiffness equation is of the form
[k] {} = {F}
96 × 96 96 × 1 96 × 1
Each element stiffness matrix is to be placed in the global stiffness matrix
appropriately. This process is called assembling global stiffness matrix. In this
problem force vector F is zero at all nodes except at nodes 45, 46, 47 and 48 in x
direction. For the given loading nodal equivalent forces are found and the force
vector F is assembled.
Step 6: In this problem, due to symmetry transverse displacements along AB and
BC are zero. The system equation [k] {}  {F} is modified to see that the solution for
{} comes out with the above values. This modification of system equation is
called imposing the boundary conditions.
Step 7: The above 96 simultaneous equations are solved using the standard
numerical procedures like Gauss- elimination or Choleski’s decomposition
techniques to get the 96 nodal displacements.
Step 8: Now the interest of the analyst is to study the stresses at various points. In
solid mechanics the relationship between the displacements and stresses are well
established. The stresses at various points of interest may be found by using shape
functions and the nodal displacements and then stresses calculated. The stress
concentrations may be studies by comparing the values obtained at various points
in the fillet zone with the values at uniform zone, far away from the fillet (which is
equal to P/b2t).
FINITE ELEMENT METHOD VS CLASSICAL METHODS

1. In classical methods exact equations are formed and exact solutions are
obtained where as in finite element analysis exact equations are formed
but approximate solutions are obtained.
2. Solutions have been obtained for few standard cases by classical methods,
where as solutions can be obtained for all problems by finite element
analysis.
3. Whenever the following complexities are faced, classical method makes
the drastic assumptions’ and looks for the solutions:
(a) Shape
(b) Boundary conditions

(c) Loading
Fig. 1.4 shows such cases in the analysis of slabs (plates).
To get the solution in the above cases, rectangular shapes, same boundary
condition along a side and regular equivalent loads are to be assumed. In
FEM no such assumptions are made. The problem is treated as it is.
4. When material property is not isotropic, solutions for the problems
become very difficult in classical method. Only few simple cases have
been tried successfully by researchers. FEM can handle structures with
anisotropic properties also without any difficulty.
5. If structure consists of more than one material, it is difficult to use
classical method, but finite element can be used without any difficulty.
6. Problems with material and geometric non-linearities can not be handled
by classical methods. There is no difficulty in FEM.
Hence FEM is superior to the classical methods only for the problems
involving a number of complexities which cannot be handled by classical methods
without making drastic assumptions. For all regular problems, the solutions by
classical methods are the best solutions. Infact, to check the validity of the FEM
programs developed, the FEM solutions are compared with the solutions by
classical methods for standard problems.
(a) Irregular shaper (b) Irregular boundary condition
(c) Irregular loading
Fig. 1.4
FEM VS FINITE DIFFERENCE METHOD (FDM)

1. FDM makes pointwise approximation to the governing equations i.e. it
ensures continuity only at the node points. Continuity along the sides of
grid lines are not ensured.
FEM make piecewise approximation i.e. it ensures the continuity at
node points as well as along the sides of the element.
2. FDM do not give the values at any point except at node points. It do not
give any approximating function to evaluate the basic values (deflections,
in case of solid mechanics) using the nodal values.
FEM can give the values at any point. However the values obtained at
points other than nodes are by using suitable interpolation formulae.
3. FDM makes stair type approximation to sloping and curved boundaries as shown in Fig.
1.5.
FEM can consider the sloping boundaries exactly. If curved elements
are used, even the curved boundaries can be handled exactly.
4. FDM needs larger number of nodes to get good results while FEM needs fewer nodes.
5. With FDM fairly complicated problems can be handled where as FEM can
handle all complicated problems.
6 Finite Element Analysis
Fig. 1.5 VDW approximation of shape
A BRIEF HISTORY OF FEM
Engineers, physicists and mathematicians have developed finite element method independently. In
1943 Courant
[1] made an effort to use piecewise continuous functions defined over triangular domain.
After that it took nearly a decade to use this distribution idea. In fifties renewed
interest in this field was shown by Polya [2], Hersh [3] and Weinberger [4].
Argyris and Kelsey [5] introduced the concept of applying energy principles to the
formation of structural analysis problems in 1960. In the same year Clough [6]
introduced the word ‘Finite Element Method’.
In sixties convergence aspect of the finite element method was pursued more
rigorously. One such study by Melesh [7] led to the formulation of the finite
element method based on the principles of minimum potential energy. Soon after
that de Veubeke [8] introduced equilibrium elements based on the principles of
minimum potential energy, Pion [9] introduced the concept of hybrid element
using the duel principle of minimum potential energy and minimum
complementary energy.
In Late 1960’s and 1970’s, considerable progress was made in the field of
finite element analysis. The improvements in the speed and memory capacity of
computers largely contributed to the progress and success of this method. In the
field of solid mechanics from the initial attention focused on the elastic analysis of
plane stress and plane strain problems, the method has been successfully extended
to the cases of the analysis of three dimensional problems, stability and vibration
problems, non-linear analysis. A number of books [10 – 20] have appeared and
made this field interesting.
NEED FOR STUDYING FEM
Now, a number of users friendly packages are available in the market. Hence one
may ask the question ‘What is the need to study FEA?’.
The above argument is not sound. The finite element knowledge makes a good
engineer better while just user without the knowledge of FEA may produce more
dangerous results. To use the FEA packages properly, the user must know the
following points clearly:
1. Which elements are to be used for solving the problem in hand.
2. How to discritise to get good results.
3. How to introduce boundary conditions properly.
4. How the element properties are developed and what are their limitations.
5. How the displays are developed in pre and post processor to understand their limitations.
6. To understand the difficulties involved in the development of FEA
programs and hence the need for checking the commercially available
packages with the results of standard cases.
Unless user has the background of FEA, he may produce worst results and may go with overconfidence.
Hence it is necessary that the users of FEA package should have sound knowledge of FEA.
WARNING TO FEA PACKAGE USERS
When hand calculations are made, the designer always gets the feel of the
structure and get rough idea about the expected results. This aspect cannot be
ignored by any designer, whatever be the reliability of the program, a complex
problem may be simplified with drastic assumptions and FEA results obtained.
Check whether expected trend of the result is obtained. Then avoid drastic
assumptions and get more refined results with FEA package. User must remember
that structural behaviour is not dictated by the computer programs. Hence the
designer should develop feel of the structure and make use of the programs to get
numerical results which are close to structural behaviour.
QUESTIONS
1. Explain the concept of FEM briefly and outline the procedure.

2. Discuss the advantages and disadvantages of FEM over
(i) Classical method
(ii) Finite difference method.
3. Clearly point out the situations in which FEM is preferred over other methods.
4. When there are several FEM packages are available is there need to study this method?
Discuss.
References
1. R. Courant, “Variational Methods for the Solutions of Problems of Equilibrium and Vibrations”,
Bulletin of American Mathematical Society, Vol. 49, 1943.
2. G. Polya, Estimates for Eigen Values, Studies presented to Richard Von Mises, Academic Press,
New York, 1954.
3. J. Hersch, “Equations Differentielles et Functions de cellules”, C.R. Acad. Science, Vol. 240, 1955.
4. H.F. Weinberger, “Upper and Lower Bounds for Eigen Values by Finite Difference Method”,
Pure Applied Mathematics, Vol. 9, 1956.
5. J.H. Argyris and S. Kelsey, “Energy Theorems and Structural Analysis”, Aircraft Engineering,
Vol. 27, 1955.
6. R.W. Clough, “The Finite Element Method in Plane Stress Analysis”, Proceeding of 2nd ASCE
Conference on Electronic Computation, Pittsburg, PA, September, 1960.
7. R.J. Melosh, “Basis for the Derivation for the Direct Stiffness Method”, AIAA Journal, Vol. 1,
1963.
8. B. Fraeijs de Veubeke, “Upper and Lower Bounds in Matrix Structural Analysis”, AGARD ograph
72, B.F. de Veubeke (ed). Pergaman Press, New York, 1964.
9. T.H.H. Pian, “Derivation of Element Stiffness Matrices”, AIAA Journal, Vol. 2, 1964. pp. 556–57.
10. O.C. Zienkiewicz, The Finite Element Method in Engineering Science, McGraw-Hill, London 1971.
11. K.H. Huebner, The Finite Element Methods for Engineers, John Wiley and Sons, 1971.
12. Desai and Abel, Introduction to the Finite Element Method, CBS Publishers & Distributors, 1972.
13. H.C. Martin and G. F. Carey, Introduction to Finite Element Analysis- Theory and Applications,
Tata McGraw-Hill Publishing Company Ltd., New Delhi, 1975.
14. K.L. Bathe and E.L. Wilson, Finite Element Methods, Prentice Hall, 1976.
15. Y.K. Cheuny and M.F. Yeo, A Practical Introduction to Finite Element Analysis, Pitman
Publishers, 1979.
16. R.D. Cook, D.S. Makus and M.F. Plesha, Concept and Applications of Finite Element Analysis,
John Wiley and Sons, 1981.
17. J.N Reddy, An introduction to the Finite Element Method, McGraw-Hill International Edition,
1984.
18. C.S. Krishnamoorthy, Finite Element Analysis, Theory and Programming, Tata McGraw-Hill
Publishing Company Ltd., New Delhi, 1987.
19. T.R. Chandrapatla and A.D. Belegundu, Introduction to Finite Elements in Engineering, Prentice
Hall, 1991.
20. S. Rajasekharan, Finite Element Analysis in Engineering Design, Wheeler Publisher, 1993.
INTRODUCTION
This chapter summarizes the results from theory of elasticity which are useful in
solving the problems in structural and continuum mechanics by the finite element
method.
STRESSES IN A TYPICAL ELEMENT
In theory of elasticity, usually right hand rule is used for selecting the coordinate
system. Fig. 2.1 shows various orientations of right hand rule of the coordinate
systems. Equations derived for any one such orientation hold good for all other
orientations of
z
z y
y x x
(b)
y
x
z
(a) (c)
Fig. 2.1
coordinate system with right hand rule. In this Chapter orientation shown in Fig.
2.1(a) is used for the explanation. Fig. 2.2 shows a typical three dimensional
element of size dx × dy × dz. Face abcd may be called as negative face of x and the
face efgh as the positive face of x since the x value for face abcd is less than that
for the face efgh. Similarly the face aehd is negative face of y and bfgc is positive
face of y. Negative and positive faces of z are dhgc and aefb.
The direct stresses  and shearing stresses  acting on the negative faces are
shown in the Fig. 2.3 with suitable subscript. It may be noted that the first
subscript of shearing stress is the plane and the second subscript is the direction.
Thus the  xy means shearing stress on the plane where x value is constant and y is
the direction.
z z
dy
a e
dx x
dz xy
yx
b f
d y xz
h y
y yz zx
zy
z
c
g
x x
Fig. 2.2 Fig. 2.3
In a stressed body, the values of stresses change from face to face of an

element. Hence on positive face the various stresses acting are shown in Fig. 2.4
with superscript ‘+’.
All these forces are listed in table 2.1
Note the sign convention: A stress is positive when it is on positive face in positive
direction or on negative face in negative direction. In other words the stress is + ve
when it is as shown in Figs 2.3 and 2.4.
z
zy
zx yz
y
xz
y
zx
xy
x
Fig. 2.4
Table 2.1 Stresses on a typical element
Face Stress on –ve Face Stresses on +ve Face

 x
 dx
x x x x
x
 xy
 xy  
  dx
xy xy
x
 xz
 xz  dx
xz xz
x
 y
 dy
y y
y y y
 yx
 yx    dy
yx yx y
 yz
 yz  

dy
yz yz y
 z
 dz
z z
z z z
 zx
 zx  

dz
zx zx z
 zy
 zy  

dz
zy zy z
Y
y
dz
X Fig. 2.5
Let the intensity of body forces acting on the element in x, y, z directions be X,

Y and Z respectively as shown in Fig 2.5. The intensity of body forces are uniform
over entire body. Hence the total body force in x, y, z direction on the element
shown are given by
(i) X dx dy dz in x – direction
(ii) Y dx dy dz in y – direction and
(iii) Z dx dy dz in z – direction
Matrix Displacement Formulation
INTRODUCTION
Though mathematicians, physicists and stress analysts worked independently in

the field of FEM, it is the matrix displacement formulation of the stress analysts
which lead to fast development of FEM. Infact till the word FEM became popular,
stress analyst worked in this field in the name of matrix displacement method. In
matrix displacement method stiffness matrix of an element is assembled by direct
approach while in FEM though direct stiffness matrix may be treated as an
approach for assembling element properties (stiffness matrix as far as stress
analysis is concerned), it is the energy approached which has revolutionized entire
FEM.
Hence in this chapter, a brief explanation of matrix displacement method is
presented and solution techniques for simultaneous equations are discussed briefly.
MATRIV DISPLACEMENT EQUATIONS

The standard form of matrix displacement equation is,
[k] {}  {F}
where [k] is stiffness matrix
{} is displacement vector and
{F} is force vector in the coordinate directions
The element kij of stiffness matrix maybe defined as the force at coordinate i
due to unit displacement in coordinate direction j.
The direct method of assembling stiffness matrix for few standard cases is briefly given in this
article.
1. Bar Element
Common problems in this category are the bars and columns with varying cross
section subjected to axial forces as shown in Fig. 3.1.
For such bar with cross section A, Young’s Modulus E and length L (Fig. 3.2 (a))
extension/shortening 
is given by
  PL
EA
P1
L1
P2
A3
P3 A2 P2 A1 x
L2
P 1
L3 L2 L1
P3
Fig. 3.2
 P  EA 
L
 If   1, P  EA
L
By giving unit displacement in coordinate direction 1, the forces development
in the coordinate direction 1 and 2 can be found (Fig. 3.2 (b)). Hence from the
definition of stiffness matrix,
   EA
k11 E and L
A k21
L
Similarly giving unit displacement in coordinate direction 2 (refer Fig. 3.2 (c)), we get
k12   
EA and E
L k22 A
L
EA µ 1 1y
Thus, k
 L ¡µ1 …
(3.5)
1
¡j
2. Truss Element
Members of the trusses are subjected to axial forces only, but their orientation in
the plane may be at any angle to the coordinate directions selected. Figure 3.3
shows a typical case in a plane truss. Figure 3.4 (a) shows a typical member of the
truss with Young’s Modulus E, cross sectional area A, length L and at angle  to
x-axis
L
k  P sin  EA cos sin

2 L
1
  P cos   EA cos2 
k3
1 L
  P sin   EA cos sin
L
k
4
1
4
P P
3
2 1
P P
1 (a)
1
P 1 P
(b) (c
)
  P sin   EA sin2 
L
k42
(iii) Unit displacement in coordinate direction 3,
Extension along the axis is 1  sin and hence the forces developed
are as shown in the Fig. 3.4 (d)
 P  EA cos
L
k13   P cosL  EA cos2 

k   P sin   EA cos sin

2 L
3
 P cos  EA cos2 
k3
3 L
 P sin  EA cos sin
L
k
4
3
(vi) Due to unit displacement in coordinate direction 4,
Extension of the bar is equal to 1 and hence the forces developed are as shown in
 sin,
Fig. 3.4 (e).
 P  EA sin 
L
  P cos   EA sin cos
L
  P sin   EA sin2 
L
 P cos  EA sin cos
L
 P sin  EA sin2 
L
¡µ cos2  cos sin y¡

cos co
sin s2 
EA ¡ cos sin sin2  cos –sin2  ¡
 sin
k
L ¡ –cos2  cos cos2 cos sin ¡
sin 
¡µcos sin – cos sin2  ¡j
sin2 sin

µ¡ l2
lm lm y¡
l2
 EA ¡ lm m2 lm m2 ¡
L ¡l2 lm l2 lm ¡ …
(3.6)
¡µlm m2 lm m2 ¡j
Where l and m are the direction cosines of the member i.e. l = cos  and m = cos (90 –  ) = sin  .
(v) Beam Element
In the analysis of continuous beams normally axial deformation is
negligible (small deflection theory) and hence only two unknowns may be
taken at each end of a element (Fig. 3.5). Typical element and the
coordinates of displacements selected are shown in Fig. 3.5 (b). The end
forces
developed due to unit displacement in all the four coordinate directions are
shown in Fig. 3.6 (a, b, c, d).
y 1 3 5 7 9
E1, I1 E2, I2 E3, I3 E4, I4
x 2 L1 4 L2 6 L3 8 L4 10
z
(a)
1 3
E,I, L
From the definition of stiffness matrix and looking at positive senses indicated, we can write
(a) Due to unit displacement in coordinate direction 1,
     6EI
k11 12EI k21 6EI k31 12EI k L2
L
3 L
2 L
3 4
1
(b) Due to unit displacement in coordinate direction 2,
2
  4EI 
k12 6EI k22 L k32 6EI k42
L L2
 2EI L
(c) Due to unit displacement in coordinate direction 3,
      6EI
k13 12EI k 6EI k33 12EI k43 L2
L
3 L
2 L
3
2
3
(d) Due to unit displacement in coordinate direction 4,

    4EI
k14 6EI k24 2EI k34 6EI k44 L
L
2 L L
2
µ¡ 12 6L 12 6L y¡
 k  FI ¡ 6L
4L 6 2L2¡
µ3 6L 2 L …(3.7)
6L
L
6 12 j
4L2¡
L 6
2L L
2
If axial deformations in the beam elements are to be considered as in
case of columns of frames, etc. (Fig. 3.7), it may be observed that axial
force do not affect values of bending moment and shear force and vice
versa is also true. Hence stiffness matrix for the element shown in Fig.
3.8 is obtained by combining the stiffness matrices of bar element and
beam element and arranging in proper locations. For this case
µ¡ EA 0 0  0 0 y¡
EA
¡ L
0 6E
L
0  6E
I 12 I
¡ 12 EI
k  ¡¡ EI L L L
L 2 3 2
3 4EI 0  2EI ¡
1 6E L 6 L¡
I E
L2 I
L
2
¡ EA 0 0 0 ¡
E
¡ 0 A
L L
12 6EI …
¡ 0  12EI
 L3
6EI
EI ¡ (3.8)
¡ 0L2 L 
L2 ¡
3
¡ 0 6EI 2EI 0  4EI ¡ ¡
µ¡ L2 L 6EI
L2 L j
(a) (b)
Fig. 3.7
2 5
1 4
3 6
Fig. 3.8
The following special features of matrix displacement equations are worth noting:
(i) The matrix is having diagonal dominance and is positive definite. Hence
in the solution process there is no need to rearrange the equations to get
diagonal dominance.
(ii) The matrix is symmetric. It is obvious from Maxwell’s reciprocal theorem.
Hence only upper or lower triangular elements may be formed and others
obtained using symmetry.
(iii) The matrix is having banded nature i.e. the nonzero elements of stiffness
matrix are concentrated near the diagonal of the matrix. The elements
away from the diagonal are zero. Considerable saving is effected in
storage requirement of stiffness matrix in the memory of computers by
avoiding storage of zero values of stiffness matrices. The banded nature of
matrix is shown in Fig. 3.9.
I
II
III
IV
Fig. 3.9 V
VI
In this case instead of storing N ×N size matrix only N × B size matrix can be stored.
V
II VII
I
IX
no force in the z-direction and no variation of any forces in z-direction. Hence

 z   xz   yz  0
The conditions  xz   yz  0 give  xz   yz  0 and the condition  z  0 gives,
 z   x   y  1    z  0
Basic Equations in Elasticity 17

i.e.   

z 1 x y
If this is substituted in equation 2.13 the constitutive law reduces to

j¡ x y¡ µ¡1  0 y¡ j¡ x y¡
E
¡  ¡y
¡ ¡ ¡
2 ¡0 01y ¡ 1 0¡ ¡ …
 ¡
¡ 1 ¡ (2.14)
t xy 1 µ¡ ¡j xy 1
t
2
(c) (d)
no force in the z-direction and no variation of any forces in z-direction. Hence

 z   xz   yz  0
The conditions  xz   yz  0 give  xz   yz  0 and the condition  z  0 gives,
 z   x   y  1    z  0
i.e.   

z 1 x y
If this is substituted in equation 2.13 the constitutive law reduces to

j¡ x y¡ µ¡1  0 y¡ j¡ x y¡
E
¡  ¡y
¡ ¡ ¡
2 ¡0 01y ¡ 1 0¡ ¡ …
 ¡
¡ 1 ¡ (2.14)
t xy 1 µ¡ ¡j xy 1
t
2
x
x
z z
(c) (d)
Fig. 2.7
z
(e)
Axi-Symmetric Problems
Axi-symmetric structures are those which can be generated by rotating a line or

curve about an axis. Cylinders (refer Fig. 2.8) are the common examples of
axisymmetric structures. If such structures are subjected to axisymmetric loadings
like uniform internal or external pressures, uniform self weight or live load
uniform over the surface,
there exist symmetry about any axis. The advantage of symmetry may be made
use to simplify the analysis. In these problems cylindrical coordinates can be used
advantageously. Because of symmetry, the stress components are independent of
the angular (  ) coordinate. Hence all derivatives with respect to  vanish i.e. in
these cases.
v   r    z   r    z  0
z, w
z, w
In these cases stress-strain

relation is j1    0 y j r y

j r y
¡ z ¡ E ¡ 1 0 ¡ ¡ z ¡

  (1  )(1 1 0 …
(2.17)
 2) ¡
¡
¡ ¡ 1 2 ¡ ¡  ¡
¡
t rz 1 y 2 j t rz 1
Fig. 2.7
z
(e)
Axi-Symmetric Problems
Axi-symmetric structures are those which can be generated by rotating a line or

curve about an axis. Cylinders (refer Fig. 2.8) are the common examples of
axisymmetric structures. If such structures are subjected to axisymmetric loadings
like uniform internal or external pressures, uniform self weight or live load
uniform over the surface,
there exist symmetry about any axis. The advantage of symmetry may be made
use to simplify the analysis. In these problems cylindrical coordinates can be used
advantageously. Because of symmetry, the stress components are independent of
the angular (  ) coordinate. Hence all derivatives with respect to  vanish i.e. in
these cases.
v   r    z   r    z  0
z, w
z, w
In these cases stress-strain

relation is
J1    0 y J r y
Jr y 
¡ ¡ z E ¡ 1 
0
¡ ¡ ¡
z
 
 (1  )(1  2) 1 0 ¡ …
¡  (2.17)
¡ ¡ 1 2

¡ ¡¡ ¡
Sem-6th ,
t 1
rz y 2 J t 1
rz

Sa-Iii 6TH Sem Sukrit

Uploaded by

Copyright:

Available Formats

Sa-Iii 6TH Sem Sukrit

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Sa-Iii 6TH Sem Sukrit

Uploaded by

Copyright:

Available Formats

GATEWAY EDUCATION, SONEPAT

∆b = Upward Deflection of point B on primary structure due to all causes

 ∆b = ∆bo + ∆bb Substituting for ∆bb –

∆b = 0 we get Xb = – ∆bo / δbb

FORCE METHOD DISPLACEMENT METHOD

To study matrix methods there are some pre-requisites :

M – B M due to applied loads & mj – B M due to unit load at j

When unit load is applied at i and is called flexibility coefficient.

Roller Support : No. of reactions,r = 1

Hinged Support : No. of reactions, r = 2

Fixed Support : No. of reactions, r = 3

No. of reaction components r = 5 (as shown)

Fixed beam : Kinematically determinate :

Simply supported beam Kinematically indeterminate

Any joint – Moves in three directions in a plane structure

i.e. reaction components prevent the displacements  no. of restraints = no. of

Reaction components r=6

∆1 = ∆1L + ∆1R Where ∆1 - - - - ∆n = Displ. At Co-ord.at 1,2 - -n

The above equations may be return as [∆] = [∆L] + [∆R]---------(1)

∆n = ∆nL + δn1 P1 + δn2 P2 +---------δnnPn

If the net displacements at the redundants are zero then

P1 = P’1 +K11 ∆1+ K12, ∆2+ K13∆3 +-----------K1n ∆n

TWO CO-ORDINATE SYSTEM

δ11=L3 / 3EI δ21=L2 / 2EI

Unit Force At Co-ord.(2) δ12=δ21 =L2 / 2EI δ22=L / EI

The Flexibility matrix can be developed by

δij = ∫ mi mj / EI x ds δij =displacement at I due to unit load at j

INVERSING THE FLEXIBILITY MATRIX [ δ ]

 Condition of Compatibility – connected ends (nodes) of adjacent

Model spring system

Derivation of the Stiffness Matrix for a Single Spring Element

The ‘global’ stiffness relation is written in Eqn.16, which we

split it into its component parts in the following way

(Note that, to create the global stiffness matrix by assembling the

We want to populate the cells to generate the global stiffness matrix.

 The term in location ii consists of the sum of the direct stiffnesses of

 Add a negative for reaction terms (–kij)

This type of assembly process is handled automatically by

Solving for (u)

To find {u} solve

Beam Stiffness Matrix.

The stiffness matrix is symmetrical. The stiffness matrix is partitioned to separate

The forces are equal to,

q11  k11u'1 ; q21  k21u'1 ; q31  k31u'1 ; q41  k41u'1 (27.3)

This may also be written in compact form as, q  k  u'

Beam (global) Stiffness Matrix.

order 6  6 . From Fig. 27.6b, K  may be written as,

Formation of load vector.

Solution of equilibrium equations

⎧⎪pk ⎫⎪ ⎡k11  k12 ⎤⎧⎪uu ⎫⎪

Le R1 , R3 and R be the reactions along the constrained degrees of freedom as

From the method of superposition,

Equation (27.14a) may be written as,

Member end actions q , q , q , q are calculated as follows. For example consider

[k]e {}e  {F}e

where [k]e is element stiffness matrix, { }e is nodal displacement vector of the

Fig. 1.1 †ypical tension flat

Fig. 1.2 Discretisation of quarter of tension flat

Step 2: The portion to be analysed is to be discretised. Fig. 1.2 shows discretised