Finite Element Method (MSTR 515)

Module 2

The Stiffness Method and The Plane Truss

Course Instructure
Shyam Sundar Khadka, PhD
Assistant Professor
Department of Civil Engineering, Kathmandu University
Contents
Introduction
Structure Stiffness Equations
Common Problems
• Bar Element
• Truss Element
• Beam Element
Nature of Stiffness (K)
Solution of Matrix Displacement Equations
Flowchart of matrix method
• Steps in solving a problem with examples

Department of Civil Engineering 2

 Introduction
In this module the stiffness method introduced in Module 1 is explained in detail
• How to assemble elements to form a structure
• How to impose the boundary conditions and
• How to solve the equations

• Plane truss is used to illustrate these concepts, which apply

with equal forces to all finite element structure

• The bars of a truss are its finite elements

• Truss is a “natural” finite element structure, since no

conceptual division into element is needed

Department of Civil Engineering 3

Introduction (contd…..)

In this module, each bar is assumed to be uniform, linearly elastic, pin-connected at

its ends and axially loaded

Actual displacement are assumed to be small enough that if θ is the angle of rotation
of any bar, Sin θ= θ and Cos θ= θ

Consider only static problems. Within these restriction, the analysis is exact, not
approximate.

Department of Civil Engineering 4

Structure Stiffness Equations
The standard form of matrix displacement equation is
[K] {d}= {F}
[K]= stiffness matrix, {d} = displacement vector ,and {F} = force vector in coordinate direction
𝐹𝑖 𝐿𝑖
𝑒𝑖 =
𝐴𝑖 𝐸𝑖
The element kij of stiffness matrix may be defined as the force at coordinate i due to
due to unit displacement in coordinate direction j

Nodes and elements (bars) are numbered arbitrarily

Ai= cross-sectional area,Ei = elastic modulus , Li= length

From the elementary mechanics of material,
𝐹𝐿 Fig. 3 elements plane truss. DoF
axial force Fi and change in length ei have the relation 𝑒𝑖 = 𝐴𝑖𝐸𝑖 u1,v1 and v3 active while DoF
𝑖 𝑖
u2,v2 and u3 are restrained
Since, Stiffness is defined as the ratio of force to displacement and
𝐹 𝐴𝐸
is given the symbol k. Thus for bar i , 𝑘𝑖 = 𝑒𝑖 = 𝐿𝑖 𝑖
𝑖 𝑖

Department of Civil Engineering 5

Common Problems

Bar Element Bars and Columns with varying Cross-Section, subjected

to axial force

For bar element having cross-section area A, elastic

Modulus E and length L, extension/shortening is given by

𝑃𝐿 𝐴𝐸∆
∆= and 𝑃 =
𝐴𝐸 𝐿

𝐀𝐄
if ∆= 1; 𝐏 =
𝐋

Department of Civil Engineering 6

 Fig. Beam/Bar element

Unit displacement is applied at node 1, along x-direction, fig. b, the force development at
nodes 1 and 2 can be found

Hence, from definition of stiffness matrix

Similarly, unit displacement at node 2, fig. c above

Department of Civil Engineering 7

Truss Element

Members of truss element are subjected to axial force only, but their
orientation in the plane may be at any angle

First, apply unit displacement at node 1, coordinate 1, along x-

direction fig. b, displacement along x axis =1* Cos θ
𝐀𝐄
Force developed at the ends are P = Cos θ
𝐋

From the definition of element of stiffness matrix, we get

AE AE
𝑘11 = PCosθ = 𝐶𝑜𝑠 2 θ;𝑘21 = PSinθ = CosθSinθ
L L

AE AE
𝑘31 = −PCosθ = − 𝐶𝑜𝑠 2 θ;𝑘41 = −PSinθ = CosθSinθ
L L

Department of Civil Engineering 8

