Victor Saouma - Matrix Structural Analysis PDF
Victor Saouma - Matrix Structural Analysis PDF
Victor Saouma - Matrix Structural Analysis PDF
DRAFT
Lecture Notes in:
CVEN4525/5525
c VICTOR E. SAOUMA,
Fall 1999
Blank page
Contents
1 INTRODUCTION 1{1
1.1 Why Matrix Structural Analysis? . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{1
1.2 Overview of Structural Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{2
1.3 Structural Idealization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{4
1.3.1 Structural Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{5
1.3.2 Coordinate Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{6
1.3.3 Sign Convention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{6
1.4 Degrees of Freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{9
1.5 Course Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{11
List of Figures
1.1 Global Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{7
1.2 Local Coordinate Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{8
1.3 Sign Convention, Design and Analysis . . . . . . . . . . . . . . . . . . . . . . . . 1{9
1.4 Total Degrees of Freedom for various Type of Elements . . . . . . . . . . . . . . 1{10
1.5 Independent Displacements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{11
1.6 Examples of Global Degrees of Freedom . . . . . . . . . . . . . . . . . . . . . . . 1{13
1.7 Organization of the Course . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{14
2.1 Example for Flexibility Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2{3
2.2 Denition of Element Stiness Coecients . . . . . . . . . . . . . . . . . . . . . . 2{5
2.3 Stiness Coecients for One Dimensional Elements . . . . . . . . . . . . . . . . . 2{6
2.4 Flexural Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2{7
2.5 Torsion Rotation Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2{9
2.6 Deformation of an Innitesimal Element Due to Shear . . . . . . . . . . . . . . . 2{11
2.7 Eect of Flexure and Shear Deformation on Translation at One End . . . . . . . 2{13
2.8 Eect of Flexure and Shear Deformation on Rotation at One End . . . . . . . . . 2{14
3.1 Problem with 2 Global d.o.f. 1 and 2 . . . . . . . . . . . . . . . . . . . . . . . . 3{3
3.2 *Frame Example (correct K23 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3{6
3.3 Grid Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3{10
4.1 Arbitrary 3D Vector Transformation . . . . . . . . . . . . . . . . . . . . . . . . . 4{2
4.2 3D Vector Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{4
4.3 2D Frame Element Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{7
4.4 Grid Element Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{7
4.5 2D Truss Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{8
4.6 Simple 3D Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{10
4.7 Arbitrary 3D Rotation; Rotation with respect to . . . . . . . . . . . . . . . . . 4{10
4.8 Arbitrary 3D Rotation; Rotation with respect to
. . . . . . . . . . . . . . . . . 4{11
4.9 Special Case of 3D Transformation for Vertical Members . . . . . . . . . . . . . . 4{12
4.10 Arbitrary 3D Rotation; Rotation with respect to . . . . . . . . . . . . . . . . . 4{13
4.11 Rotation of Cross-Section by . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4{14
4.12 Arbitrary 3D Element Transformation . . . . . . . . . . . . . . . . . . . . . . . . 4{15
Draft
0{2 LIST OF FIGURES
5.1 Example for [ID] Matrix Determination . . . . . . . . . . . . . . . . . . . . . . . 5{5
5.2 Flowchart for Assembling Global Stiness Matrix . . . . . . . . . . . . . . . . . . 5{5
5.3 Example of Bandwidth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{7
5.4 Numbering Schemes for Simple Structure . . . . . . . . . . . . . . . . . . . . . . 5{7
5.5 Beam Element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{10
5.6 ID Values for Simple Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{20
5.7 Simple Frame Anlysed with the MATLAB Code . . . . . . . . . . . . . . . . . . 5{20
5.8 Simple Frame Anlysed with the MATLAB Code . . . . . . . . . . . . . . . . . . 5{22
5.9 Program Flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{27
5.10 Program's Tree Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{28
5.11 Flowchart for the Skyline Height Determination . . . . . . . . . . . . . . . . . . . 5{30
5.12 Flowchart for the Global Stiness Matrix Assembly . . . . . . . . . . . . . . . . . 5{31
5.13 Flowchart for the Load Vector Assembly . . . . . . . . . . . . . . . . . . . . . . . 5{33
5.14 Flowchart for the Internal Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{34
5.15 Flowchart for the Reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{35
5.16 Structure Plotted with CASAP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5{56
6.1 Example of [B] Matrix for a Statically Determinate Truss . . . . . . . . . . . . . 6{3
6.2 Example of [B] Matrix for a Statically Determinate Beam . . . . . . . . . . . . . 6{5
6.3 Example of [B] Matrix for a Statically Indeterminate Truss . . . . . . . . . . . . 6{7
6.4 *Examples of Kinematic Instability . . . . . . . . . . . . . . . . . . . . . . . . . . 6{13
6.5 Example 1, Congruent Transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6{19
6.6 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6{21
7.1 Stable and Statically Determinate Element . . . . . . . . . . . . . . . . . . . . . 7{6
7.2 Example 1, [k] ! [d] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7{7
8.1 Stress Components on an Innitesimal Element . . . . . . . . . . . . . . . . . . . 8{1
8.2 Stress Traction Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8{2
8.3 Equilibrium of Stresses, Cartesian Coordinates . . . . . . . . . . . . . . . . . . . 8{4
8.4 Fundamental Equations in Solid Mechanics . . . . . . . . . . . . . . . . . . . . . 8{8
9.1 Variational and Dierential Operators . . . . . . . . . . . . . . . . . . . . . . . . 9{2
9.2 *Strain Energy and Complementary Strain Energy . . . . . . . . . . . . . . . . . 9{9
9.3 Eects of Load Histories on U and Wi . . . . . . . . . . . . . . . . . . . . . . . . 9{11
9.4 Torsion Rotation Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9{14
9.5 Flexural Member . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9{15
9.6 Tapered Cantilivered Beam Analysed by the Vitual Displacement Method . . . . 9{23
9.7 Tapered Cantilevered Beam Analysed by the Virtual Force Method . . . . . . . . 9{28
9.8 Three Hinge Semi-Circular Arch . . . . . . . . . . . . . . . . . . . . . . . . . . . 9{30
9.9 Semi-Circular Cantilevered Box Girder . . . . . . . . . . . . . . . . . . . . . . . . 9{31
9.10 Single DOF Example for Potential Energy . . . . . . . . . . . . . . . . . . . . . . 9{34
9.11 Graphical Representation of the Potential Energy . . . . . . . . . . . . . . . . . . 9{35
List of Tables
1.1 Example of Nodal Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{5
1.2 Example of Element Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{5
1.3 Example of Group Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1{6
1.4 Degrees of Freedom of Dierent Structure Types Systems . . . . . . . . . . . . . 1{12
2.1 Examples of In
uence Coecients . . . . . . . . . . . . . . . . . . . . . . . . . . 2{2
6.1 Internal Element Force Denition for the Statics Matrix . . . . . . . . . . . . . . 6{2
6.2 Conditions for Static Determinacy, and Kinematic Instability . . . . . . . . . . . 6{13
9.1 Essential and Natural Boundary Conditions . . . . . . . . . . . . . . . . . . . . . 9{5
9.2 Possible Combinations of Real and Hypothetical Formulations . . . . . . . . . . . 9{20
9.3 Comparison of 2 Alternative Approximate Solutions . . . . . . . . . . . . . . . . 9{49
9.4 Summary of Variational Terms Associated with One Dimensional Elements . . . 9{52
10.1 Characteristics of Beam Element Shape Functions . . . . . . . . . . . . . . . . . 10{6
10.2 Interpretation of Shape Functions in Terms of Polynomial Series (1D & 2D) . . . 10{12
10.3 Polynomial Terms in Various Element Formulations (1D & 2D) . . . . . . . . . . 10{12
Draft
0{2 LIST OF TABLES
Volume of body
Chapter 1
INTRODUCTION
1.1 Why Matrix Structural Analysis?
1 In most Civil engineering curriculum, students are required to take courses in: Statics,
Strength of Materials, Basic Structural Analysis. This last course is a fundamental one which
introduces basic structural analysis (determination of reactions, de
ections, and internal forces)
of both statically determinate and indeterminate structures.
2 Also Energy methods are introduced, and most if not all examples are two dimensional. Since
the emphasis is on hand solution, very seldom are three dimensional structures analyzed. The
methods covered, for the most part lend themselves for \back of the envelope" solutions and
not necessarily for computer implementation.
3 Those students who want to pursue a specialization in structural engineering/mechanics, do
take more advanced courses such as Matrix Structural Analysis and/or Finite Element Analysis.
4 Matrix Structural Analysis, or Advanced Structural Analysis, or Introduction to Structural
Engineering Finite Element, builds on the introductory analysis course to focus on those meth-
ods which lend themselves to computer implementation. In doing so, we will place equal
emphasis on both two and three dimensional structures, and develop a thorough understanding
of computer aided analysis of structures.
5 This is essential, as in practice most, if not all, structural analysis are done by the computer
and it is imperative that as structural engineers you understand what is inside those \black
boxes", develop enough self assurance to be capable of opening them and modify them to
perform certain specic tasks, and most importantly to understand their limitations.
6 With the recently placed emphasis on the nite element method in most graduate schools,
many students have been tempted to skip a course such as this one and rush into a nite element
one. Hence it is important that you understand the connection and role of those two courses.
The Finite Element Method addresses the analysis of two or three dimensional continuum. As
such, the primary unknowns is u the nodal displacements, and internal \forces" are usually
Draft
1{2 INTRODUCTION
restricted to stress . The only analogous one dimensional structure is the truss.
7 Whereas two and three dimensional continuum are essential in civil engineering to model
structures such as dams, shells, and foundation, the majority of Civil engineering structures
are constituted by \rod" one-dimensional elements such as beams, girders, or columns. For
those elements, \displacements" and internal \forces" are somehow more complex than those
encountered in continuum nite elements.
8 Hence, contrarily to continuum nite element where displacement is mostly synonymous with
translation, in one dimensional elements, and depending on the type of structure, generalized
displacements may include translation, and/or
exural and/or torsional rotation. Similarly,
\internal forces" are not stresses, but rather axial and shear forces, and/or
exural or torsional
moments. Those concepts are far more relevant in the analysis/design of most civil engineering
structures.
9 Hence, Matrix Structural Analysis, is truly a bridge course between introductory analysis
and nite element courses. The element stiness matrix [k] will rst be derived using methods
introduced in basic structural analysis, and later using energy based concepts. This later
approach is the one exclusively used in the nite element method.
10 An important component of this course is computer programing. Once the theory and the
algorithms are thoroughly explained, you will be expected to program them in either Fortran
(preferably 90) or C (sorry, but no Basic) on the computer of your choice. The program
(typically about 3,500 lines) will perform the analysis of 2 and 3 dimensional truss and frame
structures, and many students have subsequently used it in their professional activities.
11 There will be one computer assignment in which you will be expected to perform sim-
ple symbolic manipulations using Mathematica. For those of you unfamiliar with the Bechtel
Laboratory, there will be a special session to introduce you to the operation of Unix on Sun
workstations.
19 Group number will then dene both element type, and elastic/geometric properties. The
last one is a pointer to a separate array, Table 1.3. In this example element 1 has element code
1 (such as beam element), while element 2 has a code 2 (such as a truss element). Material
group 1 would have dierent elastic/geometric properties than material group 2.
20 From the analysis, we rst obtain the nodal displacements, and then the element internal
forces. Those internal forces vary according to the element type. For a two dimensional frame,
those are the axial and shear forces, and moment at each node.
21 Hence, the need to dene two coordinate systems (one for the entire structure, and one for
28 Fig. 1.4 also shows the geometric (upper left) and elastic material (upper right) properties
associated with each type of element.
34 This table shows the degree of freedoms and the corresponding generalized forces.
35 We should distinguish between local and global d.o.f.'s. The numbering scheme follows the
following simple rules:
Local: d.o.f. for a given element: Start with the rst node, number the local d.o.f. in the same
order as the subscripts of the relevant local coordinate system, and repeat for the second
node.
Global: d.o.f. for the entire structure: Starting with the 1st node, number all the unrestrained
global d.o.f.'s, and then move to the next one until all global d.o.f have been numbered,
Fig. 1.6.
Part I
Matrix Structural Analysis of
Framed Structures
Draft
Draft
Chapter 2
ELEMENT STIFFNESS MATRIX
2.1 Introduction
1 In this chapter, we shall derive the element stiness matrix [k] of various one dimensional
elements. Only after this important step is well understood, we could expand the theory and
introduce the structure stiness matrix [K] in its global coordinate system.
2 As will be seen later, there are two fundamentally dierent approaches to derive the stiness
matrix of one dimensional element. The rst one, which will be used in this chapter, is based on
classical methods of structural analysis (such as moment area or virtual force method). Thus,
in deriving the element stiness matrix, we will be reviewing concepts earlier seen.
3 The other approach, based on energy consideration through the use of assumed shape func-
tions, will be examined in chapter 11. This second approach, exclusively used in the nite
element method, will also be extended to two and three dimensional continuum elements.
11 Hence kij will correspond to the reaction at dof i due to a unit deformation (translation or
rotation) at dof j , Fig. 2.2.
12 The actual stiness coecients are shown in Fig. 2.3 for truss, beam, and grid elements in
terms of elastic and geometric properties.
13 In the next sections, we shall derive those stiness coecients.
16We start from the dierential equation of a beam, Fig. 2.4 in which we have all positive
known displacements, we have from strength of materials
d2 v = M ; V x
M = ;EI dx (2.11)
2 1 1
19 Applying the boundary conditions at x = L and combining with the expressions for C1 and
C2 ) (
v 0 = 2 ) ;;EIv
EI2 = M1 L ; 21 V1 L2 ; EI1 (2.15)
v = v2 2 = 12 M1 L2 ; 16 V1 L3 ; EI1 L ; EIv1
V1 = 6EI z ( + ) + 12EIz (v ; v )
L 1 2
2 L3 1 2 (2.21)
V2 = ; 6EI z 12EIz
L2 (1 + 2 ) ; L3 (v1 ; v2) (2.22)
k = 0:3 (2.27-b)
2
1 + db
For other sections, J is often tabulated.
27 Note that J corresponds to Ix where x is the axis along the element.
28 Having developed a relation between torsion and shear stress, we now seek a relation between
torsion and torsional rotation. In Fig. 2.5-b, we consider the arc length BD
9 9
max dx = dc >
= >
= d = T
) ddx =
max d = max (2.28)
c
max > dx Gc > dx GJ
max = G
;
max = TC
J
;
Z Z
29 Finally, we can rewrite this last equation as Tdx = Gjd and obtain:
T = GJ
L
(2.29)
τ
γ dvs
γ=
y dx
τ τ
γ
τ
x
to be determined, I is the moment of inertia of the cross sectional area about the neutral axis,
and b is the width of the rectangular beam.
33 The preceding equation can be simplied as
= AV (2.32)
s
where As is the eective cross section for shear (which is the ratio of the cross sectional area
to the area shear factor)
34 Let us derive the expression of As for rectangular sections. The exact expression for the
shear stress is
= VIbQ (2.33)
where Q is the moment of the area from the external bers to y with respect to the neutral
axis; For a rectangular section, this yields
= VIbQ (2.34-a)
!
V Z h=2
V h 2
= Ib by0dy0 = 2I 4 ; y2 (2.34-b)
y
!
6 V h 2
= bh3 4 ; 4 (2.34-c)
and we observe that the shear stress is zero for y = h=2 and maximum at the neutral axis
where it is equal to 1:5 bhV .
dvs = V (2.36)
dx GAs
Assuming V to be constant, we integrate
V x+C
vs = GA (2.37)
1
s
36If the displacement vs is zero at the opposite end of the beam, then we solve for C1 and
obtain
V (x ; L)
vs = ; GA (2.38)
s
37 We dene
def
= GA12EI (2.39)
sL2
2
= 24(1 + ) AA Lr (2.40)
s
Hence for small slenderness ratio Lr compared to unity, we can neglect .
38 Next, we shall consider the eect of shear deformations on both translations and rotations
Eect on Translation Due to a unit vertical translation, the end shear force is obtained from
Eq. 2.21 and setting v1 = 1 and 1 = 2 = v2 = 0, or V = 12LEI3 z . At x = 0 we have, Fig.
2.7 VL 9
vs = GA s >
=
12 EI
V = L3 > vs =
z (2.41)
= GA 12EI2 ;
sL
Hence, the shear deformation has increased the total translation from 1 to 1+. Similar
arguments apply to the translation at the other end.
1
6EI
L2
12EI
L3
Figure 2.7: Eect of Flexure and Shear Deformation on Translation at One End
Eect on Rotation Considering the beam shown in Fig. 2.8, even when a rotation 1 is
applied, an internal shear force is induced, and this in turn is going to give rise to shear
deformations (translation) which must be accounted for. The shear force is obtained from
Eq. 2.21 and setting 1 = 1 and 2 = v1 = v2 = 0, or V = 6EI L2 . As before
z
9
vs = VL >
GAs =
V = 6EIz vs = 0:5L (2.42)
L2EI
12 >
= GAs L2
;
in other words, the shear deformation has moved the end of the beam (which was supposed
to have zero translation) by 0:5L.
0.5βLθ1
0.5βLθ1
Figure 2.8: Eect of Flexure and Shear Deformation on Rotation at One End
2 u1 u2 3
[kt ] = AE
p 14 1 ; 15 (2.43)
L p 2 ; 1 1
2
v1 1 v2 2 3
V166 L312(1+EIz
y)
6EIz
L2 (1+y ) ; L312(1+
EIz
y)
6EIz 7
L2 (1+y ) 7
6 6 EI
M166 L2 (1+y )z (4+y )EIz ; L2(1+y )
6 EIz (2;y )EIz 7
(1+y )L L(1+y ) 7 (2.49)
[kbV ] = V2666 ; L312(1+
EIz ; L26(1+
EIz 12EIz ; L26(1+
EIz 7
7
y) y) L3 (1+y ) y) 7
EIz
M264 L26(1+ (2;y )EIz ; L26(1+
EIz (4+y )EIz 7 7
y) L(1+y ) y) L(1+y ) 5
2 u1 v1 1 u2 v2 2 3
P16 EA
L 0 0 ; EA
L 0 0 7
V1666 0 12EIz
L3
6EIz
L2 0 ; 12LEI3 z 6EIz
L2
7
7
7
M 1 6 0 6EI z 4EI
L
z 0 ; 6EI z 2EI
L
z 7 (2.51)
[k2dfr ] = P2666 ; EA L
0
2
0 EA L2
0 0
7
7
L L 7
V266 0 ; L3 ; EI
6 12 EI z 6
L2
z 0 12EIz
L3 ;6EI 77
4EILz 7
2
M24 0 6 EI
L2
z 2 EI
L
z 0 ; 6EI
L2
z
L 5
48Note that if shear deformations must be accounted for, the entries corresponding to shear
and
exure must be modied in accordance with Eq. 2.49
2 1 1 v1 2 2 v2 3
T166 GIx
L 0 0 ; GILx 0 0 7
M16 06 4 EI
L
z ; EI
6
L2
z 0 2EIz
L
6EIz
L2
7
7
7
V 0 ; L2
6 6 EI z 12 EI z 0 ; 6EI z ; 12LEI3 z 7 (2.53)
[kg ] = T21666 ; Gix 0
L3
0 GIx L2
0 0
7
7
L L 7
M266 0
6 2 EI
L
z ; EI
6
LEI
2
z 0 4EIz
L
6EIz 77
12LEIz 7
2
V24 0 6 EI
L2
z ; L3 z
12 0 6EI
L2
z
L3 5
50Note that if shear deformations must be accounted for, the entries corresponding to shear
and
exure must be modied in accordance with Eq. 2.49
2.6.5 3D Frame Element
2
u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z 2 3
Px1 t
k11 0 0 0 0 t
0 k21 0 0 0 0 0
Vy1666 0 b
k11 0 0 0 b b
k12 0 k13 0 0 0 kb14 777
V z16 0 0 b
k11 0
g
;k12b 0 0 0 b
k13 0
g
;k14b 0 77
Tx1666 0 0 0 k11 0 0 0 0 0 k12 0 0 77
My166 0 0 b
k32 0 b
k22 0 0 0 b
k12 0 b
k24 0 77
[k3dfr ] = M z16
6
0 b
k21 0 0 0 b
k22 0 ;k12b 0 0 0 b
k24 7
7
t
Px2666 k21 0 0 0 0 0 t
k22 0 0 0 0 0 7
7
b b b b 7
Vy266 0 k13 0 0 0 k14 0 k33 0 0 0 k34 7
7
V z266 0 0 b
k13 0
g
;k14b 0 0 0 b
k33 0
g
;k34b 0 7
7
Tx2666 0 0 0 k12 0 0 0 0 0 k22 0 0 7
7
My24 0 0 b
k12 0 b
k24 0 0 0 ;k43b 0 b
k44 0 7
5
Mz2 0 ;k12b 0 0 0 b
k24 0 b
k43 0 0 0 b
k44
(2.54)
Victor Saouma Matrix Structural Analysis
Draft
2.7 Remarks on Element Stiness Matrices 2{19
For [k311D ] and with we obtain:
2 u1 v1 w1 x1 y 1 z 1 u2 v2 w2 x2 y 2 z2 3
Px1 EA l 0 0 0 0 0 ; EA
L 0 0 0 0 0
Vy16 0 12EIz 0 0 0 6EIz 0 ; 12EIz 0 0 0 6EIz 7
6 L 3
12EIy
L 2 L 3 L 2 7
Vz16
6
0 0 L3 0 ; 6EI
L 2
y 0 0 0 ; 12LEI3 y 0 ;6 EI
L
y
2 0 7 7
Tx16 0 0 0 GIx 0 0 0 0 0 ; GILx 0 0 7
6 6EIy L 4EIy 6EIy 2EIy
My61 0 0 ; L2 0 L 0 0 0 L2 0 L 0 7 7
6
M1 0 6 EI z 0 0 0 4 EI z 0 ; L2
6 EI z 0 0 0 2EI z 7
[k3dfr ] = P z6 L 2 L L 7
x26 ; EA
L 0
EI
0 0 0 0
EI
EA
L 0
EI
0 0 0 0 7
EI
6
Vy26 0 ; L312 z 0 0 0 ; L2
6 z 0 12 z 0 0 0 ; L2 7
6 z
12EIy EIy L3 12EIy 6EIy
7
Vz26
6 0 0 ; L3 0 6 L 2 0 0 0 L 3 0 L 2 0 7 7
Tx26 0 0 0 ; GILx 0 0 0 0 0 GIx
L 0 0 7
6 6EIy 2EIy 6EIy 4EIy
My42 0 0 ; L2 0 L 0 0 0 L2 0 L 0 7 5
EI z EI z EI z EI
Mz2 0 ; L2 Lz
6 0 0 0 2 0 6 0 0 0 4
L2 L
(2.55)
51Note that if shear deformations must be accounted for, the entries corresponding to shear
and
exure must be modied in accordance with Eq. 2.49
Chapter 3
STIFFNESS METHOD; Part I:
ORTHOGONAL STRUCTURES
3.1 Introduction
1 In the previous chapter we have rst derived displacement force relations for dierent types
of rod elements, and then used those relations to dene element stiness matrices in local
coordinates.
2 In this chapter, we seek to perform similar operations, but for an orthogonal structure in
global coordinates.
3 In the previous chapter our starting point was basic displacement-force relations resulting in
element stiness matrices [k].
4 In this chapter, our starting point are those same element stiness matrices [k], and our
objective is to determine the structure stiness matrix [K], which when inverted, would yield
the nodal displacements.
5 The element stiness matrices were derived for fully restrained elements.
6 This chapter will be restricted to orthogonal structures, and generalization will be discussed
later.
The stiness matrices will be restricted to the unrestrained degrees of freedom.
7 From these examples, the interrelationships between structure stiness matrix, nodal dis-
placements, and xed end actions will become apparent. Then the method will be generalized
in chapter 5 to describe an algorithm which can automate the assembly of the structure global
stiness matrix in terms of the one of its individual elements.
Draft
3{2 STIFFNESS METHOD; Part I: ORTHOGONAL STRUCTURES
3.2 The Stiness Method
8 As a \vehicle" for the introduction to the stiness method let us consider the problem in Fig
3.1-a, and recognize that there are only two unknown displacements, or more precisely, two
global d.o.f: 1 and 2 .
9 If we were to analyse this problem by the force (or
exibility) method, then
1. We make the structure statically determinate by removing arbitrarily two reactions (as
long as the structure remains stable), and the beam is now statically determinate.
2. Assuming that we remove the two roller supports, then we determine the corresponding
de
ections due to the actual laod (B and C ).
3. Apply a unit load at point B , and then C , and compute the de
ections fij at note i due
to a unit force at node j .
4. Write the compatibility of displacement equation
" #( ) ( ) ( )
fBB fBC RB ; 1 = 0 (3.1)
fCB fCC RC 2 0
5. Invert the matrix, and solve for the reactions
10 We will analyze this simple problem by the stiness method.
1. The rst step consists in making it kinematically determinate (as opposed to statically
determinate in the
exibility method). Kinematically determinate in this case simply
means restraining all the d.o.f. and thus prevent joint rotation, Fig 3.1-b.
2. We then determine the xed end actions caused by the element load, and sum them for
each d.o.f., Fig 3.1-c: FEM1 and FEM2 .
3. In the third step, we will apply a unit displacement (rotation in this case) at each degree
of freedom at a time, and in each case we shall determine the reaction forces, K11 , K21 ,
and K12 , K22 respectively. Note that we use [K], rather than k since those are forces
in the global coordinate system, Fig 3.1-d. Again note that we are focusing only on the
reaction forces corresponding to a global degree of freedom. Hence, we are not attempting
to determine the reaction at node A.
4. Finally, we write the equation of equilibrium at each node:
( ) ( ) " #( )
M1 = FEM1 + K 11 K12 1 (3.2)
M2 FEM2 K21 K22 2
3.3 Examples
11. For element BC, the local and global coordinates do not match, hence we will need to
transform the displacements from their global to their local coordinate components. But since,
vector (displacement and load), and matrix transformation have not yet been covered, we note
by inspection that the relationship between global and local coordinates for element BC is
Local 1 2 3 4 5 6
Global 0 0 0 2 ;1 3
and we observe that there are no local or global displacements associated with dof 1-3; Hence
Mathematica:
Ic=Ib
M= 0
w= 0
H= L
alpha=
K={
{E A /L + 12 E Ic /H^3, 0, 6 E Ic/H^2},
{0, 12 E Ib/L^3 + E A/H, -6 E Ib/L^2},
{6 E Ic/H^2, -6 E Ib/L^2, 4 E Ib/L + 4 E Ic/H}
}
d=Inverse[K]
load={-P/2 - w H/2, -P/2, M+P L/8 -w H^2/12}
displacement=Simplify[d . load]
6. Determine the element internal forces. This will be accomplished by multiplying each element
stiness matrix [k] with the vector of nodal displacement f g. Note these operations should
be accomplished in local coordinate system, and great care should be exercized in writing the
nodal displacements in the same local coordinate system as the one used for the derivation of
the element stiness matrix, Eq. 2.53.
7. For element AB and BC, the vector of nodal displacements are
8
>
> 1 9
>
>
8
>
> 0 9
>
>
8
>
> ;3 9
>
>
>
>
>
>
>
2 >
>
>
>
>
>
>
>
>
>
0 >
>
>
>
>
>
>
>
>
>
;1 >
>
>
>
>
<
3 =
=>
<
0 =
= > 02
< =
(3.26)
>
>
>
4 >
>
> >
> 1 >
>
> >
>
>
>
>
>
>
>
>
:
5 >
>
>
>
;
>
>
>
>
:
;3 >
>
>
>
;
>
>
>
>
:
0 >
>
>
>
;
6 2 0
| {z } | {z }
AB BC
8. Hence, for element AB we have
8
> p1 9
>
2 GIx
L 0 0 ; GILx 0 0 38
0 9
> >
0 4EIy ; 6EI y 0 2EIy 6EIy 7 > > >
>
>
>
>
> p2 >
>
>
>
6
6 L L2 L L 2 7>
>
>
> 0 >
>
>
>
>
<
p3 >
= 6
6
= ; Gix
6 0 ; 6EI
L2
y 12EI
L3
y 0 ; 6EI
L2
y ; 12EIy 7
L 73 7
>
<
0
>
=
= (3.27)
> p4 > 6 0 0 GIx 0 0 7> 1 >
>
> >
> 6 L L 7>> >
>
>
>
>
>
:
p5 >
>
>
>
;
6
4 0 2EIy
L ; 6EI
LEI
2
y 0 4EIy
L
6EIy 7
L 2 5
>
>
>
>
:
;3 >
>
>
>
;
p6 0 6EI
L2
y 12
; L3 y 0 6EI
L2
y 12EI
L3
y 2
3.4 Observations
13 On the basis of these two illustrative examples we note that the global structure equilibrium
equation can be written as
fPg = fFEAg + [K]fg (3.29)
where [K] is the global structure stiness matrix (in terms of the unrestrained d.o.f.) fPg the
vector containing both the nodal load and the nodal equivalent load caused by element loading,
fg is the vector of generalized nodal displacements.
14 Whereas the preceding two examples were quite simple to analyze, we seek to generalize the
method to handle any arbitrary structure. As such, some of the questions which arise are:
1. How do we determine the element stiness matrix in global coordinate systems, [Ke ],
from the element stiness matrix in local coordinate system [ke ]?
2. How to assemble the structure [KS ] from each element [KE ]?
3. How to determine the fFEAg or the nodal equivalent load for an element load?
4. How to determine the local nodal displacements from the global ones?
5. How do we compute reactions in the restrained d.o.f?
6. How can we determine the internal element forces (P , V , M , and T )?
7. How do we account for temperature, initial displacements or prestrain?
Those questions, and others, will be addressed in the next chapters which will outline the
general algorithm for the direct stiness method.
Chapter 4
TRANSFORMATION MATRICES
4.1 Derivations
4.1.1 [ke ] [Ke] Relation
1 In the previous chapter, in which we focused on orthogonal structures, the assembly of the
structure's stiness matrix [Ke ] in terms of the element stiness matrices was relatively straight-
forward.
2 The determination of the element stiness matrix in global coordinates, from the element
stiness matrix in local coordinates requires the introduction of a transformation.
3 This chapter will examine the 2D and 3D transformations required to obtain an element
stiness matrix in global coordinate system prior to assembly (as discussed in the next chapter).
4 Recalling that
premultiplying by [;];1
fPg = [;];1 [ke][;]fg (4.6)
7 But since the rotation matrix is orthogonal, we have [;];1 = [;]T and
fPg = [|;]T [{zke][;}]fg (4.7)
[Ke]
which is the general relationship between element stiness matrix in local and global coordi-
nates.
4.1.2 Direction Cosines
8 The problem confronting us is the general transfoormation of a vector V from (X; Y; Z)
coordinate system to (X; Y; Z), Fig. 4.1: where:
8 9 2 38 9
>
< Vx >
= lxX lxY lxZ >
< VX >
=
>
Vy >
= 4 lyX lyY lyZ
6 7
5
>
VY >
(4.9)
: Vz ; lzX lzY lzZ : VZ ;
| {z }
[
]
13 If we use indecis instead of cartesian system, then direction cosines can be expressed as
Vx = VX l11 + VY l12 + VZ l13 (4.16)
or by extension: 8 9 2 38 9
>
< Vx >
= l11 l12 l13 >
< VX >
=
>
Vy >
= 64 l21 l22 l23 7
5
>
VY >
(4.17)
: Vz ; l31 l32 l33 : VZ ;
| {z }
[
]
Alternatively, [
] is the matrix whose columns are the direction cosines of x; y; z with respect
to X; Y; Z : 2 3
l11 l21 l31
[
]T = 64 l12 l22 l32 7
5 (4.18)
l13 l23 l33
The transformation of V can be written as:
fvg = [
] fVg (4.19)
where: fvg is the rotated coordinate system and fVg is in the original one.
Victor Saouma Matrix Structural Analysis
Draft
4.1 Derivations 4{5
14 Direction cosines are unit orthogonal vectors
3
X
lij lij = 1 i = 1; 2; 3 (4.20)
j =1
i.e:
2 + l2 + l2 = 1
l11 (4.21)
12 13
cos2 + cos2 + cos2
= 1 = 11 (4.22)
and
8
3
X
>
< i = 1; 2; 3
lij lkj = 0 k = 1; 2; 3 (4.23)
j =1 >
: i 6= k
l11 l21 + l12 l22 + l13 l23 = 0 = 12 (4.24)
4.2.1 2 D cases
4.2.1.1 2D Frame, and Grid Element
20 The vector rotation matrix [
] is identical for both 2D frame and grid elements, Fig. 4.3,
and 4.4 respectively.
21 From Eq. 4.10 the vector rotation matrix is dened in terms of 9 direction cosines of
9 dierent angles. However for the 2D case, we will note that four angles are interrelated
(lxX ; lxY ; lyX ; lyY ) and can all be expressed in terms of a single one , where is the direction
of the local x axis (along the member from the rst to the second node) with respect to the
global X axis. The remaining 5 terms are related to another angle, , which is between the
Z axis and the x-y plane. This angle is zero because we select an orthogonal right handed
coordinate system. Thus, the rotation matrix can be written as:
2 3 2 3 2 3
l xX lxY lxZ cos cos( 2 ; ) 0 cos sin 0
[
] = 64 lyX lyY lyZ 75 = 64 cos( 2 + ) cos 0 75 = 64 ; sin cos 0 75 (4.28)
lzX lzY lzZ 0 0 1 0 0 1
and we observe that the angles are dened from the second subscript to the rst, and that
counterclockwise angles are positive.
Victor Saouma Matrix Structural Analysis
Draft
4.2 Transformation Matrices For Framework Elements 4{7
4.2.1.2 2D Truss
23For the 2D truss element, the global coordinate system is two dimensional, whereas the local
one is only one dimensional, hence the vector transformation matrix is, Fig. 4.5.
h i h i h i
[
] = lxX lxY = cos cos( 2 ; ) = cos sin (4.30)
24 The element rotation matrix [;] will then be assembled from the vector rotation matrix [
].
8 9
( ) " # " #>
> P1 >
>
= [
0 ] [
0 ] = cos0 sin0 cos0 sin0
> >
p1 <
P2 =
(4.31)
p2 >
> P3 >
>
| {z }>
:
P4
>
;
[;]
4.2.2 3D Frame
25GIven that rod elements, are dened in such a way to have their local x axis aligned with
their major axis, and that the element is dened by the two end nodes (of known coordinates),
Victor Saouma Matrix Structural Analysis
Draft
4.2 Transformation Matrices For Framework Elements 4{9
then recalling the denition of the direction cosines it should be apparent that the evaluation
of the rst row, only, is quite simple. However evaluation of the other two is more complex.
26 This generalized transformation from X; Y ; Z to x; y; z was accomplished in one step in the
two dimensional case, but intermediary ones will have to be dened in the 3D case.
27 Starting with a reference (X1 ; Y1 ; Z1 ) coordinate system which corresponds to the global
coordinate system, we can dene another one, X2 ; Y2 ; Z2 , such that X2 is aligned along the
element, Fig. ??.
28 In the 2D case this was accomplished through one single rotation , and all other angles
where dened in terms of it.
29 In the 3D case, it will take a minimum of two rotations and
, and possibly a third one
(dierent than the one in 2D) to achieve this transformation.
30 We can start with the rst row of the transformation matrix which corresponds to the
direction cosines of the reference axis (X1 ; Y1 ; Z1 ) with respect to X2 . This will dene the
rst row of the vector rotation matrix [
]:
2 3
CX CY CZ
[
] = 4 l21 l22 l23
6 7
5 (4.32)
l31 l32 l33
q
where CX = xj L;xi , CY = yj ;L yi , CZ = zj ;L zi , L = (xj ; xi )2 + (yj ; yi )2 + (zj ; zi )2 .
31 Note that this does not uniquely dene the new coordinate system. This will be done in two
ways: a simple and a general one.
4.2.2.1 Simple 3D Case
32 We start by looking at a simplied case, Fig. 4.6, one in which Z2 is assumed to be horizontal
in the X1 ; Z1 plane, this will also dene Y2 . We note that there will be no ambiguity unless
the member is vertical.
33 This transformation can be used if:
1. The principal axis of the cross section lie in the horizontal and vertical plane (i.e the web
of an I Beam in the vertical plane).
2. If the member has 2 axis of symmetry in the cross section and same moment of inertia
about each one of them (i.e circular or square cross section).
34The last two rows of Eq. 4.32 can be determined through two successive rotations (assuming
that (X1 ; Y1 ; Z1 and X2 ; Y2 ; Z2 are originally coincident):
1. Rotation by about the Y1 axis, 4.7 this will place the X1 axis along X . This rotation
[R ] is made of the direction cosines of the axis (X ; Y ; Z ) with respect to (X1 ; Y1 ; Z1 ):
2 3
cos 0 sin
[R ] = 64 0 1 0 75 (4.33)
; sin 0 cos
q
we note that: cos = CCXZ
X , sin = CZ , and CXZ =
CXZ CX2 + CZ2 .
2. Rotation by
about the Z2 axis 4.8
2 3
cos
sin
0
[R
] = 4 ; sin
cos
6
0 7
5 (4.34)
0 0 1
where cos
= CXZ , and sin
= CY .
35 Combining Eq. 4.33 and 4.34 yields:
2 3
CX CY CZ
[
] = [R
][R ] = 64 ;CCXXZCY ; C
CXZ CXZ Y CZ 7
5 (4.35)
;CZ 0 CX
CXZ CXZ
Victor Saouma Matrix Structural Analysis
Draft
4{12 TRANSFORMATION MATRICES
36 For vertical member the preceding matrix is no longer valid as CXZ is undened. However
we can obtain the matrix by simple inspection, Fig. 4.9 as we note that:
1. X2 axis aligned with Y1
2. Y2 axis aligned with -X1
3. Z2 axis aligned with Z1
hence the rotation matrix with respect to the y axis, is similar to the one previously derived
for rotation with respect to the z axis, except for the reordering of terms:
2 3
0 CY 0
[
] = 64 ;CY 0 0 7
5 (4.36)
0 0 1
which is valid for both cases (CY = 1 for
= 90 deg, and CY = ;1 for
= 270 deg).
4.2.2.2 General Case
37 In the most general case, we need to dene an additional rotation to the preceding trans-
formation of an angle about the X
axis, Fig. 4.10. This rotation is dened such that:
Chapter 5
STIFFNESS METHOD; Part II
5.1 Introduction
1 The direct stiness method, covered in Advanced Structural Analysis is brie
y reviewed in
this lecture.
2 A slightly dierent algorithm will be used for the assembly of the global stiness matrix.
3 Some of the prescribed steps are further discussed in the next sections.
1. At input stage read ID(idof,inod) of each degree of freedom for every node such that:
(
;
ID(idof inod) =
0 if unrestrained d.o.f. (5.1)
1 if restrained d.o.f.
2. After all the node boundary conditions have been read, assign incrementally equation
numbers
(a) First to all the active dof
(b) Then to the other (restrained) dof.
(c) Multiply by -1 all the passive dof.
Note that the total number of dof will be equal to the number of nodes times the number
of dof/node NEQA.
3. The largest positive global degree of freedom number will be equal to NEQ (Number Of
Equations), which is the size of the square matrix which will have to be decomposed.
6 For example, for the frame shown in Fig. 5.1:
1. The input data le may contain:
5.3 LM Vector
7 The LM vector of a given element gives the global degree of freedom of each one of the element
degree of freedom's. For the structure shown in Fig. 5.1, we would have:
bLMc = b ;10 ;11 4 5 6 7 c element 1 (2 ! 3)
bLMc = b 5 6 7 1 2 3 c element 2 (3 ! 1)
bLMc = b 1 2 3 ;12 8 9 c element 3 (1 ! 4)
5.4 Assembly of Global Stiness Matrix
8 As for the element stiness matrix, the global stiness matrix [K] is such that Kij is the force
in degree of freedom i caused by a unit displacement at degree of freedom j .
9 Whereas this relationship was derived from basic analysis at the element level, at the structure
level, this term can be obtained from the contribution of the element stiness matrices [Ke ]
(written in global coordinate system).
10 For each Kij term, we shall add the contribution of all the elements which can connect degree
of freedom i to degree of freedom j , assuming that those forces are readily available from the
individual element stiness matrices written in global coordinate system.
11 Kij is non-zero if and only if degree of freedom i and degree of freedom j are connected by
an element or share a node.
12 There are usually more than one element connected to a dof. Hence, individual element
stiness matrices terms must be added up.
Victor Saouma Matrix Structural Analysis
Draft
5{4 STIFFNESS METHOD; Part II
13 Because each term of all the element stiness matrices must nd its position inside the global
stiness matrix [K], it is found computationally most eective to initialize the global stiness
matrix [KS ]NEQANEQA to zero, and then loop through all the elements, and then through
each entry of the respective element stiness matrix Kije .
14 The assignment of the element stiness matrix term Kij e (note that e, i, and j are all known
since we are looping on e from 1 to the number of elements, and then looping on the rows and
columns of the element stiness matrix i; j ) into the global stiness matrix KklS is made through
the LM vector (note that it is k and l which must be determined).
15 Since the global stiness matrix is also symmetric, we would need to only assemble one side
of it, usually the upper one.
16 Contrarily to Matrix Structural Analysis, we will assemble the full augmented stiness matrix.
2 3 8 9
1 3 11 17 24 32 41 68 >
>
>
1 >
>
>
6
6 2 5 10 16 23 31 40 67 7
7 >
>
>
> 2 >
>
>
>
6 7 > >
6
6
4 9 15 22 30 39 66 7
7
>
>
>
>
4 >
>
>
>
6
6 6 8 14 21 29
48 5638 65 7
7
>
>
>
>
>
6 >
>
>
>
>
6
6 7 13 20 28
47 5537 64 7
7
>
>
>
>
7 >
>
>
>
K = 66
6
6 12 19 27
46 5436 63 7
7
7 MAXA =
<
12 =
6
18 45 53 62 77
26 35 >
>
>
18 >
>
>
6
6 25
44 5234 61 7
7
>
>
>
>
>
25 >
>
>
>
>
6
6 43 5133 60 7
7
>
>
>
>
33 >
>
>
>
6
6
6
42 50 59 7
7
7
>
>
>
>
>
42 >
>
>
>
>
4 49 58 5 >
>
>
>
49 >
>
>
>
57 :
57 ;
Thus, to locate an element within the stiness matrix, we use the following formula:
30For internal book-keeping purpose, since we are assembling the augmented stiness matrix,
we proceed in two stages:
1. First number all the unrestrained degrees of freedom, i.e. those on ;t .
2. Then number all the degrees of freedom with known displacements, on ;u , and multiply
by -1.
31 Considering a simple beam, Fig. 5.5 the full stiness matrix is equal to
2
v1 1 v2 2 3
V1 12EI=L 3 6EI=L2 ;12EI=L 6EI=L2 7
3
[K] = M
6
16 6EI=L2 4EI=L ;6EI=L2 2EI=L 77 (5.9)
V2 4 ;12EI=L3 ;6EI=L2
6
12EI=L3 ;6EI=L2 5
M2 6EI=L2 2EI=L ;6EI=L2 4EI=L
This matrix is singular, it has a rank 2 and order 4 (as it embodies also 2 rigid body motions).
32 We shall consider 3 dierent cases, Fig. 5.6
>
> ;ML=6EI >
>
>
< ML=3EI >
=
= M=L
>
> >
>
> >
: ;M=L ;
Solution:
1. Determine the structure ID matrix and the LM vector for each element Initial ID matrix
Node" 1 2 3 4 5 #
ID = 01 0 1 0 0
0 1 0 0
Final ID matrix
Node" 1 2 3 4 5#
ID = 1 2 ;9 4 6
;8 3 ;10 5 7
LM vectors for each element
2 3
Element 1 1 ;8 4 5
6 1
Element 26 ;8 2 3 777
Element 36
6 2
6
3 4 5 77
6 4
Element 46 5 6 7 77
6 ;9
Element 56 ;10 4 5 777
6 2
Element 66 3 6 7 77
6
Element 74 2 3 ;9 ;10 5
Element 8 ;9 ;10 6 7
Victor Saouma Matrix Structural Analysis
Draft
5{16 STIFFNESS METHOD; Part II
2. Derive the element stiness matrix in global coordinates
2 3
c2 cs ;c2 ;cs 7
[K ] = EA 6 cs
6
s2 ;cs ;s2 77
L 64 ;c2 ;cs c2 cs 5
;cs ;s2 cs s2
where c = cos = x2 ;L x1 ; s = sin = Y2 ;L Y1
Element 1 L = 200 , c = 1620;0 = 0:8, s = 1220;0 = 0:6,
EA
L = (30;000)(10)
20 = 15; 000
2
1 ;8 4 5 3
1 9600 7200 ;9600 ;7200
[K1 ] = ;8666 7200 5400 ;7200 ;5400 777
4 4 ;9600 ;7200 9600 7200 5
5 ;7200 ;5400 7200 5400
Element 2 L = 160 ; c = 1; s = 0; EA
L = 18; 750.
1
2
;8 2 33
1 18; 750 0 ;18; 750 0
[K2 ] = ; 8666 0 0 0
2 4 ;18; 750 0 18; 750
0 777
05
3 0 0 0 0
Element 3 L = 120 ; c = 0; s = 1; EA
L = 25; 000
22
3 4 5 3
2 0 0 0 0
3 60
6
[K3 ] = 464 0 25; 000 0 ;25; 000 777
0 0 0 5
5 0 ;25; 000 0 25; 000
Element 4 L = 160 ; c = 1; s = 0; EA
L = 18; 750
2
4 5 6 73
4 18; 750 0 ;18; 750 0
[K4 ] = 56664 ;180; 750 0 0 0 777
6
0 18; 750 05
7 0 0 0 0
2
;9 ;10 4 5 3
;9 6 9600 ;7200 ;9600 7200
[K5 ] = ; 1066 ;7200 5400 7200
4 4 ;9600 7200 9600
;5400 777
;7200 5
5 7200 ;5400 ;7200 5400
Element 6 L = 200 ; c = 0:8; s = 0:6; EAL = 15; 000
2
2 3 6 7 3
2 9600 7200 ;9600 ;7200
[K6 ] = 36664 ;7200
6
5400 ;7200 ;5400 777
9600 ;7200 9600 7200 5
7 ;7200 ;5400 7200 5400
Element 7 L = 160 ; c = 1; s = 0; EA
L = 18; 750
2
2 3 ;9 ;10 3
2 18; 750 0 ;18; 750 0
[K7 ] = 3;9 664 ;180; 750 0 0 0 777
6
0 18; 750 0 5
;10 0 0 0 0
Element 8 L = 120 ; c = 0; s = 1; EA
L = 25; 000
2
;9 ;10 6 7 3
;9 6 0 0 0 0
[K8 ] = ;1066 0 25; 000 0 ;25; 000 77
7
6 4 0 0 0 0 5
7 0 ;25; 000 0 25; 000
3. Assemble the augmented global stiness matrix in kips/ft.
1 2 3 4 5 6 7
12 9; 600 + 18; 750 ;18; 750 0 ;9; 600 ;7; 200 0 0 3
2 9; 600 + (2)18; 750 7; 200 0 0 ;9; 600 ;7; 200
36 5; 400 + 25; 000 0 ;25; 000 ;7; 200 ;5; 400 7
ktt = 46
4
18; 750 + (2)9; 600 7; 200 ; 7; 200 ;18; 750 0 7
5
5 SYMMETRIC 25; 000 + 5; 400(2) 0 0
6 18; 750 + 9600 7200
7 25; 000 + 5; 400
;8 ;9 ;10 3
12 0 + 7; 200
2 0 ;18; 750 ; 18; 750 0
36 0 0 0 7
ktu = 44 ;7; 200
6 ;9; 600 7; 200 7
5 ;5; 400 7; 200 ;5; 400 5
6 0
7 ;25; 000
Victor Saouma Matrix Structural Analysis
Draft
5{18 STIFFNESS METHOD; Part II
;8 ;9 ;10
;8 0 + 5; 400
kuu = ;9 9; 600 + 18; 750 0 ; 7; 200 + 0
;10 ;7; 200 + 0 0 + 5; 400 + 25; 000
v2 >
>
>
u2 >
>
>
:
v2 ;
( )
n o1 h i ;0:0223 n o
v1 = 15;000 0:8 0:6 ;0:8 ;0:6 0 = 52:1 Compression
v2 12 ;0:8 ;0:6 0:8 0:6 ;0:0102 ;52:1
( ; :
) 0856
0
n o2 h i ;0:0233 n o
v1 = 18;750 1 0 ;1 0 0 = ;43:2 Tension
v2 12 ;1 0 1 0 0:00433 43:2
" ;0:116 #
n o3 h i 0:00433 n o
v1 = 25;000 0 1 0 ;1 ;0:116 = ;63:3 Tension
v2 12 0 ;1 0 1 ;0:0102 63:3
( ;0:0856 )
n o4 h i ;0:0102 n o
v1 = 18;750 1 0 ;1 0 ;0:0856 = ;1:58 Tension
v2 12 ;1 0 1 0 ;0:00919 1:58
n o h ;0: i n o n o
v1 5 = 15;000 ;0:8 0:6 0:8 ;0:6 ;0:0102 = 54:0
v2 12 0:8 ;0:6 ;0:8 0:6 ;0:0856 ;54:0 Compression
n o h i ;0:116 n o
v1 6 = 15;000 0:8 0:6 ;0:8 ;0:6 ;0:00919 = ;60:43 Tension
v2 12 ;0:8 ;0:6 0:8 0:6 ;0:0174 60:43
n o7 h
i ;0:116 n o
v1 = 18;750 1 0 ;1 0 0 = 6:72 Compression
v2 12 ;1 0 1 0 ;6:72
( 0 )
n o8 h i 0 n o
v1 = 25;000 0 1 0 ;1 0 = 36:3 Compression
v2 12 0 ;1 0 1 ;0:00919 ;36:3
;0:0174
50kN
4 kN/m
0
1
1
0
0
1
0
1
8m
3m
000
111
111
000
000
111
000
111
000
111 7.416 m 8m
;4 ;5 ;6 1 2 3
;4 A11 A12 A13 A14 A15 A16 3
2
We note that some terms are equal to zero because we do not have a connection between
the corresponding degrees of freedom (i.e. node 1 is not connected to node 3).
00
11
11
00
00
11
00
11
8m
3m
000
111
111
000
000
111
000
111
000
111 7.416 m 8m
function [k,K,Gamma]=stiff(EE,II,A,i,j)
% Determine the length
L=sqrt((j(2)-i(2))^2+(j(1)-i(1))^2);
% Compute the angle theta (carefull with vertical members!)
if(j(1)-i(1))~=0
alpha=atan((j(2)-i(2))/(j(1)-i(1)));
else
alpha=-pi/2;
end
% form rotation matrix Gamma
Gamma=[
cos(alpha) sin(alpha) 0 0 0 0;
-sin(alpha) cos(alpha) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(alpha) sin(alpha) 0;
0 0 0 -sin(alpha) cos(alpha) 0;
0.0010
-0.0050
-0.0005
Reactions =
130.4973
55.6766
13.3742
-149.2473
22.6734
-45.3557
int_forces = int_forces =
141.8530 149.2473
2.6758 9.3266
13.3742 -8.0315
-141.8530 -149.2473
-2.6758 22.6734
8.0315 -45.3557
We note that the internal forces are consistent with the reactions (specially for the second node
of element 2), and amongst themselves, i.e. the moment at node 2 is the same for both elements
(8.0315).
1. Initialize the vector storing the compacted stiness matrix to zero.
2. Loop through each element, e, and for each element:
(a) Retrieve its stiness matrix (in local coordinates) [ke ], and direction cosines.
(b) Determine the rotation matrix [;] of the element.
(c) Compute the element stiness matrix in global coordinates from [bK e ] = [;]T [ke ][;].
(d) Dene the fLMg array of the element
(e) Loop through each row and column of the element stiness matrix, and for those
degree of freedom not equal to zero, add the contributions of the element to the
structure's stiness matrix (note that we assemble only the upper half).
Victor SaoumaFigure 5.11: Flowchart for the Skyline Height Matrix Structural Analysis
Determination
Draft
5.7 Computer Program Flow Charts 5{31
Figure 5.12: Flowchart for the Global Stiness Matrix Assembly
Victor Saouma Matrix Structural Analysis
Draft
5{32 STIFFNESS METHOD; Part II
5.7.4 Decomposition
38 Decompose the global stiness matrix. Since the matrix is both symmetric and positive
denite, the matrix can be decomposed using Cholesky's method into: [K] = [L][L]T . Should a
division by zero occur, or an attempt to extract the square root of a negative number happen,
then this would be an indication that either the global stiness matrix is not properly assembled,
or that there are not enough restraint to prevent rigid body translation or rotation of the
structure.
5.7.5 Load
39 Once the stiness matrix has been decomposed, than the main program should loop through
each load case and, Fig. 5.13
1. Initialize the load vector (of length NEQ) to zero.
2. Read number of loaded nodes. For each loaded node store the non-zero values inside the
load vector (using the [ID] matrix for determining storage location).
3. Loop on all loaded elements:
(a) Read element number, and load value
(b) Compute the xed end actions and rotate them from local to global coordinates.
(c) Using the LM vector, add the xed end actions to the nodal load vector (unless the
corresponding equation number is zero, ie. restrained degree of freedom).
(d) Store the xed end actions for future use.
5.7.6 Backsubstitution
40 Backsubstitution is achieved by multiplying the decomposed stiness matrix with the load
vector. The resulting vector stores the nodal displacements, in global coordinate system, cor-
responding to the unrestrained degree of freedom.
5.7.7 Internal Forces and Reactions
41The internal forces for each element, and reactions at each restrained degree of freedom, are
determined by, Fig. 5.15
1. Initialize reactions to zero
2. For each element retrieve:
(a) nodal coordinates
44 It is essential that the structure be idealized such that it can be discretized. This discretiza-
tion should dene each node and element uniquely. In order to decrease the required amount
of computer storage and computation it is best to number the nodes in a manner that mini-
mizes the numerical separation of the node numbers on each element. For instance, an element
connecting nodes 1 and 4, could be better dened by nodes 1 and 2, and so on. As it was
noted previously, the user is required to have a decent understanding of structural analysis and
structural mechanics. As such, it will be necessary for the user to generate or modify an input
le input.m using the following directions. Open the le called input.m and set the existing
variables in the le to the appropriate values. The input le has additional helpful directions
given as comments for each variable. After setting the variables to the correct values, be sure
to save the le. Please note that the program is case-sensitive.
npoin=3;
nelem=2;
nload=1;
% INPUT THE ID MATRIX CONTAINING THE NODAL BOUNDARY CONDITIONS (ROW # = NODE #)
ID=[1 1 1;
0 0 0;
1 1 1];
nodecoor=[
0 0;
7416 3000;
15416 3000
];
lnods=[
1 2;
2 3
];
E=[200 200];
A=[6000 6000];
Iz=[200000000 200000000];
% INPUT THE LOAD DATA. NODAL LOADS, PNODS SHOULD BE IN MATRIX FORM. THE COLUMNS CORRESPOND
% TO THE GLOBAL DEGREE OF FREEDOM IN WHICH THE LOAD IS ACTING AND THE THE ROW NUMBER CORRESPONDS
% WITH THE LOAD CASE NUMBER. PELEM IS THE ELEMENT LOAD, GIVEN IN A MATRIX, WITH COLUMNS
% CORRESPONDING TO THE ELEMENT NUMBER AND ROW THE LOAD CASE. ARRAY "A" IS THE DISTANCE FROM
% THE LEFT END OF THE ELEMENT TO THE LOAD, IN ARRAY FORM. THE DISTRIBUTED LOAD, W SHOULD BE
% IN MATRIX FORM ALSO WITH COLUMNS = ELEMENT NUMBER UPON WHICH W IS ACTING AND ROWS = LOAD CASE.
% ZEROS SHOULD BE USED IN THE MATRICES WHEN THERE IS NO LOAD PRESENT. NODAL LOADS SHOULD
% BE GIVEN IN GLOBAL COORDINATES, WHEREAS THE ELEMENT LOADS AND DISTRIBUTED LOADS SHOULD BE
% GIVEN IN LOCAL COORDINATES.
% IF YOU WANT THE PROGRAM TO DRAW THE STUCTURE SET DRAWFLAG=1, IF NOT SET IT EQUAL TO 0.
% THIS IS USEFUL FOR CHECKING THE INPUT DATA.
drawflag=1;
format short e
clear
istrtp=3;
indat
idrasmbl
elmcoord
% 2DFRAME CALCULATIONS
% CALCULATE THE LENGTH AND ORIENTATION ANGLE, ALPHA FOR EACH ELEMENT
% CALL SCRIPTFILE LENGTH3.M
length3
stiffl3
trans3
assembl3
print_general_info
print_loads
loads3
disp3
react3
intern3
end
draw
st=fclose('all');
% END OF MAIN PROGRAM (CASAP.M)
Pnods=Pnods.';
count=1;
negcount=-1;
if istrtp==3
ndofpn=3;
nterm=6;
else
error('Incorrect structure type specified')
end
orig_ID=ID;
for inode=1:npoin
for icoord=1:ndofpn
if ID(inode,icoord)==0
for ielem=1:nelem
LM(ielem,1:ndofpn)=ID(lnods(ielem,1),1:ndofpn);
LM(ielem,(ndofpn+1):(2*ndofpn))=ID(lnods(ielem,2),1:ndofpn);
end
% ASSEMBLE THE ELEMENT COORDINATE MATRIX, ELEMCOOR FROM NODECOOR AND LNODS
for ielem=1:nelem
elemcoor(ielem,1)=nodecoor(lnods(ielem,1),1);
elemcoor(ielem,2)=nodecoor(lnods(ielem,1),2);
%elemcoor(ielem,3)=nodecoor(lnods(ielem,1),3);
elemcoor(ielem,3)=nodecoor(lnods(ielem,2),1);
elemcoor(ielem,4)=nodecoor(lnods(ielem,2),2);
%elemcoor(ielem,6)=nodecoor(lnods(ielem,2),3);
end
% END OF ELMCOORD.M SCRIPTFILE
% COMPUTE THE LENGTH AND ANGLE BETWEEN LOCAL AND GLOBAL X-AXES FOR EACH ELEMENT
for ielem=1:nelem
L(ielem)=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
alpha(ielem)=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
for ielem=1:nelem
k(1:6,1:6,ielem)=...
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
end
for ielem=1:nelem
rotation(1:6,1:6,ielem)=...
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
ktemp=k(1:6,1:6,ielem);
% CALCULATE THE ELEMENT STIFFNESS MATRIX IN GLOBAL COORDINATES
Rot=rotation(1:6,1:6,ielem);
K(1:6,1:6,ielem)=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
end
% END OF TRANS3.M SCRIPTFILE
% RENUMBER DOF INCLUDE ALL DOF, FREE DOF FIRST, RESTRAINED NEXT
new_LM=LM;
number_gdofs=max(LM(:));
new_LM(find(LM<0))=number_gdofs-LM(find(LM<0));
aug_total_dofs=max(new_LM(:));
Ktt=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Ktu=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Kut=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Kuu=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% END OF ASSEMBL3.M SCRIPTFILE
fprintf(fid,'\nNode Info:\n');
for inode=1:npoin
fprintf(fid,' Node %d (%d,%d)\n',inode,nodecoor(inode,1),nodecoor(inode,2));
freedof=' ';
if(ID(inode,1))>0
freedof=strcat(freedof,' X ');
end
if(ID(inode,2))>0
freedof=strcat(freedof,' Y ');
end
if(ID(inode,3))>0
freedof=strcat(freedof,' Rot');
end
if freedof==' '
freedof=' none; node is fixed';
end
fprintf(fid,' Free dofs:%s\n',freedof);
end
fprintf(fid,'\nElement Info:\n');
for ielem=1:nelem
fprintf(fid,' Element %d (%d->%d)',ielem,lnods(ielem,1),lnods(ielem,2));
fprintf(fid,' E=%d A=%d Iz=%d \n',E(ielem),A(ielem),Iz(ielem));
end
Load_case=iload
if iload==1
fprintf(fid,'\n_________________________________________________________________________\n\n');
end
% CALCULATE THE FIXED END ACTIONS AND INSERT INTO A MATRIX IN WHICH THE ROWS CORRESPOND
% WITH THE ELEMENT NUMBER AND THE COLUMNS CORRESPOND WITH THE ELEMENT LOCAL DEGREES
% OF FREEDOM
for ielem=1:nelem
b(ielem)=L(ielem)-a(ielem);
Ffl=((w(ielem)*L(ielem))/2)+((Pelem(ielem)*(b(ielem))^2)/(L(ielem))^3)*(3*a(ielem)+b(ielem));
Mfl=((w(ielem)*(L(ielem))^2))/12+(Pelem(ielem)*a(ielem)*(b(ielem))^2)/(L(ielem))^2;
Ffr=((w(ielem)*L(ielem))/2)+((Pelem(ielem)*(a(ielem))^2)/(L(ielem))^3)*(a(ielem)+3*b(ielem));
Mfr=-((w(ielem)*(L(ielem))^2))/12+(Pelem(ielem)*a(ielem)*(b(ielem))^2)/(L(ielem))^2;
feamatrix_global=...
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
end
for idofpn=1:ndofpn
fea_vector(idofpn,1)=0;
end
for ielem=1:nelem
for idof=1:6
if ielem==1
if LM(ielem,idof)>0
fea_vector(LM(ielem,idof),1)=feamatrix_global(idof,ielem);
end
elseif ielem>1
if LM(ielem,idof)>0
fea_vector(LM(ielem,idof),1)=fea_vector(LM(ielem,1))+feamatrix_global(idof,ielem);
end
for ielem=1:nelem
for iterm=1:nterm
if feamatrix_global(iterm,ielem)==0
else
if new_LM(ielem,iterm)>number_gdofs
fea_vector_react(iterm,1)=feamatrix_global(iterm,ielem);
end
end
end
end
% CREATE A TEMPORARY VARIABLE EQUAL TO THE INVERSE OF THE STRUCTURAL STIFFNESS MATRIX
Ksinv=inv(Ktt);
Delta=
XXXXXXXXXXXXXXXXXXXXXX COMPLETE XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
fprintf(fid,' Displacements:\n');
for k=1:size(Delta,1)
%WORK BACKWARDS WITH LM MATRIX TO FIND NODE# AND DOF
LM_spot=find(LM'==k);
elem=fix(LM_spot(1)/(nterm+1))+1;
dof=mod(LM_spot(1)-1,nterm)+1;
node=lnods(elem,fix(dof/4)+1);
switch(dof)
case {1,4}, dof='delta X';
case {2,5}, dof='delta Y';
otherwise, dof='rotate ';
end
%PRINT THE DISPLACEMENTS
fprintf(fid,' (Node: %2d %s) %14d\n',node, dof, Delta(k));
end
fprintf(fid,'\n');
5.8.2.12 Reactions
%**********************************************************************************************
% Scriptfile name : react3.m (for 2d-frame structures)
%
% Main program : casap.m
%
% When this file is called, it calculates the reactions at the restrained degrees of
% freedom.
%
% Variable Descriptions:
%
% Reactions = Reactions at restrained degrees of freedom
% Kut = Upper left part of aug stiffness matrix, normal structure stiff matrix
% Delta = vector of displacements
% fea_vector_react = vector of fea's in restrained dofs
%
%
% By Dean A. Frank
% CVEN 5525 - Term Project
% Fall 1995
%
% Edited by Pawel Smolarkiewicz, 3/16/99
% Simplified for 2D Frame Case only
%
%**********************************************************************************************
fprintf(fid,' Reactions:\n');
for k=1:size(Reactions,1)
%WORK BACKWARDS WITH LM MATRIX TO FIND NODE# AND DOF
LM_spot=find(LM'==-k);
elem=fix(LM_spot(1)/(nterm+1))+1;
dof=mod(LM_spot(1)-1,nterm)+1;
node=lnods(elem,fix(dof/4)+1);
switch(dof)
case {1,4}, dof='Fx';
case {2,5}, dof='Fy';
otherwise, dof='M ';
end
%PRINT THE REACTIONS
fprintf(fid,' (Node: %2d %s) %14d\n',node, dof, Reactions(k));
end
fprintf(fid,'\n');
Pglobe=zeros(6,nelem);
Plocal=Pglobe;
%PRINT RESULTS
fprintf(fid,'\n Element: %2d\n',ielem);
for idof=1:6
if idof==1
fprintf(fid,' At Node: %d\n',lnods(ielem,1));
end
if idof==4
fprintf(fid,' At Node: %d\n',lnods(ielem,2));
end
switch(idof)
case {1,4}, dof='Fx';
case {2,5}, dof='Fy';
otherwise, dof='M ';
end
fprintf(fid,' (Global : %s ) %14d',dof, Pglobe(idof,ielem));
fprintf(fid,' (Local : %s ) %14d\n',dof, Plocal(idof,ielem));
end
end
fprintf(fid,'\n_________________________________________________________________________\n\n');
Node Info:
Node 1 (0,0)
Free dofs: none; node is fixed
Node 2 (7416,3000)
Free dofs: X Y Rot
Node 3 (15416,3000)
Free dofs: none; node is fixed
Element Info:
Element 1 (1->2) E=200 A=6000 Iz=200000000
8000
6000
4000
2 2 3
2000 1
1
0
−2000
−4000
−6000
0 2000 4000 6000 8000 10000 12000 14000 16000
_________________________________________________________________________
Load Case: 1
Nodal Loads:
Node: 2 Fx = 1.875000e+001
Node: 2 Fy = -4.635000e+001
Elemental Loads:
Element: 1 Point load = 0 at 0 from left
Distributed load = 0
Element: 2 Point load = 0 at 0 from left
Distributed load = 4.000000e-003
Displacements:
(Node: 2 delta X) 9.949820e-001
(Node: 2 delta Y) -4.981310e+000
(Node: 2 rotate ) -5.342485e-004
Reactions:
(Node: 1 Fx) 1.304973e+002
(Node: 1 Fy) 5.567659e+001
(Node: 1 M ) 1.337416e+004
(Node: 3 Fx) -1.492473e+002
(Node: 3 Fy) 2.267341e+001
(Node: 3 M ) -4.535573e+004
Element: 2
At Node: 2
(Global : Fx ) 1.492473e+002 (Local : Fx ) 1.492473e+002
(Global : Fy ) 9.326590e+000 (Local : Fy ) 9.326590e+000
(Global : M ) -8.031549e+003 (Local : M ) -8.031549e+003
At Node: 3
(Global : Fx ) -1.492473e+002 (Local : Fx ) -1.492473e+002
(Global : Fy ) 2.267341e+001 (Local : Fy ) 2.267341e+001
(Global : M ) -4.535573e+004 (Local : M ) -4.535573e+004
_________________________________________________________________________
Chapter 6
EQUATIONS OF STATICS and
KINEMATICS
Note: This section is largely based on chapter 6 of Mc-Guire and Gallagher, Matrix Structural
Analysis, John Wiley
1 Having developed the stiness method in great details, and prior to the introduction of energy
based methods (which will culminate with the nite element formulation), we ought to revisit
the
exibility method. This will be done by rst introducing some basic statics and kinematics
relationship.
2 Those relations will eventually enable us not only to formulate the
exibility/stiness relation,
but also other \by-products" such as algorithms for: 1) the extraction of a statically determinate
structure from a statically indeterminate one; 2) checking prior to analysis whether a structure
is kinematically unstable; 3) providing an alternative method of assembling the global stiness
matrix.
Table 6.1: Internal Element Force Denition for the Statics Matrix
of independent element internal forces which can be selected.
6 Matrix [B ] will be a square matrix for a statically determinate structure, and rectangular
(more columns than rows) otherwise.
7 For the case of a statically indeterminate structure, Eq. 6.1 can be generalized as:
( )
fPg2n1 = [ [B0]2n2n [Bx]2nr F0
] F (6.6)
x (2n+r)1
where [B0 ] is a square matrix, fF0 g the vector of unknown internal element forces or external
reactions, and fFx g the vector of unknown redundant internal forces or reactions.
8 Hence, we can determine fF0 g from
( )
fPg2n1 = [ [B0]2n2n [Bx]2nr F0
] F (6.11-a)
x (2n+r)1
8
fP9g = [2B0] fF0g + [Bx] fFxg 38 9 8 9
>
> Px1 >
> 0 0 ; C 0 1 0 0 0 > F1
> >
> >
> 0 >
>
6 ;1 0 ;S 0 0 1 0 0 77 >
> > 7> > > >
>
>
>
> Py1 >
>
>
>
6 >
>
> F2 >
>
>
>
>
>
>
> 0 >
>
>
>
>
>
>
>
>
Px2 >
>
>
>
>
6
6
6
0 ; 1 0 0 0 0 0 0 7
7>
>
>
>
>
F3 >
>
>
>
>
>
>
>
>
>
;C >
>
>
>
>
<
Py2 =
=
6
6 1 0 0 <
0 0 0 0 0 7
7 F4 =
+ > S0
< =
fF5 g(6.11-b)
>
>
>
P
x3 >
>
>
6
6 0 1 S 0 0 0 0 0 7>
7> Rx1 >
> > >
> | {z }
7> > > >
f xg F
>
>
>
>
>
P
y3
>
>
>
>
>
6
6
6
0 0 S 1 0 0 0 0 7>
7>
>
>
>
Ry1 >
>
>
>
>
>
>
>
>
>
0 >
>
>
>
>
>
>
>
> xP4 >
>
>
>
4 0 0 0 0 0 0 1 0 5>>
>
>
Rx4 >
>
>
>
>
>
>
>
C >
>
>
>
:
P
y4
| {z
;
} |
0 0 0
{z
;1 0 0 0 1 :
}|
Ry4
{z
;
}
:
| {z
;S ;
}
fPg [B0 ] fF0 g [Bx ]
We can solve for the internal forces in terms of the (still unknown) redundant force
fF0g = [B0 ];1 fPg ; [B0];1 [Bx] fFxg 2 3
(6.12-a)
8 9 2 3 68 9 8 9 7
>
> F1 >
> 0 0 0 1 0 0 0 0 6>
6> 0 >
> >
> 0 >
>
7
7
>
>
>
>
> F2 >
>
>
>
>
6
6 0 0 ;1 0 0 0 0 0 7 6>
7 6>
>
> 0
>
>
>
>
>
>
>
> 0
>
>
>
>
7
7
7 6> > > >
>
>
>
> F3 >
>
>
>
6
6 0 0 C1 0 C1 0 0 0 7 6>
7 6>
>
> x2 P >
>
>
>
>
>
>
> ;C >
>
>
>
7
7
0 ; CS 0 ; CS 1 0 0
> > 6 > > > > 7
<
F4 =
=
6
6 0 7 6< y 2
76 P =
; > S0
< = 7
fF g (6.12-b)
7
> Rx1 > 6 1 0 1 0 1 0 0 0 7 6> 0 > >
5
| {z }7
>
> >
> 6 7 6>> >
> >
> >
> 7
>
>
>
> Ry 1 >
>
>
>
6
6 0 1 CS 1 CS 0 0 0 7 6>
7 6>> 0 >
>
>
>
>
> 0 >
>
> F
f x g 77
> > 6 7 6>> >
> >
> >
> 7
>
>
>
>
Rx4 >
>
>
>
4 0 0 0 0 0 0 1 0 5 6>
6>
>
>
0 >
>
>
>
>
>
>
>
C >
>
>
>
7
7
:
|
Ry 4
{z
;
} |
0 0 ; CS 0 ; CS 1 0 1
{z
6:
4
} | {z
0 ;
}
:
| {z
;S ;
}
7
5
fF0 g [B0 ];1 [C1 ] fPg [Bx ]
Or using the following relations [B0 ];1 [C1 ] and ;[B0 ];1 [Bx ] [C2 ] we obtain
8 9 2 38 9 8 9
>
> F1 >
> 0 0 0 1 0 0 0 0 >
> 0 >
> >
> ;S >
>
>
>
>
>
> F2 >
>
>
>
>
6 0
6 0 ;1 0 0 0 0 0 7>
7>
>
>
> 0
>
>
>
>
>
>
>
>
>
> ;C >
>
>
>
>
>
>
>
> F3 >
>
>
>
6
6 0 0 C1 0 C1 0 0 0 7>
7>>
> Px2 >
>
>
>
>
>
>
> 1 >
>
>
>
0 ; CS 0 ; CS ;S fF5 g
> > 6 7> > > >
<
F4 =
=
6 0
6 1 0 0 7<
7 Py2 +
= < =
(6.13)
>
>
>
Rx1 >
>
>
6 1
6 0 1 0 1 0 0 0 7>
7> 0 >
> > > C >
> | {z }
7> >
> > >
fFx g
>
>
>
> Ry 1 >
>
>
> 6 0
6
1 CS 1 CS 0 0 0 7>
7>
>
> 0 >
>
>
> > >
>
> 0 >
>
>
>
>
>
>
>
Rx4 >
>
>
>
>
6
4 0 0 0 0 0 0 1 0 5>
>
>
>
>
0 >
>
>
>
>
>
>
>
>
>
;C >
>
>
>
>
:
Ry 4 ;
0 0 ; CS 0 ; CS 1 0 1 :
0 ; :
0 ;
| {z } | {z }| {z } | {z }
fF0 g [C1 ] fPg [C2 ]
2
F1 F2 F3 F4 F5 Rx1 Ry1 Rx4 Ry4 Px2 Py2 3
A 0 0 ;C 0 0 1 0 0 0 0 0
B666 ;1 0 ;S 0 0 0 1 0 0 0 0 777
C66 0 ;1 0 0 ;C 0 0 0 0 ;1 0 77
D66 1 0 0 0 S 0 0 0 0 0 ;1 77 (6.18)
E66 0 1 C 0 0 0 0 0 0 0 0 77
F666 0 0 S 1 0 0 0 0 0 0 0 777
G4 0 0 0 0 C 0 0 1 0 0 0 5
H 0 0 0 ;1 ;S 0 0 0 1 0 0
2. Interchange columns
2
Rx1 Ry1 F2 F1 F3 F4 F5 Rx4 Ry4 Px2 Py2 3
A 1 0 0 0 ;C 0 0 0 0 0 0
B666 0 1 0 ;1 ;S 0 0 0 0 0 0 777
C66 0 0 ;1 0 0 0 ;C 0 0 ;1 0 77
D66 0 0 0 1 0 0 S 0 0 0 ;1 77 (6.19)
E66 0 0 1 0 C 0 0 0 0 0 0 77
F666 0 0 0 0 S 1 0 0 0 0 0 777
G4 0 0 0 0 0 0 C 1 0 0 0 5
H 0 0 0 0 0 ;1 ;S 0 1 0 0
2
Rx1 Ry1 F2 F1 F3 F4 F5 Rx4 Ry4 Px 2 Py 2 3
A0 = A + CE 0 1 0 0 0 0 0 ;C 0 0 ;1 0
B 00 = B 0 + SE6660 0 1 0 0 0 0 0 0 0 ;S=C ;1 777
C 0 = ;C 6 0 0 1 0 0 0 C 0 0 1 0 77
D 6
6
6 0 0 0 1 0 0 S 0 0 0 ;1 77 (6.21)
E0 6
6 0 0 0 0 1 0 ;1 0 0 ;1=C 0 777
F 0 = F ; SE 0 66
6
0 0 0 0 0 1 S 0 0 S=C 0 77
G 4 0 0 0 0 0 0 C 1 0 0 0 5
H0 = H + F 0 0 0 0 0 0 0 0 0 1 S=C 0
5. Interchange columns and observe that F5 is the selected redundant.
2
Rx1 Ry1 F2 F1 F3 F4 Rx4 Ry4 F5 Px2 Py2 3
A0 1 0 0 0 0 0 0 0 ;C ;1 0
B 00666 0 1 0 0 0 0 0 0 0 ;S=C ;1 777
C 06 0 0 1 0 0 0 0 0 C 1 0 7
D666 0 0 0 1 0 0 0 0 S 0 ;1 777 (6.22)
E 066 0 0 0 0 1 0 0 0 ;1 ;1=C 0 777
F 0666 0 0 0 0 0 1 0 0 S S=C 0 77
G4 0 0 0 0 0 0 1 0 C 0 0 5
H0 0 0 0 0 0 0 0 1 0 S=C 0
| {z }
[A]
| {z }
[A]
We should observe that [A] is indeed the transpose of the [B] matrix in Eq. 6.17
where Fx correspond to the redundant forces. The external work will then be
" #
t
Wext = 21 b F0 Fx c [[BB0]]t fg (6.35)
x
27 As before, equating the external to the internal work Wext = Wint and simplifying, we obtain:
[B0 ]T = [A0 ] (6.38)
[Bx]T = [Ax] (6.39)
This equation relates the unknown relative displacements to the relative known ones.
and the congruent transformation approach: [K]neqneq = [A]Tneq6n [Ke ]6n6n [A]6nneq .
6. Congruent approach appears to be less ecient than the direct stiness method as both
[Ke ] and [A] are larger than [K].
The stiness matrices of elements AB and BC in local coordinate system are given by:
2
:75 0: 0: ;:75 0: 0: 3
6
6 :00469 18:75 0 ;:0048 18:75 77
6
[k]AB = [k]BC = 200 666 1 105 0 ;18:75 :5 105 777 (6.64)
6
:75 0: 0: 77
4 sym :00469 ;18:75 5
1 105
while the rotation matrix is given by:
2
:9272 :375 0: 0: 0 : 0: 3
6 ;:375 :927 0: 0: 0: 0: 77
6
[;]AB = 666 00:: 0: 1 0: 0: 0: 777
6
8
f9g = [A] fg (6.72)
2 3
>
> u1 >
> 1 0 0
>
>
>
>
>
v1 >
>
>
>
>
6
6 0 1 0 78
7> u
9
>
<
1 =
=
6
6 0 0 1 7<
7
v
=
(6.73)
>
>
>
u2 >
>
>
6
6
6
1 0 0 7
7>
7: >
;
>
>
>
>
:
v2 >
>
>
>
;
4 0 1 0 5
2 0 0 1
3. Finally, if we take the product |[A{z]T} [|K{ze}] |{z}
[A] we obtain the structure global stiness
36 66 63
matrix [K] in Eq. 6.69
Congruent Transformation (local axis):
1. Unconnected stiness matrix in local coordinates:
fpe g = [ke ] fe g (6.74)
8 9 2 38 u1x 9
> PX1 > :75 0: 0: > >
>
>
< PY1 >
>
= 6 :00469 ;18:75 7>
>
< vy1 >
>
=
MZ1 = 200 6 sym 1 105 7 z1 (6.75)
PX2 6 :75 0: 0: 7 u2x
>
> >
> 4 5>> >
>
>
: PY2 >
; :00469 18:75 >
: vy2 >
;
M2
Z sym 1 105 z2
Chapter 7
FLEXIBILITY METHOD
7.1 Introduction
1 Recall the denition of the
exibility matrix
where fg, [d], and fpg are the element relative displacements, element
exibility matrix, and
forces at the element degrees of freedom free to displace.
2 As with the congruent approach for the stiness matrix, we dene:
9 This equation should be compared with Eq. 7.1 and will be referred to as the unsolved global
assembled
exibility equation.
Victor Saouma Matrix Structural Analysis
Draft
7.2 Flexibility Matrix 7{3
7.2.1 Solution of Redundant Forces
10 We can solve for the redundant forces (recall that in the
exibility method, redundant forces
are the primary unknowns as opposed to displacements in the stiness method) by solving the
lower partition of Eq. 7.10:
fFxg = ; [Dxx];1 [Dxp] fPg (7.14)
0
7
7 1 >
(7.19)
6 7 6 7 >
> >
> >
> >
>
6 ;S=C 0 7 6 ;0:75 0 ; S ;0:6
6 7 6 7 >
> >
> >
> >
>
7 >
> >
> >
> >
>
6
4 0 0 5 64 0
7
0 7
5
>
>
>
>
>
; C >
>
>
>
>
>
>
>
>
>
;0:8 >
>
>
>
>
;S=C 0 ;0:75 0 : ; :
0 0 ;
>
f3 >
= 66
66 1:25 0 77 66 0:790 ;0:185
77
77 Py 2 (7.24-b)
>
>
>
>
f4 >
>
>
>
66
66
;0:75 0 777 666 ;0:474 0:111
77
77
>
>
>
>
>
Rx4 >
>
>
>
>
66
66
44 0 0 5 4 ;0:632 0:148
77
77
55
>
>
>
>
Ry 4 >
>
>
> ;0:75 0 0 0
:
f5 ;
2 3
0:368 0:148
6
6 0:75 1:000 7
7
6
6
6
;0:368 ;0:148 7
7(
7 )
=
6
6 0:474 0:889 7
7 Px2 (7.24-c)
6
6 0:460 0:185 7
7 Py2
6
6 ;0:276 ;0:111 7
7
6
4 0:632 ;0:148 7
5
;0:750 0
Finally, the displacements are obtained from Eq. 7.17
h i
fpg = [Dpp] ; [Dpx] [Dxx];1 [Dxp] fPg (7.25-a)
| {z }
[D]
( ) "" # " ## ( )
u2 L
= AE 4:797 0 ; 1 14:730 ;3:454 Px2 (7.25-b)
v2 0 1:500 4:860 ;3:454 0:810 Py 2
" #( )
L 1 : 766 0
= AE 0:711 1:333: 711 P x 2 (7.25-c)
Py2
It should be noted that whereas we have used the
exibility method in its algorithmic
implementation (as it would lead itself to computer implementation) to analyse this simple
problem, the solution requires a formidable amount of matrix operations in comparaison with
the \classical" (hand based)
exibility method.
those d.o.f.'s which are supported, and 2) subscript `f ' for those dof which are free.
( ) " #( )
Pf = kff kfs f (7.26)
Ps ksf kss s
7.3.1 From Stiness to Flexibility
15 To obtain [d] the structure must be supported in a stable and statically determinate way, as
for the beam in Fig. 7.1. for which we would have:
( )
ff g = 1 (7.27)
2
( )
fsg = v1 (7.28)
v2
( )
fPf g = M1 (7.29)
M2
( )
fPsg = V1 (7.30)
V2
Since fs g = f0g the above equation reduces to:
( ) " #
Pf = kff ff g (7.31)
Ps ksf
and we would have:
fPf g = [kff ] ff g (7.32)
20 In summary we have:
" #
[k] = [B[d][d] ];1 [B[d][d] ];[1B[B] ]T
;1 ;1 T
(7.45)
21 A very important observation, is that the stiness matrix is obviously singular, since the
second \row" is linearly dependent on the rst one (through [B]) and thus, its determinent is
equal to zero.
2. The statics matrix [B] relating external to internal forces is given by:
( ) " #( )
V1 = L1 ;11 ;11 M1 (7.47)
V2 M2
| {z }
[B]
6. [kss]:
[kss ] = [B] [d];1 [B]T = [ksf ] [B]T (7.51)
" #" # " #
EI 1
= L2 L ;6 ;6 6 6 1 ; 1 EI 12 ;12
1 ;1 = L3 ;12 12 (7.52)
and
M = M3 + N (R ; R cos ) + V (R2 sin ) (7.57-a)
M + 13 = R(1 ; cos1 ) + R sin2 (7.57-b)
fij = Disp:inDOFicausedbyunitloadinDOFj (7.57-c)
f11 = R Z M M d (7.57-d)
EI Zo p1 D1
= R R2 (1 ; cos )2 d (7.57-e)
EI o
= R3 Z (1 ; 2 cos + cos2 )d (7.57-f)
EI o
= R3 [ ; 2 sin + =2 = 1=4 sin 2] (7.57-g)
EI o
h i
R3 3 ; 2 sin + 1 sin 2
= EI (7.57-h)
2 4
f12 = R Z
f21 = EI R2 sin (1 ; cos )d (7.57-i)
o
= R3 Z (sin ; cossin)d (7.57-j)
EI o
= R3 ; cos ; 1 sin2 (7.57-k)
EI 2 o
= R3 (; cos ; 1 sin2 ) ; (;1 ; 0) (7.57-l)
EI 2
h i
R3 1 ; cos ; 1 sin2
= f21 = EI (7.57-m)
2
f13 R Z
= f31 = EI R(1 ; cos )d (7.57-n)
o
R 2
= EI [ ; sin ]o (7.57-o)
R2 [ ; sin ]
= EI (7.57-p)
R
f22 = EI
Z 2 2
R sin d (7.57-q)
o
Victor Saouma Matrix Structural Analysis
Draft
7.5 Duality between the Flexibility and the Stiness Methods 7{11
7.5 Duality between the Flexibility and the Stiness Methods
FLEXIBILITY STIFFNESS
Indeterminancy Static Kinematic
Primary Unknows Nodal Forces Nodal Displacements
Variational Principle Virtual Force Virtual Displacement
fpg = [B0] fPg fg = [;] fg
[d] [k]
[ D ] = [ B ] T [d][B] [ k ] = [ ;]T [k][;#]( )
( ) " #( ) ( ) "
P = D11 D12 P P kff Dfr f
R D21 D22 R R = Drf Drr r
;
fRg = [D22 ] (fR g ; [D21 ] fpg) ff g = [Kff ] (fPg ; [kfr ] fr g)
1 ;1
fP g = [D11 ];1 (fPg + [D12 ] fRg) fRg = [Krf ] ff g + [krr ] fr g
fpg = [B] fPg fpg = [k][;] fg
Part II
Introduction to Finite Elements
Draft
Draft
Chapter 8
REVIEW OF ELASTICITY
8.1 Stress
1 A stress, Fig 8.1 is a second order cartesian tensor, ij where the 1st subscript (i) refers to
ti = nj ij (8.2)
2 3
11 12 13
b t1 t2 t3 c = b| n1 n{z2 n3 c} 64 21 22 23 7
5 (8.3)
direction cosines 31 32 33
| {z }
stress tensor
where n is a unit outward vector normal to the plane.
4 Note that the stress is dened at a point, and a traction is dened at a point and with respect
to a given plane orientation.
5 When expanded in cartesian coordinates,, the previous equation yields
= Lu (8.9)
or 2 @ 3
8
> "xx 9
> @x 0 0
> >
0 @ 0
>
>
>
> "yy >
>
>
>
6
6 @y
78
7
ux
9
>
<
"zz >
=
=
6
6
6
0 0 @z@ 7>
7<
7 uy
>
=
(8.10)
"xy @ @
@x 0
> > 6 7> >
>
>
"xz
>
> 6 @y
@
7: uz ;
>
>
>
>
:
>
>
>
>
;
6
4 @z 0 @x@ 7
5
"yz 0 @z@ @y@
LT + b = 0 (8.13)
or 8 9
>
> xx >
>
2 @ 0 @y@ @z@ 0
0 3>>
>
> yy >
>
>
>
8
bx
9
@x >
<
zz >
= >
< >
=
6
4 0 @y@ 0 @x@ 0 @z@ 7
5
> xy >
+> by >
=0 (8.14)
0 0 @z@ 0 @x@ @y@ >
>
>
>
> xz
>
>
>
>
>
: bz ;
>
: >
;
yz
13 Expanding
@xx + @xy + @xz + bx = 0
@x @y @z
@yx + @yy + @yz + b = 0 (8.15)
@x @y @z y
@zx + @zy + @zz + bz = 0
@x @y @z
16 In 3D, this would yield 9 equations in total, however only six are distinct. In 2D, this results
in (by setting i = 2, j = 1 and l = 2):
@ 2 "11 + @ 2 "22 = @ 2
12 (8.17)
@x22 @x21 @x1 @x2
(recall that 2"12 =
12 .
17 When he compatibility equation is written in term of the stresses, it yields:
19 The (fourth order) tensor of elastic constants Cijklhas 81(34 ) components however, due to
the symmetry of both and , there are at most 36 9(92;1) distinct elastic terms.
20 For the purpose of writing Hooke's Law, the double indexed system is often replaced by a
simple indexed system with a range of six:
62 =36
z }| {
k = Dkm "m k; m = 1; 2; 3; 4; 5; 6 (8.20)
21 For isotropic bodies (elastic properties independent of reference system used to describe it),
it can be shown that the number of independent elastic constants is two. The stress strain
relations can be written in terms of E and as:
" = 1 + ;
ij E ij E ij kk (8.21)
E
ij = 1 + "ij + 1 ; 2 ij "kk (8.22)
where ij is the kroneker delta and is equal to 1 if i = j and to 0 if i 6= j . When the strain
equation is expanded in 3D cartesian coordinates it would yield:
"xy = E1 [xx ; (yy + zz )]
"yy = E1 [yy ; (zz + xx )]
"zz = E1 [zz ; (xx + yy )] (8.23)
"xy = 1+E xy
"yz = 1+E yz
"zx = 1+E zx
8.6 Summary
26 Fig. 8.4 illustrates the fundamental equations in solid mechanics.
Essential B.C.
ui : ;u
Body Forces
?
Displacements
bi ui
?
Equilibrium
?
Kinematics
ij;j + bi = 0 @ui + @uj
"ij = 12 @xj @xi
? ?
Stresses Constitutive Rel. Strain
ij
- ij = Dijkl "kl
"ij
6
Nonessential B.C.
ti : ;t
Chapter 9
VARIATIONAL AND ENERGY
METHODS
9.1 y Variational Calculus; Preliminaries
From Pilkey and Wunderlich book
u, u
u(x)
C
B
u(x)
du
A dx
x
x=a x=c x=b
5 Letting u~ to be a family of neighbouring paths of the extremizing function u(x) and we assume
that at the end points x = a; b they coincide. We dene u~ as the sum of the extremizing path
and some arbitrary variation, Fig. 9.1.
u~(x; ") = u(x) + "(x) = u(x) + u(x) (9.3)
where " is a small parameter, and u(x) is the variation of u(x)
u = u~(x; ") ; u(x) (9.4-a)
= "(x) (9.4-b)
and (x) is twice dierentiable, has undened amplitude, and (a) = (b) = 0. We note that
u~ coincides with u if " = 0
6 It can be shown that the variation and derivation operators are commutative
)
d (u) = u~0 (x; ") ; u0 (x)
dx 0 d (u) = du (9.5)
0 0
u = u~ (x; ") ; u (x) dx dx
7 Furthermore, the variational operator and the dierential calculus operator d can be simi-
larly used, i.e.
(u0 )2 = 2u0 u0 (9.6-a)
(u + v) = u
Z Z
+ v (9.6-b)
udx = (u)dx (9.6-c)
@u x + @u y
u = @x (9.6-d)
@y
Victor Saouma Matrix Structural Analysis
Draft
9.1 y Variational Calculus; Preliminaries 9{3
however, they have clearly dierent meanings. du is associated with a neighbouring point at
a distance dx, however u is a small arbitrary change in u for a given x (there is no associated
x).
8 For boundaries where u is specied, its variation must be zero, and it is arbitrary elsewhere.
The variation u of u is said to undergo a virtual change.
9 To solve the variational problem of extremizing , we consider
Z b
(u + ") = (") = F (x; u + "; u0 + "0 )dx (9.7)
a
11 From Eq. 9.3 u~ = u + ", and u~0 = u0 + "0 , and applying the chain rule
d(") = Z b @F du~ + @F du~0 dx = Z b @F + 0 @F dx (9.9)
d" a @ u~ d" @ u~0 d" a @ u~ @ u~0
for " = 0, u~ = u, thus
d(") = Z b @F + 0 @F dx = 0 (9.10)
d" "=0 a @u @u0
Integration by part (Eq. 5.1 and 5.1) of the second term leads to
b 0 @F b
Z
@u0 dx = @F
Z b
d @F
@u0 ; (x) dx @u0 dx
(9.11)
a a a
Using the end conditions (a) = (b) = 0, Eq. 9.10 leads to
Z b @F d @F
(x) @u ; dx @u0 dx = 0 (9.12)
a
12The fundamental lemma of the calculus of variation states that for continuous (x) in a
x b, and with arbitrary continuous function (x) which vanishes at a and b, then
Z b
(x) (x)dx = 0 , (x) = 0 (9.13)
a
Thus,
@F ; d @F = 0 (9.14)
@u dx @u0
Victor Saouma Matrix Structural Analysis
Draft
9{4 VARIATIONAL AND ENERGY METHODS
13 This dierential equation is called the Euler equation associated with and is a necessary
condition for u(x) to extremize .
14 Generalizing for a functional which depends on two eld variables, u = u(x; y ) and v =
v(x; y) Z Z
= F (x; y; u; v; u;x ; u;y ; v;x ; v;y ; ; v;yy )dxdy (9.15)
There would be as many Euler equations as dependent eld variables
8
< @F ; @ @F ; @ @F + @ 2 @F + @ 2 @F + @ 2 @F
@u @x @u;x @y @u;y @x22 @u;xx @x@y @u;xy @y22 @u;yy = 0 (9.16)
: @F ; @ @F ; @ @F + @ 2 @F + @
2 @F @ @F
@v @x @v;x @y @v;y @x @v;xx @x@y @v;xy + @y2 @v;yy = 0
15 We note that the Functional and the corresponding Euler Equations, Eq. 9.1 and 9.14, or
Eq. 9.15 and 9.16 describe the same problem.
16 The Euler equations usually correspond to the governing dierential equation and are referred
to as the strong form (or classical form).
17 The functional is referred to as the weak form (or generalized solution). This classication
stems from the fact that equilibrium is enforced in an average sense over the body (and the
eld variable is dierentiated m times in the weak form, and 2m times in the strong form).
18 It can be shown that in the principle of virtual displacements, the Euler equations are the
equilibrium equations, whereas in the principle of virtual forces, they are the compatibility
equations.
19 Euler equations are dierential equations which can not always be solved by exact meth-
ods. An alternative method consists in bypassing the Euler equations and go directly to the
variational statement of the problem to the solution of the Euler equations.
20 Finite Element formulation are based on the weak form, whereas the formulation of Finite
Dierences are based on the strong form.
21 We still have to dene . The rst variation of a functional expression is
F = R@F @F 0 ) Z b
@F u + @F u0 dx
@ub u + @u0 u = (9.17)
= a Fdx a @u @u0
As above, integration by parts of the second term yields
=
Z b @F d @F
u @u ; dx @u0 dx (9.18)
a
22 We have just shown that nding the stationary value of by setting = 0 is equivalent to
nding the extremal value of by setting d(d"") "=0 equal to zero.
23 Similarly, it can be shown that as with second derivatives in calculus, the second variation
2 can be used to characterize the extremum as either a minimum or maximum.
Victor Saouma Matrix Structural Analysis
Draft
9.1 y Variational Calculus; Preliminaries 9{5
9.1.2 Boundary Conditions
24 Revisiting the integration by parts of the second term in Eq. 9.10, we had
Z b 0 @F
dx = @F b ; Z b d @F dx (9.19)
a @u0 @u0 a a dx @u0
We note that
1. Derivation of the Euler equation required (a) = (b) = 0, thus this equation is a state-
ment of the essential (or forced) boundary conditions, where u(a) = u(b) = 0.
@F0 = 0 at x = a and b.
2. If we left arbitrary, then it would have been necessary to use @u
These are the natural boundary conditions.
25 For a problem with, one eld variable, in which the highest derivative in the governing
dierential equation is of order 2m (or simply m in the corresponding functional), then we have
Essential (or Forced, or geometric) boundary conditions, involve derivatives of order zero
(the eld variable itself) through m-1. Trial displacement functions are explicitely required
to satisfy this B.C. Mathematically, this corresponds to Dirichlet boundary-value problems.
Nonessential (or Natural, or static) boundary conditions, involve derivatives of order m
and up. This B.C. is implied by the satisfaction of the variational statement but not
explicitly stated in the functional itself. Mathematically, this corresponds to Neuman
boundary-value problems.
26 Table 9.1 illustrates the boundary conditions associated with some problems
Problem Axial Member Flexural Member
Distributed load Distributed load
Dierential Equation AE ddx2 u2 + q = 0 EI ddx4 w4 ; q = 0
m 1 2
Essential B.C. [0; m ; 1] u w; dw
dx
Natural B.C. [m; 2m ; 1] du
dx
d2 w2
dxand ddx3 w3
or x = Eu;x or M = EIw;xx and V = EIw;xxx
σ σ
* *
U0 U0
A A
U0 U0
A A
ε ε
Nonlinear Linear
U def
= U0 d
(9.36)
35 To obtain a general form of the internal strain energy, we rst dene a stress-strain relation-
ship accounting for both initial strains and stresses
= D( ; 0 ) + 0 (9.37)
where D is the constitutive matrix; is the strain vector due to the displacements u; 0 is the
initial strain vector; 0 is the initial stress vector; and is the stress vector.
36 The initial strains and stresses are the result of conditions such as heating or cooling of a
system or the presence of pore pressures in a system.
37 The strain energy U for a linear elastic system is obtained by substituting
= D (9.38)
with Eq. 9.33 and 9.37
Z Z Z
U = 21 T Dd
; T D0 d
+ T 0 d
(9.39)
where
is the volume of the system.
38 Considering uniaxial stresses, in the absence of initial strains and stresses, and for linear
elastic systems, Eq. 9.39 reduces to
Z
U = 12 " |{z}
E" d
(9.40)
Victor Saouma Matrix Structural Analysis
Draft
9{10 VARIATIONAL AND ENERGY METHODS
39 When this relation is applied to various one dimensional structural elements it leads to
Axial Members:
U = "
Z 9
d
>
>
P
2 >
>
= Z L P2
=A U = 12 dx (9.41)
P
" = AE >
>
> 0 AE
>
;
d
= Adx
Torsional Members:
Z 9
U = 12 " |{z}
E" d
>
>
>
>
Z
>
>
>
>
>
1
U = 2 xy G xy d
>
>
>
>
| {z } >
>
= Z L T2
xy
xy = JTr U = 21 dx (9.42)
>
> 0 GJ
xy = Gxy >
>
>
>
>
dZ
Z= rddrdx >
>
>
>
r 2 2 >
>
>
r ddr = J >
;
o 0
Flexural Members:
Z 9
U = 21 " |{z}E" >
>
>
>
>
>
>
x = MIzz y >
>
=
1
Z L
M 2 dx
" = MEIzzy >
U= 2 0 EIz (9.43)
>
dZ
= dAdx >
>
>
>
>
>
y2 dA = Iz >
;
A
9.2.2.1 Internal Work versus Strain Energy
40During strain increment, the work done by internal forces in a dierential element will be
the negative of that performed by the stresses acting upon it.
Z
Wi = ; dd
(9.44)
41 If the strained elastic solid were permitted to slowly return to their unstrained state, then
the solid would return the work performed by the external forces. This is due to the release of
strain energy stored in the solid.
Victor Saouma Matrix Structural Analysis
Draft
9.2 Work, Energy & Potentials; Denitions 9{11
σx σx
C,D C D
σ x dε x σ x dε x
A B εx A,B εx
ε0 εx ε0
εx
(εx)F
(εx)F
(a) (b)
U = ;Wi (9.45)
43 The internal work depends on the load history, this is illustrated by considering an axial
member subjected to two load cases, Fig. 9.3: a) Initial thermal strains (with no corresponding
stress increase), followed by an external force; and b) External force, followed by thermal strain.
In both cases the internal work is equal to the area under the curve ABCD.
Z !
L Z ("x )F
Ui = x d"x Adx (9.46-a)
0 0
Wia = ;Ui !
(9.46-b)
Z L Z "0
Wib = ;Ui ; xd"x Adx (9.46-c)
0 0
Hence, Wi is not always equal to ;Ui .
9.2.3 External Work
44 External work W performed by the applied loads on an arbitrary system is dened as
Z Z
We def
= uT bd
+ uT ^td; (9.47)
;t
Z f Z f
We = Pd + Md (9.48)
0 0
We =
Z f
Pd ; e W = K (9.49)
0
0
When this last equation is combined with Pf = K f we obtain
We = 21 Pf f (9.50)
We = 12 Mf f (9.51)
dW = Fxdx + Fy dy (9.52)
50 From calculus, a necessary and sucient condition for dW to be an exact dierential is that
@Fx = @Fy (9.53)
@y @x
51If the force were to move along a closed contour (or from A to B and then back to A along
any arbitrary path), corresponding to ;, then from Green's theorem (Eq. 5.3) we have
Z
I
(Rdx + Sdy) = @S ; @R dxdy (9.54)
; @x @y
Victor Saouma Matrix Structural Analysis
Draft
9.2 Work, Energy & Potentials; Denitions 9{13
If we let R = Fx and S = Fy , then
Z
I
W = (Fx dx + Fy dy) = @Fy ; @Fx dxdy (9.55)
; | @x {z @y }
0
where all the terms have been previously dened and b is the body force vector.
9.2.4.1 Internal Virtual Work
57 Next we shall derive a displacement based expression of U for each type of one dimensional
structural member. It should be noted that the Virtual Force method would yield analogous
ones but based on forces rather than displacements.
58 Two sets of solutions will be given, the rst one is independent of the material stress strain
relations, and the other assumes a linear elastic stress strain relation.
Victor Saouma Matrix Structural Analysis
Draft
9{14 VARIATIONAL AND ENERGY METHODS
Shear Members:
Z 9
U = xy
xy d
>
>
>
> Z L Z Z L
Z
=
V = xy dA U = ( xy dA)
xy dx ) U = V
xy dx (9.62)
A
>
>
> 0 | A {z } 0
>
d
= dAdx ;
V
d
= Adx >
; \00 \"00
Torsional Members: With reference to Fig. 9.4
Z Z
U = xy
xy d
= r d
(9.65)
or: 9
Eq.
Z 9.68 = Z L d2 v d2 (v) dx
y2 dA = Iz U = EIz dx 2 dx2 (9.69)
; 0 | {z } | {z }
A
\00 \"00
Shear Members:
Z 9
U = xy
xy d
>
>
>
> Z L Z Z L
Z
=
U = ( xy dA)
xy dx ) U = V
xy dx (9.75)
V = xy dA >
> 0 | A {z 0
A >
> }
d
= dAdx ;
V
uT bd
+ ;t
uT ^td; + uP (9.83)
where u are the displacements, b is the body force vector; ^t is the applied surface traction
vector; ;t is that portion of the boundary where ^t is applied, and P are the applied nodal
forces.
70 Note that the potential of the external work is dierent from the external work itself (usually
by a factor of 1/2)
9.2.6.3 Potential Energy
71 The potential energy of a system is dened as
def
= UZ ; We Z Z
(9.84)
= U0 d
; ubd
+ u^td; + uP (9.85)
;t
72 Note that in the potential the full load is always acting, and through the displacements of
its points of application it does work but loses an equivalent amount of potential, this explains
the negative sign.
1. The second approach is the same one on which the method of virtual or unit load is based.
It is simpler to use than the third as a internal force distribution compatible with the
assumed virtual force can be easily obtained for statically determinate structures. This
approach will yield exact solutions for staticaally determinate structures.
2. The third approach is favoured for kinematically indeterminate problems or in conjunction
with approximate solution. It requires a proper \guess" of a displacement shape and is
the basis of the stifness method.
9.3.1 Principle of Virtual Work
9.3.1.1 y Derivation
75 Derivation of the principle of virtual work starts with the assumption of that forces are in
equilibrium and satisfaction of the static boundary conditions.
76 The Equation of equilibrium (Eq. 8.11) which is rewritten as
or
LT + b = 0 (9.89)
Z
Z
; uT td; ; uT ^td; = 0 (9.97-a)
Z Z
;
Z u
;
Z t
; uT bd
+ divuT d
; uT td; ; uT ^td; = 0 (9.97-b)
;u ;t
84 Virtual displacement must be kinematically admissible, i.e. u must satisfy the essential
boundary conditions u = 0 on ;u , (note that the exact solution had to satisfy the natural B.C.
instead), hence the previous equation reduces to
Z Z Z
T d
; uT bd
; uT ^td; = 0 (9.100)
|
{z }|
{z
;t }
;Wi =Ui ;We
Each of the preceding equations is a work expression, (Eq. 9.59). The rst one corresponds to
the internal virtual work, and the last two are expressions of the work done by the body forces
and the surface tractions through the corresponding virtual displacement u, hence
or
Ui = We (9.102)
which is the expression of the principle of virtual work (or more specically of virtual
displacement) which can be stated as
A deformable system is in equilibrium if the sum of the external virtual work and
the internal virtual work is zero for virtual displacements u which are kinematically
admissible.
The major governing equations are summarized
Victor Saouma Matrix Structural Analysis
Draft
9.3 Principle of Virtual Work and Complementary Virtual Work 9{23
Figure 9.6: Tapered Cantilivered Beam Analysed by the Vitual Displacement Method
Z Z Z
T d
; uT bd
; uT ^td; = 0 (9.103)
|
{z }|
{z
;t }
;Wi ;We
= Lu in
(9.104)
u = 0 on ;u (9.105)
85 Note that the principle is independent of material properties, and that the primary unknowns
are the displacements.
86 For one dimensional elements, with no initial strains (U = ;Wi )
Z
"d
= Pv
|{z} (9.106)
| {z } W
U
Note that both Eqn. 9.107 and Eqn. 9.108 satisfy the essential (geometric) B.C.
Z L
U = v00 EIz v00 dx (9.109)
0
W = P2v2 (9.110)
Solution 1:
U =
Z L 2
cos x v 6 ; 12x v EI 1 ; x dx
0 4L2 2l 2 L2 L3 2 1 2L
= 3EI1 1 ; 10 + 16 v v
2L3 2 2 2
= P2 v2 (9.111)
which yields:
P2 L3
v2 = 2:648 (9.112)
EI 1
Solution 2:
U =
Z L4 cos2 x v v EI 1 ; x dx
0 16L4 2l 2 2 1 2l
4
= 32EI 1 3 1
L3 4 + 2 v2 v2
= P2 v2 (9.113)
which yields:
3
v2 = 2:P572 LEI (9.114)
1
96 Note that the principle is independent of material properties, and that the primary unknowns
are the stresses.
97 The principle of virtual forces leads to the
exibility matrix.
Figure 9.7: Tapered Cantilevered Beam Analysed by the Virtual Force Method
From Mathematica we note that:
Z 0
x2 = 1 1 (a + bx)2 ; 2a(a + bx) + a2 ln(a + bx) (9.135)
0 a + bx b3 2
Thus substituting a = L and b = 1 into Eqn. 9.135, we obtain:
= 2EI P2 L 1 (L + x)2 ; 2L(L + x) + L2 ln(L + x) jL
0
1 "2 #
2 P 2 L L 2
= EI 2L ; 4L + L ln 2L ; 2 + 2L + L log L
2 2 2 2 2
1
2 P2 L
= EI L (ln 2 ; 2 )
2 1
1
= 2:5887P 2 L3 (9.136)
1EI
This exact value should be compared with the approximate one obtained with the3 Virtual
Displacement method in which a displacement eld was assumed in Eq. 9.215 of 2:PL
55EI1 .
Similarly:
=
Z L
M; (1)
0 EI1 :5 + Lx
= EI2 ML Z L
1
1 0 L+x
2 ML
= EI ln(L + x) jL0
1
= 2 ML (ln 2L ; ln L)
EI1
ML ln 2
= 2EI
1
ML
= :721 (9.137)
EI1
4
= :0337 wR
EI (9.146)
103 Hence, comparing this last equation, with Eq. 9.85 we obtain
= 0 (9.159)
def
= UZ ; We Z (9.160)
Z
= U0d
; ubd
+ u^td; + uP (9.161)
;t
104 We have thus derived the principle of stationary value of the potential energy:
Of all kinematically admissible deformations (displacements satisfying the essential
boundary conditions), the actual deformations (those which correspond to stresses
which satisfy equilibrium) are the ones for which the total potential energy assumes
a stationalry value.
2 Note that the variation of strain energy density is, U0 = ij "ij , and the variation of the strain energy itself
R
is U =
U0 d
.
k= 500 lbf/in
106 It can be shown that the minimum potential energy yields a lower bound prediction of
displacements.
107 As an illustrative example (adapted from Willam, 1987), let us consider the single dof system
shown in Fig. 9.10. The strain energy U and potential of the external work W are given by
U = 21 u(Ku) = 250u2 (9.163-a)
We = mgu = 100u (9.163-b)
Thus the total potential energy is given by
= 250u2 ; 100u (9.164)
and will be stationary for
@ = d = 0 ) 500u ; 100 = 0 ) u = 0:2 in
du (9.165)
Substituting, this would yield
U = 250(0:2)2 = 10 lbf-in
W = 100(0:2) = 20 lbf-in (9.166)
= 10 ; 20 = ;10 lbf-in
Fig. 9.11 illustrates the two components of the potential energy.
Victor Saouma Matrix Structural Analysis
Draft
9.4 Potential Energy 9{35
0.0
−20.0
−40.0
0.00 0.10 0.20 0.30
Displacement [in]
Z Z
; uT bd
; uT ^td; = 0 (9.170)
;t
is better suited for obtaining the discrete system of equations.
112 To obtain the Euler equations for the general form of the potential energy variational prin-
ciple the volume integrals dening the virtual strain energy U in Equation 9.169 must be
integrated by parts in order to convert the variation of the strains (Lu) into a variation of the
displacements u.
113 Integration by parts of these integrals using Green's theorem (Kreyszig 1988) yields
Z I Z
(Lu)T Dd
= uT G(D)d; ; uT LT (D)d
Z
I@
Z
(Lu)T D0 d
= u G(D0 )d; ; uT LT (D0 )d
T (9.171)
Z I@
Z
(Lu)T 0 d
= uT G0 d; ; uT LT 0 d
where G is a transformation matrix containing the direction cosines for a unit normal vector
such that the surface tractions t are dened as t = G and the surface integrals are over the
entire surface of the body @
.
114 Substituting Equation 9.171 into Equation 9.169, the variational statement becomes
Z
= ; uT fLT [D( ; 0 ) + 0 ] + bgd
Z
+ uT fG[D( ; 0 ) + 0 ] + ^tgd; = 0 (9.172)
@
115 Since u is arbitrary the expressions in the integrands within the braces must both be equal
to zero for to be equal to zero. Recognizing that the stress-strain relationship (i.e. Equation
9.37) appears in both the volume and surface integrals, the Euler equations are
LT + b = 0 on
(9.173)
G ; ^t = 0 on ;t
where the rst Euler equation is the equilibrium equation and the second Euler equation denes
the natural boundary conditions. The natural boundary conditions are dened on ;t rather
than @
because both the applied surface tractions ^t and the matrix-vector product G are
identically zero outside ;t .
= ; @@U
1 1 ; @U ; ; @U
@ 2 2 @ n n
+P^1 1 + P^2 2 + P^n n = 0 (9.176-a)
or
; @@U ^ @U ^ @U ^
1 + P1 1 + ; @ 2 + P2 2 + + ; @ n + Pn n = 0 (9.177)
but since the variation i is arbitrary, then each factor within the parenthesis must be equal
to zero. Thus
@U ^ (9.178)
@ k = Pk
which is Castigliano's rst theorem:
If the strain energy of a body is expressed in terms of displacement components in
the direction of the prescribed forces, then the rst partial derivative of the strain
energy with respect to a displacement, is equal to the corresponding force.
v = a1 x3 + a2 x2 + a3 x + a4 (9.179)
1. First, this solution must satisfy the essential B.C.: v = v0 = 0 at x = 0; and v = vmax and
v0 = 0 at x = L2 . This will be enforced by determining the four parameters in terms of a
single unknown quantity (4 equations and 4 B.C.'s):
@x = 0 v = 0 ) a4 = 0
@x = 0 dx dv = 0 ) a3 = 0
(9.180)
@x = L2 v = vmax ) vmax = a1 L83 + a2 L42
@x = L2 dxdv = 0 ) 43 a1L2 + a2L = 0 ) a2 = ; 34 a1 L
upon substitution, we obtain:
!
3 2
v = ; 16Lx3 + 12Lx2 vmax (9.181)
Hence, in this problem the solution is in terms of only one unknown variable vmax.
2. In order to apply the principle of Minimum Potential Energy we should evaluate:
2
Internal Strain Energy U : for
exural members is given by U = 2MEI dx (Eq. 9.43);
Z
z
M = d2 v2 , thus we must evaluate d2 v2 from above:
recalling that EI z dx dx
!
dv 2
dx = ; 48Lx3 + 24L2x vmax (9.182)
Victor Saouma Matrix Structural Analysis
Draft
9.4 Potential Energy 9{39
d2 v = ; 24 1 ; 4x v (9.183)
dx2 L2 L max
which yields 2 !2 3
1 Z
d2 v
U = 2 2 E dx
4 Iz dx5 (9.184)
2
or:
U = E Z L4 242 1 ; 4x 2 v2 Iz dx
2 2 0 L4 L max 2
Z L 2
E 2 242 4 x
+ 2 L L4 1 ; L vmax 2 I dx
z
4
U = 72LEI z 2
3 vmax (9.185-a)
Figure 9.13: Uniformly Loaded Simply Supported Beam Analysed by the Rayleigh-Ritz Method
=U ;W =
Z L M 2 dx ; Z Lwv(x)dx (9.193)
o 2EIz 0
d2 v2 , the above simplies to:
Recalling that: EIMz = dx
2 ! 3
=
Z L
4
EIz d2 v 2 ; wv(x)5 dx (9.194)
0 2 dx2
Z L EIz
= (;2a1 )2 ; a1 wx(L ; x) dx
0 2
3 3
= EI2 z 4a21 L ; a1 w L2 + a1 w L3
3
= 2a21 EIz L ; a1 wL
6 (9.195)
@ = 0, we would obtain:
If we now take @a 1
3
4a1 EIz l ; wL
6 = 0
2
a1 = 24wL EIz (9.196)
Having solved the displacement eld in terms of a1 , we now determine vmax at L2 :
!
4 x x2
v = 24wL
EI L ; L2
| {z z}
a1
Victor Saouma Matrix Structural Analysis
Draft
9.4 Potential Energy 9{43
4
= 96wL
EIz (9.197)
exact = 5 wL4 = wL4 which constitutes
This is to be compared with the exact value of vmax 384 EIz 76:8EIz
17% error.
Note: If two terms were retained, then we would have obtained: a1 = 24wLEI2z & a2 = 24wEIz
exact . (Why?)
and vmax would be equal to vmax
2l cos 2l L m=n
(9.209)
0 :
2
8
0 m 6= n
x cos mx nx dx =
Z L <
2l cos 2l L2 ; L22 2 m = n
(9.210)
0 :
8 2n
Thus combining Eqns. 9.208, 9.209, and 9.210, we obtain:
4 1X
U = 64EI 3 + 1 n4 a2 (9.211)
L3 1;3;5 4 n22 n
n(;1) n;2 1 an
X
Similarly we solve for 2 = 2L
8
< PL2 + Ml n=1
1:65EI 1:05EI1
2 : PL2 + ML n=3 (9.216)
1:65EIz 1:04EI1
129 The complementary strain energy can also be expressed in terms of the forces Pi thus the
complementary potential energy will be dened in terms of generalized coordinates or generalized
forces.
130 For the solid to be in equilibrium, = 0 or
= @W
@P1
i P1 + @Wi P2 + + @Wi Pn
@P2 @Pn
;^1P1 ; ^2 P2 ; ^nPn = 0 (9.221-a)
or
@Wi ; ^ P + @Wi ; ^ P + + @Wi ; ^ P (9.222)
@P1 1 1 @P2 2 2 @Pn n n
= @U
@R R=0
R L2 2 R L2
h i
2 @M
= EI 0 M @R dx = EI 0
wL + R
2 2 x ; w x22 x2 dx
4 R=0 R=0
5wL
= 384 EI
9.7 Summary
135 A summary of the various methods introduced in this chapter is shown in Fig. 9.15.
136 The duality between the two variational principles is highlighted by Fig. 9.16, where be-
ginning with kinematically admissible displacements, the principle of virtual work provides
Victor Saouma Matrix Structural Analysis
Draft
9{50 VARIATIONAL AND ENERGY METHODS
Natural B.C.
Essential B.C.
? ? ? ?
div + b = 0 ; Du = 0
"ij ; 1 (ui;j
2 + uj;i ) = 0 ij;j = 0
t ; bt = 0 ;t u = 0 ;u ; ui ; ub = 0 ;u ti = 0 ;t
U0 def
= 0" d" U0 def
= 0 "d
R R
6 Gauss 6Gauss
? ? ? ?
Principle of Virtual Work Principle of Complementary
R VirtualR Work
ij ij d
; ;u ubi ti d; = 0
R R R
d
;
u bd
; ;t u btd; = 0 "
T T T
Wi ; We = 0 Wi ; We = 0
?
Principle of Stationary
?
Principle of Complementary
Potential Energy Stationary Potential Energy
= 0 = 0
def def We
R R U ; We R
= R = Wi + R
=
U0 d
; (
ui bi d
+ ;t ui tbi d;) =
U0 d
+ ;u ubi ti d;
?
Castigliano's Castigliano's
?
First Theorem Second Theorem
@Wi = Pc @Wi = c
@ k k @Pk k
?
Rayleigh-Ritz
n
X
uj cji ji + j0
i=1
@ = 0 i = 1; 2; ; n; j = 1; 2 ; 3
@cji
Figure 9.15: Summary of Variational Methods
?
Principle of Stationary Principle of Virtual Work
Complementary Energy
Principle of Complementary Principle of Stationary
Virtual Work Potential Energy
6
?
Statically Admissible Stresses
Stresses satisfy the equilibrium conditions
and the static boundary conditions
statically admissible solutions. Similarly, for statically admissible stresses, the principle of
complementary virtual work leads to kinematically admissible solutions.
137 Finally, Table 9.4 summarizes some of the major equations associated with one dimensional
rod elements.
Table 9.4: Summary of Variational Terms Associated with One Dimensional Elements
Chapter 10
INTERPOLATION FUNCTIONS
10.1 Introduction
27Application of the Principle of Virtual Displacement requires an assumed displacement eld.
This displacement eld can be approximated by interpolation functions written in terms of:
1. Unknown polynomial coecients, most appropriate for continuous systems, and the Rayleigh-
Ritz method
y = a1 + a2 x + a3x2 + a4 x3 (10.1)
A major drawback of this approach, is that the coecients have no physical meaning.
2. Unknown nodal deformations, most appropriate for discrete systems and Potential Energy
based formulations
y = = N1 1 + N2 2 + : : : + Nn n (10.2)
28For simple problems both Eqn. 10.1 and Eqn. 4y
10.2 can readily provide the exact solutions of
d
the governing dierential equation (such as dx4 = EIq for
exure), but for more complex ones,
one must use an approximate one.
10.2.2 Generalization
36 The previous derivation can be generalized by writing:
( )
u = a1 x + a2 = b x 1 c aa1 (10.14)
| {z } 2
[p] | {z }
fag
where [p] corresponds to the polynomial approximation, and fag is the coecient vector.
37 We next apply the boundary conditions:
( ) " #( )
u1 = L0 11 a1 (10.15)
u2 a2
| {z } | {z }| {z }
fg [L] fag
following inversion of [L], this leads to
( ) " #( )
a1
a2 = L1 ;L1 10 u1
u2 (10.16)
| {z } | {z }| {z }
fag [L];1 f g
10.2.3 Flexural
40 With reference to Fig. 10.2. We have 4 d.o.f.'s, fg41 : and hence will need 4 shape
functions, N1 to N4 , and those will be obtained through 4 boundary conditions. Therefore we
need to assume a polynomial approximation for displacements of degree 3.
v = a1 x3 + a2 x2 + a3 x + a4 (10.19)
= dx dv = 3a x2 + 2a x + a (10.20)
1 2 3
0.8
N1
0.6
N3
N2
N4
0.4
N
0.2
0.0
−0.2
0.0 0.2 0.4 0.6 0.8 1.0
ξ(x/L)
=0 =1
Function Ni Ni;x Ni Ni;x
N1 = (1 + 23 ; 32 ) 1 0 0 0
N2 = x(1 ; )2 0 1 0 0
N3 = (32 ; 2 3 ) 0 0 1 0
N4 = x(2 ; ) 0 0 0 1
48 As before, we rst seek the shape functions, and hence we apply the boundary conditions at
the nodes for the u displacements rst:
8 9 2 38 9
>
< u1 >
= 1 0 0 >
< 1 a >
=
>
u2 >
= 4 1 x2 0
6 7
5
> 2 a >
(10.33)
: u3 ; 1 x3 y3 :
3 a ;
| {z } | {z } | {z }
fg [L] fag
49 We then multiply the inverse of [L] in Eq. 10.33 by [p] and obtain:
u = N1 u1 + N2 u2 + N3 u3 (10.34)
where
N1 = x 1y (x2 y3 ; xy3 ; x2 y + x3 y)
2 3
Victor Saouma Matrix Structural Analysis
Draft
10{8 INTERPOLATION FUNCTIONS
v = N1 v1 + N2 v2 + N3 v3 (10.36)
58 For the axial member, m = 1, x1 = 0, and x2 = L, the above equations will result in:
= (x ; L) + x = (1 ; x ) + x
;L 1 L 2 L} 1 L (10.40)
2
| {z |{z}
N1 N2
which is identical to Eq. 10.12.
10.3.1.1 Constant Strain Quadrilateral Element
59 Next we consider a quadrilateral element, Fig. 10.5 with bi-linear displacement eld (in both
x and y).
60 Using the Lagrangian interpolation function of Eq. 10.38, and starting with the u displace-
ment, we perform two interpolations: the rst one along the bottom edge (1-2) and along the
top one (4-3).
61 From Eq. 10.38 with m = 1 we obtain:
62 Similarly
u43 = xx2 ;;xx u4 + xx1 ;;xx u3
2 1 1 2
a ; x x
= 2a u4 + 2a u3 + a (10.42)
Victor Saouma Matrix Structural Analysis
Draft
10{10 INTERPOLATION FUNCTIONS
63Next, we interpolate in the y direction along 1-4 and 2-3 between u12 and u43 . Again, we
use Eq. 10.38 however this time we replace x by y:
u = yy2 ;;yy u12 + yy1 ;;yy u43 (10.43)
2 1 1 2
= b 2;b y a 2;a x u1 + b 2;b y x 2+a a u2 + y 2+b b a 2;a x u4 + y 2+b b x 2+a a u3
= (a ; x4)(abb ; y) u1 + (a + x4)(abb ; y) u2 + (a + x4)(abb + y) u3 + (a ; x4)(abb + y) u4
| {z } | {z } | {z } | {z }
N1 N2 N3 N1
64 One can easily check that at each node i the corresponding Ni is equal to 1, and all others
to zero, and that at any point N1 + N2 + N3 + N4 = 1. Hence, the displacement eld will be
given by: 8 9
>
>
>
u1 >
>
>
>
>
>
> v1 >
>
>
>
( ) "
>
>
#>
>
>
u2 >
>
>
>
>
u = N01 N0 N02 N0 N03 N0 N04 N0
<
v2 =
(10.44)
v 1 2 3 4 >
>
>
u3 >
>
>
>
>
>
>
>
v3 >
>
>
>
>
>
>
>
>
u4 >
>
>
>
:
v4 ;
where
Ni = (a x)(b8abc
y)(c z) (10.46)
Table 10.2: Interpretation of Shape Functions in Terms of Polynomial Series (1D & 2D)
Element Terms # of Nodes
(terms)
Linear a1 , a2 2
Quadratic a1 , a2 , a4 3
Bi-Linear (triangle) a1 , a2 , a3 , 3
Bi-Linear (quadrilateral) a1 , a2 , a3 , a5 4
Bi-Quadratic (Serendipity) a1 , a2 , a3 , a4 , a5 , a6 , a8 , a9 8
Bi-Quadratic (Lagrangian) a1 , a2 , a3 , a4 , a5 , a6 , a8 , a9 , a13 9
Table 10.3: Polynomial Terms in Various Element Formulations (1D & 2D)
67 Hermitian interpolation functions are piecewise cubic functions which satisfy the conditions
of displacement and slope (C 0 , C 1 ) continuities. They are exensively used in CAD as Bezier
curves.
Chapter 11
FINITE ELEMENT
FORMULATION
27 Having introduced the virtual displacement method in chapter 9, the shape functions in
chapter 10, and nally having reviewed the basic equations of elasticity in chapter 8, we shall
present a general energy based formulation of the element stiness matrix in this chapter.
28 Whereas chapter 2 derived the stiness matrices of one dimensional rod elements, the ap-
proach used could not be generalized to general nite element. Alternatively, the derivation of
this chapter will be applicable to both one dimensional rod elements or contnuum (2D or 3D)
elements.
35 Let us now apply the principle of virtual displacement and restate some known relaations:
U = W
Z
(11.14)
U = bcf gd
(11.15)
36 Combining Eqns. 11.14, 11.15, 11.16, 11.19, and 11.17, the internal virtual strain energy
is given by:
Z Z
U = b{zc[B]T} [|D][B{z]fg} d
; b|{zc[B]T} [|D]{zf0g} d
fg f g fg f 0 g
Z Z
= bc [B]T [D][B] d
fg ; bc [B]T [D]f"i gd
(11.20)
Chapter 12
SOME FINITE ELEMENTS
12.1 Introduction
27 Having rst introduced the method of virtual displacements in Chapter 9, than the shape
functions [N] (Chapter 10) which relate internal to external nodal displacements, than the
basic equations of elasticity (Chapter 8) which dened the [D] matrix, and nally having
applied the virtual displacement method to nite element in chapter 11, we now revisit some
one dimensional element whose stiness matrix was earlier derived, and derive the stiness
matrices of additional two dimensional nite elements.
Draft
12{2 SOME FINITE ELEMENTS
sectional area we obtain: ( )
[k] = A
Z L ; L1 E b ; 1 1 cdx
1 L L
0 L
" #
[k] = AE
Z L
1 ;1
L 2 0 ;1 1 dx
" # (12.2)
= AE 1 ;1
L ;1 1
31 We observe that this stiness matrix is identical to the one earlier derived in Eq. 2.45.
Which is identical to the beam stiness matrix derived in Eq. 2.45 from equilibrium relations.
Victor Saouma Matrix Structural Analysis
Draft
12.4 Triangular Element 12{3
12.4 Triangular Element
35 Having retrieved the stiness matrices of simple one dimensional elements using the principle
of virtual displacement, we next consider two dimensional continuum elements starting with
the triangular element of constant thickness t made out of isotropic linear elastic material. The
element will have two d.o.f's at each node:
n o
= b u1 u2 u3 v1 v2 v3 ct (12.8)
2 3
; x12 0 x3 ;x2
x2 y3
6 1
x2 0 ; x;2xy33 7 " #
Z 6
6 0 0 1 7
7 E 1 0
= 6 x3 ;x2 y3 7 1 0
6
6 0 x2 y3 ; x12 7 1 ; 2
7 0 0 1; )
4 0 ; x;2xy33 1
x2 5| {z 2 }
0 1 0 [D]
|
y{z
3
}
[B]T
2 3
; x12 1
x2 0 0 0 0
0 0 0 x3 ;x2 ; x;2xy33 1 (12.12)
4
x3 ;x2 x2 y3 y3 5 tdxdy
x2 y3 ; ;x3
x2 y3
1 ;
y3 {z x2
1 1
x2 0 |dvol{z }
| }
[B]
2 3
y32 + x23;2 ;y32 2; x3 x23;2 x2 x3;2 ;y3 x3;2 x3 y3 + y3 x3;2 ;x2 y3
6 ;y32 ; x3 x3;2 y 3 + x 3 ; x 2 x 3 y3 x 3; 2 + x 3 y 3 ;x3 y3 x2 y3 7
6 7
6 x
=
6 ;y3 x3;2 2 x 3; 2 ; x 2 x3 x 2
2 ; x 2 y 3 x 2 y 3 0 7
6 y3 x3;2 + x3 y3 ;x2 y3 y32 + x23;2 ;y32 ; x3 x3;2 x2 x3;2 7
7
4 x3 y3 + y3 x3;2 ;x3 y3 x2 y3 ;y32 ; x3 x3;2 y32 + x23 ;x22x3 5
;x2y3 x2 y3 0 x2 x3;2 ;x2 x3 x2
where = 1;2 , = 1+2 ,
= 2(1;ET 2 )x2 y3 , x3;2 = x3 ; x2 , and y3;2 = y3 ; y2 .
Chapter 13
GEOMETRIC NONLINEARITY
13.1 Strong Form
27 Column buckling theory originated with Leonhard Euler in 1744.
28 An initially straight member is concentrically loaded, and all bers remain elastic until
buckling occur.
29 For buckling to occur, it must be assumed that the column is slightly bent as shown in Fig.
13.1. Note, in reality no column is either perfectly straight, and in all cases a minor imperfection
P P
x
x and y are
principal axes
x
35The fundamental buckling mode, i.e. a single curvature de
ection, will occur for n = 1; Thus
Euler critical load for a pinned column is
2
Pcr = LEI
2 (13.7)
41 Neglecting the terms in dx2 which are small, and then dierentiating each term with respect
to x, we obtain
d2 M ; dV ; P d2 v = 0 (13.10)
dx2 dx dx2
42 However, considering equilibrium in the y direction gives
dV = ;w (13.11)
dx
43 From beam theory, neglecting axial and shear deformations, we have
d2 v
M = ;EI dx (13.12)
2
44Substituting Eq. 13.11 and 13.12 into 13.10, and assuming a beam of uniform cross section,
we obtain
4 2
EI d v4 ; P d v2 = w
dx dx
(13.13)
45 P , the general solution of this fourth order dierential equation to any set
Introdcing k2 = EI
of boundary conditions is
v = C1 sin kx + C2 cos kx + C3 x + C4 (13.14)
46 If we consider again the stability of a hinged-hinged column, the boundary conditions are
v = 0; v;xx = 0 at x = 0 (13.15)
v = 0; v;xx = 0 at x = L
Victor Saouma Matrix Structural Analysis
Draft
13{4 GEOMETRIC NONLINEARITY
w(x)
P P
x
dx
y,u
w
M
δV
P i V+ δx dx
V
δv P
δx
dx
M+ δM
δx dx
θi i P
P
j
P
θj
dx
Figure 13.2: Simply Supported Beam Column; Dierential Segment; Eect of Axial Force P
8.0
6.0
4.0
2.0
0.0
-2.0
-4.0
-6.0
-8.0
-10.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0
Figure 13.3: Solution of the Tanscendental Equation for the Buckling Load of a Fixed-Hinged
Column
13.2 Weak Form
13.2.1 Strain Energy
51 Considering a uniform section prismatic element, Fig. ??, subjected to axial and
exural
deformation (no shear), the Lagrangian nite strain-displacement relation is given by ??
"xx = u;x + 21 (u2;x + v;x2 + w;x2 ) (13.22)
thus, the total strain would be
du d2v ! 1
dv 2
"xx = dx ; y dx2 + 2 {zdx } (13.23)
|{z} |
Axial |Flexure
{z }
Large Deformation
52 we note that the rst and second terms are the familiar components of axial and
exural
strains respectively, and the third one (which is nonlinear) is obtained from large-de
ection
strain-displacement.
53 The Strain energy of the element is given by
Z
e 1
U = 2 E"2xx d
(13.24)
i.p.
l/4 i.p.
i.p.
<kl<l
kl=l kl= l l l
2
2
l
i.p.
l/4
i.p.
i.p.
l l<kl<
8
i.p. l
l
kl=21
kl=1
Pcr
Pcr Pcr
i.p.
l<kl<
8
i.p.
k=2 k=1
Figure 13.4: Critical lengths of columns
55 Noting that Z Z Z
dA = A; ydA = 0; y2 dA = I (13.26)
A A A
for y measured from the centroid, U e reduces to
2 !2 3
2
Z
U e = 21 E 4A du + I d2 v
dv 4 + A du dv 2 5 dx
+ A4 dx (13.27)
L dx dx2 dx dx
4
We discard the highest order term A4 dx
dv in order to transform the above equation into a
linear instability formulation.
56 Under the assumption of an independent prebuckling analysis for axial loading, the axial
load Px is
du
Px = EA dx (13.28)
Thus Eq. 13.27 reduces to
2 ! 3
U e = 21
Z
4
du + EI
EA dx
2
d2 v 2 + P dv 25 dx (13.29)
dx2 x dx
L
57We can thus decouple the strain energy into two components, one associated with axial and
the other with
exural deformations
U e = Uae + Ufe (13.30-a)
2
Uae = 21 du dx
EA dx
Z
(13.30-b)
L2 3
1 Z
d 2 v !2
dv 2
Ufe = 2 4EI dx2 + Px dx 5 dx (13.30-c)
L
Substituting into the Euler equation, and assuming constant Px , and EI , we obtain
d4 v ; P d2 v = 0
EI dx (13.36)
4 x dx2
61 Substituting this last equation into Eq. 13.30-c, the element potential energy is given by
e = Ufe + W e (13.38-a)
1 1
= 2 bve c[ke ] fve g + 2 bve c[kg ] fveg ; bvc fPg (13.38-b)
where Z
[ke ] = EI fN;xxg bN;xxcdx (13.39)
L
and
Z
[kg ] = P fN;xg bN;xcdx (13.40)
L
where [ke ] is the conventional element
exural stiness matrix.
62 [kg ] introduces the considerations related to elastic instability. We note that its terms solely
depend on geometric parameters (length), therefore this matrix is often referred to as the
geometric stiness matrix.
63 Using the shape functions for
exural elements, Eq. 10.29, and substituting into Eq. 13.39
and Eq. 13.40 we obtain
2
u1 v1 1 u2 v2 2 3
EA
L 0 0 ; EA
L 0 0
6
6 0 12EI
L3
6EI
L2 0 ; 12LEI3 6EI
L2
7
7
6
6
ke = ; EA
6
0 6EI
L2
4EI
L 0 ; 6LEI2 2EI
L
7
7
7 (13.41)
6 0 0 EA 0 0 7
6 L L 7
6
40 ; 12LEI 3 ; 6LEI2 0 12LEI3 ; 6EI 7
L2 5
0 6EI
L2
2EI
L 0 ; 6LEI2 4EI
L
which is the same element stiness matrix derived earlier in Eq. 12.7.
64 The geometric stiness matrix is given by
2 u1 v1 1 u2 v2 2 3
0 06 L0 0 0 0
6
6
6
0 5 10 0 ; 65 10L 2 777
0 10L 152 L2 0 ; 10L ; L30 77
kg = PL (13.42)
6
6
6 0 06 0L 0 0 0 7
6
6
6
0 ; 5 ; 10 0 65 ; 10L 777
4 0 10L ; L302 0 ; 10L 152 L2 5
K = Ke + Kg (13.45)
67 We assume that conservative loading is applied, that is the direction of the load does not
\follow" the de
ected direction of the member upon which it acts.
73 Alternatively, it can simply be argued that there is no unique solution (bifurcation condition)
to v.
74 The lowest value of , crit will give the buckling load for the structure and the buckling
loads will be given by
Pcrit = critP (13.50)
75 The corresponding deformed shape is directly obtained from the corresponding eigenvector.
2
3
(1) l
5
6 4
l
(2)
8
9
Solution:
2
1 2 3 4 5 6 3
EA
L 0 0 ; EA
L 0 0
6
6 0 12EI
L3
6EI
L2 0 ; 12LEI3 6EI
L2
7
7
k1e =
6
6
6
0 6EI
L2
4EI
L 0 ; 6LEI2 2EI
L
7
7
7 (13.51-a)
6
6 ; EA
L 0 0 EA
L 0 0 7
7
6
4 0 ; 12LEI3 ; 6LEI2 0 12EI
L3 ; 6EI 7
L2 5
0 6EI
L2
2EI
L 0 ; 6LEI2 4EI
L
2
4 5 6 7 8 9 3
EA
L 0 0 ; EA
L 0 0
6
6 0 12EI
L3
6EI
L2 0 ; 12LEI3 6EI
L2
7
7
k2e =
6
6
6
0 6EI
L2
4EI
L 0 ; 6LEI2 2EI
L
7
7
7 (13.51-b)
6
6 ; EA
L 0 0 EA
L 0 0 7
7
6
4 0 ; 12LEI3 ; 6LEI2 0 12EI
L3 ; 6EI 7
L2 5
0 6EI
L2
2EI
L 0 ; 6LEI2 4EI
L
Similarly, the geometric stiness matrices are given by
2
1 2 3 4 5 6 3
0 0 0 0 0 0
60
6
6
5
L
10 0 ; 65 10L 2 777
k1g = ;LP ; 10L ; L30 77
6 L 2 2
60
6 10 15 L 0 (13.52-a)
60 0 0 0 0 0 7
6 7
4 0 ;5 ; 10L2 ; 10L
6 6 0 6 7
5 5
0 10L ; L30 0 ; 10L 2 2
15 L
2
4 5 6 7 8 9 3
0 65 L
10 0 ; 65 L
10 2
6
60 L 2 2
15 L 0 ; 10L ; L30 7
k2g = ;LP
7
6 10 7
60 0 0 0 0 0 7 (13.52-b)
6 7
6 0 ;6 ; 10L2 0 6 ; 10L 7
4 5 5 5
0 10L ; L30 0 ; 10L 2 2
15 L
The structure's stiness matrices Ke and Kg can now be assembled from the element sti-
nesses. Eliminating rows and columns 2, 7, 8, 9 corresponding to zero displacements in the
column, we obtain
1 4 3 5 6
I 2 ; I2
2 AL2 AL 2
0 0 0 3
6; I
6 AL
EI 2 ALI 0 0 0 777
Ke = L3 66 0
6
0 4L2 ;6L 2L2 77 (13.53)
4 0 0 ;6L 24 0 5
0 0 2L2 0 8L2
Victor Saouma Matrix Structural Analysis
Draft
13{14 GEOMETRIC NONLINEARITY
and
1 4 3 5 6
2
0 0 0 0 0 3
60 0 0 0 0 2 77
; P 6
; ;
Kg = L 66 0 0 15;LL 1012 30L 777
6 2 2 L (13.54)
40 0 10 5 0 5
0 0 ;30L2 0 154 L2
noting that in this case Kg = Kg for P = 1, the determinant jKe + Kg j = 0 leads to
1 4 3 5 6
1
AL2
I ; ALI 22 0 0 0
4
; ALI 2 2 ALI 0 0 0
4
3
0 0 4L ; 15 EI ;6L + 10 EI 2L + 30 EI = 0
2 2 L 4 1 L 3 2 1 L (13.55)
5
0 0 ;6L + 101 L
EI4
3
24 ; 125 L
EI
2
0
6 0 0 1 L
2L2 + 30 EI 0 8L2 ; 15 EI
4 L4
1 4 5 3 6
1 ;
0 0 0
4 ; 2 ;
0 0 0
3 0 0 2 2 ; 15 ;6;L + 10 2 + 30 = 0
(13.56)
5 0 0 ;6L + 10 12 2 ; 5
;
0
6 0 0
2 + 30 0 4 2 ; 15
Expanding the determinant, we obtain the cubic equation in
33 ; 2202 + 3; 840 ; 14; 400 = 0 (13.57)
and the lowest root of this equation is = 5:1772 .
We note that from Eq. 13.21, the exact solution for a column of length L was
2 (4:4934)2 EI = 5:0477 EI
Pcr = (4:4934)
l2 EI = (2L)2 L2 (13.58)
and thus, the numerical value is about 2.6 percent higher than the exact one. The mathematica
code for this operation is:
(* Define elastic stiffness matrices *)
ke[e_,a_,l_,i_]:={
{e a/l , 0 , 0 , -e a/l , 0 , 0 },
{0 , 12 e i/l^3 , 6 e i/l^2 , 0 , -12 e i/l^3 , 6 e i/l^2 },
{0 , 6 e i/l^2 , 4 e i/l , 0 , -6 e i/l^2 , 2 e i/l },
{-e a/l , 0 , 0 , e a/l , 0 , 0 },
{ 0 , -12 e i/l^3 , -6 e i/l^2 , 0 , 12 e i/l^3 , -6 e i/l^2 },
P P
15’
I=200
I=50 6’
10’
I=100
Solution:
The element stiness matrices are given by
2
u1 2 0 0 3
20 1; 208
k1e = 664 1;208
6
96; 667 777 (13.59-a)
5
u1 2
2
0 0 3
0:01 0:10 7
k1g = ;P 64 0:10 16 :00
6
6 77 (13.59-b)
5
2
0 2 0 3 3
7
k2e = 664 128 ; 890 64; 440 77
6
2
u1
3 0 0 3
0:01667 0:1
k3g = ;P 64 0 :1 9 :6 775
6 7
6
(13.59-e)
The global equilibrium relation can now be written as
(Ke ; P Kg ) = 0 (13.60)
u1 2 3
(66:75) ; P (0:026666) (1; 208:33) ; P (0:1) (1; 678:24) ; P (0:1)
(1; 208:33) ; P (0:1) (225; 556:) ; P (16:) (64; 444:4) ; P (0) = 0 (13.61)
(1; 678:24) ; P (0:1) (64; 444:) ; P (0) (209; 444:) ; P (9:6)
The smallest buckling load amplication factor is thus equal to 2; 017 kips.
(* Initialize constants *)
a1=0
a2=0
a3=0
i1=100
i2=200
i3=50
l1=10 12
l2=15 12
l3=6 12
e1=29000
e2=e1
e3=e1
(* Define elastic stiffness matrices *)
ke[e_,a_,l_,i_]:={
{e a/l , 0 , 0 , -e a/l , 0 , 0 },
{0 , 12 e i/l^3 , 6 e i/l^2 , 0 , -12 e i/l^3 , 6 e i/l^2 },
{0 , 6 e i/l^2 , 4 e i/l , 0 , -6 e i/l^2 , 2 e i/l },
{-e a/l , 0 , 0 , e a/l , 0 , 0 },
{ 0 , -12 e i/l^3 , -6 e i/l^2 , 0 , 12 e i/l^3 , -6 e i/l^2 },
{ 0 , 6e i/l^2 , 2 e i/l , 0 , -6 e i/l^2 , 4 e i/l }
}
ke1=ke[e1,a1,l1,i1]
ke2=ke[e2,a2,l2,i2]
ke3=ke[e3,a3,l3,i3]
(* Define geometric stiffness matrices *)
kg[l_,p_]:=p/l{
{0 , 0 , 0 , 0 , 0 , 0 },
{0 , 6/5 , l/10 , 0 , - 6/5 , l/10 },
{0 , l/10 , 2 l^2/15 , 0 , - l/10 , - l^2/30 },
{0 , 0 , 0 , 0 , 0 , 0 },
{0 , -6/5 , - l/10 , 0 , 6/5 , - l/10 },
{0 , l/10 , - l^2/30 , 0 , - l/10 , 2 l^2/15 }
50
80,000
6m 6m
Solution:
Using two elements for the beam column, the only degrees of freedom are the de
ection and
rotation at midspan (we neglect the axial deformation).
Victor Saouma Matrix Structural Analysis
Draft
13.4 Geometric Non-Linearity 13{19
The element stiness and geometric matrices are given by
0 0 0 0 v1 2
2 3
0 0 0 0 0 0
60
6 222; 222: 666 ; 666: 0 ;222; 222: 666; 666: 77
[K1e ] = 666 00 666;0666: 2; 6660; 666
6
0: ;666; 666: 1; 333; 333 777 (13.62)
6
0 0 0 7
7
4 0 ;222; 222: ;666; 666: 0 222; 222: ;666; 666: 5
0 666; 666: 1; 333; 333 0: ;666; 666: 2; 666; 666
0 v1 2 0 0 0
2 3
0 0 0 0 0 0
60
6 222; 222: 666; 666: 0 ;222; 222: 666; 666: 77
[K2e ] = 666 00 666;0666: 2; 6660; 666 00: ;6660; 666: 1; 3330; 333 777
6 7
(13.63)
6 7
4 0 ;222; 222: ;666; 666: 0 222; 222: ;666; 666: 5
0 666; 666: 1; 333; 333 0: ;666; 666: 2; 666; 666
0 0 0 0 v1 2
2
0 0 0 0 0 0 3
6 0 ;16; 000 ;8; 000 0 16; 000 ;8; 000 777
6
[K1g ] = 666 00 ;8;0000 ;640; 000 00 8; 000 16; 000 77
6
0 0 77 (13.64)
6
40 16; 000 8; 000 0 ;16; 000 8; 000 5
0 ;8; 000 16; 000 0: 8; 000 ;64; 000
0 v1 2 0 0 0
2
0 0 0 0 0 0 3
6 0 ;16; 000 ;8; 000 0 16; 000 ;8; 000 777
6
[K2g ] = 666 00 ;8;0000 ;640; 000 00 8; 000 16; 000 77
6
0 0 77 (13.65)
6
40 16; 000 8; 000 0 ;16; 000 8; 000 5
0 ;8; 000 16; 000 0: 8; 000 ;64; 000
Assembling the stiness and geometric matrices we get
v1
"
2 #
[K] = 4120; 444
:
: 0:
5; 205; 330 (13.66)
6
0 0 0 0 0 7>
7>>
0 >
>
>
40 ;206; 222 ;658; 667: 0 206; 222: ;658; 667: > 5 >
>
> ; 0 : 00012123 >
>
>
>
0 658; 667: 1; 349; 330 0 ;658; 667: 2; 602; 670 :
0 ;
8 9
>
> 0 > >
>
>
>
>
>
25: > >
>
>
>
=
<
79:8491 = (13.68-a)
>
>
>
0: > >
>
>
>
>
>
:
;25: > >
>
>
79:8491 ;
Note that had we not accounted for the axial forces, then
( ) ( )
v1 = ;0:0001125 (13.69-a)
2 0
8 9 8 9
>
> Plft >
> >
> 0 >
>
>
>
>
>
>
Vlft >
>
>
>
>
>
>
>
>
>
25: >
>
>
>
>
<
Mlft =
=
<
75: =
(13.69-b)
>
>
>
Prgt >
>
>
>
>
>
0: >
>
>
>
>
>
>
:
Vrgt >
>
>
>
;
>
>
>
>
:
;25: >
>
>
>
;
Mrgt 75:
Alternatively, if instead of having a compressive force, we had a tensile force, then
( ) ( )
v1 = ;0:000104944 (13.70-a)
2 0
8 9 8 9
>
> Plft >
> >
> 0 >
>
>
>
>
>
>
Vlft >
>
>
>
>
>
>
>
>
>
25: >
>
>
>
>
<
Mlft =
=
<
70:8022 =
(13.70-b)
>
>
>
Prgt >
>
>
>
>
>
0: >
>
>
>
>
>
>
:
Vrgt >
>
>
>
;
>
>
>
>
:
;25: >
>
>
>
;
Mrgt 70:8022
Victor Saouma Matrix Structural Analysis
Draft
13.4 Geometric Non-Linearity 13{21
We observe that the compressive force increased the displacements and the end moments,
whereas a tensile one stiens the structure by reducing them.
The Mathematica to solve this problem follows
(* Initialize constants *)
OpenWrite["mat.out"]
a1=0
a2=0
e=2 10^9
i=2 10^(-3)
i1=i
i2=i1
l=6
l1=l
l2=6
p=-80000 (* negative compression *)
load={-50,0}
(* Define elastic stiffness matrices *)
ke[e_,a_,l_,i_]:={
{e a/l , 0 , 0 , -e a/l , 0 , 0 },
{0 , 12 e i/l^3 , 6 e i/l^2 , 0 , -12 e i/l^3 , 6 e i/l^2 },
{0 , 6 e i/l^2 , 4 e i/l , 0 , -6 e i/l^2 , 2 e i/l },
{-e a/l , 0 , 0 , e a/l , 0 , 0 },
{ 0 , -12 e i/l^3 , -6 e i/l^2 , 0 , 12 e i/l^3 , -6 e i/l^2 },
{ 0 , 6e i/l^2 , 2 e i/l , 0 , -6 e i/l^2 , 4 e i/l }
}
ke1=N[ke[e,a1,l1,i1]]
ke2=N[ke[e,a2,l2,i2]]
(* Assemble structure elastic stiffness matrices *)
ke=N[{
{ ke1[[5,5]]+ke2[[2,2]], ke1[[5,6]]+ke2[[2,3]]},
{ ke1[[6,5]]+ke2[[3,2]], ke1[[6,6]]+ke2[[3,3]]}
}]
WriteString["mat.out",MatrixForm[ke1]]
WriteString["mat.out",MatrixForm[ke2]]
WriteString["mat.out",MatrixForm[ke]]
(* Define geometric stiffness matrices *)
kg[p_,l_]:=p/l {
{0 , 0 , 0 , 0 , 0 , 0 },
{0 , 6/5 , l/10 , 0 , - 6/5 , l/10 },
{0 , l/10 , 2 l^2/15 , 0 , - l/10 , - l^2/30 },
{0 , 0 , 0 , 0 , 0 , 0 },
{0 , -6/5 , - l/10 , 0 , 6/5 , - l/10 },
{0 , l/10 , - l^2/30 , 0 , - l/10 , 2 l^2/15 }
}
kg1=N[kg[p,l1]]
kg2=N[kg[p,l2]]
(* Assemble structure geometric stiffness matrices *)
kg=N[{
{ kg1[[5,5]]+kg2[[2,2]], kg1[[5,6]]+kg2[[2,3]]},
{ kg1[[6,5]]+kg2[[3,2]], kg1[[6,6]]+kg2[[3,3]]}
}]
1,000
12 12
θ1 θ2
13.5 Summary
? ? ?
2 2
dx ; y
"x = du
2nd Order D.E. 4th Order D.E. dv dv
+ 21 dx
dx2
2 B.C. 4 B.C.
? ?
d2 y
; EIP = 0 U= 1R 2
dx2 2
E" d
v = ;A sin kx ; B cos kx
? ?
EI ; P = w
d4 v
dx4
d2 v
dx2 K= Ke Kg
v = C1 sin kx + C2 cos kx + C3 x + C4
?
P = (Ke + Kg )v
?
jKe + Kg j = 0
Appendix A
REFERENCES
Basic Structural Analysis :
1. Arbabi, F., Structural Analysis and Behavior, McGraw-Hill, Inc., 1991
2. Beaufait, F.W., Basic Concepts of Structural Analysis, Prentice-Hall Inc., Englewood
Clis, N.J., 1977
3. Chajes, A., Structural Analysis, Prentice-Hall, Inc., Englewood Clis, N.J., 1983
4. Gerstle, K.H., Basic Structural Analysis, (Local Reprint 1984.
5. Ghali, A., and Neville, A.M., Structural Analysis, Chapan and Hall, London, 1978
6. Gutowski, R.M., Structures: Fundamental Theory and Behavior, Van Nostrand Rein-
hold Co., N.Y., 1984
7. Hsieh, Y.Y., Elementary Theory of Structures, 2nd Ed., Prentice Hall, Inc., Englewood
Clis, N.J., 1982
8. Laursen, H.I., Structural Analysis, 2nd Ed., McGraw-Hill, N.Y., 1978
9. Morris, J.C., Wilbur, S., and Utku, S., Elementary Structural Analysis, McGraw-Hill,
N.Y., 1976
10. Wang, C.K., Intermediate Structural Analysis, McGraw-Hill, N.Y., 1983
Matrix Analysis :
1. Argyris, J.H., Recent Advances in Matrix Methods of Structural Analysis, Pergamon
Press, Oxford, 1964
2. Beaufait, F.W., Rowan Jr., W.H., Hoadley, P.G., and Hackett, R.M., Computer Meth-
ods of Structural Analysis, 4th Edition, 1982
3. Bhatt, P., Programming the Matrix Analysis of Skeletal Structures, Halsted Press,
1986
4. Elias, Z.M., Theory and Methods of Structural Analysis, John Wiley & Sons, 1986
5. Holzer, S.M., Computer Analysis of Structures. Elsevier, 1985
6. Livesley, R., Matrix Methods of Structural Analysis, Pergamon Press, Oxford, 964
7. Martin, H.C., Introduction to Matrix Methods of Structural Analysis, McGraw-Hill,
N.Y., 1966
Draft
A{2 REFERENCES
8. McGuire, W., and Gallagher, R.H., Matrix Structural Analysis, John Wiley and Sons
Inc., N.Y., 1979
9. Meek, J.L., Matrix Structural Analysis, McGraw-Hill, N.Y., 1971
10. Meyers, V.J., Matrix Analysis of Structures, Harper and Row, Publ., N.Y., 1983
11. Przemieniecki, J.S., Theory of Matrix Structural Analaysis, McGraw-Hill, N.Y., 1968
12. Weaver Jr, W., and Gere, J.M., Matrix Analysis of Framed Structures, 2nd Ed., Van
Nostrand Co., N.Y., 1980
Introduction to Finite Element and Programming :
1. Bathe, K.J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Inc.,
Englewood Clis, N.J., 1982
2. Cook, Malkus, and Plesha, Concepts and Applications of Finite Element Analysis,
John Wiley & Sons, 1989 (Third Edition)
3. Gallagher, R.H., Finite Element Analysis Fundamentals, Prentice Hall, Inc., Engle-
wood Clis, N.J., 1979
4. Hinton and Owen, An Introduction to Finite Element Computation, Pineridge Press,
Swansea U.K., 1978
5. Hughes, T.R., The Finite Element Method, Linear Static and Dynamic Finite Element
Analysis, Prentice Hall, 1987
6. Zienkiewicz, O., and Taylor, R., The Finite Element Method, Vol. 1 Basic Formulation
and Linear Problems, 4th Ed., McGraw-Hill, 1989
Energy Methods :
1. Pilkey and Wunderlich, Mechanics of Structures, Variational and Computational Meth-
ods, CRC Press, 1994
2. Langhaar, H., Energy Methods in Applied Mechanics, John Wiley and Sons, N.Y.,
1962
3. Reddy, J.N., Energy and Variational Methods in Applied Mechanics, John Wiley and
Sons, 1984.
Numerical Techniques :
1. Jennings, A., Matrix Computations for Engineers and Scientists, John Wiley and
Sons, N.Y., 1977.
2. Hilderbrand, F.B., An Introduction to Numerical Analysis, McGraw-Hill, N.Y., 1974
3. Press, W., et. al., Numerical Recipes, The Art of Scientic Computing, Cambridge
University Press, 1987
Journals :
1. Journal of Structural Engineering, American Society of Civil Engineering
2. Computers and Structures
3. International Journal for Numerical Methods in Engineering
Appendix B
REVIEW of MATRIX ALGEBRA
Because of the discretization of the structure into a nite number of nodes, its solution will
always lead to a matrix formulation. This matrix representation will be exploited by the
computer ability to operate on vectors and matrices. Hence, it is essential that we do get a
thorough understanding of basic concepts of matrix algebra.
B.1 Denitions
Matrix: 2 3
A11 A12 : : : A1j : : : A1n
6
6 A21 A22 : : : A2j : : : A2n 7
7
6 .. .. ... .. ... .. 7
[A] =
6
6 . . . . 7
7
(2.1)
6
6 Ai1 Ai2 : : : Aij : : : Ain 7
7
6
6 .. .. ... .. ... .. 7
7
4 . . . . 5
Am1 Am2 : : : Amj : : : Amn
We would indicate the size of the matrix as [A]mn , and refer to an individual term of
the matrix as Aij . Note that matrices, and vectors are usually boldfaced when typeset, or
with a tilde when handwritten A~.
Vectors: are one column matrices: 8 9
>
>
>
B1 >
>
>
>
>
>
>
B2 >
>
>
>
>
> .. >
>
fXg = > B. i
< =
>
(2.2)
> >
>
>
> .. >
>
>
>
>
>
> . >
>
>
>
:
Bm ;
Orthogonal matrices: [A]mn and [B]mn are said to be orthogonal if [A]T [B] = [B]T [A] =
[I]
A square matrix [C]mm is orthogonal if [C]T [C] = [C] [C]T = [I]
Trace of a matrix: tr(A) = Pni=1 Aii
Submatrices: are matrices within a matrix, for example
2 3
5 3 1 "
[ A ] [ A ]
#
[ A] = 4 4 6 2 5=
6 7 11
[A21 ] [A22 ]
12 (2.8)
10 3 4
2 3
1 5 "
[ B ]
#
[B] = 4 2 4 5=
6 7 1
[B2 ] (2.9)
3 2
" #
[A] [B] = [A11 ] [B1 ] + [A12 ] [B2 ] (2.10)
[A21 ] [B1 ] + [A22 ] [B2 ]
" #" # " #
[A11 ] [B1 ] = 5 3 1 5 = 11 37 (2.11)
4 6 2 4 16 44
Victor Saouma Matrix Structural Analysis
Draft
B.2 Elementary Matrix Operations B{3
[A22 ] [B2 ] = [4]
2
[3 2] = 3[12 8] (2.12)
14 34
[A] [B] = 4 22
6
48 75 (2.13)
28 70
Operate on submatrices just as if we were operating on individual matrix elements.
B.5 Inversion
The inverse of a square (nonsingular) matrix [A] is denoted by [A];1 and is such that
[A] [A];1 = [A];1 [A] = [I] (2.31)
Some observations
1. The inverse of the transpose of a matrix is equal to the transpose of the inverse
h i;1 h iT
AT = A;1 (2.32)
2. The inverse of a matrix product is the reverse product of the inverses
([A] [B]);1 = [B];1 [A];1 (2.33)
3. The inverse of a symmetric matrix is also symmetric
4. The inverse of a diagonal matrix is another diagonal one with entries equal to the inverse
of the entries of the original matrix.
5. The inverse of a triangular matrix is a triangular matrix.
6. It is computationally more ecient to decompose a matrix ([A] = [L] [D] [U]) using upper
and lower decomposition or Gauss elimination) than to invert a matrix.
When the determinant is expanded, we obtain an nth order polynomial in terms of which is
known as the characteristic equation of [A]. The n solutions (which can be real or complex)
are the eigenvalues of [A], and each one of them i satises
[A] fxi g = i fxi g (2.39)
where fxi g is a corresponding eigenvector.
It can be shown that:
1. The n eigenvalues of real symmetric matrices of rank n are all real.
2. The eigenvectors are orthogonal and form an orthogonal basis in En .
Eigenvalues and eigenvectors are used in stability (buckling) analysis, dynamic analysis, and to
assess the performance of nite element formulations.
Appendix C
SOLUTIONS OF LINEAR
EQUATIONS
Note this chapter is incomplete
C.1 Introduction
76 Given a system of linear equations [A]nn fxg = fbg (which may result from the direct
stiness method), we seek to solve for fxg. Symbolically this operation is represented by:
fxg = [A];1 fbg
77 There are two approaches for this operation:
Direct inversion using Cramer's rule where [A];1 = [adj[AA] ] . However, this approach is compu-
tationally very inecient for n 3 as it requires evaluation of n high order determinants.
Decomposition: where in the most general case we seek to decompose [A] into [A] = [L][D][U]
and where:
[L] lower triangle matrix
[D] diagonal matrix
[U] upper triangle matrix
There are two classes of solutions
Direct Method: characterized by known, nite number of operations required to achieve
the decomposition yielding exact results.
Indirect methods: or iterative decomposition technique, with no a-priori knowledge of
the number of operations required yielding an aapproximate solution with user dened
level of accuracy.
Draft
C{2 SOLUTIONS OF LINEAR EQUATIONS
C.2 Direct Methods
C.2.1 Gauss, and Gaus-Jordan Elimination
78 Given [A]fxg = fbg, we seek to transform this equation into
1. Gaus Elimination: [U]fxg = fyg where [U]is an upper triangle, and then backsubstitute
from the bottom up to solve for the unknowns. Note that in this case we operate on both
[A] & fbg, yielding fxg.
2. Gauss-Jordan Elimination: is similar to the Gaus Elimination, however tather than
transforming the [A] matrix into an upper diagonal one, we transform [AjI] into [IjA;1 ].
Thus no backsubstitution is needed and the matrix inverse can be explicitely obtained.
C.2.1.1 Algorithm
79 Based on the preceding numerical examples, we dene a two step algorithm for the Gaussian
ellimination.
80 Dening akij to be the coecient of the ith row & j th column at the k th reduction step with
i k & j k:
Reduction:
akik+1 = 0 k<in
ak ak
aij = akij ; akkk k < i n; k < j n
k +1 ik kj
(3.9)
k k
bkij+1 = bkij ; aikakkkbkj k < i n; 1 < j m
Backsubstitution: Xn
i
xij =
biij ; k=ii+1 aik xkj (3.10)
aii
Note that Gauss-Jordan produces both the solution of the equations as well as the inverse of
the original matrix. However, if the inverse is 3not desired it requires three times (N 3 ) more
operations than Gauss or LU decomposition ( N3 ).
C.2.2 LU Decomposition
81 In the previous decomposition method, the right hand side (fbg must have been known before
decomposition (unless we want to detemine the inverse of the matrix which is computationaly
more expensive).
82 In some applications it may be desirable to decompose the matrix without having the RHS
completed. For instance, in the direct stiness method we may have multiple load cases yet we
would like to invert only once the stiness matrix.
Victor Saouma Matrix Structural Analysis
Draft
C.2 Direct Methods C{5
This will be achieved through the following decomposition:
[A] = [L][U] (3.11)
It can be shown that:
1. Both decompositions are equivalent.
2. Count on number of operation show that the 2 methods yield the same number of opera-
tions. Number of operations in LU decomposition is equal to the one in Gauss elimination.
83The solution consists in:
Decomposition: of the matrix independently of the right hand side vector
[A] = [L] [U] (3.12)
[L] [|U]{zfxg} = fbg (3.13)
fy g
Backsubstitution: for each right hand side vector
1. Solve for fyg from [L]fyg = fbg starting from top
2. Solve for fxg from [U]fxg = fyg starting from bottom
84 The vector fyg is the same as the one to which fbg was reduced to in the Gauss Elimination.
C.2.2.1 Algorithm
1. Given:
2 3 2 32 3
a11 a12 a1n 1 u11 u12 u1n
6
6 a21 a22 a2n 7
7
6
6 l21 1 76
76 u22 u2n 7
7
6
6 .. .. .. .. 7
7 = ..
6
6 .. ... 76
76 ... 7
7 (3.14)
4 . . . . 5 4 . . 54 5
an1 an2 ann ln1 ln2 1 unn
2. solve:
a11 = u11 a12 = u12 a1n = u1n
a21 = l21 u11 a22 = l21 u12 + u22 a2n = l21 u1n + u2n
.. (3.15)
. Xn;1
an1 = ln1 u11 an2 = ln1 u12 + ln2 u22 ann = k=1 lnk ukn + unn
3. let:
2 3
u11 u12 u1n
6
6 l21 u22 u2n 7
7
[A]F = ..
6
6 .. 7
7 (3.16)
4. . 5
ln1 ln2 unn
Victor Saouma Matrix Structural Analysis
Draft
C{6 SOLUTIONS OF LINEAR EQUATIONS
4. Take row by row or column by column
Xj ;1
lij =
aij ; k=1 lik ukj i>j
u jj
X i;1 (3.17)
uij = aij ; k=1 lik ukj i j
lii = 1
Note:
1. Computed elements lij or uij may always overwrite corresponding element aij
2. If [A] is symmetric [L]T 6= [U], symmetry is destroyed in [A]F
For symmetric matrices, LU decomposition reduces to:
Xi;1
uij = aij ; k=1 lik ukj ij
lii = 1u (3.18)
lij = ujjji
|
2 0
{z
;8 2
}|
0 0 0
{z
2
}
[L] [U]
| {z }
[A]
j x ;xlkx j "
k k;1
(3.31)
2. The convergence can be accelerated by relaxation
xki = xki + (1 ; )xki ;1 (3.32)
where is a weight factor between 0. and 2. For values below 1 we have underrelaxation,
and for values greater than 1 we have overrelaxation. The former is used for nonconvergent
systems, whereas the later is used to accelerate convergence of converging ones. optimum
for frame analysis is around 1.8.
Appendix D
TENSOR NOTATION
NEEDS SOME EDITING
Ax Ay Az
(4.1-e)
Bx By Bz
grad A = rA = i @A +
@x @y j @A + k @A
@z (4.1-f)
div A = rA = @ + j @ + k @ (iA + jA + kA )
i @x @y @z x y z
1. If there is one letter index, that index goes from i to n. For instance:
8 9
>
< a1 >
=
ai = ai = b a1 a2 a3 c = > a2 >
i = 1; 3 (4.2)
: a3 ;
assuming that n = 3.
A11 = B11 C11 D11 + B11 C12 D12 + B12 C11 D21 + B12 C12 D22
A12 = B11 C11 D11 + B11 C12 D12 + B12 C11 D21 + B12 C12 D22
A21 = B21 C11 D11 + B21 C12 D12 + B22 C11 D21 + B22 C12 D22
A22 = B21 C21 D11 + B21 C22 D12 + B22 C21 D21 + B22 C22 D22 (4.10-a)
Appendix E
INTEGRAL THEOREMS
76Some useful integral theorems are presented here without proofs. Schey's textbook div grad
curl and all that provides an excellent informal presentation of related material.
Draft
E{2 INTEGRAL THEOREMS
or
Z Z
vini d; = vi;id
(5.5)
;
79 Alternatively
Z Z Z Z Z Z Z Z
div qdV = qT :ndS ; (r)T qdV (5.8)
V S V