Introduction to Matrix Methods
of Structural Analysis
CE 222: Structural Analysis
Durgesh C Rai
Department of Civil Engineering
Indian Institute of Technology Kanpur
Kanpur 208016
Historical Context
Airplanes and Digital Computers
revolutionize structural analysis
In 40s -50s, structural engineers
confronted with two highly
indeterminate structural systems
Swept-wing aircrafts in 1947
analyzed using the force
method (compatibility)
simultaneous linear
algebraic equations with
flexibility as coefficients
were manipulated using
matrix algebra
CE222: Structural Analysis/Dr Durgesh Rai/2008/2
Historical Context
Airplanes and Digital Computers
revolutionize structural analysis
In 40s -50s, structural engineers
confronted with two highly
indeterminate structural systems
Delta-wing aircraft in 1953
analyzed using the
displacement method
(equilibrium)
simultaneous linear
algebraic equations with
stiffness as coefficients were
manipulated using matrix
algebra
CE222: Structural Analysis/Dr Durgesh Rai/2008/3
Historical Context
Using Energy Principles and Finite Element
Method
Structures must be in equilibrium, with their
displacement in a compatible state and
material laws satisfied.
Their behaviour can be studied by
Solving fundamental equations
Employing energy principles
Argyris and Kelsey in 1954 formulated matrix
structural analysis using energy principles
In 1950s, a group of Boeing engineers showed
that the method can be extended to continuum
(i.e., Finite Element Method)
CE222: Structural Analysis/Dr Durgesh Rai/2008/4
Historical Context
Computer-oriented versus
Matrix
Matrix methods were suitable
for use with digital computers
ENIAC (1946) contained
17,468 vacuum tubes, weighed
30 tons and was roughly 2.6 m
by 0.9 m by 26 m and
consumed 150 kW of power.
Transistors in 1947 and silicon
chip in 1959 pushed
development of computer and
structural analysis
CE222: Structural Analysis/Dr Durgesh Rai/2008/5
Historical Context
Computer-oriented Methods
Matrix methods are no new structural
theory
They are simple system of notations and
book-keeping amenable to computer
Computer employ on engineering
judgment
Engineers must use their judgment and
understanding of structural behaviour for
proper use of computer and interpreting
computer output.
Garbage in, garbage out, applies more often
than not!
CE222: Structural Analysis/Dr Durgesh Rai/2008/6
STIFFNESS METHOD USING BASIC EQNS
Truss Element
Truss is composed of individual discrete
elements connected at node (joints).
Basic equations has to be satisfied within each
element and also by the structure.
Prismatic, cross-sectional area A
Force is collinear with centroidal axis
Element is straight and plane sections remain L
p p
normal.
Small displacements meaning strain-
displacement eqn are linear and equilibrium
equations relate to undeformed state
Material is homogeneous, isotropic and
linearly elastic
CE222: Structural Analysis/Dr Durgesh Rai/2008/7
TRUSS ELEMENT USING BASIC EQNS
P =x A equilibrium is satisfied
' u ' is continuous function of x geometrically compatible
strains, no rips or tears associated with displn. strain-displn.
relation
L
p p
Strain-displn. relation x x
u
x + x x
o ' a ' oa x u
x = lim = lim = p A x
x 0 oa x 0 x x
u
Constitutive law
u 0 u+ x
0
a x
x = E x a a
x x
P u du
=E =E Pdx = AE du
A x dx
CE222: Structural Analysis/Dr Durgesh Rai/2008/8
TRUSS ELEMENT USING BASIC EQNS
Integrate
px
u= + c1
AE
u(0) = 0 c1 = 0 and x = k , u = L
AE
p= L , L is total displn
L
If the element is a part of truss, then both joints move,
u ( 0 ) = ui and u ( L ) = u j then c1 = ui
p =
AE
L
(
u j ui )
Invoking conditions of equilibrium, compatibility
and relations, force displn of axial element is obtained.
CE222: Structural Analysis/Dr Durgesh Rai/2008/9
ELEMENT STIFFNESS MATRIX
p yj , v j
Global Co-ordinates
Local co-ordinates
j
pxj u
j
y
Degrees of Freedom
No of displn. co-ordinate to describe
pxi x
uniquely its position in space ui i
It has infinite dof. But all other points p , v
yi i
can be described uniquely if end axial
deflections are known.
p = pxi pyi pxj pyi
Four independent dof. 2 for each
joint. u = ui vi u j v j
CE222: Structural Analysis/Dr Durgesh Rai/2008/10
ELEMENT STIFFNESS MATRIX
First displn. case
ui 0; u j = v j = ui = 0 AE
( ui cos )
L
Shortening of element = ui cos
(Based on assumption that the ui y
is small compared to L)
AE
piI = ui cos ui cos
L
Its components: pxiI x
AE AE ui
pxI i = ui cos cos = ui cos 2 (IA)
L L piI pyiI
AE (IB)
pyI i = ui cos sin
L
AE (IC)
px Ij = p x iI = ui cos 2
L
AE (ID)
py Ij = p y iI = ui sin 2 cos
L
CE222: Structural Analysis/Dr Durgesh Rai/2008/11
ELEMENT STIFFNESS MATRIX
Second displn. case
AE
vi 0; ui = u j = v j = 0 ( ui cos )
L
Shortening of element = vi sin y
Its components
AE
pxiII = piII cos = vi sin cos vi
L
pxiII x
AE
pyiII = piII sin = vi sin 2 vi sin
L AE
piII = vi sin
L pyiII
Equilibrium requires
AE
px IIj = p II xi = vi sin cos
L
AE
py IIj = pyi
II
= vi sin 2
L
CE222: Structural Analysis/Dr Durgesh Rai/2008/12
ELEMENT STIFFNESS MATRIX
Third element displn. case
u j 0; ui = vi = v j = 0
AE
px III
j = p III
xi = u j cos 2
L
AE
py III
j = p III
yi = u j sin cos
L
Fourth displn. case
u j 0; ui = vi = u j = 0
AE
pxj IV = p IV
x
= v j sin cos
i L
AE
pyj IV = pyi IV = v j sin 2
L
CE222: Structural Analysis/Dr Durgesh Rai/2008/13
ELEMENT STIFFNESS MATRIX
General Displacement Case:
vi 0 vi 0 u j 0 v j 0
Superposition is valid, for linear elastic system
undergoing small displn.
px i = px Ii + px IIi + px III
i
+ p x
IV
i
=
AE 2
L
( )
cos ui + ( sin cos ) vi cos 2 u j ( sin cos ) v j
AE
py i = ( sin cos ) ui + sin 2
vi ( sin cos ) u j sin 2
v j
L
px j =
AE
L
( )
cos 2 ui + ( sin cos ) vi + cos 2 u j + ( sin cos ) v j
py j =
AE
L
( ) (
( sin cos ) ui sin 2 vi + ( sin cos ) u j + sin 2 v j
)
CE222: Structural Analysis/Dr Durgesh Rai/2008/14
ELEMENT STIFFNESS MATRIX
In matrix form:
px i cos 2 sin cos cos 2 sin cos ui
p
yi sin cos sin 2 sin cos 2
sin vi
px =
j cos
2
sin cos cos 2 sin cos u j
p v
y j sin cos sin 2 sin cos sin 2 j
p=k u
Important properties
Column of k contains forces associated with specified displaced config.
e.g. for u = [ 0 vi 0 0 ] gives forces corresponding 2nd displaced condn.
Sum of column element = 0, (inherent of truss element only).
(
k is symmetric kij = k ji
)
k has no inverse. (because element is not constrained
and is free to move to undergo rigid body motion)
CE222: Structural Analysis/Dr Durgesh Rai/2008/15
NODAL EQUILIBRIUM OF THE STRUCTURE
Basic Eqns. has to be satisfied at Pxa, Ua Pxb, Ub Pxc, Uc
structure level.
Interelement compatibility a A ,L ,E b A ,L ,E c
1 1 1 2 2 2
A1 E1 A E
Equilibrium k1 = ; k2 = 2 2 ;
L1 L2
Ist Assemblage displn. case PxaI PxbI PxcI
U a 0; U b = Uc = 0
k1Ua k1Ua 0 0
IInd Assemblage Displn. case
PxaII PxbII PxcII
U b 0; U a = Uc = 0
k1Ub k1Ub k2Ub k2Ub
IIIrd Assemblage Displn. case
PxaIII PxbIII PxcIII
Uc 0; U a = U b = 0
0 0 k2Uc k2Uc
CE222: Structural Analysis/Dr Durgesh Rai/2008/16
NODAL EQUILIBRIUM OF THE STRUCTURE
General Assemblage Case Pxa, Ua Pxb, Ub Pxc, Uc
U a 0; Ub Uc 0
I II III a A ,L ,E b A ,L ,E c
Pxa =Pxa +Pxa + Pxa 1 1 1 2 2 2
= + k1 U a k1U b + 0 Uc k1 =
A1 E1 A E
; k2 = 2 2 ;
L1 L2
Pxb = k1 U a + k1 U b + k2 Ub k2 Uc
Pxc = 0 U a k2 U b + k2 Uc
In Matrix Form
Pxa k1 k1 0 U a K : column of K = equilibrium force
P = k k1 + k 2 k2 U b for a prescribed
xb 1 displaced shape.
Pxc 0 k2 k2 Uc
P=K U
K11 = Pxa due to unit displn. of U a with Ub = Uc = 0
K 12 = Pxa due to unit displn. U (with U a = Uc = 0 )
b
CE222: Structural Analysis/Dr Durgesh Rai/2008/17
NODAL EQUILIBRIUM OF THE STRUCTURE
Get element stiffness matrix
= 0 and ignore d.o.f. in y-direction ( Pyj , pyj , vi u j = 0 )
pxi AE 1 1 ui Pxa, Ua Pxb, Ub Pxc, Uc
p = 1 1 u
xj L j
a k1 b k2 c
p=k u
k1 k1 k2 k2
k ab
= ; k bc
= k ab
k1 k1 k2 k2 Ua Ub Uc
Node equilibrium of assemblage is a
convenient method of assembling structural
stiffness K k=
Nodal equilibrium is assured by combining
stiffness coefficient where they overlap.
Start with a null matrix of order = structure
DOF and add element stiffness matrix. k bc
CE222: Structural Analysis/Dr Durgesh Rai/2008/18
NODAL EQUILIBRIUM OF THE STRUCTURE
Example
1 2 3 4 5 6
6 Ua Va Ub Vb Uc Vc
c 1
5
2
3 2 3
4
4
a 1
3 5
2 1 b
6
Addition is different from matrix addition.
Sub matrix addition.
Enlarge element stiffness matrix to be conformable with
structure stiffness matrix by filling zero at appropriate places.
CE222: Structural Analysis/Dr Durgesh Rai/2008/19
ASSEMBLAGE OF ELEMENT STIFFNESS MATRICES
Example y
2
Aab = 600 mm
Abc = 1000 mm2 c
Aac = 800 mm2
E = 200 GPa 4m
AE kN
= 40,000
L m
x
a b
For ab , = 0 0 cos = 1 sin = 0 3m
U a Va Ub Vb
1 0 1 0
k ab = 40000 0 0 0 0
1 0 1 0
0 0 0 0
CE222: Structural Analysis/Dr Durgesh Rai/2008/20
ASSEMBLAGE OF ELEMENT STIFFNESS MATRICES
y
For bc = 126.9
0
cos = .60 and sin = 0.8
c
Ub Vb Uc Vc
0.36 .48 .36 0.48
4m
k bc = 40,000 .48 .64 0.48 .64
.36 .48 0.36 .48
.48 .64 .48 0.64 x
a b
3m
For ac , = 90 , cos = 0 and sin = 1
0
U a Va Uc Vc
0 0 0 0
k ac = 40,000 0 1 0 1
0 0 0 0
0 1 0 1
CE222: Structural Analysis/Dr Durgesh Rai/2008/21
ASSEMBLAGE OF ELEMENT STIFFNESS MATRICES
Pxa 1 0 1 0 0 0 U a
P 0 1
ya 0 0 0 1 Va
Pxb 1 0 1.36 .48 .36 0.48 Ub
= 40,000
Pyb 0 0 .48 0.64 0.48 .64 Vb
0 0 .36 0.48 0.36 .48 Uc
Pxc
Pyc 0 1 0.48 .64 0.48 1.64 Vc
y
4m
x
a b
3m
CE222: Structural Analysis/Dr Durgesh Rai/2008/22
NODAL EQUILIBRIUM OF THE STRUCTURE
2 6
Computer Implementation 4 c
1 1 5
ID Table
3
a b 3 2
No. of rows = no. of global DOF 4
No. of col. = no. of elements a 1
3
2 1 b
Element number
1 2 3
Element number
1 2 3
1 1 0 1
1 1 3 1
Global 2 2 0 2
Dof
3 3 1 0 Element 2 2 4 2
DOF
4 4 2 0 3 3 5 5
5 0 3 3
4 4 6 6
6 0 4 4
More capacity Used in CAL 89
CE222: Structural Analysis/Dr Durgesh Rai/2008/23
BOUNDARY CONDITIONS
Truss element k = p = k u
Nodal equilibrium K = P = K U Pxa Pxb Pxc
Displn. can only be solved if system is
constrained against rigid body motion. a k1 b k2 c
Pxa k1 k1 0 U a
For k2 U b
Pxb = k1 k1 + k 2
Pxc 0 k2 k2 Uc
K can not be inverted because both
P and U contain known and unknown quantities.
If an arbitrary force is applied which is not in equilibrium,
Pxa + Pxb + Pxc 0 , the elements will exhibit elastic deformation
+ rigid body translations.
P=KU
1 1
K P=K k U = 1 U but K 1 = 0 , no unique solution unless system is
constrained.
CE222: Structural Analysis/Dr Durgesh Rai/2008/24
BOUNDARY CONDITIONS
Now constrain node c.
k1 c
Pf Pxa k1 0 Ua
Uf
P = k k1 + k 2 k2 U b Pxa Pxb
xb 1
Ps Pxc 0 k2 k2 Uc = 0 U s
Pxa k1 k1 U a
P = k k1 + k2 Ub
xb 1
1 1 1 1 1 Pxb
+ + P
xa +
U a 1 k1 + k 2 k1 Pxa k1 k2 k2 Pxa k k k
U = = = 1 2 2
b k1 k2 k1 k1 Pxb 1 1 Pxb P + P
xa xb
k k2
2 k2
Therefore,
U a
Pxc = [ 0 k2 ] = Pxa Pxb
U b
CE222: Structural Analysis/Dr Durgesh Rai/2008/25
BOUNDARY CONDITIONS
If U f are free-dof and Us are support nodes, then
Pf K ff K fs U f
= =
Ps K sf K ss Us
Pf = K ff U f + K fs Us
Ps = K sf U f + K ss Us
If U s = 0 Pf = K 1 ff Pf
Ps = K sf U f
CE222: Structural Analysis/Dr Durgesh Rai/2008/26
BOUNDARY CONDITIONS
You might have to re arrange equilibrium
equations, before partitioning is performed.
Pxa k1 k1 0 U a Pxb
P = k k1 + k 2 k2 U b
xb 1 c
Pxc 0 k2 k2 Uc
Pxb k1 k1 + k 2 k2 U a
P = k k1 0 U b
xa 1
Pxc 0 k2 k2 Uc
Kff Kfs
Pf Pxb k1 + k2 k1 k2 Ub Uf
P = k k 0 U
Ps xa 1 1 a U
Pxc k2 0 k2 Uc s
Ksf Kss
CE222: Structural Analysis/Dr Durgesh Rai/2008/27
CO-ORDINATE TRANSFORMATION
x
Global co-ordinate system for descending nodal displn.
k
Local co-ordinate system for defining axial force in each
element. y
p 'x ' g = pxg cos + pyg sin
y g
p 'y ' g = pxg sin + pyg cos
Similarly, p 'x ' k = pxk cos + pyk sin
pyk , vk
p 'y ' k = pxk sin + pyk cos
x
pyk (vk) pxk, uk
p 'x ' g c s 0 0 pxg pxk, uk
k
p y'
p '
y'g s c 0 0
= yg y
p' 0 0 c s p
pyg (vg)
x'k
xk
0 0 s c pyk
p 'y ' k pxg (ug)
g
x
p' T u pxg (ug)
pxg (vg)
CE222: Structural Analysis/Dr Durgesh Rai/2008/28
CO-ORDINATE TRANSFORMATION
p' = T p T= co-ordinate Transformation matrix
u '
g c s 00 ug
v ' g = s c 0 0 vg
u' 0 0 c s uk
k
v 'k 0 0 s c vk
u' T u u' = T u
Since T is an orthogonal matrix, T T = T 1
T T
which means p = T p' and u = T u'
CE222: Structural Analysis/Dr Durgesh Rai/2008/29
CO-ORDINATE TRANSFORMATION
p' = T p T= co-ordinate Transformation matrix
u '
g c s 00 ug
v ' g = s c 0 0 vg
u' 0 0 c s uk
k
v 'k 0 0 s c vk
u' T u u' = T u
Since T is an orthogonal matrix, T T = T 1
T
T
which means p = T p ' and u = T u '
we have p = k u
T
T
T p' = k T u'
CE222: Structural Analysis/Dr Durgesh Rai/2008/30
CO-ORDINATE TRANSFORMATION
we have
p=ku
TT p' = k TT u'
Pre multiply both sides by T ,
T T T p ' = T k T T u
p = T k T T u
p = k u k = T k T T
k = Element Stiffness Matrix in Local Co-ord. (4x4)
k = Element Stiffness Matrix in Global Co-ord. (4x4)
CE222: Structural Analysis/Dr Durgesh Rai/2008/31
CO-ORDINATE TRANSFORMATION
0 c sc c s
2
c s 0 sc c 2 0 0 1 0 1 0
s c
0 sc
0 +s2 sc 2
s s c 0 0 AE 0 0 0 0
k = T k T =
T =
0 0 c s c 2 sc c2 sc 0 0 c s L 1 0 1 0
0 0 s c 2 0 0 s c 0 0 0 0
sc s 2 sc s
Similarly,
p = k ' u '
T p = k T u
T
T It is easier to store stiffness matrix
T T p = T k T u ' in local co-ordinate and then
transfer to get global.
p = T T k T u '
It is probably more efficient to
T implement directly and cheaper.
k = T k T
CE222: Structural Analysis/Dr Durgesh Rai/2008/32
ELEMENT FORCES
Solving Eqns. we get U And construct u for each element.
p = k u
ui
v
u=
i
and
ij
u = T u pxi 1 0 1 0 ui
p
u j 0
0 0 0 v j
AE
yi
v j pxj =
L 1 0 1 0 uj
ui pyj 0 0 0 v
v 0 j
AE i
=
pxj
L
[ 1 0 1 0 ]
u
j
v j
s ij = e u' ij e = Local element force transformation matrix
CE222: Structural Analysis/Dr Durgesh Rai/2008/33
ELEMENT FORCES
p = k u
=k Tu
1 0 1 0 c s 0 0 ui
AE 0 0 0 0 s c 0 0 vi
=
L 1 0 1 0 0 0 c s u j
0 0 0 0 0 0 s c v j
ui
v
AE i
ij
s =
L
[ cos sin cos sin ]
u j
v j
ij
=eu e= global element force deformation matrix.
CE222: Structural Analysis/Dr Durgesh Rai/2008/34
SUMMARY
k Element stiffness matrix in local co-ord
For all
elements
T
k = T k T Element stiffness matrix in Global co-ord
K Nodal Eqm.
Pf = K ff U f + K fs Us Apply boundary conditions
Pf = K ff U f Solve to get U f
ij Extract from U f
U
ij ij
s = eU Member forces for each elements
CE222: Structural Analysis/Dr Durgesh Rai/2008/35
Example
Pxb = 72 1.36 .48 0.48 1 0 .36
P = 90 .48 0.64
yb .64 0 0 0.48 c
Pyc = 0 .48 .64 1.64 0 1 .48
= 40,000
Pxa 1 0 0 0 0
0 0 1 0 1 0
Pya
P .36 .48 .48 0 0 .36
xc 90kN
U f = K ff1 Pf
a b 72kN
Ub 0.64 0.48 0 72 3.49
1 3
V =
25600
b 0.48 2.0 0.64 90
= 8.38 10 m Nodal Displacement
Vc 0 0.64 0.64 0 2.25
139.5
Ps = K sf U f = 90 kN Support Reactions
67.5
CE222: Structural Analysis/Dr Durgesh Rai/2008/36
Example
Element force in ab y
T=I
= {U a Va U b Vb } = 10 3 {0 0 3.49 8.381}
U ab c
0
0
(
s ab = e u ab = 40,000 [ 1 0 1 0 ] 10 3 ) = 139.5kN ( t )
3.49
4m
Element force in bc 8.381
= T ubc
ubc a b
x
3m
.6 .8 0 0 3.49 4.61
.8 .6 0 0 3 8.38 7.82
= (10 ) = (10 )
3
0 0 .6 .8 0 1.80
0 0 .8 .6 2.25 1.35 4.61
7.82
sbc = e u 'bc = 40,000 [ 1 0 1 0 ] 10 3( )
1.8
= 112.5 kN ( c )
1.35
CE222: Structural Analysis/Dr Durgesh Rai/2008/37
Example
Element force in ac y
uac = T u ac c
0 1 00 0.00 0.00
1 0 0 0 3 0.00 0.00
= (10 ) = (10 )
3
0 0 0 1 0.00 2.25 4m
0 0 1 0 2.25 0.00
0.00 x
0.00 a b
(
s ac = e u 'ac = 40,000 [ 1 0 1 0 ] 10 3 )
0.00
= 90.0 kN ( t ) 3m
2.25
CE222: Structural Analysis/Dr Durgesh Rai/2008/38