Truss Element (Contd…..)
Second, apply unit displacement at node 1, Coordinate 2, case fig. c.
axial deformation =1*Sin θ and developed force at each ends are shown
𝐀𝐄
in Fig. c. 𝐏 = 𝐋 Sin θ
Therefore,
AE AE
𝑘12 = PCosθ = L SinθCosθ; 𝑘22 = PSinθ = L Sin2θ
AE AE
𝑘32 = −PCosθ = − SinθCosθ;𝑘42 = −PSinθ = − Sin2θ
L L

Third, unit displacement at node 2, coordinate 3, fig. d, extension along

the axis= 1*Cos θ
𝐀𝐄
Force developed, 𝐏 = 𝐋 Cos θ
AE 2 θ;𝑘 AE
𝑘13 = −PCosθ = − 𝐶𝑜𝑠 23 = −PSinθ = − CosθSinθ
L L

AE AE
𝑘33 = PCosθ = 𝐶𝑜𝑠 2 θ; 𝑘43 = PSinθ = CosθSinθ
L L

Department of Civil Engineering 9

 Truss Element (Contd…..)
Fourth, Fig. e, apply unit load 2, coordinate 4, bar extension= 1* Sin θ, force Direction Cosines
developed,
𝐀𝐄
𝐏 = 𝐋 Sin θ

Therefore,
AE AE
𝑘14 = PCosθ = L SinθCosθ; 𝑘24 = −PSinθ = − L Sin2θ
AE AE
𝑘34 = PCosθ = SinθCosθ;𝑘44 = PSinθ = Sin2θ
L L

The Stiffness matrix is

𝐶𝑜𝑠 2 θ SinθCosθ −𝐶𝑜𝑠 2 θ −SinθCosθ

𝐾 =𝐿
𝐴𝐸 SinθCosθ Sin2θ −SinθCosθ −Sin2θ
−𝐶𝑜𝑠 2 θ −SinθCosθ 𝐶𝑜𝑠 2 θ SinθCosθ
−SinθCosθ −Sin2θ SinθCosθ Sin2θ
𝑙2 𝑙𝑚 −𝑙2 −𝑙𝑚
𝐴𝐸 𝑙𝑚 𝑚2 −𝑙𝑚 −𝑚2 where, l and m are the Direction Cosines of the
=
𝐿 −𝑙2 −𝑙𝑚 𝑙2 𝑙𝑚 member le l=Cosθ and m =Cos(90−θ)= Sin θ
−𝑙𝑚 −𝑚2 𝑙𝑚 𝑚2

Department of Civil Engineering 10

Beam Element

Analysis of continuous beams normally axial deformation is negligible and hence only
two unknowns may be taken at each end of the element

Department of Civil Engineering 11

Beam Element (Contd….)
If axial deformation in the beam element are
to be considers: in case of frames

𝐸𝐴 𝐸𝐴
0 0 − 0 0
𝐿 𝐿
12𝐸𝐼 6𝐸𝐼 12𝐸𝐼 6𝐸𝐼
0 0 −
𝐿3 𝐿2 𝐿3 𝐿2
6𝐸𝐼 4𝐸𝐼 6𝐸𝐼 2𝐸𝐼
0 0 − 𝐿2
𝐿2 𝐿2 𝐿
𝐾= 𝐸𝐴 𝐸𝐴
− 0 0 0 0
𝐿 𝐿
12𝐸𝐼 6𝐸𝐼 12𝐸𝐼 6𝐸𝐼
0 − − 0 −
𝐿3 𝐿2 𝐿3 𝐿2
The Stiffness matrix is 6𝐸𝐼 2𝐸𝐼 6𝐸𝐼 4𝐸𝐼
0 0 − 𝐿2
12 6𝐿 −12 6L 𝐿2 𝐿 𝐿
𝐸𝐼 6𝐿 4𝐿2 −6L 2L2
𝐾= 3
𝐿 −12 −6L 12 −6L
6L 2L2 −6L 4L2

Department of Civil Engineering 12

Nature of Stiffness (K)
• Each diagonal coefficient of 𝑘𝑖𝑗 is +ve. If 𝑘𝑖𝑖 were zero, displacement 𝐷𝑖 would generate no resisting force, which
implies that the structure is unstable
• If 𝑘𝑖𝑖 were negative, the force and its corresponding displacement would be oppositely directed, which is physically
unreasonable
• The matrix is symmetric. Only upper or lower triangular elements may be formed and others obtained using symmetry.
[structure that has a linear force displacement relationship, Betti-Maxwell reciprocal theorem]
• Each column of [K] sums to zero because each column represent a set of nodal forces that satisfies static equilibrium.
Caution: If {D} contains linear and rotational DoF, for a beam or a frame, column do not sum to zero.
• The matrix is having banded nature i.e. the non-zero elements of matrix are concentracted near the diagonal of the
matrix. The elements away from the diagonal are zero

• It is considerable save the storage of data in computer by avoiding

storage of zero values of stiffness matrix.

Department of Civil Engineering 13

 Solution of Matrix Displacement Equations

• Matrix displacement equations are linear simultaneous equations

• These equations can be solved using Gaussian Elimination Method

The method is based on the idea of reducing the given system of equations Ax = B
to an upper triangular system of equations Ux = z, using elementary row operations.
That is, the solutions of both the systems are identical. Here 𝑥 = 𝑢𝑛𝑘𝑛𝑜𝑤𝑛 𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒′𝑠 𝑚𝑎𝑡𝑟𝑖𝑥

We illustrate the method using the 3 × 3 system

𝑎11𝑥 + 𝑎12𝑦 + 𝑎13𝑧 = 𝑏1 ∶ R1
𝑎21𝑥 + 𝑎22𝑦 + 𝑎23𝑧 = 𝑏2 ∶ 𝑅2
𝑎31𝑥 + 𝑎32𝑦 + 𝑎33𝑧 = 𝑏3 ∶ 𝑅3
To make the pattern of equation upper triangular matrix we have to use the operation :
𝑅𝑖 = 𝑅𝑖 − (𝑎𝑖𝑘 /𝑎𝑘𝑘)∗ 𝑅𝑘 Where i=(k+1), (k+2),……n ; k= coloumn number

Department of Civil Engineering 14

Solution of Matrix Displacement Equations (contd…..)
So after the operation, the equation becomes :
𝑎11𝑥 + 𝑎12𝑦 + 𝑎13𝑧 = 𝑏1
0+𝑎′22𝑦 + 𝑎′23𝑧 = 𝑏′2 [𝑅′2 = 𝑅2 −(𝑎21/𝑎11)∗ 𝑅1 ]
0+𝑎′32𝑦 + 𝑎′33𝑧 = 𝑏′3 [ 𝑅′3 = 𝑅3 −(𝑎31/𝑎11)∗ 𝑅1 ]

After that we have to apply the same process to make that pattern to a upper triangular
matrix.
𝑎11𝑥 + 𝑎12𝑦 + 𝑎13𝑧 = 𝑏1
0+𝑎′22𝑦 + 𝑎′23𝑧 = 𝑏′2
0 + 0 + 𝑎′′33𝑧 = 𝑏′′3 [ 𝑅′′3 = 𝑅′3 −(𝑎′32/𝑎′22)∗ 𝑅′2 ]
We can then successively find the values z, y, and x

Department of Civil Engineering 15

Solution of Matrix Displacement Equations (contd…..)
For example we take three equations:
𝑥 − 2𝑦 + 9𝑧 = 8
3𝑥 + 𝑦 − 𝑧 = 3
2x − 8𝑦 + 𝑧 = −5

Now we have to make it a upper triangular matrix using the operation:

𝑅2 = 𝑅2 −(𝑎21/𝑎11)∗ 𝑅1 and 𝑅3 = 𝑅3 − (𝑎31/𝑎11) ∗ 𝑅1

They become :
𝑥 − 2𝑦 + 9𝑧 = 8
0 + 7𝑦 − 28𝑧 = −21
0 − 4𝑦 − 17𝑧 = −21

To get the pattern of upper triangular matrix we again do the operation:

Department of Civil Engineering 16

Solution of Matrix Displacement Equations (contd…..)

The operation is : 𝑅3 = 𝑅3 − (−4/7) ∗ 𝑅2

Now, they become :

𝑥 − 2𝑦 + 9𝑧 = 8
0 + 7𝑦 − 28𝑧 = −21
0 + 0 − 33𝑧 = −33

So that :
𝑧 = 1; 𝑦 =(−21+28𝑧)/7= 1;
𝑥 = (8 − 9𝑧 + 2𝑦) = 1 .

The set of solution is : x=1; y=1; z=1.

Department of Civil Engineering 17

Flowchart of matrix method
Classification
of members Stiffness matrices are composed according to
member models

Stiffness matrices for

members
Stiffness matrices are transformed from local
to global coordinates
Transformed stiffness
matrices
Stiffness matrices of separate members are
Final equation assembled into a single stiffness matrix K
F=K·Z

Unknown displacements and reaction forces

Stress-strain state of are calculated
structure

Department of Civil Engineering 18

Steps in solving a problem

Step 1: Write down the node-element connectivity table inking local and global nodes;
also form the table of direction cosines (l, m)

Step 2: Write down the stiffness matrix of each element in global coordinate system
with global numbering

Step 3: Assemble the element stiffness matrices to form the global stiffness matrix for
the entire structure using the node element connectivity table

Step 4: Incorporate appropriate boundary conditions

Step 5: Solve resulting set of reduced equations for the unknown displacements

Step 6: Compute the unknown nodal forces

Department of Civil Engineering 19
Node element connectivity table

ELEMENT Node 1 Node 2

1 1 2
2 2 3

3 3 1

1 2 (x2,y2)
L
El 1 60 El 3
θ
2 60 60 3 1 (x ,y )
1 1
El 2
Department of Civil Engineering 20
Stiffness matrix of element 1
Stiffness matrix of element 2
d1x d1y d2x d2y
d2x d2y d3x d3y
  d1x
    d2x
 
=  d1y
(1)
=  d2y
k ( 2)
  d2x k
    d3x
  d2y  
  d3y
Stiffness matrix of element 3
There are 4 degrees of
d3x d3y d1x d1y
freedom (dof) per element
  d3x (2 per node)
 
=  d3y
( 3)
k
  d1x
 
  d1y
Department of Civil Engineering 21
 (1)
k
Global stiffness matrix
d1x d1y d2x d2y d3x d3y
  d1x
  ( 2)
  d1y k
  d2x
K=  d2y
 
  d3x
( 3)
k
  d3y
  66

How do you incorporate boundary conditions?

Department of Civil Engineering 22

 Example 1 The length of bars 12 and 23 are equal (L)
y E: Young’s modulus
3 A: Cross sectional area of each bar
El#2 P2
Solve for
P1 (1) d and d
2x 2y

El#1 (2) Stresses in each bar

2
45o
x
1 Solution

Step 1: Node element connectivity table

ELEMENT Node 1 Node 2

1 1 2
2 2 3

Department of Civil Engineering 23

Table of nodal coordinates
Node x y
1 0 0
2 Lcos45 Lsin45

3 0 2Lsin45

Table of direction cosines

ELEMENT Length x2 − x1 y −y
l= m= 2 1
length length
1 L cos45 sin45
2 L -cos45 sin45

Department of Civil Engineering 24

Step 2: Stiffness matrix of each element in global coordinates with global numbering
Stiffness matrix of element 1
 l2 lm −l 2 −lm 
 
EA  lm m2 −lm −m2 
k =
(1)

L  −l 2 −lm l2 lm 
 
 −lm − m 2
lm m2 

d1x d1y d2x d2y

 1 1 −1 −1  d1x
 1 1 −1 −1 
EA   d1y
=
2L  −1 −1 1 1  d2x
 
 −1 −1 1 1  d2y

Department of Civil Engineering 25

 Stiffness matrix of element 2

d2x d2y d3x d3y

 1 −1 −1 1  d2x
 −1 1 1 −1 
EA   d2y
=
(2)
k
2L  −1 1 1 −1  d3x
 
 1 −1 −1 1  d3y

Department of Civil Engineering 26

Step 3: Assemble the global stiffness matrix

 1 1 −1 −1 0 0 
 1 1 −1 −1 0 0 
 
EA  −1 −1 2 0 −1 1 
K=  
2L  −1 −1 0 2 1 −1
 0 0 −1 1 1 −1
 
 0 0 1 −1 −1 1 

The final set of equations is Kd = F

Department of Civil Engineering 27

Step 4: Incorporate boundary conditions
 0 
 0 
 
d2 x 
d = 
d 2 y 
 0 
 
 0 

Hence reduced set of equations to solve for unknown displacements at node 2

EA  2 0   d 2 x   P1 
 = 
2L 
0

2   d 2 y   P2 

Department of Civil Engineering 28

Step 5: Solve for unknown displacements

 P1L 
d2 x  
 EA  
 = 
d
 2 y   P2 L 

 EA  
Step 6: Obtain stresses in the elements 0

For element #1:  d1x  0

d 
E 1 1 1 1   1y 
 = −
(1)
−   
L 2 2 2 2  d2 x 
d 2 y 
E P +P
= (d 2 x + d 2 y ) = 1 2
2L A 2

Department of Civil Engineering 29

For element #2: d2 x 
d 
E 1 1 1 1   2y 
 = 
(2)
− −    0
L 2 2 2 2   d3 x  0
 d3 y 
E P −P
= (d 2 x − d 2 y ) = 1 2
2L A 2

Department of Civil Engineering 30

Problem 1: For the plane truss
P=1000 kN,
y L=length of elements 1 and 2 = 1m
P El#2 3 E=210 GPa
A = 6×10-4m2 for elements 1 and 2
2 = 6 2 ×10-4 m2 for element 3
El#1
El#3
Determine the unknown displacements and reaction
45o forces.
x
1

Solution
Step 1: Node element connectivity table ELEMENT Node 1 Node 2
1 1 2
2 2 3
3 1 3

Department of Civil Engineering 31

 Table of nodal coordinates
Node x y
1 0 0
2 0 L

3 L L

Table of direction cosines

ELEMENT Length x2 − x1 y2 − y1
l= m=
length length
1 L 0 1
2 L 1 0

3 L 2 1/ 2 1/ 2

Department of Civil Engineering 32

Step 2: Stiffness matrix of each element in global coordinates with global numbering
Stiffness matrix of element 1
 l2 lm −l 2 −lm 
 
EA  lm m2 −lm −m2 
=
(1)
k
L  −l 2 −lm l2 lm 
 
 −lm − m 2
lm m2 

d1x d1y d2x d2y

0 0 0 0  d1x
-4  
(210 10 )(6 10 ) 0 1
9
0 −1  d1y
=
1 0 0 0 0  d2x
 
0 −1 0 1  d2y

Department of Civil Engineering 33

Stiffness matrix of element 2 d2x d2y d3x d3y
1 0 −1 0  d2x
-4  
(210 10 )(6 10 )  0
9
0 0 0 
= d2y
(2)
k
1  −1 0 1 0  d3x
 
0 0 0 0  d3y
Stiffness matrix of element 3
d1x d1y d3x d3y
 0.5 0.5 −0.5 −0.5  d1x
 
(210 109 )(6 2 10-4 )  0.5 0.5 −0.5 −0.5  d1y
=
(3)
k
2  −0.5 −0.5 0.5 0.5  d3x
 
 −0.5 −0.5 0.5 0.5  d3y

Department of Civil Engineering 34

Step 3: Assemble the global stiffness matrix

 0.5 0.5 0 0 −0.5 −0.5

 0.5 1.5 0 −1 −0.5 −0.5
 
 0 0 1 0 −1 0 
K = 1260 10 
5
 N/m
 0 −1 0 1 0 0 
 −0.5 −0.5 −1 0 1.5 0.5 
 
 −0.5 −0.5 0 0 0.5 0.5 

The final set of equations is Kd = F Eq(1)

Department of Civil Engineering 35

Step 4: Incorporate boundary conditions y
x
 0  y
 0  3
  P El#2

d 2 x 

d =  2
 0  El#1
 d3 x  El#3
  o
 3y 
 d  45
x
1

Also, d 3y = 0 in the local coordinate system of element 3

How do I convert this to a boundary condition in the global (x,y)

coordinates?
Department of Civil Engineering 36
 y
x
 F1 x  y
F  3
 1y  P El#2

 P  
F =  2
F
 2y  El#1
 F3 x  El#3
 
 F3 y 
  45o
x
1

Also, F 3x = 0 in the local coordinate system of element 3

How do I convert this to a boundary condition in the global (x,y)

coordinates?
Department of Civil Engineering 37
 Using coordinate transformations

d 3x   l m   d3 x  1
 =  −m    l=m=

d 3 y 
  l  d3 y  2

 1 1   1 

d 3x   2 2   d  
 2
( d 3x + d 3 y )

 =    3x  =  
  − 1
d 3 y  1  d3 y   1
( d − d ) 

 2 2 
 2
3 y 3 x



d 3y = 0 (Multi-point constraint)

 d 3y =
1
2
( d3 y − d3 x ) = 0

 d3 y − d3 x = 0 Eq (2)

Department of Civil Engineering 38

 Similarly for the forces at node 3

 F 3x 
  l m   F3 x  1
 =  −m    l=m=

F 3 y 
  n   F3 y  2

 1 1   1 
  
 F 3x  2 2   F  
 2
( F3x + F3 y )

 =    3x  =  
  − 1
F 3 y  1   F3 y   1
( F − F ) 

 2 2  2

3 y 3 x



F 3x = 0
 F 3x =
1
2
( F3 y + F3 x ) = 0

 F3 y + F3 x = 0 Eq (3)

Department of Civil Engineering 39

Therefore we need to solve the following equations simultaneously

Kd = F Eq(1)
d3 y − d3 x = 0 Eq(2)

F3 y + F3 x = 0 Eq(3)

Incorporate boundary conditions and reduce Eq(1) to

1 −1 0  d 2 x   P 
   
1260 105 
 −1 1.5 0.5 =
  3x   3x 
d F
 0.5 d   F 
0 0.5   3y   3y 

Department of Civil Engineering 40

 Write these equations out explicitly

1260 105 (d 2 x − d3 x ) = P Eq(4)

1260 105 (−d 2 x + 1.5d3 x + 0.5d3 y ) = F3 x Eq(5)
1260 105 (0.5d3 x + 0.5d3 y ) = F3 y Eq(6)

Add Eq (5) and (6)

1260 105 ( −d 2 x + 2d3 x + d3 y ) = F3 x + F3 y = 0 using Eq(3)

 1260 105 (−d2 x + 3d3 x ) = 0 using Eq(2)

 d 2 x = 3d3 x Eq(7)
 1260  105 (3d 3 x − d 3 x ) = P
Plug this into Eq(4)
 2520  105 d 3 x = 106
Department of Civil Engineering 41
  d 3 x = 0.003968m
d 2 x = 3d 3 x = 0.0119m

Compute the reaction forces

 F1x  0 −0.5 −0.5 
F  0
 1y 
 −0.5 −0.5 

d 2 x 
  
 F2 y  = 1260  10  0 
5
0 0  d3 x 
F   
 −1 1.5 0.5   d 3 y 

 3x 

 F3 y 
 
0 0.5 0.5 
 −500 
 −500 

 

=  0  kN
 −500 
 

 500 
Department of Civil Engineering 42
Thank You.

Department of Civil Engineering 43