LEHIGH U IVERSITY

Residual Stresses in Thick Welded Plates

APPLICATIONS OF THE FINITE
ELEMENT METHOD TO
BEAMCOLUMN PROBLEMS
FRITZ ENC1NEERING
i:Ar:JOR/\TOr,'( UBRAR'f
by
Negussie Tebedge
September, 1972
Fritz Laboratory Report No. 337.32
APPLICATIONS OF THE FINITE ELEMENT METHOD
TO BEAM-COLUMN PROBLEMS
by
Negussie Tebedge
A Dissertation
Presented to the Graduate Committee
of Lehigh University
in Candidacy of the Degree of
Doctor of Philosophy
in
Civil Engineering
Lehigh University
1972
CERTIFICATE OF APPROVAL
Approved and recommended for acceptance as a dissertation
in partial fulfillment of the requirements for the degree of Doctor
of Philosophy.
1(72
(Date)
Accepted
Lambert Tall
Professor in Charge
Special committee directing the
doctoral work of Mr. Negussie
Tebedge
e (h/.
George C. Drisco 11, Jt .

Professor Lambert Tall
tJ(7
Professor Wai-Fah Chen
Professor Robert A. Lucas

Professor David A. VanHorn
Ex Officio
ACKNOWLEDGMENTS
The author is deeply indepted to Professor Lambert Tall,
Professor in Charge of this dissertation, for his guidance and con-
tinuous encouragement in the preparation of this thesis and through-
out the author's graduate studies. The guidance of Professors George
C. Driscoll (Chairman), Wai-Fah Chen, Robert A. Lucas, and David A.
VanHorn, members of the special committee directing the author's
doctoral work, is gratefully acknowledged.
Sincere appreciation is expressed to all the author's
associates in Fritz Laboratory, in particular Suresh Desai, for
fruitful discussions throughout this investigation.
The investigation leading to this dissertation is part of a
research project "Residual Stresses in Thick We lded Plates", which is
jointly sponsored by the National Science Foundation and the Column
Research Council. Dr. L. Tall and Dr. W. F. Chen are the directors
of the project.
The author owes a special debt of gratitude to the African-
American Institute for their support during the first year of his
graduate study.
iii
Acknowledgments are also due to the Computer Center at Lehigh
University which provided facilities for the computer work. The
assistance of Shirley Matlock in typing the manuscript and of John
Gera and his staff for preparing the drawings, is gratefully acknowl-
edged.
iv
LIST OF TABLES
Page
TABLE 1 CRITICAL LOADS OF COLUMNS WITH DISTRIBUTED AXIAL LOADS
114
TABLE 2 CRITICAL LOADS OF UNIFORMLY TAPERED COLUMNS
115
TABLE 3 CR ITICAL LOADS OF TWO-COMPONENT COLUMNS
116
TABLE 4 CRITICAL LOADS OF PIECEWISE PRISMATIC COLUMNS
117
TABLE 5 EXPERIMENT ON PRE TWIS TED COLUMNS
118
TABLE 6 CRITICAL LOADS OF PRETWISTED COLUMNS
119
TABLE 7 LATERAL TORSIONAL BUCKLING LOADS
120
TABLE 8 SPACE FRAME BUCKLING
121
TABLE 9 BUCKLING OF SHAFTS UNDER PlffiE TORSION
122
Fig. 2.1
Fig. 2.2
Fig. 2.3
Fig. 2.4
Fig. 2.5
Fig. 2.6
Fig. 2.7
Fig. 2.8
Fig. 2.9
Fig. 3.1
Fig. 3.2
Fig. 3.3
Fig. 3.4
Fig. 3.5
Fig. 3.6
Fig. 3.7
Fig. 3.8
v
LIS T OF FIGURES
Page
Reference Axes System for the Beam Element 123
Notations for Generalized Stresses and Displacements of the 123
Beam Element in Flexure and Extension
The Centroida1-Principa1 Axes and the Generalized Coordinate 124
System
Notations for a Beam Component Subjected to Torsion and 124
Warping
A Cantilever Beam Loaded by Torque
125
A Cantilever Beam Loaded by a Concentrated Bimoment
126
Comparison of Beam Element Stiffness Matrices using a Canti- 127
lever Beam under Torsion and Bimoment
Comparison of Beam Stiffness Matrices 128
Continuous Beam l,oaded by Concentrated Bimoment 129
Beam Element Under Axial and Torsional Loading 130
Kinematics of Wide Flange Shape Under Torsion 130
Schematic Description of the Kinematics of a Shaft Under 131
Large Torsional Displacement
Convergence of Finite Element Solution for Axially Loaded 132
Columns
Comparison of Convergence Between Finite Element and Finite 133
Differences Solutions for the Euler Column
A Column Under Distributed Axial Loads 134
Comparison of Finite Element and Analytical Solutions of 135
Columns with Distributed Axial Loads
The First Buckling Modes of Columns with Distributed Axial 136
Loads
Fig. 3.9 Comparison of Finite Element and Analytical Solutions of ' 137
Tapered Columns
Fig. 3.10 Finite Element Solution of Critical Loads of Tapered Columns 138
Fig. 3.11 Convergence of Finite Element Solution of Tapered Columns 139
Fig. 3.12 Comparison of Finite Element and Analytical Solutions of 140
Columns on Elastic Foundations
vi
Page
Fig. 3.13 The Generalized Forces and Displacements of a Two-Component 141
Column
Fig. 3.14 Buckling Loads of Two-Component Columns with Pin Ends 142
Fig. 3.15 Buckling Loads of Two-Component Columns with Fixed Ends 143
Fig. 3.16 Buckling Loads of Piecewise Prismatic Columns with Pin Ends 144
Fig. 3.17 Buckling Loads of Piecewise Prismatic Columns with Fixed 145
Ends
Fig. 3.18 Preparation of a pretwisted Column
146
Fig. 3.19 Buckling Test of a Pretwisted Column with Knife Edge Condi-
147
tions
Fig. 3.20 Comparison of Finite Element Solution and Experimental 148
Results of pretwisted Columns
Fig. 3.21 Comparison of Theoretical and Experimental Buckling Modes 149
of pretwisted Columns
Fig. 3.22 Finite Element Solution of pretwisted Columns with Knife 150
Edge Conditions
Fig. 3.23 Buckling Loads of Prismatic Columns with Crossed Pins 151
(Pin Ends)
Fig. 3.24 Buckling Loads of Prismatic Columns with Crossed Pins 152
(Fixed-Pinned Ends)
Fig. 3.25 Convergence of Finite Element Solution of Torsional and 153
Lateral-Torsional Buckling
Fig. 3.26 Comparison of Finite Element and Analytical Solutions of 154
Lateral Torsional Buckling
Fig. 3.27 Buckling Loads for Two-Span Continuous Beams 155
Fig. 3.28 Simply Supported Tapered I-Beams 156
Fig. 3.29 Comparison of Finite Element Solution and Experimental 157
Results of Lateral Torsional Buckling of Tapered Beams
Fig. 3.30 Comparison of Finite Element Solution and Finite Integral 158
Solution of Tapered Beams
Fig. 3.31 Finite Element Solutions of Beams of Combined Tapers 159
Fig. 3.32 Convergence Characteristics of Knee-Bent Frames 160
vii
Fig. 3.33 Buckling of Portal Frames 161
Fig. ;1,34 Space Frame BuckUng by Finite Element Method 162
Fig. 3,35 Buckling Modes of Shafts under Pure Torsion 163
Fig. 4.1 Basic Incremental Procedure 164
165
Fig. Effect of Number of Load Increments when Us ing the Incrementa1166
Method
Fig. 4,4 Use of the Incrementa,l Method to Determine 167
Curves of Columns
Flg. l},5 Use of the Incremental Nethod for the Analysis of Large 168
Deflections
TABLE OF CONTENTS
ABSTRACT
1. INTRODUCTION
2. FORMULATION FOR THE STIFFNESS MATRIX OF THE BEAM ElEMENT
2.1 INTRODUCTION
2.2 A VARIATIONAL FORMULATION FOR THE BEAM ELEMENT
A Generalized Variational Principle
A Variational Principle for a Beam MOdel
A Finite Element Model
2.3 EVALUATION OF THE BEAM STIFFNESS MATRIX
Stiffness Matrix in Extension and Flexure
Stiffness Matrix in Torsion
Hyperbolic Functions
Polynomials
2.4 APPLICATIONS AND SAMPLE SOLUTIONS
Numerical Examples
3. LINEAR STABILITY ANALYSIS OF BEAM-COLUMNS
viii
Page
1
3
9
9
10
10
11
15
17
19
31
35
36
38
40
43
3.1 INTRODUCTION 43
3.2 A FINITE ELEMENT MODEL FOR DERIVING GEOMETRIC STIFFNESS 45
MATRICES
3.3 LARGE AXIAL DISPLACEMENTS 51
Flexural Geometric Stiffness Matrix
Torsional Geometric Stiffness Matrix
3.4 LARGE LATERAL DISPLACEMENTS
Flexural-Torsional Geometric Stiffness Matrix
3.5 LARGE TORSIONAL DISPLACEMENTS
Torsional-Lateral Geometric Stiffness Matrix
3.6 LINEAR EIGENVALUE PROBLEMS
Solution Technique
Numerical Examples and Results
Columns with Distributed Axial Loads
Tapered Columns
Columns on Elastic Foundations
Piecewise Prismatic Columns
Pretwisted Columns
Lateral-Torsional Buckling Problems
Space Frames
Buckling of Shaft under Torsion
51
53
55
55
59
59
62
62
64
65
67
69
71
74
78
81
83
4 NONLINEAR ANALYS IS OF BEAM-COLUMNS
4.1 INTRODUCTION
4.2 SOLUTION TECHNIQUES
Incremental Methods
Iterative Methods
4.3 GEOMETRIC NONLINEARITY
4.4 MATERIAL NONLINEARITY
4.5 SAMPLE PROBLEMS AND RESULTS
5. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE DEVELOPMENTS
6 APPENDICES
APPENDIX I
ELASTIC STIFFNESS MATRIX FOR THE BEAM ELEMENT
APPENDIX II
FLEXURAL STIFFNESS MATRIX OF TAPERED BEAMS
ix
85
85
87
88
90
91
93
95
96
100
100
102
APPENDIX III 105
CONSISTENT STIFFNESS MATRIX FOR BEAMS ON ELASTIC FOUNDATION
APPENDIX IV 107
FLEXURAL STIFFNESS MATRIX OF PRETWISTED BEAMS
7. NOTATIONS 109
8. TABLES AND. FIGURES 113
9. REFERENCES 169
10. VITA 176
-1
ABSTRACT
This dissertation is concerned with developing a finite element
formulation of the analysis of beam-columns, and with demonstrating the
applicability of the method to general beam-column problems as a practical
tool. The treatment includes linear static, linear stability and non-
linear analyses of beam-columns.
In the linear static analysis, a formulation is presented to
develop a one-dimensional finite element using variational principles.
The formulation is based on a functional consisting of two independent
fields: a polynomial approximation of a strain field in the domain, and
displacements at the boundary. The beam element has two nodes and
seven degrees of freedom at each node.
The linear stability analysis, which reduces to an eigenvalue
determination, utilizes a displacement formulation based on a finite
element idealization. A systematic procedure is developed to evaluate
the geometric stiffness matrices for beam-columns. The matrices derived
correspond to large displacements in axial and transverse directions and
also in twist. The f.inite element solutions are compared to analytical
solutions and the convergence characteristics are studied for a variety
of problems which include: columns with distributed axial loads, tapered
columns, columns on elastic foundations, pretwisted columns, space frames
and the lateral buckling of beams.
Finally, a finite element formulation of nonlinear analysis is
given to study general instability problems of beam-columns. A solution
procedure, using a direct incremental approach, is applied to numerical
examples to demonstrate the validity of the procedure.
The contributions achieved in this dissertation are:
-A one-dimensional finite element model is developed to
analyze general linear static beam-column problems.
-A systematic procedure is presented to evaluate geometric
stiffness matrices for beam-columns which are required to
perform a finite element analysis of stability problems.
-The geometric stiffness matrices are derived Which corres-
pond to large lateral and torsional displacements.
-2
-The advantages of the finite element method are demonstrated
in the solution of a few stability problems, such as the
buckling of pretwisted columns and the lateral buckling
of tapered beams. tIle analytical solutions of which are not
yet available.
-3
1. INTRODUCTION
A beam-column, also known as a rod in classical terms, is
defined in this dissertation as a three dimensional body having one
dimension significantly greater than the other two. Examples of bodies
which may be so regarded are numerous, such as members of framed struc-
tures, arches and curved beams. In the case of framed structures,
for example, the members are identified further depending on the
loading conditions, where columns represent the one limiting case where
the bending moments become zero, and beams the case in which the axial
force vanishes.
The purpose of a theory of beams is to provide appropriate
one-dimensional equations applicable to beam-type bodies. A one-
dimensional analysis, referred to as a beam theory, is necessarily
approximate and furnishes only partial or limited information. Indeed,
the desire for such limited information is the basic motivation for
the construction of a one-dimensional theory with the aim of providing
a simpler theory for the limited information sought. While the three-
dimensional viewpoint is certainly the most fundamental, the possibility
of employing a one-dimensional model for beam-type bodies presents it-
self in a natural way because of the considerable difficulties associated
with the derivation of the beam theory from the three-dimensional equa-
tions. The model, however, must be capable of supplying a substantial
portion of information the three-dimensional theory would furnish. The
notion of employing a model for an idealized body is frequently used
in classical continuum mechanics, in fact, the continuum itself is a
model representing an idealized body in some sense.
-4
An approximate system of equations for beam-type bodies may
be developed by converting formally the three-dimensional field rela-
tionships to their one-dimensional analogue. Historically, interest
in the construction of more elaborate theories of beams arose from the
desire to treat wave propagation and vibrations of elastic rods. After
the three-dimensional theories were accepted in certain domains of
mechanics, Cauchy and Poisson sought to obtain theories by averaging
over a cross-section the results from a three-dimensional theory and
then letting the cross-sectional area approach zero(1,2,3). Recently,
the use of polynomial approximations has been adopted extensively
to develop analytical beam theories. For example, the three-dimensional
field relations may be converted to their one dimensional analogues by
replacing the field variables by series expansions in products of
Legendre polynomials (4,5,6) . In these efforts, the Legendre expansions
lead to reducing the governing partial differential equations to either
ordinary differential equations, or to more tractable partial differen-
tial equations. The exact analysis of beam behavior, when treated
in this fashion, will be intrinsically more complex, necessitating the
satisfaction of the boundary conditions on numerous planes as compared
to one pair of surfaces for plates or shells.
Inasmuch as considerable difficulties remain in the deriva-
tion of a system of equations from the three-dimensional theory, the
alternative development is to utilize a direct approach if a simplified
formulation is sought. The tIc lass ica1 beam theory", for example, is
based on a direct approach, In the development of the classical theory,
Bernoulli (1705) was the first to make kinematical assumptions to solve
, ...
,
-5
flexural problems, and hypotheses regarding the constitution of the
material were first given by Euler (1771). Saint-Venant (1855) was the
first to remark that six equations are needed to express the equilibruim
of rods which are twisted as well as bent, based on special simplifying
hypotheses. The general equations were given in principle, but obscurely,
by Kirchhoff (1859). The process by which Kirchhoff developed his
theory was, to a great extent, kinematical. I ~ b s c h (1862) modified
the theory and gave explicit general equations which were confirmed
by later writers(1,7),
Recently, it has been established that the Euler-Bernoulli
theory of beams was not applicable to thin-walled beams because of the
inherent distortion of the cross section that occurs during bending.
Wagner (1929) was the first to introduce the concept of "warping" in
the analysis of thin-wa11ed beams (8) , Comprehensive reviews on the
bending and torsion of open sections including the buckling character-
istics were made by Goodier(9) and Timoshenko(lO). A general treat-
ment of beams is fully described in Vlasov's treatise of thin-
walled beams (11) . A common feature in these investigations is that
each formulation results in the development of the governing differ-
ential equations from consideration of equilibrium conditions. Many
particular problems based on these general formulations have been
solved either exactly or approximately by seeking the analytical solu-
tions of the differential equations or by employing different numerical
techniques. A historical review on this subject, particularly on the
development and utilization of the various numerical approaches that
have been used to solve the governing equations is given in Ref. 12.

More recently, the calculation of complex structural problems
by means of the concept of piecewise approximations has received a
great impetus. A significant intermediary step in the evolution of
modern structural mechanics is the discrete element method. Here,
in the context of beam-type bodies, the structural beam is physically
replaced as a combination of elastic blocks, rigid bars, torsional springs
and flexural springs. This is equivalent to the early works of
Hrenikoff(13), representing a plane solid as an assembly of discrete
systems, which is regarded as a forerunner to the development of general
discrete methods of structural mechanics. It has been shown that
the discrete element approach is mathematically equivalent to the finite
difference method and thus it may be considered as a physical
pretation of the finite difference method (14) .
During the past decade, great strides have been made on the
development and utilization of the finite element method. This method
can be considered as the most powerful and versatile technique presently
available for the numerical solution of complex structural problems.
Moreover, it can be formulated in terms of simple physical concepts
without recourse to complex differential equations. The method was
developed originally as an application of the standard structural
analysis procedure to a physically discretized approximation of the
actual system. The concept has been extensively described in Refs.
lS and 16. Study of the mathematical foundations of the method(16,l7)
as well as its application to a wider class of field prob1ems(lS,l6)
has clarified the basic requirements for its effective formulation.
-7
Two parallel developments were responsible for the widespread
acceptance of the finite element method; the formulation of the matrix
transformation theory of structures and the introduction of high-speed
digital computers. The use of matrices allows a very efficient, sys-
tematic and simplified calculation superior to any other currently
available scheme. Once the initial matrices are assembled, the sub-
sequent operations involve merely elementary matrix algebra, which are
ideally suitable for automatic computations using the computer.
While the advantages of the finite element method have been
widely recognized and its applications extensively demonstrated
particularly to a variety of problems in solid mechanics and in struc-
tural mechanics, more specifically to plate and shell structures, the
application of the method to beam-column analyses has not been explored
to an equivalent degree. Analysis of beam-columns, and in particular,
beam-columns under general loading and support conditions, is a subject
of wide interest in current research. Most of the previous develop
ments in beam-column analysis by finite elements are found in Refs.
18 to 22.
In this study, a direct approach is employed to analyse the
beam-column problem where a numerical solution is sought by utilizing
the finite element concept and its applications. The beam is ficti-
tiously subdivided by imaginary planes into an assembly of elements
and is regarded as a one-dimensional problem. This notion of sub-
division, which is mathematical and not physical does not consider the
beam to be divided into separate physical elements that are assembled
in the analysis procedure. Using this concept, a variational prin-
ciple is employed in constructing a finite element formulation to
-8
evaluate the properties of the elements and finally to solve the com-
plete system. In order to demonstrate the best balance in practical
usage, the study takes into consideration factors such as Slllplicity
of formulation, versatility of application, reliability, computational
efforts and accuracy of results.
2.1 INTRODUCTION
2. FORMULATION FOR THE STIFFNESS
MATRIX OF THE BEAM ELEMENT
-9
The finite element technique may be applied to analyse the
beam problem in which the beam can be treated as a three-dimensional
body with the use of three-dimensional elements, or as a general plane
stress problem when two-dimensional elements are used. Here the sole
interest is to construct a one-dimensional model capable of furnishing
a substantial portion of the information a three-dimensional theory
would furnish. Indeed, the evaluation of the element is one of the
most important aspects of the finite element analysis. A fundamental
property of finite element models is that typical elements can be
isolated from the idealized system and their behavior can be studied
independently. The process of connecting the elements to form the final
system is mainly a topological one and is independent of the physical
nature of the problem.
In evaluating the element properties, either a direct method
or a variational method may be used. The direct approach, in which
direct consideration is given to the conditions of equilibrium and
compatibility, is not used in this study. A thorough treatment of
the direct approach is given in Ref. 23.
The formulation presented herein is based on an appropriately
constructed functional where variational principles are applied to
develop the finite element model. The functional is established based
on two independent fields: a polynomial approximation of the strain
i '
-}O
field in the domain, and displacements at the boundaries. The use of
polynomials is advantageous since it permits differentiation and in-
tegration, with relative ease. The main purpose of the choice of
such a formulation is its ability to incorporate the underlying hypo-
theses given by the beam theories in a rather simple manner. Further
discussions and justifications for the choice of the formulation are
given in a later section. Of course, in most structural problems,
assumption of the displacement field alone usually will furnish good
results. This is because the strain field can be derived in a straight-
forward manner, by using the relationships given by the deformation
theory, as the derivatives of the displacement. However, the same
logic does not hold true for the case of a one-dimension.al analysis
of the beam problem due to the complex: nature of the problem"
2.2 A VARIATIONAL FORMULATION FOR THE BEAM ELEMENT
A Generalized Variational Principle
The variational principle may be regarded as one of the most
important bases for the finite element method(24,25). It has con-
tributed to the development of structural analysis in leading to finite
.. 11
element formulations. Numerous finite element models may be derived,
based on variational principles, by introducing different constraining
conditions within the element or at the interelement boundaries. Since
the problem usually cannot be solved exactly, the variational method
provides an approximate formulation of the problem which yields a
solution compatible with the assumed degree of approximation(26).
The variational formulation has several other advantages
once the existence of a functional is assured .. The functional which
is subject to variation may usually be given a physical interpreta-
tion (such as strain energy or complementary energy) and is invariant
under coordinate transformation. The original problem may be trans-
formed into an equivalent one that can be solved more easily, for
example, by applying the method of the Lagrangian multiplier for problems
having subsidiary conditions. Another advantageous aspect of the
variational principles is that they may lead to es'tablishing upper and
lower bounds of the exact solution; also, they may provide convergence
(24 26)
proofs ' .
The generalized variational principle, often referred to as
the Washizu-Hu principle, involves several free and independent fields
(26,27)
The general principle is based on three independent fields
in the domain, namely, the displacement field u
i
' strain field eij
and stress field crij; and two fields on the boundary, the displacement
field U. on S and boundary traction Pi on S. The generalized func-
u
tional may be expressed by
U()
(*)Standard tensor notation and the summation convention is used. A
comma denotes partial derivation with respect to the variable that
follows.
-1Z
0'1'J' - -Zl(U .. + u .. )] - Fi u. dV
1,J J,1 . 1
- J
(Z .1)
In
8
0'
this expression
D
ijkt
ij
O'ij
u
i
F.
1
T.
1
Pi
u.
1
8 ,8
0' u
=
genera lized Hookean constant
=
strain tensor
=
stress tensor
=
displacement
=
prescribed body force
= prescribed traction on boundary S
0'
=
prescribed traction on boundary 8
I:: prescribed displacement on S
u
portions of boundaries over which T. and u. are prescribed,
1. 1
respectively
s = s + S = whole surface
(J u
There are eighteen independent variables subject to variation
in the functional TI , with no constra.ining or subsidiary conditions,
w
these are: three displacements u., six strains e'j' six stresses 0'
1 1 ij
and three boundary tractions P . Taking variations with respect to
1
these quantities leads to the following(Z8)
1
- J [ei' - -2(U' . + u. ,)] 8 O"j dv
V J J.J. 1,J 1
(constitutive equations) ,
(strain-displacement
relations)
- J
V
- J
S
0'
Is
u
- J
S
u
(O'ij . + F.)
& u
j
dV (equations of
,l. J
equilibrium)
(T
j - n
i
O"ij) 0
u, ciS (static boundary
J
conditions on S )
a
(u. - ~ .
6
P. dS (kinematic boundary
l. l. I.
conditions on S )
u
(P. - n.
crij)
I) u. dS (continuity of
J
1.
J
tractions on S )
u
The vanishing of 811 will establish the relations between
w
the fields and impose on them the appropriate field equations, and
(2.2)
boundary and continuity conditions as expressed by the following Euler
equations
(i)
aij
0= Dijk{, E:k{,
in V
(ii)
1
in V
e:ij
= I(u j i + u
i
.)
, ,J
(iii) 0' + F. 0 in V
1.J ,J 1-
(2.3)
(iv) T. = n.
O"ij
= T
j
on S
J
l.
(v) u, = u. on S
I. I. t!
(vi) P
j
= n
i o'ij
=
T
j
on S
u
Based qn the generalized variational principle given by Eq,
2.1, different forms of variational principles may be derived by making
a priori assumptions on one or more subsidiary conditions. For example,
by stipulating that the stress and the strain fields are related by
the constitutive equations, the variat ional principle will be involved
with one less independent field, The functional TT thus will be
w
-14
reduced to another functional TTR' equivalent to the Hellinger-Reissner
principle(29). Stipulating further that the strain fields are remted
to the displacement field, and by satisfying the kinematic boundary
conditions, the functional TTR reduces to another functional T T p ~ which
is equivalent to the principle of minimum potential energy. Similarly)
the principle of minimum complementary energy also can be derived by
introducing the appropria.te subsidiary conditions (26).
A Variational Principle for a Beam Model
At this stage, a variational principle can be established.
from the generalized principle to evalua.te the properties of the finite
element modeL The functions that will be assumed in this variational
principle are
a)
b)
strain fields e .. in the domain V
~
boundary displacement fields u. in S .
~ u
If, in addition, it is stipulated that the stress and the strain fields
are related by the constitutive equations and the static boundary c o n ~
ditions are satisfied, then incorporating these constraint conditions
in the functional given by Eq. 2.1, and introducing the Lagrangia.n
multiplier technique, the functional will be reduced to:
TI
;:,
J
(1
Dijke, ij
eke, - F i u
i
)
dV
2
V
(2.4)
'"
- J
T. u. dS +
Is
T. u. ds
S
~ ~ ~ ~
u
where u. is the interelement boundary displacement and is the same for
~
two adjacent elements on their common boundary.
-15
Taking the variations with respect to the independent quanti-
ties results in the following,
+J
(D
ijkt
:U
n.
- T ) 8
u. dS
V
J
i
(2.5)
,....
+J
(u
i
- u ) B T. ds
i
S
,..,
+J
T.
5
u. dS
S
1. 1.
U
The vanishing of Bn for an arbitrary au. in V and an arbitrary

au. on the interelement boundaries S , will give the following Euler
u
equations
(Di.U :U) .
-+ F. = 0 in V
J d
1.
D
ijU i3kt nj
- Ti
0 on S (2.6)
.....
u. - u =
0 on S
u
A finite element model that satisfies the conditions inEq.
2.6 is developed in the following section. In the development of this
model, the functional given by Eq. 2.4 is applied directly.
A Finite Element Model
In this analysis, the body forces F i are ignored and the matrix
notation following pian
i
s(30) notations is employed, The functions fe}
and [u} are simply chosen as polynomials with unknown coefficients. The
strain field [s} is expressed in terms of polynomial functions of the
coordinates [PJ and undetermined strain coefficients fa}. The displacements
along the intere1ement boundaries [u} are represented by the interpo-
1ating functions [LJ and the generalized displacements fa} at the
nodes. In matrix form they may be written respectively as
te} = [PJ ~ }
(2.7)
[u} = [LJ [a)
(2.8)
(2.9)
The stress field may be derived from the strain field using
the constitutive equations, thus
fa) = [D] [P] f ~ }
(2.10)
The tractions at the boundaries tab} are expressed in terms
of the stress field fa} and the undetermined coefficients ~ } as
follows
(2.11)
where [R] contains the coordinates on the surface.
The functional given in Eq. 2.4 when written in matrix form
will reduce to
where
1 T T
TI = S 2 ~ } [PJ [D] [PJ ~ } dV
V
[H] = S [p]T [D] [P] dV
V
(2.12)
and
[TJ = f [RJT [LJ dS
S
-17
Minimizing the functional TI with respect to each fs}, that
is en/d ~ = 0, yields
[H] [S} - [T] to} = 0 (2.13)
from which the undetermined coefficients are solved as
-1
fa} = [H] [T] [5}
(2.14)
Since the functional can be expressed in terms of the element
stiffness matrix [k], the generalized force vector ff}, and the general-
ized displacements [5} as
1 T T
TI = 2 [5} [k] f5} - [5} [f}
(2.15)
comparison with Eq. 2.12 yields the element stiffness matrix,
(2.16)
and the generalized force vector
(2.17)
2.3 EVALUATION OF THE BEAM STIFFNESS MATRIX
Following the outline described above, the stiffness matrix
for the beam element is derived. The beam is assumed to be a straight
bar of uniform cross section. Among the many possible, and perhaps
equally acceptable, ways of representing generalized displacements,
the chosen set consists of extensions, bending rotations, transverse
displacements, torsional and warping rotations. The corresponding
generalized stresses, sometimes referred to as stress resultants,
are determined by lumping the integrals of the boundary stresses at
the nodes of the element.
In this study, the generalized force system is reduced to
smaller uncoupled systems in order to demonstrate the application
-18
of the formulation more effectively. This is achieved, without much
loss of generality, by making an appropriate selection of the reference
axes. The flexural and torsional components, for instance, will be
uncoupled by employing as reference axes a right-handed rectangular
coordinate (cartesian) system where o ~ of the axes is oriented parallel
to the element. Of course, this is true only if the material is
linearly elastic, The reference axes system adopted for the beam
element is shown in Fig. 2.1, where the x-axis is directed parallel
to the element.
It has been indicated earlier in this study, that the major
problem associated in formulating the beam problem is in developing
a kinematical model from which the strain-displacement relationship
can be easily established. The formulation of the beam problem is
based on V1asov's hypothesis of the invariability of the beam cross
section(ll). The hypothesis implies that distances between points
on the normal plane of the beam do not change during deformation.
This reduces the strain field to fewer strain components which can be
represented directly by using a polynomial expansion. Polynomials
are used, as in many other situations, because of their simplicity in
-19
manipulation. Based on these approximate fields and utilizing the
finite element model described in Section 2.2, the beam problem can then
be formulated by treating separately the flexural and torsional problems.
Stiffness Matrix in Extension and Flexure
For a linearly elastic material the flexural problem can be
further uncoupled into three systems; two bending components about the
two axes and the extension component, provided the axes of the beam
element coincide with the centroidal-principal axes of the cross section.
However, it is found inexpedient to limit the analysis based on the use
of the centroidal-principal axes, therefore, the general flexural problem
will be analyzed.
The strain field [e} in the beam element when expressed in
terms of polynomial functions of the coordinate [P J and the undetennined
coefficients ~ } is rewritten as
[e} = [PJ ~ }
(2.18)
The assumption that the cross section remains invariant will set the
strain components e... , e ,
yy zz
zero. For flexural problems, the remain-
ing three strain components may be approximated by taking the linear tenns
in x, y, and z of a polynomial. Thus, Eq. 2.18 is written as
e
xx
1 y xy z xz
13
1
0 0 cp(y,z) 0 X (y ,z)
13
2
e
xy
-
133 (2.19)
xz
0 0 X (y,z) 0
134
135
- ~
It is noted that, when the shear strain functions ~ y , z ) , X(y,z)
and W(y,z) in Eq. 2.19 are set to zero, shear strains are eliminated
in the formulation and the strain field will be equivalent to Bern-
oulli's hypothesis: plane sections remain plane where normals to the
reference axis before bending remain normal after bending (elementary
beam theory). Similarly, prescribing the functions ~ y , z ) and *(y,z)
with constant values introduces shear strains in the formulation, and
the resulting strain field will be equivalent to that derived from
Timoshenko's kinematical mode1(31), where points on normals to a
reference axis before deformation remain on a straight line after
deformation (Timoshenko beam theory). Or, when stated differently,
plane sections remain plane but planes normal to a reference axis
before deformations do not necessarily remain normal after deforlnation.
In a similar fashion, various forms of kinematical models may be
developed by manipulating the shear strain functions. Inclusion
of higher order terms of y and z in e will also furnish mainfold
xx
forms of more sophisticated kinematical models.
In the discussion presented so far, it has been tacitly
assumed that the shearing strains e and e are linearly dependent
xy xz
functions of e . It is shown at a later state that this relationship
xx
exists as a result of a minimization process. From a different stand-
point, the established relationship can be viewed as a result of the
overall equilibrium condition of the beam element. The manner in which
these relationships are established is now presented.
The stress field corresponding to the strain field given by
Eq. 2.19 is obtained by employing the constitutive relationships,
-21
<J
xx Dl1 D12 D13
1
Y
xy z xz
Sl
D2l DZ2
D
23
a a ep(y,z) 0 x (y ,z)
Sz
(2.20)
<J
xy
<J
xz D3l D32 D33
La
a x (y ,z) 0 (y, z)
For a beam whose material is elastic, isotropic and the
Hookean constants are taken as
D11
E (Young's Modulus)
])22
=
D33
=:
G (Shear Modulus)
Dij
a (for i =f j)
For different as in the case of the inelastic beam or a
beam of anisotropic material, the appropriate Hookean constants must
be used.
Figure 2.2 shows the stress resultants which are equivalent
to the integrals of the boundary stresses lumped at nodes 1 and 2.
The stress resultants are,
P.
SA <Jxxi
dA

V
= f
<Jxy .
dA
Yi
A
V
=-J
A
<Jxz .
dA
z 0
l.
(2,21 )
M == S Z (J dA
Yi
xx.
A l.
M
= f
y
<Jxx .
dA
z.
A 1.
where i = 1,2 which are the node points of the element.
'-22
The unbalanced stress resultants in the beam element,
determined from the equilibrium condition, become
b,P == I (0'
- O"xx )
dA
A xX
2
I
b,V = b.V == 0
y z
()J;tf.y =1
z(O" - 0" )
dA
A
xx
2 xxI
(2.22)
b,M
z = I yeo-xx
- (J ) dA
A 2
xxI
Or, when expressed in terms of the strain coefficients they are written
in matrix notation as
b,P = ELI [y
A
z] {f3
3
} dA
f35
b.M = ELI [yz
Z2 ]
{::}
dA (2.23)
y
A
b.M 0: ELI [y2 yz]
{S3}
dA
z
A
f3s
where L length of the beam element.
These unbalanced forces are counteracted by
shearing forces at the nodes. The magnitude of the shearing stress
resultants [V}, which satisfy the equilibrium condition of the element,
are
,.., 1
rV} = - f()J;tf.}
. L .
(2.24)
or, in terms of the axial forces [b,P} as
""' 1 '"
[V} = L [eJ [b. P}
(2.25 )
-23
where [e] consists of the associated distances from the reference
axis to the resultants of Substitution of Eq. 2.23 in Eq.
2.25 will furnish the shearing stress resultants, which are required
for the overall equilibrium condition, in terms of and which are
written in the form
[v} == E [J] (2.26 )
where
__ [Y2 yz]
[J] --
yz z2
dy dz
and
It is noted that [J] is equivalent to the inertia matrix of the cross
section.
In the same manner, the shearing stress resultants of Eq,
2.20 may be expressed as
where
[PJ '" J
A
[V] G [PJ
[
c:p(Y'Z)
X(Y'z)
x (y,z)]
*(y,z)
(2.27)
dy dz
Obviously, the shearing stress resultants given by Eq. 2.27 are
alent to those given by Eq, thus
[V] [V] (2.28)
At this stage, a vector consisting of the average shearing
strains [v} is introduced to account for the variation in shearing
~
-24
strains over the cross section. The expression for [y} can be written
in the form
[v}
1
= AG ~ ] [V}
1 ~ ~
= A ~ ] [PJ ~ }
(2.29)
where
~ ]
z z
The matrix [a] is composed of numerical factors which are commonly
known as shear deformation coefficients. A coefficient detennines the
shear deformation, by considering an average value of the shear induced
transverse displacement, due to a transverse shearing force. For
example, a is the coefficient in the y direction due to a shear force
yy
in the same direction, a is the counterpart of a in the z direction,
zz yy
and ayz is the coefficient that determines the shear deformation in the
y direction due to a force in the z direction, and vice versa.
In order to evaluate the elements in raJ, it is necessary
to establish the general solution of the problem of bending by terminal
transverse loads. The customary approach to the solution of this
problem, based on the semi-inverse method of St. Venant, has been given
(2 32 33) . (34)
by several authors ' , . A recent contribution reduces the
elasticity solution, by introducing appropriate simplifying assumptions,
in order to derive a formula for shear coefficients which is applicable
to symmetric shapes only. Based on this formula, numerical values
of [a] are calculated for simple geometric cross sections. The latest
..,25
contributions (35,36) , furnish a numerical solution based on displacement
formulation by a finite element technique.
Comparison of Eq. 2.29 and Eq. 2.26 yields the following,
(2.30)
Substitution of Eq. 2.30 in the strain field equations given by Eq.
2.19 results in defining the strain field explicitly as
[e} = [PJ fS}
e
xx
1 y xy z xz
13
1
Er
13
2
e
xy
0 0
-Ii.
0 133
(2.31)
-
GA
..
134
Er
135
0 0
~
0
e
xz
GA
where
r =
OI
yy
J yy + OI
yz
J
yy zy
ryz
=
OI
yy
J
+ OI
yz
J
yz zz
r
zy
OI
zy
J
+ O z ~
J
yy zy
r
zz
=
OI
zy
J + 01 J
yz zz zz
or, in tensor notation it may be expressed as
r ij OIik J
kj
where the summation convention for repeated indexes is employed ..
Once the strain field is defined explicitly in terms of the
undetermined coefficients, the stiffness matrix can be determined by
following the outline described in Section 2.2. In Eq. 2.16 the stiff-
ness matrix [kJ is expressed in terms of the matrices [H] and [TJ.
These matrices are determined, by integrating over the volume of the beam
element, the nmtrices [P], [DJ, [R] and [LJ which are defined in Eq.
2.7 to 2.9.
The [II] is written in the form
[H] = S [PJT [D] [PJ dV
V
(2,12 )
where the matrix [P] is given by Eq. 2.31 for the problem at hand.
Substitution of [PJ results in the following symmetric matrix
r-
E
Ey
Ey2
[H] = S
V
Exy Exy2
Ez Eyz
Exz Exyz
-
Ex2y2+G[r 2
yy
Exyz
-
Ex2y
2
+G[r r
yy
+r 2J
yz
EZ2
SYlYlME TR Ie
EX2Z
2
+G[r 2+r 2J
zy ZZ
dV
(2.32 )
Obviously, the matrix reduces to a diagonal matrix when the centroida1-
principal axes of the cross section are used.
In order to determine [TJ, the matrices [R] and [LJ must be
evaluated first. The matrix [RJ is obtained by relating the boundary
force vector
f CJb}
in terms of strain coefficients
.
The boundary
force vector
[CJb}
consists of six elements representing the y and z
components of the boundary forces at the two nodes. Thus.
fO"b} = [RJ [ ~ }
O"xx
-E -Ey Ety -Ez Etz
~
1
O"XYI
0 0 -Er fA 0 -Er /A
13
2 yy yz
O"xz
0 0 -Er fA 0 -Er fA
~ 3
1
zy zz
=
(2.33)
O"xx
E Ey Ety Ez
Etz
134
2
C1
xy
0 0 Er /A 0 Er /A
~ 5
2
yy yz
O"xz
0 0 Er fA 0 Er /A
2
zy zz
The matrix [LJ relates the displacements at the boundary
,-w
[u}
and the generalized displacements Co}. For the problem at hand,
r>J
there are six prescribed displacements at the boundary [u} which
are related to the ten generalized displacements [oJ at the nodes in
the following form
fu} = [LJ fa}
1 0 -y 0 -z
01000 o
w
l
0 0 0 1 0
=
(2.34)
u
z
1 0 -y 0
v
2
0 010 0
w
2
000 1
Once the matrices [RJ and [LJ are determined, the matrix
[TJ is obtained by integrating the product over the boundary
-28
[TJ = J
[RJT [LJ dy dz
A
-E 0 Ey E Ez E 0 -Ey 0 -Ez
-Ey 0
Ey2
Ey Eya Ey 0 -Ey2 0 -Ey2
.=J
E.f.,y -Er /A -E.f.,y2 E.f.,y -E.f.,ya E.f.,y Er /A -E.f.,r E /A -E.f.,yz dy dz
A
yy yy zy
-Ez 0 Eyz Ez EZ2
Ez 0 -Eyz 0
_EZ2
E.f.,z -Er /A
yz
-E.f.,yz E-l,Z E-l,z2 E.f.,z Er /A E-l,yz
yz
Er /A-E.f.,Yz
zz
(2.35)
At this stage, the stiffness matrix [kJ can be determined
from Eq. 2.16 since the matrices [HJ and [TJ are known and are given
by Eq. 2.33 and Eq. 2.35, respectively. This matrix is not evaluated
here in the manner described above since it requires a rather tedious
manipulation consisting of manual integrations and matrix operations.
However, review of the derivation process discloses a complete sequence
of numerical integration and matrix operations, which can be performed
in a systematic manner, by developing a suitable computer program.
Alternatively, the manipulation required to evaluate the
general stiffness matrix is significantly reduced from the viewpoint
of manual computations, by performing transformations to the stiffness
matrix computed for the centroidal-principal axes. For this particular
set of axes the non-diagonal elements in [HJ will vanish and most of
the elements in [TJ reduce to zero, thus the required matrix operation
in Eq. 2.16 is simplified. Once the stiffness matrix for these sets of
axes has been evaluated, the corresponding stiffness matrix for an
arbitrarily assigned set of axes can be easily established following
-'
\
v'
V
-29
the standard transformation procedure of stiffness matrices used in
structural mechanics.
The computation for the element stiffness matrix [kJ corres-
ponding to the centroidal-principa1 axes results as given below:
PI
M
L
12EJ
V
01
0
-L
L3 (Hi
o
)
MYl
6EJ (4+t
z
)EJ
X
0
--L..
L" (Ht
z
> L(l-lt
z
>
V)'l
0 0 0
Mol
0 0 0
P
2
EA
0 0
-1"
v
02
12EJ 6EJ
0
.--L..
-.----"---
L3(H ) L" (H )
MY2
0

(2-'
z
}EJ
X
IP (l-lt
z
> L(I-1t
z
>
V)'2
0 0 0
M02
0 0 0
where
and
[f} = [kJ [a}
12EJ
0
L3 (Itiy)
6EI
z
- L" (l-lty)
0
0
0
12EJ '
z
- L3 (Itiy)
6EJ
z
- L" (ltiy)
(4ti
X
)EJ
z
L(ltiy)
0
EA
1"

0 0
L3 (1+-q;.)
0 0

L" (ltiy>
6EJ
z
0 0
L" (l+ty)

0 0
L(lti.y)
12Er
zz
=
<P
z
=
12Er
yy
SYMMETRIC
(4-1t
z
)EJ
X
L(Ht.)
12EJ
0
z
L3 (1+ly)
6EJ (4+' y>EJ.
0
0
IP (1+.
y
) L(l+.y)
-
In the above the terms rand r are defined in Eq.
yy zz
(2.36 )
2.31. It is observed that use of the centroida1-principa1 axes sets
the cross terms such as r ,J , ,etc. to zero thus resulting in
yz yz yz
a more simplified form of stiffness matrix.
-3U
In order to obtain the general stiffness matrix, for an
arbitrary set of axes, the stiffness matrix [k] given by Eq. 2.36
is subjected to a transformation. The new set of axes which pass
through po in t P, is shown in Fig. 2.,3, where the axes are trans lated by
v and w from the original point 0, (the centroid) and are rotated
p p
by an arbitrary angle a"
The transformation matrix, designated by [T ]. is obtained
r
by expressing the displacement field at the boundary as coordinates
of each reference axes. This relationship is,
u
l
I 0 -v 0 w
til
p
P
WI
0 -sina 0 cosa 0
WI
eyl
0 0 -sina 0 COSCy 0
e
yl
VI
0 COBa 0 Bina 0
VI
e
zi
0 0 COBa 0 sina
ez 1
-
-
u
2
1 0 -v 0 w
ti2
p
P
w
2
0 -sina 0 cosO! 0 w
2
e
y2
0 0 0 -sina 0 cosO!
9
y2
v
2
0 cosQ' 0 sinQ' 0 v
2
9
z2
0 0 COBa 0 sina 9
z2
(2.37)
Since the transformation matrix [T J is orthogonal the
r
general flexural stiffness matrix is obtained from
-31
(2.38)
Stiffness Matrix in Torsion
In this section, the stiffness matrix is derived for a beam
element having the characteristics of "thin-walled beams". A thin-
walled beam is composed of plates whi.ch are assumed to undergo in-
plane strains alone when subjected to loads. The theory of torsion of
thin-walled beams, unli.ke solid beams, has as a distinctive feature
the occurrence of considerable axial strains as a result of torsion,
The general theory of thin-walled beams as developed by Vlasov(ll)
is essentially based on the assumptions that the cross section remains
undeformed and the shearing deformation in the middle surface vanishes.
In this study, the cross section will also be assumed to be rigid but
the shearing strains at the middle surface are not neglected.
The stiffness matrix for the beam element subjected to
torsion is derived by following the same procedure adopted in the
flexural problem. Figure 2.4 shows a component plate of the beam model
used for deriving this stiffness matrix. In order to simplify the
computatrons, the z-axis is oriented parallel to the mid-surface of
the element and the x-axis passes through the center of torsion.
The strain field for a torsional problem is approximated
in a similar manner shown in Eq. 2.19 as
[ e} :: [PJ ~ }
e
xx
z zx 0
S6
S7
(2.39)
xy
0 S (z) (y-a)
Sa
It is noted that the shear strain function is assumed constant
through the thickness. But the variation along the z-axis may be
accounted for by introducing an average shearing strain -; and a shear
xz
deformation factor CY.w in a similar manner described earlier in this
section. Thus, the approximated strai.n field is defined explicitly
as follows
z xz 0

-32
.xxl _
(2.40)
xy f -
ex.;" E J
becomes
0
_. _......::a.
(y-a)
GA
The corresponding [H] matrix for the given strain field
EZ2
[H] = S EXZ2
V
o
SYMMETRIC
lGA)
-w yy
2
G(y-a)
dx dy dz
Relating the boundary force vector [crb} and the strain
coefficients [s} yields the matrix [R], thus
[crb}
= [RJ [S}
O'xxl
-Ez E,f,z 0
S6
crxzl
0
-CY.wEJ I A -G(y-a)
yy
- S7
-
O'xx2
Ez E,f,Z 0
O'xz2
0 Cl. EJ IA G(y-a)
Sa
w yy
....,
(2.41)
(2.42)
There are four prescribed displacements [u} at the boundary
which are given in terms of four generalized displacements [a} at the
-33
nodes. The re lationship is written in the form,
[tl} == [LJ fa}
o z I
I 0
(
1
ep
1
-a-2 (y-a) 0 I
tol
I
= ------r----- (2.43 )
I 0 .
0
eP2
o I
I =a=2(y-a) 0
W2
It is noted that the term (y-a) in the displacement 11] g ivan by Eq < 2.43
is multiplied by a factor of 2. The Hell known explanation given by
Lord Kelvin and Tait (32 ,37) on the ma.nner in which an applied torque
is resisted by a bar having a thin rectangular cross section, refers
to the stress distribution in the cross section. They have explained
that one-half of the applied torque is carried by a system of shearing
stresses parallel to the longer dimension of the cross section. The
other half, is carried by the transverse shearing stresses which act
normal to the plate. These transverse stresses are normally neglected
even though they become of an appreciable magnitude near the short
sides of the rectangle. However, since they act at a greater distance,
their contribution to the torque is significant and thus they constitute
the other half. In this study, the aforementioned difficulty is cir-
cumvented by establishing an appropriate kinematical mode 1 which will
yield a solution consistent with the approximate elasticity solution.
The matrix [TJ is computed as follows,
[TJ = S [RJT [LJ dy dz
A
.:.34
o EZ2 -EZ2
[TJ = J
A
EJ fA)
W yy
o
EJ fA)
(I) yy
=E.{,Z2 dy dz
-G(y-a) (a-2y) o o
(2.44)
The general stiffness matrix is obtained by substituting
Eq. 2.41 and Eq. 2.44 in Eq. 2.16. If the y-axis passes through the
centroid of the element, the stiffness matrix reduces to the following:
MT1
M(I)I
=
MT2
M(I)2
where
(f} = [kJ [5}
u
2
l+f2(l+iJ?(I)
SYMMETRIC
L
L2
I2EI 2
12 (4+gj (I)
!!.l
L3 (1+iJ? )
u
2
L u
2 (I)
-1--(1+iJ? )
2
1-+]2 (1 +ip (I) )
12 (J)
L
)
L IP
-Z
2
12 (L!.i1j )
12 ill
- ill
I = warping moment of inertia about the shear center
(I)
KT
=
St. Venant torsion constant
12O!EJ
iJ?
=
y.y.
W
GAL2
CPl
wI
Cfl2
w2
(2.45)
Of course, for profiles with zero warping rigidity (I =: 0), the factor
ill
U must be expressed such that EI does not become a denominator.
ill
The torsional stiffness of a beam element may be derived
consistently from the displacement approach also as a minimization of
-35
the total potential energy. In the following, two different forms
of displacement functions are considered, namely, hyperbolic functions
and polynomials. The stiffness matrices derived from these two forms
of displacement functions are compared to that given by Eq. 2.45.
Hyperbolic Functions
According to V1asov's beam theory(ll) the governing differen-
tial equation for the beam element in torsion is given by
(2.46)
The solution of Eq. 2.46 yields the following displacement functions
written in matrix form as
x
- [:
1
where
cosh kx
k sinh kx
(l)==CP
,x
k =lK
T
'
EI
Ul
sinh kxkx:l
k cosh J
Corresponding to the displacement function given by Eq. 2.47 the
11 ff ix b d
. d(38,39)
fo owing sti ness matr can e er1ve
(2.47)
-36
CPl
S YNI"IE TR I C
~
11.
2
L( l-coshx)
-sinhx}
-11.
3
(l-coshu) -x.L:d (l-coshx)
-XL2 (l-coshx)
-sinhx)
(2.48)
where
D = 2 (1 - coshx.) + 11 s inx.
and x = kL
Polynomials
The displacement field for the beam element may be expressed
by a polynomial of the third order as follows,
[u} = [PJ fO'}
cP
1 X
X2
x
3
0'1
=
0'2
Er
OJ
0 1 2x
( 3X2 7)
Q'3
Q'4
In order to take account of shearing deformations due to.warping,
an appropriate term is incorporated in [PJ as described earlier in
this section.
The vector fa} consists of the coefficients which are to be
determined in terms of the nodal displacements [a} from the relation-
ship
-37
{5} = [C] fa}
(2.50)
The strain field fs} corresponding to the displacement field of Eq.
2.49 may be written as follows,
sxx
-azq>
,xx
z(y-a)cp
Er
,x
+ --i'f-cp,xxx
[e} = [Q] fa}
o o
o 2(y-a)
-6xz
01
1
Er
0,12
6[x
2
(y-a)+cfl 0'3 (2.51)

Following the standard finite element procedure, the stiffness matrix
can be obtained from
(2.52)
Thus,

x2 1 2
l-+1z[S+(l+'P
w
) J
SYMMETRIC
L x
2
L2 a

I2[3+(1+<P )
w 2
MW1 ,wI
12EI
_ II)
- L3 (1+<P )2
._ W
x2 1 2
-=+( 1 +<P ) ]
12 5 w
L X

(1+<P )

4 5 3
L X
2[l+tQJ
L2 a
-[3-(1+ )
12 W a
(l+<p )
-tlf,('!:. U!
J
4 5 3
XS 1 a
1+'i2[s+(1+'P
w
) J
L x
La 2
I[l+tQJ
12 (3+(1+ )
W a
(l+<p )
ofl!:.(l, m
J
453
(2.53)
If warping shearing deformations are ignored, the parameters are
(l)
set to zero, and the stiffness matrix will reduce to the following,
CPl
SYMMETRIC
l2EI
-38
=
L3
(2.54)
1{a
-I--
I 10
The stiffness matrix g:i "en ,)y Eq. 2.54, wh :ch is a special case of
Eq. 2.53, is identical to previous derivations(40).
2.4 APPLICATIONS AND SAMPLE SOLUTIONS
Before any finite element solution can be used with confidence
some idea of its accuracy and are required.
Suitable evidence is usually provided the finite element
results with accurate results derived by alternate means. To illustrate
the validity and application of the method the finite element formu1a-
tion developed is applied to a number of the few numerical examples
whose analytical solutions are straightforward. The procedure of
analysis is based on the displacement method which is adequately covered
. bl i (18,19,23)
many pu ons .
The structural member under consideration is suitably idealized
by a set of basic beam elements with a 7-degrees-of-freedom at each node.
The stiffness matrix for such an element is given in Appendix I when
the centroidal-principal axes system is used. Once the member is
idealized as an assemblage of beam elements, the over-all unconstrained
-39
structural stiffness matrix is generated following the rules that
govern the assembly process used in the matrix analysis of framed
structures. This matrix is generated by a simple summation of the
individual stiffnesses and loads at the nodes using nodal compatibility
for this process. Alternatively. the variational concept may be used
on the entire assemblage to derive a mathematical statement of the
assembly rules. The assembly rules for the assemblage of the
ness matrix and load vector is written as
N
[KJ
=
[kiJ
i=l
(2.55)
N
[FJ
[fiJ

(2.56)
where N is the total number of elements.
It is evident that structural members are subjected to bound-
ary conditions in forms of tractions or displacements. The traction
boundary conditions are incorporated automatically into the load vector
[F}. When imposing the displacement boundary conditions, the standard
procedure, which involves eliminating the equilibrium equation at
which the particular are specified, results in reducing
the size of the master stiffness matrix and thus requires reorganization
of the computer storage. However, the same conditions can be imposed
without changing the size of the matrix, simply by modifying the
ness matrix and the load vector. This is accomplished by multiplying
with a very large number the element on the diagonal of the matrix [KJ
at the location concerned and also by replacing the corresponding element
I .
in the load vector fF} by the same large number multiplied by the
specified displacement (15) This procedure applies whether the p r ~
-40
scribed displacements are homogeneous or nonhomogeneous. For the case
of elastic restraints, the matrix [KJ is modified by adding the sup
port stiffness on the appropriate matrix element on the diagonal of
the stiffness matrix.
The resulting equilibrium equations of the complete system
are expressed in the form
(2.57)
in which the number of simultaneous equations in the preceding relation-
ship is equal to seven times the number of nodal points. The nodal
displacements f6} are unknown and are determined by solving the set
of simultaneous equations (Eq. 2.57) and then the stress resultants
are evaluated by using the relationships of the individual elements.
Numerical Examples
The first numerical example to illustrate the application
and validity of the procedure described is that of a cantilever beam
subjected to a concentrated torque at the free end. The beam repre-
sented in Fig. 2.5 allows warping and twist at one end and these dis.
placements are restrained at the other end. The twist and warping
displacements are computed at the nodal points for different values
of the parameter K = /G K
t
cross-sectional properties.
L2/EI' in order to cover a wider range in
W '
The computed values plotted in Fig. 2.5
agree very well with the charts given in Ref. 41.
-41
The second example consists of the same beam used in the first
example but loaded with a concentrated bimoment at the free end instead
of a torque. The results are plotted in Fig. 2.6 where close agreement
is observed with those given in Ref. 42.
In order to compare the differences in solutions that may
arise when using the element stiffness matrices given by Eqs. 2.45,
2.48 and 2.54, the cantilever beam used in the earlier examples was
used again. The stiffness matrices are determined based on three
different formulations, namely, the strain field formulation, the
polynomial formulation and the hyperbolic functions formulation as
described in Section 2.3. The cantilever beamwas loaded by a con-
centrated torque and bimoment at the free end as shown in Fig. 2.7,
Different values of the cross sectional parameter nwereused in order
to cover a wider spectrum in beam characteristics ranging from those
having resistances in pure warping to those in pure torsion (St. Venant
torsion). As indicated in Fig. 2.7, good correlation is observed for
the values of n normally regarded as thin-walled beams. However, for
larger values of the parameter n, or when I approaches zero, the
w
differences in the solutions increase and the computational errors
grow when using the stiffness matrix based on the strain formulation
(Eq. 2.45). Moreover, for this particular formulation, the results
oscillate for larger values of n as shown in Fig. 2.8 whereas good
agreement is observed for values of n less than 10. This problem of
numerical instability arises from the fact that the strain field
formulation necessarily assumes that the cross section has warping
resistance, I. This is not regarded as a serious practical problem,
w
since useful solutions can be obtained using ordinary analyses once
it is known that the cross section has no warping resistance.
The final numerical example is intended to demonstrate the
versatility of the method by considering a continuous beam with two
equal spans subjected to a concentrated bimoment M = 1.0 acting at
m
the right support as shown in Fig. 2.9. The shear centers of both spans
are assumed to form one straight line which is considered as the axis
of the beam. For both spans, the cross section parameter u is assumed
equal to 1.0. The variations of the torque MT and the bimoment M
W
are shown in Fig. 2.9. It is important to distinguish the two components
of the torque MT which are present in thin-walled beams. The torque
Tsv is the more familiar St. Venant's torsion, the other component Too
is the warping torsion. The latter results not from warping but from
the suppression of warping.
3. LINEAR STABILITY ANALYSIS OF BEAM-COLUMNS
3.1 INTRODUCTION
The theory of linear stability, which is based on the concept
introduced by Euler, represents closely the circumstances of failure
of beam-columns. Although the actual conditions of failure require
the inclusion of nonlinear influences for the precise determinationat
the failure condition, as in determining the complete response of an:'
imperfect column, the critical condition obtained from linear stability
analysis is useful from the design standpoint. Linear stability analysis ,-
is defined here as the calculation of the bifurcation of equilibrium,
the point of bifurcation occuring at the critical load which is charac-
terized by the existence of a fundamental state of equilibrium.
In the conventional linear stability analysis of structural
problems, two approaches are normally used. The first approach is to
determine the lowest eigenvalue of the governing differential equations
of the structural system for a given set of boundary conditions.
Alternatively, if the governing differential equations are too diffi-
cult to prescribe, numerical methods are utilized by establishing a
strain energy expression for large deflections which is subsequently
minimized, leading to roots representing instability conditions.
The use of matrix methods in solving problems of stability
based on the concept of discrete element idealization has recently
received considerable attention. In particular, displacement formu-
lations based on finite element idealization are found to be more
suitable. The adoption of the matrix force method has been
-44
(19) but is employed to a lesser extent. From the viewpoint of sta-
bility and finite displacement analyses, the possibility of incorpora-
ting geometric nonlinearities within a displacement method offers a
suitable means of utilizing the finite element concept and its appli-
cations.
i
Moreover, the finite element approach plays an important
j:
role through its ability to lead to solutions to problems with irre-
I
gularities in loading and geometry which defy adequate treatment by
the classical means. Another important characteristic for the use of
finite elements has been their intrinsic simplicity. Useful r v i ~ w s
I
I
of the accomplishments in finite element stability analyses are found
in Ref. 43 to 47.
I
, 1
As in other aspects of the finite element method, the treat-
ment of problems dealing with linear stability consists of two component
parts: the formulation of the element relationship and the solution of
the complete system. Furthermore, the formulation of the element
relationship involves the calculation of corrective terms to the linear-
ized equations. Consequently, application of the conventional matrix
displacement methods to problems in elastic stability has been concerned
with the derivation of so-called geometric stiffness matrices to account
for the instability effects. The inclusion of the geometric stiffness
matrices into the formulation is performed by adding them directly to
the elastic stiffnesses to form a resultant stiffness matrix. The
derivation of the elastic stiffness matrix for the beam element is
given in Sectwn 2.3. In the following section, a formulation is
presented for deriving the element geometric stiffness matrices.
-45
3.2 A FINITE ELEMENT MODEL FOR DERIVING GEOMETRIC STIFFNESS MATRICES
In this section a general method for evaluating geometric
stiffness matrices for the stability analysis of general discrete
structural systems is presented. The formulation provides a means for
a direct determination of geometric stiffness matrices that are con-
sistent with any kinematically admissible displacement field assumed
for the element. The derivation of stiffness matrices often is based
on the approximate displacement field, as defined by suitable interpo-
1ation polynomials or shape functions [N] of coordinates, and a set
of nodal parameters [oJ, element by element, as
[u} = [N] [5}
(3.1)
The displacement field may also be defined more conveniently by poly-
nomial functions of the coordinates [P] and a set of generalized
coordinates or generalized displacement amplitudes fa} expressed by
[u} = [P] fa} (3.2)
A relationship between the vectors to} and fa} is established using a
displacement transformation matrix [C]. This matrix is determined by
substituting the coordinates of the nodes in Eq. 3.2 in the following
form
to} = [C] fcd
(3.3)
When performing a stability analysis, use of the non-linear
-,-
total strain-displacement equations is essential and leads to relevent
(48 49)
solutions ' . The strain-displacement equations as given by the
general theory of deformations or known as the Lagrangian
. (35) .. .
tensor are tensor as
-46
(3.4)
In Eq. 3.4 the summation convention for repeated indexes is employed and
a comma denotes partial derivation with respect to the variable that
follows. A typical strain component, frequently used in structural
analysis is the axial strain, e , and can be expanded from Eq. 3.4
xx
which is written in standard mathematical notation using the engineer-
ing definition as
(3.5)
In developing the strain-displacement relationships for beam-
type bodies it is necessary to convert formally the relationships
given by the three-dimensional theory of deformations into their one-
dimensional analogue. An attempt is made here to formulate the
stability problem of beam-columns by making use of the strain-displace-
ment relationships provided by the classical theory of thin-walled
b
(11,50,51)
earns
The displacements of a thin-walled beam of rigid cross-section
is adequately described by the lateral displacements v and wand by
the The displacements at any point on the beam are functions
of the coordinate x and are given by
u = - z
oX
Q!.+ (
y oX t y,z 'ox
v = vex) + z
w = w(x) - y
(3.6)
-47
Where u, v and are arbitrary functions of the coordinate x and
is the warping function. Based on the assumption that the elongations
are negligible compared to unity, Novozhilov(48) has developed the
strain-displacement relationships for thin rods by expanding the dis-
placements u, v and w (Eq. 3.6) in series in the coordinates of y
and z of the points of the cross section.
The expressions for the nonvanishing components of strains
which become adequate for the problem at hand are:
= e + z x (x) + y x (x) + ,I,(Y z)
xx zz yy '1" dx
(3,7)
e = e + - z) rn
xy xy oY 't'
where e
ij
are the strain components given by Eq. 3.4, and xzz(x) and
x (x) are the curvatures of the deformed axis of the beam which take
yy
into account large displacements in u, v The expressions for
the curvatures for large displacements are found in Ref. 48.
The strain vector given above may be resolved into two com-
ponents and is written in matrix notation
(3.8)
where fe } is the usual linear, infinitesimal strain vector while
o .
reL} represents the non-linear strain contribution. It is well known
that the mere presence of [eL}' without regard to magnitude, has a
decisive influence on the behavior predicted in stability situations.
-48
Having established the displacement field it would be logical
to define the strain field in terms of the same parameters ~ } used in
Eq. 3.2 which may be written collectively as
(3.9)
The matrices [Q] and [QL] introduced in Eq. 3.9 represent derivatives
of the displacements corresponding to the linear and non-linear strain
contributions, respectively.
Following a similar development normally used in finite
element formulation, the strain energy in the new configuration of
equilibrium is evaluated as
u = t f [e}T [D] [e} dV (3.10)
V
The matrix [D] represents the generalized Hookean constant. On sub-
stituting Eq. 3.9 into Eq. 3.10 and neglecting the non-linear strain
product feL}T [D] [e
L
}, since it is of much higher order, the strain
energy functional reduces to
(3.11)
where a new stress matrix
has been introduced to denote the stresses corresponding to the linear
-49
From Castigliano's first theorem, which is applicable to
non-linear strains provided that the total strain energy is evaluated,
the following relationship is obtained
(3.12 )
In performing the differentiation with respect to the generalized
coordinates fa} the linear stresses [a } have been assumed to remain
o
constant. Introducing the displacement transformation matrix [C]
into Eq. 3.12 yields the force-displacement relationship
(3.13)
where
(3.14 )
is the usual stiffness matrix obtained by the linear theory, and
(3.15 )
is the geometric stiffness matrix. The geometric stiffness matrix
derives its name from the fact that it depends on the geometry of the
displaced element. It is noted that the geometric stiffness matrix
can easily be determined from an integral of simple matrix products
evaluated over the volume of the element. The approach avoids the usual
procedure of determining strain energy in terms of displacements and
its subsequent differentiation with respect to the displacement, which
in this case would be more time consuming. Furthermore, the formulation
allows for a systematic investigation of the effects of higher order
terms in the strain-displacement relationship by introducing the
appropriate matrices.
Application of the finite element model developed above is
made to derive the geometric stiffness matrix of the beam element.
This matrix may be derived in a single operation by formulating one
general strain expression which includes all the strain components
corresponding to the generalized displacements and substituting it
into the strain energy functional given by Eq. 3.11. Although this
approach may seem to offer simplicity it suffers from the drawback
-50
that the required computations become cumbersome. Alternatively, the
derivation may be carried out more conveniently by treating separately
each of the large displacements that introduce geometric nonlinearity.
The stiffness matrices corresponding to each form of large displacement
are finally aggregated together to constitute the general geometric
stiffness matrix for the beam element.
The displacements that introduce geometric nonlinearity as
they become 'large' may be put under three categories:
a) Axial displacements
b) Transverse displacements
c) Twist
In the following sections, the constitutive geometric stiffness matrices
[kGJ
i
are derived corresponding to each large displacement given above.
The derivation of the usual linear stiffness matrices [kEJ
i
is not given
here; it is presented in Section 2.3.
-51
3.3 LARGE AXIAL DISPIACEMENTS
When a beanrco1umn is subjected to an axial force, it is well
known that the stiffness of the beam in flexure and torsion become
dependent on the magnitude of the applied axial force. In the following,
the influence of the axial force on the flexural and torsional stiff-
nesses of the beam element are derived separately.
Flexural Geometric Stiffness Matrix
Cons ider a uniform beam element shown in Fig. 3. la. The element
sustains only flexural displacement in the x-y plane under the generalized
forces P, V and M applied at the nodes. A simplified model is used
y z
here merely for convenience; however, a model capable of sustaining
flexural displacements in both x-y and x-z planes simultaneously could
also be used without introducing significant complications.
For the beam element shown in Fig.3.la, the axial displacement,
u, is taken to be adequately represented by a linear polynomial and a
cubic polynomial is assumed for the transverse displacement, v. The
assumed displacement field is written in matrix notation as, (52)
[u}
=
[PI fQ'}
u I
r;
0 0 0 0
Q'l
=
0 0 1 r;2 r;3
Q'2
(3.16 ) v
r;
Q'3
e
0 0 0 1 2r;
3r;2
Q'4
Q'5
Q'6
where r; = x/L
-52
The approximate strain-displacement relationship based on the
ordinary beam theory (Bernoulli's hypothesis) and which, In addition,
included the predominant component of Eq. 3.5 to account for the large
transverse displacement is given by
e . + e = [u - yv ] + [.! v2 ]
o L .,x ,xx 2 ,x
(3.17)
Note that e is a nonlinear transformation under the assumption that
xx
the strain due to midline rotation is not small when compared to the
midline axial strain.
The matrix [QLJ is determined by expressing the nonlinear
strain component [eLl in terms of the parameter [a} as expressed in
Eq. 3.9, thus,
0 0 0 0 0 0
0 0 0 0 0 0
[QLJ
1 0 0 0 0 0 0
(3.18) =
2"
0 0 0 1
2S
3
S
2
0 0 0
2S
4;2 6;3
0 0 0 3S;2 6;3 9g
4
The transformation matrix [CJ is determined by substituting
the coordinates of the nodes (Eq. 3.2) into the equations of the
displacements given by Eq. 3.16. At this stage, the matrices [QLJ
and [C] are known. Finally, by substituting these matrices into Eq.
3.15, and integrating over the whole volume of the beam element the
geometric stiffness matrix is evaluated as
-53
0
SYMMETRIC
0 6/5
0 L/10 2If3 /15
(3.19)
P
[kG]
= -
L
0 0
0 -6/5
0 L/10
0
-L/lO
-L2/30
0
0
0
6/5
-L/10 2L2/1S
The result agrees with that found in the (52)
Torsional Geometric Stiffness Matrix
A uniform beam element is considered here (Fig. 3.2b) which
sustains only torsional and the associated warping displacements under
the. generalized forces P, MT and Mw at the nodes. For convenience,
the x-axis is chosen such that it passes through the shear center of
the beam section. As in the flexural analysis, the axial displacement,
u, is taken to be adequately represented by a linear polynomial; the
twist, is represented by.a cubic polynomial and a quadratic polynomial
is assumed for the warping displacement, w. The assumed displacement
field when written in matrix notation is
fU} = [PJ [a}
u 1 0 0 0 0 0
al
0 0 1
a2
(3.20)
qJ
=
S
a3
ill
0 0 0 -1
-2s -3
S
2 a4
Q/S
Q/6
where S = x/L
-54
It is noted that the displacement field given above is identical to
Eq. 3.16, likewise, the transformation matrix [C] will be the same.
For open thin-walled beams, which are composed of plates
assumed to undergo in-plane strains only, a strain-displacement rela-
tionship can be established based on the assumption that the cross
section remains undeformed. The approximate axial strain, e ,based
xx
on Vlasov's beam theory, and which includes, in addition, the second
order effects of large twists is given by
n n
= eo + e
L
= [u,x - I: p. z cp J +[-21 I: (Pi
2
0)2 + ri
2
0)2)J (3.21)
i=l 1. ,xx i=l
where n = number of component plates of the shape
p = the perpendicular distance between the. shear center
of the section and the middle line of a plate
r = the distance from the center of twist of a component
plate (due to St. Venant torsion) to a general point
on the middle of the plate
The kinematics of the beam section under torsion, from which the
strain-displacement relationship was established (Eq. 3.21), and the
dimensions defined above are shown in Fig. 3.2.
By expressing the nonlinear component of the axial strain
(Eq. 3.9), e ,in terms of the parameters fa} the matrix [QLJ is
xx
determined, thus
0 0 0 0 0 0
0 0 0 0 0 0
(p. + r.)
0 0 0 0 0 0
[QLJ
=
1. 1.
(3.22)
2
0 0 0 1 2g
3 ~ ?
0 0 0 2g
4g2 6g
3
0 0 0
31='2
';)
6g
B
9S
4
-55
, '
At this level, the matrices [QLJ and [CJ, required for the computation
of [kGJ, are known; on substituting these into Eq. 3.15 and performing
the required integration over the volume of the beam element, the
torsional geometric stiffness matrix of a beam element under axial
load is obtained as
0
SYMMETRIC
0 6/5
PI 0 -L/10 2L2/15
[kGJ
0
= --
LA
0 0 0 0
0 -6/5 L/10 o . 6/5
0 -L/10 L2/30 0 ,L/lO 2L2/15
where I = polar moment of inertia of the section about the shear
o
center.
3.4 LARGE LATERAL DISPLACEMENTS
(3.23)
When a beam is subjected to a major axis bending, the minor
axis flexural stiffness and the torsional stiffness become dependent
on the applied moment. For a perfect member, there is a critical load oc
which the beam buckles in a combined mode involving twist and lateral
deflection. In the following, the influence of the major axis bending
on the stiffnesses in torsion and in flexure about the weak axis are
investigated.
Flexural-Torsional Geometric Stiffness Matrix
Consider a uniform beam element that undergoes only trans-
lationa1 displacements v, in the y direction, when subjected to unequal
moments Mil and Mz2 at the two nodes (Fig. 2.2). Just prior to
buckling, the initial axial stress at any point in the beam, accord-
ing to the elementary beam theory, is given by
where
M Y V xy
= ~ ~
0"0 I I
M
oz
z z
1
= -(M + M
z2
)
2 zl
(3.24 )
Obviously, the initial stress 0" must be modified for the case when,
o
in addition, the beam element carries distributed loads between the
nodes. For instance, the additional term becomes quadratic in X for
a uniformly distributed load. However, such additional terms may not
be required whenever a finer discretization is used and a proper
lumping of the nodal forces is performed.
At the critical state, the adjacent equilibrium configuration
assumes lateral displacements associated with twisting. For a finite
element formulation, these displacements may be represented adequately
by polynomial functions which are written in matrix form as
[u} =
[PJ f Q/}
v 1
S
S2
g3
0 0 0 0 0 0 0 0
Q/1
0 1 2g 3g
2 0 0 0 0 0 0 0 0
Q/2
9
z Q/3
0 '0 0 0 1
ga ga 0 0 0 0
Q/4
w g
Q/5
-
(3.25 )
-
Q/6
9
y
0 0 0 0 0 1 2g 3
S
2
0 0 0 0
Q/7
If'
0 0 0 0 0 0 0 0 1
g
S2
g3
Q/S
Q/9
0 0 0 0 0 0 0 0 0 -1 - 2g _3g
2
Q/10
tl)
Q/1l
Q/12
where ~ = x/t
In deriving the geometric stiffness matrix, it is important
to identify the nonlinear strain components when establishing the
strain-displacement relationship. The total axial strain e may be
xx
obtained from the general strain-displacement relationship for beams
(Eq. 3.7) by .substituting the appropriate curvatures due to large
displacements. The relevent terms, in the case of lateral torsional
displacements, are those which are products of ~ and 00 and their
derivatives. It is believed that inclusion of all the terms in the
analysis will result in furnishing a more complete geometric stiffness
matrix. In the study presented herein, only a limited number of these
terms are adopted merely to constitute a formulation which is compa-
tible with the classical analysis of lateral-torsional buckling. In
the classical approach, for example, it is assumed that the flexural
rigidity about the major axis is very much greater than about the minor
axis. This is equivalent to the assumption that the deflections prior
to buckling are small and can be neglected. However, if the flexural
rigidities about both principal axes are of the same order of magni-
tude, the deflections may be of importance and should be considered. (37)
In the energy formulation for. lateral buckling of a beam
as given by Timoshenko(53) , or in the kinematical model demonstrated
by Bleich(54), the equivalent strain-displacement relationship is
written as
= [- yv xx - p ~ xx] + [ y ~ w xx]
, , ,
(3.26)
-.5.8
Once the nonlinear strain component is defined, the matrix
[QLJ
is determined by establishing the strain-displacement relation-
ship in terms of the parameters [a}, as expressed in Eq. 3.9, thus
0 0 0 .0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 2
2S
2S2
2
S
3
[QLJ
0 0 0 0 3
S 3S
2 3
S
3
~
(3.27)
0 0 2 3
S
0 0 0 0
0 0
2S 3S
2 0 0 0 0
0 0 2S2
3S
3 0 0 0 0
0 0
2
S
3
3
S
4
0 0 0 0
The transformation matrix
[CJ
is evaluated by substituting the coor-
dinates of the nodes into the equations of the displacements (Eq. 3.3).
On substituting the matrices [QLJ, [CJ and [a
o
} into Eq.
3.15 and performing the required integration over the volume of the
beam element, the geometric stiffness is obtained. The resulting
matrix, where the columns and rows corresponding to the torsional
displacements are eliminated since they are all zero, is written as
0
0 0 SYMMETRIC
-36 -33L 0
M 3L 4L2 0 0
[kGJ
oz
= --
30L
0 0 36 -3L 0
0 0 -3L -L
2
0 0
36 3L 0 0 -36 33L 0
3L -L2 0 0 -3L 4L2 0 0
0
0 0 SYMMETRIC
30 2lL 0
V
-3L -2L2 0 0
+60
0 0 -30 3L 0
0 0 9L ~ 0 0
30 9L 0 0 -30 2lL 0
3L ~ 0 0 -3L 2 ~ 0 0
(3.28)
3.5 LARGE TORSIONAL DISPLACEMENTS
A straight shaft may become unstable under the action of a
torque. Similar to the case of Euler buckling in flexure when subjected
to an action of axial compressive force, the bending moment in the shaft
remains zero so long as the axis remains straight. However, as
soon as a deflection occurs, bending moments are introduced about
I
both principal axes at various sections of the shaft. The deformed
configuration in this case is not a plane curve, and the bending
moments vary accordingly as components of the applied axial torque.
A comprehensive treatment of buckling of shafts, based on the conven-
tional approach of establishing the governing differential equations,
is given by ziegler(55) for different loading and support conditions.
Torsional-Lateral Geometric Stiffness Matrix
In the finite element formulation, the displacement of the
beam is represented adequately by a system of third order polynomial
as given by Eq. 3.25. Just prior to buckling, the initial torsional
-60
stress at any point in the beam is given by
Mx P
0'0
=
I
(3.29)
0
where p = radial distance from the center of twist
I = polar moment of inertia of the section about the shear center
o
In deriving the geometric stiffness matrix, it is important
to identify the linear and nonlinear displacement components. For a
beam subjected to a large torsional displacement, Fig. 3.3 shows
schematically the displacement components. The total unit torsional
displacement is written as
(3.30)
where a comma denotes a differentiation.
When these generalized displacements are expressed in terms of the
displacement functions of the beam, the relationship becomes
(3.31)
Note that the nonlinear strain components are associated with the curva-
ture of the beam.
The matrix [QLJ is determined by expressing the nonlinear
strain component [ell in terms of the parameter fa} as expressed in
Eq. 3.9, thus
o 0 o 000
o 0 o 000
o 0 o o 0 2 ~ }
o 0 o o 0 6xf,f}
o 0 o o 0 o
o 0 o
o
o
o 0
o 0
4X/t4 6x
2
1t5
12x2 / t
6
l8x
3
/ t
6
o o
o o
o o
o o
-61
(3.32 )
The displacement transformation matrix [C] is determined by
substituting the coordinates of the nodes into the displacement
function given by Eq. 3.25. At this level, the matrices [QL] and
[C], required for the evaluation of the geometric stiffness matrix
[kG]' are known; on substituting these into Eq. 3.15 and performing
the required integration over the volume of the beam element, the
geometric stiffness matrix of a beam element under pure tension is
obtained as
0
0 0 SYMMETRIC
0 1.0 0
M -1.0 0 0 a
KG
x
(3.33) =-
L
a 0 0 1.0 a
a 0 -1.0 -L/2 0 a
0 -1.0 0 a 0 1.0 a
1.0 L/2 a a -1.0 a 0 0
6 ~
3.6 LINEAR EIGENVALUE PROBLEMS
Solution Technique
Once the elastic and geometric stiffness matrices are deter-
mined for each structural element, the element stiffnesses can be
transformed into a common displacement system and the assembled stiff-
ness matrix is obtained using conventional procedures from the summa-
tion of the element stiffnesses. The resulting equations for the com-
plete system, which take into account the nonlinear effects of large
displacements, are expressed in general form
[F} = [ ~ J [8} + Fl [KGJ
l
{8} + F2 [K
G
J
2
[8}
+ ... + Fn [KGJ
n
[8}
(3.34)
where the matrices [KGJ
i
(i = 1,2, . n) represent various components
of geometric stiffness matrices which are conveniently resolved such
that the generalized forces become the scalar multipliers. Since
the geometric stiffness matrix depends on the nodal displacement vector
[6}, the system of equation given by Eq. 3.34 becomes nonlinear.
In dealing with linear stability problems, it is tacitly
assumed that the buckling deformations are independent of the defor-
mations prior to instability. This leads to the possiblility of
expressing the load vectors of each element as ratios of those parti-
cular loads that introduce instability to the whole structure. Since
the critical load is unknown, a factor A and an arbitrary measure
(scale factor) af the normalized load vector [F} is introduced to
represent the relative magnitude of the applied loads only
! .
! i
-63
fF} = A fr} (3.35)
The factor A is the instability factor or eigenvalue. Since the
geometric stiffness matrix is proportional to the internal forces,
it follows that
(3.36 )
where [KG] is the geometric stiffness matrix for the reference value
of the applied loading.
In performing the numerical analysis, the stiffness matrices
are resolved into two components: the effective elastic stiffness
matrix, which includes the effect of prestress in the element, and the
geometric stiffness matrices whose instability factors are unknown.
The effective elastic stiffness matrix [KE] is composed of the initial
elastic stiffness matrix and those geometric stiffness matrice whose
load parameters A are known,
o
For small displacements f ~ } the effective elastic stiffness matrix
may be taken as a constant. Hence, Eq. 3.37 reduces to
'"
~ ] ~ } + Acr [KG] [il} = [OJ
,
The requirement for a non-trivial solution is
'"
Det I [K
E
] + Acr [KG] I = 0
or
(3.37)
(3.38)
(3.39)
(3.40)
-64
From Eq. 3.40 the eigenvalues A ,which represent the critical loads;
cr
and the associated eigenvectors ~ } , representing the buckling modes,
. cr
are determined using standard matrix methods. The eigenvalue routine
used in this study is based on the Jacobi reduction scheme. The routine
yields the complete eigenvalues and eigenvectors corresponding to the
order of the matrix [KG} Consequently, the higher buckling modes and
the associated critical loads are obtained without introducing extra
provision.
Numerical Examples and Results
At this stage it is appropriate to examine the
numerical accuracy that may be attained in the solution
of actual problems. In general, the adequacy and validity of a
numerical formulation is measured by its performance on problem cases
for which accurate solutions have been derived by alternate means. A
simple and fundamental measure of accuracy may be furnished by referring
to the case of the flexural buckling of the Euler column. In Fig. 3.4
the results of the finite element solutions are shown for centrally
loaded columns with different forms of end conditions, where the
percent error of the solutions are plotted versus the number of
elements in the idealization. The convergence of the numerically
computed critical loads toward the exact solution is remarkably good
as indicated in Fig. 3.4. A similar plot is made in Fig. 3.5 of the
results obtained from the application of the finite differences(56) .
to the governing differential equations for the case of a centrally
loaded pinned-end column, where also a comparison is made with a
finite element solution. It is noted that the error in a1two element
..
-65
idealization, which is regarded as the coarsest possible grid, is less
than 0.8% for the finite element solution.
Of course, the finite element method is not intended as a
procedure for calculation of problems given above which can be treated
rather efficiently by the classical method. However, there are
numerous problems whose solutions are not straightforward when con-
ventional approaches are used especially when there are irregularities
in loading and geometry. Such problems lead to mathematical problems
which are usually intractable. The formulation given here will enable
the solution of a wide variety of problems. In order to demonstrate
the efficacy of the finite element method in the solution of linear
stability problems, several basic problem cases are selected whose
solutions by classical means are not straightforward. Another aspect
taken into consideration in the selection of the demonstrative problem
cases is that each case possesses a distinctive feature from a finite
element standpoint. In the following, several basic examples, which
are encountered in many engineering situations are studied and the
numerical results are presented.
a) Columns with Distributed Axial Loads
Variations in axial loads in columns is a condition encountered
in many engineering situations. A vertical column having significant
self-weight, the decrease in axial load with depth in a pile embedded
in a frictional medium, guyed towers, industrial racks and library
stacks are but a few examples. solutions to such problems by
conventional means are not straightforward especially when there are
-66
irregularities in loading and geometry. For instance, the solution
of the governing differential equation for a pinned-end prismatic
column under a uniformly distributed axial load requires the applica-
tion of Bessel functions (53) For such problems, the finite element
method furnishes solutions in a simple manner.
In formulating a finite element relationship consider a
column having general boundary conditions and carrying some arbitrarily
distributed axial load of intensity q(x) per unit length together with
. an end load P. Both the distributed load and the end load may be
compressive as shown in Fig. 3.6, or either one may be tensile while
the other is compressive. In all cases, however, the conditions for
buckling are represented by the critical combination of the two sets
of loads. The influence of the initial load can readily be intro-
duced by modifying the elastic stiffness matrix ekE] of each element
by adding a scalar mUltiple of the associated geometric stiffness
matrix.
For the case when the end load P is the initial load (pre-
stress load), the modified elastic stiffness matrix is
(3.41)
Or, when the initial load is a distributed axial load, the corres-
ponding modified stiffness matrix becomes
(3.42)
where qi = resultant axial load at element i
consistent load factor of element i
-67
Once the modified matrix for each element is evaluated, the critical
loads are determined by finding the non-trivial solution for the homo-
geneous equations of the complete system given by Eq. 3.40.
To illustrate the advantages of the finite element method,
examples are selected whose closed-form solution are available. Very
few analytical solutions are, found in the literature since the solutions
involve integrals which are difficult, and usually impossible, to
evaluate(37). The examples selected are prismatic columns, having the
four conventional boundary conditions, and loaded by end loads together
with uniformly distributed axial loads q(x) = q. The analytical
o
solutions for the critical combinations of buckling loads, as evaluated
by Dalal(57) are compared to the results obtained through the applica-
tion of the finite element method. The results are compared in graphi-
cal form in Fig. 3.7. A very good agreement is observed for all ranges
of combinations of loads. The results from the finite element solutions
are also given in tabular form in Table 1. For two of the numerical
examples solved, the first buckling modes, under different combinations
of loads, are shown in Fig. 3.8.
Unlike the conventional approaches, which require additional
and usually tedious computatiDnal scheme to determine the higher buck-
ling loads and the associated modes, the finite element approach readily
furnishes these values as part of the original operational scheme.
b) Tapered Columns
Tapered columns of different cross-sectional shapes are used
for structural purposes in a variety of applications. Their use is
-6b
attractive in many situations where the applied load can closely be
specified and a saving in weight is encountered. Analysis of tapered
columns of different shapes and end conditions by classical means is
difficult and has been the subject of a number of investigations
(53,58,59,60)
Such problems, however, can be treated in a simple
manner by the use of the finite element method.
In formulating the finite element relationship, the column
is idealized as an assemblage of discrete elements, where either
piecewise prismatic elements or uniformly tapered elements may be
used. For the stepped representation, the section properties at the
mid-length of the element sufficiently describe the element. However,
such idealization may furnish less accurate results when coarse
discretication is employed and when the member has a high gradient
of taper. For such members, use of tapered elements will yield
results with a better accuracy. The derivation of the stiffness matrix
for beams with a uniform taper in either one or both principal axes
is given in Appendix I.
The assembled stiffness matrix of the complete system is
obtained from the summation of the geometric stiffness matrices, which
are independent of the cross-sectional properties, and the elastic
stiffness matrices of each element. The critical loads are then
determined by finding the non-trivial solution of the homogeneous
equation (Eq. 3.40).
To illustrate the advantages of the method, the critical
loads of tapered columns with one end fixed and free at the other end
~ 69
are computed for uniformly tapered columns. For these columns, a
taper parameter is defined as the ratio of moment of inertia at the
two ends ~ = IT/IB' where IT is the moment of inertia at the free end
(top) and IB represents for the fixed end (bottom). For columns with
~ ~ 1.0, the results obtained by the finite element method are compared
to the theoretical solutions given by Timoshenko(37). Thecomparison
is shown both in graphical and tabular form in Fig. 3.9 and Table 2,
respectively. A very close correlation is observed for all values of
~ The critical loads for ~ ~ 1.0 is shown in Fig. 3.10 and a study
on convergence is demonstrated in Fig. 3.11.
c) Columns on Elastic Foundations
The behavior of axially loaded columns with continuous elastic
support can be considered as the idealized form of a number of related
problems in engineering. Embedded piles receive lateral support
from the surrounding soil, compression flanges of beams and girders
are laterally supported through the web system, and railway tracks or
continuous crane rails subjected to axial loads, such as during tem-
perature changes, also receive lateral elastic support.
A considerable' amount of literature exists regarding the
analysis of beams supported on elastic foundations by seeking solu-
i h i d
'ff . 1 (53,61,62) Th' h
tons to t e govern ng 1 erent1a equat10ns 1S approac
becomes more difficult to solve those problems with variations in
loading, the supporting medium and on the geometry of the member.
In formulating the finite element relationship, a column
supported on a Winkler-type foundation(6l) is assumed, that is, the
-70
elastic foundation can be replaced by a continuous set of springs each
of which can deflect independently. (*) The supporting medium may have
a variation in the foundation modulus k(x) along the column length
and may be capable of developing lateral, rotational and torsional
restraints.
The discrete element stiffness formulation results in a
simple matrix relationship of the form
(3.43)
where a new stiffness matrix [kFJ is introduced to represent the
consistent stiffness matrix of the foundation. The derivation of [kFJ
for a Winkler-type foundation is given in Appendix II. Once the
evaluation of the element stiffness matrices is completed, they are
assembled to obtain the equations of the complete system. The critical
loads are determined by finding the non-trivial solution of the homo-
geneous equation (Eq. 3.40).
To demonstrate the usefulness of the method, the critical
loads of axially loaded pinned-end columns on elastic foundations
(lateral springs) are evaluated for various values of foundation
modu1ii k(x) = k. In Fig. 3.12, the results are compared, in graphical
form, to the analytical solutions obtained by Hetenyi(6l). A very
good agreement of results is observed even when a coarse discretiza-
tion is used (N = 4).
(*)There are also other types of foundation models which have been
suggested by Wieghart(63) , Fi1enko-Bordich (64), V1asov (65) and
Biot (66). The use of such models may also be incorporated by
developing the appropriate consistent stiffness matrices of the
foundation as demonstrated in Appendix II for the case ofa Winkler-
type foundation.
-11
d) Piecewise prismatic Columns
piecewise prismatic columns are occasionally encountered
in special situations in engineering. The main differing feature which
is characteristic of such columns, when compared to the example cases
studied above, is the variation of the direction of the principal axes
of the cross section along the length. During buckling of a piecewise
prismatic column, the deflected configuration does not necessarily
become perpendicular to the axis of the least moment of inertia and
thus will exhibit a non-planar configuration. This feature causes the
governing differential equations to be nonlinear, consequently the
solution by the classical approach becomes difficult. For the parti-
cular case of elastic buckling of a two-component column composed of
identical components, the solution was given by Hsu(67) using classical
methods.
In order to demonstrate the application of the finite element
method to stability problems of piecewise prismatic columns,
consider a two-component column composed of two prismatic columns of
the same cross section. The components are assumed to be rigidly
connected, with one component on top of the other (Fig. 3.13), in such
a manner that their centoidal axes are coincident but the principal
axes are offset by an arbitrary angle a. The shear center of each
column cross section is assumed to coincide with its centroid.
In formulating the finite element relationship, two local
coordinate systems, namely the x-y-z system and the x-y'-z' system (Fig. 3.13)
are used such that y-z and y'-z' coincide with the principal axes
of the cross-section of the lower and upper components of the column,
-72
respectively. The x-axis is identical to both segments and coincides
with the centroidal axis of the column. The stiffness expressions,
when written with reference to the local coordinate system, are
(3.44 )
for the lower component, and
(3.45)
for the upper component. In Eq. 3.45, [F'} and f5'} are the force
and displacement vectors corresponding to the y'-z' axes system.
To obtain the equations of the complete system, a common
displacement system is used by choosing a global system of axes which
coincides with the y-z system of the lower component. Following conven-
tional procedures, the generalized displacements of the upper component
[5'} are transformed to the global system, through the relations
u'
w'
e'
y
y'
e'
z
,
cp
,
OJ
1
o
o
= 0
o
o
o
o o
cosa o
o cosa
sina o
o -sina
o o
o o
o o o o u
-sina o o o w
o sina o o
cosa o o o
o cosa .0, 0
o o 1 0
o 001
OJ
(3.46)
where a is the angle of offset shown in Fig. 3.13. In a similar manner,
the force vector [F'} is determined using the same transformation matrix,
thus
-73
fF'} = [TJ fF} (3.47)
Since the matrix [TJ is orthogonal, Eq. 3.47 may be written as
Finally, the evaluation of the element stiffness matrices is performed
with reference to the global system and the element stiffness matrices
are assembled to form the equation of the complete sY,stem., The critical
loads, P , are obtained from the non-trivial solution of Eq. 3.15.
cr
In order to demonstrate the application of the finite element
solution procedure described above, the critical loads of two-component
columns with arbitrary offset angle ~ are determined for various values
of the stiffness factor ~ = 1
1
/1
2
, The analytical solutions for such
columns are given by Hsu(67) for limited values of the factor ~ .
For the available analytical results, the finite element counterparts
are compared. Figure 3.14 compares the critical loads of columns sup-
ported on spherical pins, and in Fig. 3.15 for columns with fixed end
conditions. Additional results, which cover a wider range in values
of the factor ~ , are given in graphical form (Fig. 3.14 and Fig. 3.15)
and also in tabular form (Table 3).
As a second example, the critical loads of multi-segment
o
columns, where the segments are offset orthogonally ~ = 90 ) in
consecutive order, are determined for different values of the factor
The approximate analytical solutions are given by Fischer
The results are compared in Fig. 3.16 for columns with pin
ends and in Fig. 3.17 for fixed columns, where in both cases good
correlation is observed. The results are also given in Table 4.
-74
e) Pretwisted Columns
A pretwisted column is a structural member that has a natural
twist about the longitudinal axis which may vary in some arbitrary
manner along the length. While pretwisted beams are used in a variety
of applications, as in turbine blade ,and aircraft propellers,
pretwisted columns are also encountered occasionally. More than its
practical justification, the study of pretwisted columns here offers
a good example to demonstrate the efficacy of the finite element
method to beam-column problems. In the buckling of pretwisted columns,
as in the case of piecewise prismatic columns, the column assumes a
non-planar deformed configuration and the resulting deflectiOn is no
longer perpendicular to the axis of least stiffness. This character-
istic makes the equilibrium equations to be nonlinear differential
equations whose solution may be difficult to obtain in a simple manner.
Little information is available on the study of pretwisted
columns. Ziegler(69) investigated the buckling of pretwisted beams
and columns. Later, Zickel(70) developed a theory concerned with the
behavior of pretwisted beams and columns of thin-walled section. More
recently, Fischer(7l) investigated the influence of pretwisting on
the buckling load of a column for different boundary conditions and
moments of inertia. These investigations deal essentially with
approximating the nonlinear differential equations by a set of uncoupled
homogeneous linear differential equations, introducing certain
simplifying assumptions.
In formulating the finite element relationship, the pretwisted
column is idealized as an assemblage of either uniformly pretwisted
-75
beam elements or prismatic beam elements which are piecewise twisted
with respect to one another along the centroidal axis. For the piece-
wise prismatic representation, the section properties and the incre-
mental twist at the mid-length of the element may sufficiently describe
the interval modeled. However, for situattons when coarse discritiza-
tion are used and when the member has a high gradient in twist, the
use of pretwisted elements will furnish results with better accuracy.
The derivation of the stiffness matrix for a uniformly pretwisted element
is given in Appendix III.
Once the stiffness matrices for the elements are determined
with reference to the local coordinates they are subsequently trans-
formed to the global coordinates using Eq. 3.48. The transformed
matrices are finally assembled following the conventional procedure of
summation of the element stiffnesses. The critical loads are determined
by finding the non-trivial solution for the homogeneous equation of
the complete system (Eq. 3.15).
In order to obtain a verification of the finite element
solution, a short test program was carried out to compare the experi-
mental and theoretical strengths of pretwisted columns. It appears
that very little theoretical study and no experimental work is found
in the literature on the elastic buckling of pretwisted columns. The
test program conducted in this study consists of four high-strength
steel wide flange shapes (2-5/8 x 1-1/2 WF). The specimens were
prepared from beams which were originally prismatic by twisting the
beams, to simulate a natural twist along the length, until a permanent
-76
twist of the desired magnitude was attained. In order to induce a
uniform twist along the length, pure torque was applied to the specimens
in an engine lathe as shown in Fig. 3.18, by controlling the rotational
displacements. Four specimens were prepared having natural twists of
o. 0 0 0
o 90 , 180 and 360
The pretwisted columns were tested under a hydraulic testing
machine (Fig. 3.19) following a standard testing procedure for centrally
loaded prismatic columns(72). For each column test, a boundary condi-
tion consisting of a knife edge condition along the web of the cross
section (the minor principal axis) was used at each end of the column.
The test data and the experimental results of the stability tests are
summarized in Table 5. The Table also gives the theoretical critical
loads predicted by the finite element method and estimates of the
experimental critical loads which were derived by extrapolating from
the load and deflection using a Southwell plot method (37) The results
are also shown graphically in Fig. 3.20. It can be seen, in general,
that the critical loads are consistently more than the experimental
loads and there is good correlation between the theoretical and
experimental critical loads. Comparisons between the theoretical
and experimental buckling modes of the tested columns are shown in
Fig. 3.21.
Theoretical predictions of critical loads of pretwisted
columns for various angles of pretwist and different values of the
factor = I1/12 are shown graphically in Fig. 3.22 and also in tabular
form in Table 6. A knife edge condition about the minor principal
axis is imposed at each end of the pretwisted column. It is seen that,
-77
for cross sections with > 1.0, the critical load is always greater
than that of the prismatic column counterpart. The increase in strength
becomes significant for larger values of the factor and for certain
ranges in values of as shown in Fig. 3.22.
It is instructive to note the variation in strength for the
particular case = 1.0 shown in Fig. 3.22. Hypothetically, within the
context of the beam theory, pretwisting has no effect on the buckling
strength of columns whose cross sections have equal moment of inertia
about all centroidal axes = 1.0). However, the variation in strength
shown in Fig. 3.22 is attributed to the directions of the knife edge
condition imposed at the ends. For the particular case when = 90, the
strength of the pretwisted column is equivalent to that of a prismatic
column with pinned-fixed end conditions or to that of crossed pin end
conditions oriented orthogonally, The latter view leads to the classic
problem of the buckling of prismatic columns with crossed pins or
ginglymus joints, that is, the pins that permit rotations in a single
plane at the column ends are oriented at an arbitrary angle to each
other. Ashwell(73) solved the problem of buckling of prismatic columns
with crossed pins for cross sections with = 1.0. Analytical solutions
for cross sections with unequal moment of inertias about the principal
axes become difficult since, in pretwisted columns, the column deflec-
tions during buckling are non-planar.
These complications do not arise when using the,finite
element method, since it is needed only to have an additional trans-
at the node where the arbitrarily oriented pinned end
condition is imposed. At this node, transformation is performed from
-78
the local to the global system using the matrix given by Eq. 3.48.
Applying the finite element method described, the critical loads
were determined for prismatic columns of arbitrary cross sections with
crossed pins. The results are shown graphically in Fig. 3.23. For
the particlar case when ~ = 1.0 the theoretical solution given by Ash-
well(73) is compared to that of the finite element solution. Prismatic
columns fixed at one end and pinned at an arbitrary angle at the other
end are also solved and the numerical results are shown graphically
in Fig. 3.24.
f) Lateral-Torsional Buckling Problems
The importance of lateral-torsional buckling in governing
the strength of thin-walled beams has long been recognized. For
perfectly straight beams subjected to major axis bending, failure in
the lateral-torsional buckling mode may occur at loads considerably
below those necessary to cause failure in the plane of the applied
loads.
The classical procedure for determining the buckling loads
of beams involves the formation and solution of the governing differ-
ential equations. Alternatively, the solution is obtained by employ-
. (5354)
ing energy methods such as the ~ t z procedure ' These methods
are not, however, sufficiently general to deal with many situations
which occur in practice, such as in handling complex loading and support
conditions and irregularities in the beam geometry. The advantages
of the finite element approach are therefore considerable since a
general formulation once established may then be applied to a wide
variety of lateral-torsional problems.
-79
In order to demonstrate the accuracy and convergence charac-
teristics of the finite element formulation (Sec. 3.4), a few examples
are examined whose solutions are available by alternate means. The
convergence characteristics are demonstrated by examining a simp1y-
supported beam subjected by equal end moments M
cr
The results are
shown in Fig. 3.25 where the percent error is plotted versus the number
of elements used in the idealization, using the elastic stiffness
matrices based on cubic polynomial and strain field formulation (Sec.
2.3). For a four element idealization, the percent error is less than
0.1% which may be regarded as a remarkable convergence characteristic.
Note that the convergence characteristics of torsional buckling under
axial loading are identical to the case of lateral-torsional buckling as
shown in Fig. 3.25. To demonstrate the accuracy of the formulation,
a simply-supported beam subjected to a concentrated load at the midspan
of the beam is examined. In Fig. 3.26, the finite element solution
is compared graphically to the analytical solution given by Timoshenko(53)
for a wide range of cross section properties (see also Table 7).
The application of the finite element formulation to complex
problems may be demonstrated sufficiently by examining the experimental
investigations performed by Trahair on the elastic stability of
(74) . (75)
continuous beams and s ~ m p y supported tapered beams In both
cases the loads were applied at the top of the flange. In the finite
element approach, effects of applying a load P at distance a from the
shear center are taken into the formulation simply by adding a moment
term (pa) in the [kG] matrix corresponding to the degree of freedom
~ at the node where the load is applied.
The experimental results of two-span aluminum beams tested
by Trahair(74), which are presented in the form of an interaction
diagram, are compared to the finite element solution in Fig. 3.27.
-80
The agreement with the experimental results is satisfactory. The
theoretical values are obtained by taking the appropriate displacement
and loading conditions during the assembly of the element stiffnesses.
The moment diagram prior to buckling is taken as the loading condition
which consists of terms as mUltiples of the unknown load parameter P
cr
applied at one span and the known load P at the other span. The effect
of the known load P is added to the elastic stiffness matrix ekE]'
For different values of P the correspondin, PeR is determined by solving
the eigenvalue problem.
In the case of the buckling of tapered beams, the only
differing feature is the elastic stiffness matrix ekE] which varies
for each element. As in the case of tapered columns a stepped repre-
sentation may sufficiently describe each element by taking the section
properties at mid-length of each element. However, when beams with
high taper gradients are encountered, the elastic stiffness matrix
of uniformly tapered elements may be used (Appendix II). The experimental
results of simply supported aluminum beams (Fig. 3.28) tested by Trahair
(75) are compared to finite element solution in Fig. 3.29. Also, the
theoretical values given by Trahair(75) , where the differential equations
are solved by numerical methods, are compared in Fig. 3.30. In both
cases, the agreement is satisfactory. Using the finite element method,
additional numerical results are given in Fig. 3.31 for different
combinations of taper parameters.
-81
g) Space Frames
The study of the overall stability of a space frame to obtain
the actual buckling condition of the entire system has been a subject
of numerous investigations. A number of methods are now available for
obtaining exact or approximate solutions. A detailed account of the
work in this field is found in Ref. 54 and 76. The application of
the finite element method in solving problems of space frame buckling
is now presented.
In formulating the finite element relationship, the internal
forces in each member prior to buckling are first evaluated. The
frame is then idealized as an assembly of discrete beam elements. The
members which are not subjected to axial forces during the prebuckling
state are sufficiently represented by single elements since the geo-
metric stiffness matrices for such elements are null matrices. In
such situations, a one-element representation does not introduce addi-
tional numerical errors since the displacement expression of the
element, which is a third order polynomial, describes consistently
the assumed deflection of a beam with constant shear.
Once the element stiffness matrices are determined in terms
of the local coordinate systems, a displacement transformation is
performed for each element involving the direction cosines when rela-
ting the local and global systems. The equilibrium equations for an
element when expressed in terms of the global coordinate system is
(3.49)
where
[kEJ
g
[TJT [kEJ
t
[TJ
[kGJ
g
= [TJT [kGJ
t
[TJ
[F}g = [TJT [F}t
The subscript g and t represent the global and local reference axes
systems, respectively.
Applying the procedure described above, numerical results
-82
of a few, sample problems are now presented. In order to demonstrate the
convergence characteristics of the method, knee-bent frames subjected
by axial forces and having two different boundary conditions are
considered. The results are compared to analytical solutions in Fig.
3.32 where the percent errors are plotted versus the number of elements
used in discretizing the columns. The accuracy of the method is
demonstrated by comparing the results to the analytical solutions
of the four standard cases of single-story portal frames. The
results are shown graphically in Fig. 3.33.
The application of the method to more complex problems is
illustrated by solving the space frame used in the investigation by
Morino(77) based on the determinantal approach which makes use of the
concept that the determinant of the overall stiffness matrix becomes
zero at the critical load. The space frame is shown in Fig. 3.34
and is subjected to vertical loads p at each joint. All members are
made of the WlOx49 shape and the columns are oriented as indicated
in the Figure. Also shown is the idealization of the frame for the
-83
finite element analysis. The critical loads and the associated buckling
modes are obtained in .one operation as in the previous examples and
do not require iteration. In Table 8 the results from this analysis
are compared to the results in Ref. 77.
h) Buckling of Shaft under Torsion
In order to demonstrate the use of the geometric stiffness
matrix derived in Section 3.5, which corresponds to large torsional
displacements, the buckling strength of shafts under pure torsion
is now presented. For the purpose of simplicity, two shafts are con-
sidered having simple boundary conditions: pinned and fixed in flexural
rotations at both ends. In both cases concentrated torques are applied
at the ends and the shafts are assumed to have equal moment of inertias
about all axes. The buckling analysis of shafts having unequal moments
of inertia about the principal axes is not straightforward. For
such members, just prior to buckling, it is possible that the member
is excessively deformed in torsion, thus, the stiffness of the beam
at buckling is better represented by a pretwisted beam. However, the
pretwist of the beam at the instant of buckling is unknown since
,
the buckling load is yet to be determined. In such situatiDns, there-
fore, it is 'suggested to employ iterative schemes in solving the
problem.
In order to study the convergence characteristics of shafts
having equal moments of inertia about all axes, the shaft is discretized
into different numbers of elements. For a specified number of elements,
the assembled stiffness matrix of the complete system is obtained in
84
the usual manner from the summation of the individual matrices. The
critical loads are determined by finding the non-trivial solution of .
the homogeneous equations (Eq. 3.40). The results are given in Table
9 and are compared to the analytical solution. For both shafts the
eigenvalues corresponding to the first buckling modes are shown in
Fig. 3.35.
, .
4. NONLINEAR ANALYSIS OF BEAM-COLUMNS
4. 1 INTRODUCTION
Nonlinear behavior in-: structural systems arises from two
distinct sources, namely, material nonlinearities and geometric non-
linearities, the former being reflected in nonlinearities in the
constitutive equations and the latter in the nonlinear terms in the
kinematical equations of large displacements.
-85
The application of nonlinear theory to conduct the conven-
tional closed form analysis of structural responses leads to mathe-
matical problems which are usually intractable. In order to obtain
quantitative solutions, it is natural, therefore, to resort to numerical
methods. In general, numerical methods of structural analysis may be
described under two categories. In the first category are the methods
that lead to numerical solution of the governing algebraic or differ-
ential equations, by approximating the mathematical functions, which
are then solved by either finite differences or ,by numerical integra-
tion. In the second category is the finite element method which is
based on the concept of piecewise approximating continuous fields.
The adapt ion of the finite element method has been accelerated in recent
years and is employed extensively in a wide range of nonlinear problems
mainly due to its simplicity and generality. A comprehensive treatment
on the theoretical foundations of the method to nonlinear problems
is given in Ref. 16.
The problem of material nonlinearity arises from the non-
linearity of the constitutive equations as a consequence of describing
the laws of material behavior under multiaxial stress states. The
relationships do not involve only the current state of stresses as
-86
in the case of elastic behavior; rather, they depend as well upon the
past histories of these components. Despite these complexities,
substantial progress has been made in developing the finite element
method based on the principles of the general theory of plasticity.
Useful reviews of accomplishments in finite element inelastic analysis
are found in Refs. 16 and 17 and also in Refs. 78 to 80.
The problem of geometric nonlinearity may be considered as
having several levels of nonlinearity. For a more realistic evaluation
of actual behavior, additional measures of nonlinearities must be
considered. Such measures, however, will result in increasing the
complexity of the formulation and of the solution.
The problems that may be considered at the lowest level in
this hierarchy of nonlinearity comprise those that can be transformed
into linear eigenvalue problems. This is the classical stability
problem, such as in the buckling of perfectly straight columns, where
satisfactory predictions can be made for critical loads (eigenvalues)
and the corresponding buckling modes (eigenvectors) by solving a typical
eigenvalue problem. These problems involve a bifurcation of equili-
brium, the point of bifurcation occurring at the critical load, which
is characterized by the existence of a fundamental state of equilibrium.
A detailed treatment of such problems is given in Chapter 3.
The problems dealing with predicting post-buckling behavior
or those characterized by initial imperfections extend to the next
-87
level in the hierarchy of geometric nonlinearity. The investigation
of post-buckling behavior requires higher order approximations, which,
unlike the linear eigenvalue problem, depend on the type of singularity
occurring at the critical load. In a similar manner, the presence of
geometric imperfections will introduce nonlinearity, but of a different
nature when compared to the linear analysis, since the load-deflection
relationships are no longer based on the initial geometry. The critical
load for an imperfect system is characterized by the load for which
the deflections increase indefinitely. Further higher levels in the
hierarchy of nonlinearity include problems with large rotations and
axial strains as in snap-through type instability problems arising in
arches.
4.2 SOLUTION TECHNIQUES
The solution techniques for nonlinear problems are funda-
mentally the same despite the sources of nonlinearities. A variety
of solution procedures utilizing the finite element concept and its
applications have been employed extensively in recent years. Basically,
most of the numerical procedures fall into two broad classes: namely,
incremental and iterative methods. The incremental methods do not
necessarily satisfy equilibrium conditions, whereas the iterative
methods tend to stay on the true equilibrium path at all steps of the
computation.
The procedures are described by considering the nonlinear
equilibrium equations of the complete system
[KJ fa} = fF}
. .
(4.1)
The nonlinearity occurs in the stiffness matrix [KJ which itself is
a function of the load [F} and displacement fa}. In the following,
the fundamental principles in the two solution techniques are pre-
sented.
Incremental Methods
-88
The application of the piecewise linear incremental procedure
to the finite element analysis of nonlinear structural behavior was
first proposed by Turner et al(81). In this method nonlinear behavior
is determined by solving a sequence of linear problems. The load is
applied as a sequence of sufficiently small increments so that during
the application of each increment the structure is assumed to respond
linearly. consequently, the equations become linear. For each load
increment, the corresponding increment of displacement is computed;
it is accumulated to give the total displacement at any stage of
loading. A subsequent increment of load is applied and the incremental
process is continued until the desired number of load increments has
been completed. The solution procedure takes the following form,
where [KJ =
linear stiffness matrix
[KIJ
=
incremental stiffness matrix at load step (i-I)
fa}
= incremental displacement
[F}
incremental load
Essentially, the incremental procedure solves a sequence of linear
problems where the stiffness properties are recomputed based on the
(4.2)
current geometry prior to each load increment. This process is basically
the Euler-Cauchy integration scheme applied to Eq. 4.1 with the load
[F} playing the role of the dependent parameter. The method is
schematically indicated in Fig. 4.1.
One of the principal advantages of the incremental method is
its complete generality and simplicity in application to many types
of problems with nonlinear behavior. In addition, it provides a
relatively complete description of the load-displacement behavior.
Nevertheless, the method has the disadvantage -that a real estimate
of the solution accuracy is not possible since, in general, equilibrium
is not satisfied at any given load level. In some situations, the
incremental method may lead to a deviation from the true load-displace-
ment relationship especially in the neighborhood of instability
conditions.
In order to reduce the deviation from the true behavior,
an effect which is prominent with Euler type integration schemes,
more accurate schemes such as the Runga-Kutta method may be used.
The addition of a corrective term to the incremental method, for instance,
which requires only insignificant computational effort, has been
demonstrated by Haisler et a 1 ~ 8 2 to increase the accuracy considerably.
This procedure includes a load vector representing the out-of-balance
force [F
R
} determined from equilibrium considerations. Consequently,
the self-correcting incremental procedure becomes
(4.3)
-90
Iterative Methods
The iterative approach is one of the oldest procedures
used by many investigators for solving systems of nonlinear algebraic
equations. In this procedure a sequence of linearized equations is
solved to obtain an improved solution by starting with an initial
estimate to the displacement solutiion. This improved solution is back ..
substituted into the equations and the process is repeated until an
acceptable convergence measured by a prescribed tolerance is obtained.
For the nonlinear equilibrium equations given by Eq. 4.1 the iterative
procedure consists of successive corrections to the solution until
equilibrium condition under the total load [F} is satisfied. The
success of iterative methods, in general, depend
estimates of the displacements.
To obtain rapid convergence for problems exhibiting high
nonlinearities, many investigators have utilized the
iterative approach. This method is accurate and converges rapidly,
whenever a realistic initial estimate of the solution is made, and is
considered as one of the most re1iab1e(82). Based upon an initial
estimate of the displacement f5}. at a given load [F} in Eq. 4.1,
. 1
and using a first-order Taylor series expansion corresponding to
(0 + oJ. a linear incremental relationship is established
1
An increment to the displacements [oJ. is computed during the i th
. 1
cycle of iteration; it is added to the approximate displacement
(4.4)
[o}. to obtain a more nearly correct (i + l)th approximate solution, thus
1
(4.5)
This new solution fa}i+l is then substituted into Eq. 4.4 to obtain
a further correction by utilizing the tangent stiffness at the end of
the previous iterative step. This process is repeated until the
increments of displacements fa}. become sufficiently null. The method

is described schematically in Fig. 4.2(a).
The Newton-Raphson method requires updating and subsequently
inverting the stiffness matrix at each cycle. This process may become
lengthy in particular for larger systems of equations. In such
situations, it may be advantageous to employ the modified Newton-Raphson
procedure (45,82) . In this method the stiffness matrix is held constant
for several iterations or load increments after which the stiffness
matrix is updated based on the current displacements. Obviously,
the procedure necessitates a greater number of iterations; however,
itguarantees a substantial saving of computations as it does not
require an inversion of the matrix at each cycle. A schematic repre-
sentation of the method is shown in Fig. 4.2(b).
4.3 GEOMETRIC NONLINEARITY
This section deals with the general instability problems of
beam-columns in which the displacements are large but the engineering
strains remain 'small'. Geometric nonlinearities enter the finite
element fo'nnulation as a result of nonlinear strain-displacement
relationships, which consist of strain products of the same order
of magnitude as the engineering strains, and also from the effect of
large displacements on the equilibruim equations.
-92
A formulation that takes into account geometric nonlinearities
can be used to study the response of imperfect structures as well
as the post-critical behavior of perfect structures. In general,
actual structural systems are never built precisely as planned and thus
inevitably contain small imperfections associated with geometnic ~ r o r s
These imperfections can change the response of the system.
The nonlinear response of a structure is obtained by utilizing
linearized incremental methods (44) . The resulting incremental equili-
brium equation may be written in the form
(4.6)
where [RE] is the conventional linear stiffness matrix and [KG] is
the geometric stiffness matrix. The subscript (i-I) indicates the
stiffness matrices are evaluated for the state of displacement at the
beginning of the increment. At each load increment the local dis-
placement vector, the overall stiffness matrix and load vector are
related to those in the global system by the transformations
f 5
g
} = [T]T f 5
t
}
[Kg] = [T]T [KtJ [TJ
f F
g
} = [TJT fFt}
(4.7)
(4.8)
(4.9)
The subscripts t and g represent the local and global systems, repsec-
tively. In incremental formulations, the direction cosines in the
transformation matrix [T] become also functions of the current displace-
ment state in addition to the initial geometry. Thus, the transformation
matrix [TJ is modified at each load increment. The load [F}, which
is the scalar multiplier of the geometric stiffness matrix [KG]' is
also evaluated at the beginning of each load increment.
In the iterative method, a load increment [F} is applied
to the structure and the resulting displacements are used to revise
the new configuration of the structure. At each cycle of iteration,
the new geometry is used to recompute the stiffness matrix and the
-93
load vector by using a linear analyses. The solution procedure takes
the following form
(4.10)
where
4.4 MATERIAL NONLINEARITY
The stiffness matrix for problems of material nonlinearity
is computed from the relationship given by Eq. 2.16
[k] = [TJT [H]-l [T] (2.16)
where [H] = J [PJ [D ] [PJ dV
V ep
The matrix [D ], which describes the material behavior under mu1tiaxia1
ep
stress, is now a variable and depends not only on the state of stresses
but also upon the history of loading, The constitutive relationship
foran inelastic material can be expressed in terms of finite increments
(4.11)
The matrix [D ] is modified by updating the components of the matrix
ep
-94
according to the deformation laws of plasticity. For metallic structures,
in general, idealization of the material by a Prandtl-Ruess material
obeying the von Mises yield criterion is universally used.
Two general methods, which are based on essentially different
concepts, have been developed for inelastic analysis of solid bodies.
These are the method of "initial strain" and the method of "tangent
modulus". (78)
The initial strain method is based on the idea of modifying
the equations of equilibrium so that completely elastic behavior may
be assumed. This approach introduces modifications to compensate for
the fact that the inelastic strains do not cause any change in the
stresses. Details on the development of the method are found in Ref.
80.
The tangent modulus method is based on the linearity of the
incremental laws of plasticity and approaches the problem in a piece-
wise linear fashion. As the load is applied in increments, a new set
of coefficients is obtained for the equilibrium equations. This
approach is used now more extensively among many investigators due to
its consistency with the classical methods of plasticity analysis and
also due to its computational efficiency. In Ref. 79 the use of this
method is described in detail.
In this study further investigation on the material nonlinearity
aspect is not made. The problem of material nonlinearity in beam-
columns is identical to those problems of solid mechanics or plates
and shells. Substantial information is available in the literature
on the application of finite elements to problems of material non-
Ii i
(16,17,78,79,80)
near ty.
4.5 SAMpLE PROBLEMS AND RESULTS
-95
The first numerical example to illustrate the validity of the
procedure described is that of a cantilever column with an initial
out-of-straightness 8
0
at the free end. The out-of-straightness along
the length of the. column is assumed to be expressed by
= 6 (1 - cosTI x/L)
U
x
0
(4.12)
For this column the analytical solution for the complete load-dis place-
ment relationship is given in Ref. 53. In the application of the
finite method, the direct incremental procedure is used in evaluating
the load-displacement relationship. The effect of the number of load
increments in the final solution is shown in Fig. 4.3 when using the
direct incremental procedure. It is observed that the results agree
closely to the analytical solution as the load increments become small.
In Fig. 4.4 a set of load-displacement relationships is shown
for the same column for different values of 6 ranging from L/500 to
o
L/5000. The results indicate fairly close correlation to the analytical
solutions. In Fig. 4.5 the finite element and analytical solutions
are compared when the column is bent excessively to large displacements
having the order of magnitude of the length of the column.
-96
5. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE DEVELOPMENTS
The principal objective of this study is to develop a finite
element formulation for the analysis of beam-columns and to demonstrate
the applicability of the method to general beam-column problems as a
practical tool. In the formulation of the finite element model, the
beam column is regarded as a one-dimensional body.
The contributions achieved in this dissertation are:
-A one-dimensional finite element model is developed to analyze
general linear static beam-column problems.
-A systematic procedure is presented to evaluate geometric
stiffness matrices for beam-columns which are required to
perform a finite element analysis of stability problems.
-The geometric stiffness matrices are derived Which correspond
to large lateral and torsional displacements.
-The advantages of the finite element method are demonstrated
in the solution of a few stability problems, such as the
buckling of pretwisted columns and the lateral buckling of
tapered beams, the analytical solutions of which are not
yet available.
The study has served in demonstrating the use of the finite
element method in conducting linear static, linear stability and non-
linear analyses of beam-columns. Furthermore, the advantages of the
method are demonstrated in solving a wide variety of problems having
irregularities in geometry, loading and support conditions which defy
adequate treatment by the classical means.
-97
In the linear static analysis, a one-dimensional finite element
model is developed by making use of variational principles. The evalua-
tion of the element properties is performed by developing a formula-
tion based on a functional constituting of two independent fields:
a polynomial approximation of the strain field in the domain, and
displacements at the boundary. The beam element has two nodes with
seven degrees-of-freedom at each node: three linear displacements,
three rotations about the cartesian system, and a warping displace-
ment. The formulation has an advantage when compared to previous
work in that the required manipulation to evaluate the element pro-
perties is simple especially when additional kinematical assumptions
are introduced such as shearing deformations due to bending and warp-
ing. Reviews of the derivation process disclose a complete sequence.
of numerical integration and matrix operations which can be operated
in a systematic manner.
In the linear stability analysis, a systematic procedure
is developed to evaluate the so-called geometric stiffness matrices
for beam-columns by making use of the nonlinear strain-displacement
relationships. The model developed is used in deriving the geometric
stiffness matrices for beam-columns when the displacements are large
in axial and transverse directions and also in twist. The use of these
matrices is demonstrated in solving a variety of stability problems
through a direct eigenvalue determination. The examples are selected
-98
such that each has a differing feature from the finite element stand-
point. The types of examples include: columns with distributed axial
loads, tapered columns, columns on elastic foundations, pretwisted
columns, space frames and lateral-stability problems. In all cases,
a remarkable convergence is observed and satisfactory accuracy is
obtained by using relatively few elements. Moreover, the required
computational effort is found to be minimal.
In nonlinear analysis of beam columns, the use of geometric
stiffness matrices is found to play an important role in determining
nonlinear responses of structures. The use of the direct incremental
method is demonstrated on a few numerical examples whose analytical
solutions are available. For beam-columns with geometric nonlinearity,
it is found that the direct incremental method furnishes satisfactory
results, especially when the load increments are small. In addition,
the method is straightforward and requires very little computational
effort.
While the applications presented in this study are very simple
the finite element method enables the study of complex and practical
problems. The method derived here should find applications in linear
static and linear stability analyses of novel structures. The stability
analysis, in .particular, should find immediate application to ordinary
structural problems such as space frames. Other applications of the
method include performing optimization studies on the strength of
structural members such as tapered columns and piecewise prismatic
members.
-99
For a more general application of the method, the formulation
should be extended to enable analysis of thin-walled closed beams,
multi-cellular beams, curved beams and the like.
Finally, extension of the method to include material non-
linearity will encourage studies on a wider variety of practical
problems. Among the many possible applications, the validity of some
of the universally accepted assumptions used in inelastic analysis
of beam-columns can be investigated. For instance, the universally
accepted assumption that the shear response is always elastic for an
inelastic beam-column is but one of the many possible items for
investigation. Other problems of interest are related to the i n v s t i ~
gation of structures under repeated loading, cyclic loading and shake-
down.
6 APPEND ICES
APPENDIX I
ElASTIC STIFFNESS MATRIX FOR THE BEAM ELEMENT
-lOO
The linear elastic stiffness matrix for a straight beam e1elrent
of uniform cross section is now presented. A cartesian reference axes
system is used which coincides with the centroidal-principa1 axes of the
cross section. The element is represented by two nodes with seven
degrees of freedom at each node; three linear displacements, three
rotations about the reference axes and a set of warping displacements.
The stiffness matrix is derived based on the formulation given in Section
2.3 which take into account the shearing deformations due to flexural
and torsional-warping loadings. In the stiffness matrix given below the
following expressions are used:
where
0
e
qiy =
=
3 + (1 + 1.5
=
3 - (1 + 1.5
l2EJ
z
GA JP
sy
12EJ
qi = y
z GA L2
sz
2
cpz)
cpz)
+ ~ r +
4 ,.5
(1 + 1.5 cpz)2]
2
2
+ ~ [
(1+1.5cp )
J
z
4 5 3
A A = Effective shear area for V and V , respectively.
sy' sz y z
(AI. I)
(A1.2)
(A1.3)
(Al.4 )
EA
L
12EJ y
L3 (J+cPzl
6EJy (4+cPzl EJy
L
2
C1+cPzl L (l+cPzl
12 EJz
L3 (J+cPyl
6EJz
(4+cPylEJz
-
@(I+cPyl LCI+cPyl
12E1w'
L3 (I+1.5cPzJ2
6E1wA
L2(1+1.5cPzl2
[k]
EA
-T
12EJy 6EJy
- -
L
3
(J+cPzl L2 (J+cPzl
6EJy (2- cPzlEJy
L2 (J+cPzl L( l+cPz l
12EJz 6EJz
-
L3 (I+cPyl L2(I+cPyl
6EJ z
(2-cPyl EJz
L2C1+cPyl L( l+cPyl
12E1w'
-

I
6E1wA

E1w'u
L (1+1.5cPzl
2
EA
L
6E1wA
(1+1.5 cP
z
)2
El w8
L(J+1.5<P.zl
2
i
r
l
:--...... ........'8yl y
8" ""K
vI L-----
Z
r
2
,
("2-e,-e Z2
V2 u2 ........
;.> 8X2
f X
SYMMETRIC
12EJy
L3(J+cPz l
6EJy (4+cP
z
lEJy
L2 (1+ cPzl L(I+cPzl
12EJz
L3(1+cPyl
6EJz (4+cPylEJz
L2 (l+cPyl L (l+cPyl
12E1w'
t: (1+1.5 cPz)2
6E1wA E1w'u
UI+I.5cPzl
2
I
I
t-'
a
t-'
-102
APPENDIX II
FLEXURAL STIFFNESS MATRIX OF TAPERED BEAMS
Consider a beam element with the smaller end denoted as end
i and the larger as end j. The element may have uniform taper in
either one or both principal axes. Gere and Carter(58) defined the
depth, d , and the moment of inertia, I , at a distance x from end i
x x
of the member as
d =
d. [1 + ~ -
1) fJ
x ~ d
i
(A2.1)
d.
1) it
Ix
=1 [1+(--1-
i d.
~
(A2 .2)
in which ~ refers to the shape factor that depends upon the shape of
the cross section and may be obtained from
(A2.3)
The values of ~ for various types of cross sections are found in Ref. 59.
Once all values of ~ are known for various sectional proper-
ties, such as cross-sectional area and moments of interia, the flexi-
bi1ity matrix may be determined by the matrix integration
where [ T ] = translational matrix from j to x
. xj
[U ] = the basic sectional property matrix
x
(A2.4)
which are defined as
I 0 0 0 0
0 I 0 0 0
[T .J
= 0 0 I 0 0
xJ
0 0 -(L-x) I 0
0 (L-x) 0 0 I
and
1
0 0 0 0
U-
xx
0
1
0 0 0
GA
yx:
[UxJ
0 0
I
0 0
GA
zx
0 0 0
I
0
EI
yx
0 0 0 0
I
Elzx
The flexibility matrix is determined by performing the
integration of Eq. A2.4 over the length of the beam(59),

0 0 0 0
EAxi
0
+ )
0 0
f32IP
3EI i GA.
2Elzi
z
[f j j J
0 0
(f3
3L3
+ )

0 =
3El
Yi
GA
zi 2Elyi
0 0

f3IL
0
- 2Elyi EI .

0
f32L2
0 0

--
2Elzi
EI .

-103
(A2.5)
-104
The values of ~ 1 ~ 2 and ~ 3 are given in Ref. 59 for various types of
cross sections.
Finally, the stiffness submatrix is determined from
(A2.6)
the other stiffness matrices are found through the use of equilibrium
equations(19) and symmetry considerations as follows.
and
T
= [k
ij
] (A2.7)
Or, the element stiffness matrix may be determined in a sfmple operation
from
(A2.8)
where [N] = [T, I], [I] = identity matrix.
-105
APPENDIX III
CONSISTENT STIFFNESS MATRIX FOR BEAMS ON ELASTIC FOUNDATION
Consider a beam element on a Winkler-type elastic foundation
with spring constant (foundation modulus) k over the length of the
element. For a consistent formulation, the deformation of the founda-
tion must be identical to that of the beam it supports. Accordingly,
the displacement field for the foundation is expressed as a polynomial
of the third order which is written in matrix notation, as
fU} = [PJ fed
y (A3.1)
The displacement, y, can be either lateral or torsional, depending
on whether the foundation is represented by lateral or torsional springs.
The vector ~ } consists of the coefficients which are to be determined
in terms of the nodal displacements fa} from
(A3.2)
Using the principle of vietual work, the total external work
is written as
W = f5}T fF} (A3.3)
ext
where fF} is the vector of nodal forces. The total internal work is
= J [u}T [kJ [u} dx
L
(A3.4)
-101)
Assuming a homogeneous foundation modulus over the element length,
the total internal work is expressed in terms of the nodal displacements
(A3.5)
Comparison of Eq. A3.4 and Eq. A3.5 yields
(A3. 6)
substituting the values of [C-
1
] and [PJ and integrating over the length
L furnishes the consistent stiffness matrix as follows
156
22L SYMMETRIC
(A3.7)
54 13L 156
-13L -3L2 -22L 4L2
I
-lD7
...
APPENDIX IV
FLEXURAL STIFFNESS MATRIX OF PRETWlSTED BEAMS
Consider a uniform beam element of length L with natural
twist about the centroidal axis which varres linearly with the axial
coordinate x. The element has a total angle of Let the
end at x = 0 be denoted as end i and the other as end j. At distance
(A4.l)
By assigning = 0 and = , the y and z axes become aligned with
J 0
the principal axes of the cross section at x = O. The flexural stiff-
nesses about the minor and major principal axes of the beam cross sec-
tion are Ell' and El
2
, respectively, and are independent of x.
A direct evaluation of the stiffness matrix for a pretwisted
beam, in the manner performed for the case of a prismatic beam (Section
2.3) is difficult and is not attempted here. Nevertheless, the
determination is simplified by evaluating the flexibility matrix first,
such as by utilizing previous investigations which have
ing the governing differential equations for pretwisted beams.
The equations governing displacements of pretwisted beams
due to terminal loads are found in Refs. 83 and 84. For matrix appli-
cations, these governing equations are separated and integrated(85)
to estab lish explicit displacement functions describing the transverse
displacements and rotations of the beam about the y and z axes. The
-lOB:
flexibility matrix is found by imposing the displacement boundary
conditions at end j where unit terminal loads are assigned, thus,
SYMMETRIC
(A4.2)
where
1
fll = 6 (Ell + EI
2
) - (EI
2
- Ell) (a - sina)/of
1
f2l = 2 (EI
2
- Ell) [2(1 - cosa)/af - l/a]
1
f22 = 6 (Ell + EI
2
) + (EI
2
- Ell) (a - sina)/af
1
f32 = - 41 [Ell + EI2 + 2(EI
2
- Ell) (1 - cosa)/of]
1
f33 = 2L2 [Ell + EI2 + (EI
2
- Ell) (sina)/aJ
1
f4l = 41 [Ell + EI2 - 2(EI
2
- Ell) (1 - cosa)/ofJ
1
f42 = - 21 (EI
2
- Ell) (a - sina)/of
1
f43 = 2Lf (EI
2
- Ell) (1 - cosa)/a
1
f44 = 212 [Ell + EI2 - (EI
2
- Ell) (sina)/aJ
Once the flexibility matrix is determined, the evaluation of
the stiffness matrix is accomplished following the same procedure of
matrix operations tl-escrThed in Appendix II, thus
[kJ
T -1
= [N] [fjj] [N]
(A4.3)
where
[N] = [T, IJ, [IJ = identity matrix.
A
[C]
D,[D]
E
F
fF} ,fF'} ,(F"}
[F}
G
[H]
I
o
I
ill
J ,J
Y z
[J]
[L]
M
M
cr
-109
7. NOTATIONS
= cross sectional area
displacement transformation matriK
= generalized Hookean constant
Young's modulus of elasticity
= body force
= load vector of assembled system
incremental load vector
= shearing modulus
= stiffness matriK in terms of generalized displacement
amplitude
= moment of inertia
moment of inertia about the shear center
= warping moment of inertia
= moment of inertia about the y and z axes, respectively
= inertia matrix
= St. Venant torsional constant
= assembled (master) stiffness matriK
= linear stiffness matrix
geometric stiffness matrix
incremental stiffness matrix
= length of beam element, column length
displacement interpolating function matriK
= bending moment
critical moment
tota 1 torque
N
[N]
P
PeR
P
E
[P]
[Q], [QLJ
[R]
8,8 ,8
(J u
T
u
v ,V
Y z
w
a
a
d
e
f f}
h
= critical torque
= bimoment
number of elements
= shape function
axial thrust
= axial critical load
= Euler load
= polynomial functions of coordinates
polynomial functions of generalized displacement
amplitudes
suface coordinates of tractions
= surface area and portions of boundaries
= traction on boundary S
(J
= transformation matrices
= strain energy
= basic sectional property matrix
= volume
= shear force in y and z directions, respectively
= work or energy
=
distance from rotation axis to shear center
= point of load application from shear center
=
depth of beam
= Lagrangian strain tensor
=
element load vector
=
frame height
-110
-111
k =
foundation modulus, cross sectional parameter
[k]
element stiffness matrix
n = directional vector
p
=
traction on boundary S
qo,q =
intensity of distributed load
r distance from axis of twist
u displacement in x-direction
u =
displacement on boundary S
u
'"
u =
interelement boundary displacement
v =
displacement in y-direction
w = displacement in z-direction
x,y,z =
reference axes
r =
cross sectional parameter
=
displacement vector of complete system
11. =
cross sectional property

= cross sectional property
'
=
cross sectional property
0
= cross section property
cc
angle between pins, angle of pretwist
""'
r cc} ,r cc} =
generalized displacement amplitudes
[cc] =
shear coefficient matrix
(3
critical load parameter
[ (3}
= strain coefficients
y
average shearing strain
5
[5},(8"}
f 5}
[ e}
[ eo}
[ ell
fJ
p
cr
x
ill
-112
variational convention
element dis placement "'
= incremental displacement
strain fie Id
= linear strain
= nonlinear strain
= ratio of moment of inertia
= rotational displacements
= cross sectional parameter
= instability factor or eigenvalue, frame parameter
= nondimensionalized length
= functional or objective function
= radial distance of component plate from shear center
= stress
= stress field
= shear strain function, angle of twist
= shear strain function
= cross sectional property
= warping displacement
-113
8. TABLES AND FIGURES
-114
TABLE 1 CRITICAL LOADS OF COLUMNS WITH DISTRIBUTED AXIAL LOADS(*)
A. When the end loads (P) are prescribed:
P
QCR/PE
P
E
Case 1 Case 2 Case 3 Case 4
-2.0 5.1875 6.6611 8.5779 5.2296
-1.0 3.5931 4.7505 6.0030 3.6146
-0.5 2.7565 3.7088 4.6199 2.7694
0.0 1.8846 2.5849 3.1633 1.8909
0.2 1.5242 2.1072 2.5584 1. 5285
0.5 0.9692 1.3553 1.6257 0.9712
0.8 0.3949 0.5561 0.6613 0.3956
1.0 0.0000 0.0000 0.0000 0.0000
B. When the distributed loads (Q = are prescribed:
.JL
PCR/P
E
P
E
Case 1 Case 2 Case 3 Case 4
-2.0 1. 9351 1. 6471 1.5786 1. 9391
-3.0 2.3564 1. 9405 1.8542 2.3653
-4.0 2.7502 2.2171 2.1214 2.7664
(*) Refer to Fig. 3.7 for identification of cases
-H5
TABLE 2 CRITICAL LOADS OF UNIFORMLY TAPERED COLUMNS(*)
Boundary condition: Free at top
Fixed at bottom
Values of ~
IT
Timoshenko
Finite Element Solution
(N = No. of Elements)
11 =-
Solution(53) IB
N = 2
N = 4 N = 10
0.1 1.202 1.0906 1.1740 1.1985
0.2 1.505 1.1409 1.4818 1.5014
0.3 1.710 1.6309 1.6906 1. 7071
0.4 1.870 0.8036 1.8535 1.8672
0.5 2.002 1.9482 1.9890 2.0002
0.6 2.116 2.0740 2.1060 2.1149
0.7 2.217 2.1861 2.2097 2.2163
0.8 2.308 2.2879 2.3032 2.3075
0.9 2.391 2.3816 2.3886 2.3907
1.0 2.467 2.4687 2.4675 2.4674
2.0 3.1328 3.0502 3.0286
4.0 4.0098 3.7658 3.7025
10.0 5.6887 4.9892 4.8071
20.0 7.5899 6.2088 5.8447
40.0 10.3745 7.7947 7.1029
100.0 16.3334 10.7359 9.2046
(*) Compare with Figs. 3.9 and 3.10.
TABLE 3 CRITICAL LOADS OF TWO-COMPONENT COLUMNS(*)
P
CR
= ~ P E where P
E
= Euler load for a prismatic column
a = angle of offset
Values of the factor ~
A. Columns with spherical pins
11 = 11/12
a = 0 30
0
60 75
0
2 1.0000 1.0315 1.1280 1.2037
4 1.0000 1.0452 1.1917 1.3146
10 1.0000 1.0530 1.2292 1.3818
20 1.0000 1.0555 1.2415 1.4042
30 1.0000 1.0563 1.2456 1.4116
50 1.0000 1. 0570 1.2489 1.4175
1000 1.0000 1.0580 1.2536 1.4298
B. Columns with fixed ends
11 = 11/12
Ol = 15
0
30
0
45 60
0
75
0
2 1.0116 1.0446 1.0901 1.1524 1.2298
4 1.0152 1.0854 1.1790 1.3029 1.4468
10 1.0456 1.1822 1.3984 . 1.6620 1. 9429
20 1.0965 1.3536 1. 7461 2.2240 2.6740
50 1. 2238 1.8438 3.2228 3.5044 3.6792
100 1.4352 2.6498 3.5864 3.7836 3.8884
1000 3.6315 3.9380 4.0037 4.0273 4.0388
(*) Compare with Figs. 3.14 and 3.15.
90
1.2994
1.4810
1.5947
1.6327
1.6453
1.6554
1.6697
90
1.3209
1.6091
2.1876
2.9509
3.7818
3.9483
4.0442
-116
TABLE 4 CRITICAL LOADS OF PIECEWISE PRISMATIC COLUMNS ~ \ - )
P
CR
= yP
E
where P
E
= Euler Load
M = Number of segments
90

Of =
Values of the factor y:
A. Pinned columns
11 = 11/12
2
5
10
15
20
30
B. Fixed columns
11 = 11/12
2
5
10
15
20
30
M = 2 M=4
1.3121 1.3243
1.5428 1.6210
1.6246 1. 7447
1.6520 1. 7892
1.6658 1. 8121
1.6796 1.8355
M = 2 M=4
1.3209 1.4346
1. 7184 2.2937
2.i876 3.0442
2.5948 3.3870
2.9434 3.5605
3.4262 3.7399
(*) Compare with Figs. 3.16 and 3.17
-117
TABLE 5 EXPERIMENT ON PRETWISTED COLUMNS
Cross Section Properties:
Shape:
5 1
2- x 1- WF
8 2
B = 1. 838 in.
D = 2.628 in.
T = 0.202 in.
W = 0.165 in.
L = 75 in.
A = 1.110 ina
Ixx = 1.246 in
4
I = 0.210 in4
yy
I = 0.308 io6
OJ
K = o. 0131 in
4
T
Modulus of Elasticity, E = 31,374,000 psi
Summary of Test Results:
Critical
Specimen Angle of Failure Load (lb)
No. Pre twist Load (lb) (Southwell
Plot)
01 0
0
7800 8403
02 90
0
20000 25500
03 180
0
14050 21212
04 360
0
10600 12800
I
tlB=o:; {T
D X-:Jti-x
PCR/P
E
Experiment Theory
1.00 1.00
3.04 3.14
2.52 2.54
1.52 1.71
-118
-119
TABLE 6 CRITICAL LOADS OF PRETWISTED COLUMNS<*)
End Conditions: Knife edge along minor axis
P
cr
= SP
E
where P
E
= Euler load for a prismatic column
= angle of prestwist
Values of the factor s:
Ol 'Tl = I1/I2
= 1.0 'Tl = 2.0 'Tl = 5.0 'Tl = 10.0
0
0
1.0 1.0 1.0 1.0
60
0
1.5086 1. 8629 2.9037 3.6967
90
0
2.0499 2.7409 3.0990 3.2178
120
0
1.5086 2.0543 2.4849 2.6436
180
0
1.0000 1. 5126 2.3539 3.1961
210
0
1.1388 1.9968 3.6843 5.4572
270
0
2.0499 2.8252 3.5910 3.9160
360
0
1.0000 1.3391 1.6684 1. 8129
(*) Compare with Fig. 3.22
-120
TABLE 7 LATERAL TORSIONAL BUCKLING LOADS(*)
Beam: Simply supported, free warping at ends
Loading: Concentrated load at mid span of beam (at shear center)
Values of the factor A
L 2 G ~
Timoshenko Finite Element Solution
--
EI
Solution (53)
N = 4 N = 6 w
o ..
86.4/
85.6827 86.2257
.J/
/4
3 1 ~ 31.4772 31.6772
16 21.8 21.4772 21.6110
32 19.6 19.3055 19.4225
64
18.3 !
18.1135 18.2184
17.5 \
160 17.3448 17.4393
400 17.2 17.0136 17.1057
('''') Compare with Fig. 3.26
TABLE 8 SPACE FRAME BUCKLING
STRUCTURE:
one-story space frame (Fig. 3.34)
SHAPE:
DIMENS IONS:
W10x48 (all members)
Height, h = 20 r
y
Span, L = 60 r
y
BOUNDARY CONDITION: Fixed at bases
Comparison of results:
Mode
CRITICAL LOAD,
of
Buckling
Determinanta1(77)
SWAY
11. 93
TWIST
11.96
P (kips)
Finite Element
11. 73
12.21
-121
TABLE 9 BUCKLING OF SHAFTS UNDER PURE TORSION
Cross Section: I = I
x Y
EI
Critical Torque, (Mx)CR = lL
values of
Number End Condition
of
Pinned
Elements
2 6.405
3 6.322
4 6.298
5 6.289
6 6.284
Analytical(55)
6.283 (=2
TT
)
in Flexure
Fixed
---
9.862
9.305
9.182
9.094
8.988
-1221.,
Fig. 2.1 Reference Axes System for the Beam Element
Fig. 2.2 Notations for Generalized Stresses and
Displacements of the Beam Element in
Flexure and Extension
-123
"
Yo
Fig. 2.3
Axis of
y
The Centroidal-Principal Axes and the
Generalized Coordinate System
..
Fig. 2.4 Notations for a Beam Component
Subjected to Torsion and Warping
-124
-
z
, '
"
1.0
~
z
w
~ 0.8
w
()
<t
..J
a..
en
. 0 0.6
(.!)
z
0...
0::
~ 0.4
o
z
<t
t-
en
;: 0.2
t-
a
r:-x
~ B
-.-__ ==== __ ... 1-MT
1_ L .1
~ CP/CPs (Twist)
~ w/ws (Warping)
1<=25
0.2 0.4 0.6
LENGTH, X/L
Fig. 2.5 A Cantileyer Beam Loaded by Torque
5
o
0.8
-12'5'
1.0
1.0
~
z
w
~ 0.8
w
u
<t
...J
0..
(f)
00.6
CD
z
0..
a::
<t
3= 0.4
o
z
<t
f-
(f)
3= 0.2
f-
a
t:-
x
st
Mr
~ i
J
1-
L
cplcpa
(Twist)
--0-- wlwe (Warping)
0.2 0.4 0.6
LENGTH, X/L
Fig. 2.6 A Cantilever Beam Loaded by a
Concentrated Bimoment
-126
0.8 1.0
Pure Warping Mixed Torsion
I O ~ ~
Mw
MW(K=O)
a - Strain Field Formulation
o Cubic Polynomial Formulation
Hyperbolic Function Formulation
st. venant Torsion
"
r--x t
Mw
~ i - M
1_ LIT
A MT,Mw =0
8 Mw,MT=O
0.01 ' , , , I , , , :I
D I' ,=- "
0.1 1.0 10 100
BEAM PARAMETER, K=/GK
T
L2/EIw
\
1000
Fig. 2.7 Comparison of Beam Element Stiffness Matrices using a Cantilever Beam under Torsion and Bimoment
I
I-'
tv
.......
C/)
I-
Z
w
:E
w
u

c[ 1.0
C/).
o
(.!)
z
a..
0:::
0.5
3:
o
z

I-
C/)
3:
I-
o
t-
x
M
~ = : : : : : : J i T
~ L _I
I ..
0.2
~ CPlcP.B (Twist)
-0- WlwB (Warping)
0.4 0.6
LENGTH, X/L
A Strain Field
8 Cubic Polynomial
C Hyperbolic Functions
0.8 1.0
Fig. 2.8 Comparison of Beam Stiffness Matrices
I
t-'
N
CXl
-129
L L ..

-0.5
0.5l
-Q5l
I.Ol
-LOl
tOl
BIMOMENT
Mw
ST-VENANT
MOMENT
Tsv
WARPING
MOMENT
Tw
o
TORQUE
-I.Ol Mr
Fig.,2.9 Continuous Beam Loaded by
Bimoment
Node I Node 2
Mzl
r
MZ2
P-f1
x
'r
p

. tVYI tVY2 .
I.
L
-I
Node I Node 2
MWI L
Y
x MW2
Mn_ ____ P MT2
t r-
1_ L _I
Fig. 3.1 Beam Element Under Axial
and Torsional Loading
p
-
W
Fig. 3.2 Kinematics of Wide Flange
Shape Under Torsion
+
-130
(a)
(b)
w
y
Detail A
I
/
I
I
/'
.....
"-
\
\

\
\
Fig. 3.3 Schematic Description of the Kinematics of a
Shaft Under Large Torsional Displacement
-131
x
-132
3.0
9.45%
PcR
~ f ~
~
L
2.0
0 D 0 A
"/0 ERROR
0-: CRITICAL
LOAC
1.0
OIL ~ ~ ~ ~ ~ ~ = K ~ ~ ~ ~ ~ = = . I O
246 8
NUMBER OF ELEMENTS
Fig. 3.4 Convergence of Fini,te Element Solution
for Axially Loaded Columns
Cl
<t
0
...J
...J
<t
u
.....
0::
U
IJ..
0
0::
0
a::
a::
w

0
5
0
-5
-10
-15
-20
2
o Finite Differences (Wang)
Finite Element
3 4 5 6
NUMBER OF ELEMENTS, N
PcR=
17"2 EI
L2
L
Fig. 3.5 Comparison of Convergence Between Finite
Element and Finite Differences solutions
for the Euler Column
-133"
7
"j
<
y ....
x
l
x
Q=!q(X)dX
L
Distribution of
Axial Load
Fig. 3.6 A Column Under Distributed
Axial Loads
-134
-4.0
3.0
Per/Po
E
.2.0
o
-1.0
q (x)=qo
= Q c ~ L
-2.0
1 I
Cases CD @ @)
- Series Solution (Timoshenko, Dalal)
. Finite Element Solution.
( 10 Elements)
6.0 8.0
Qcr/P
E
Fig. 3.7 Comparison of Finite Element and Analytical
solutions of Columns with Distributed Axial Loads
I-'
W
\JI
0.8
0.6
L/LCOL.
0.4
Pinned -End
Column
1.0
0.8
0.6
L/L
cOL

0.4
0.2
Fixed- Free
Column
Curve PeRI P
E
a 2.0
b
c
d
1.0 (End Load Only)
0.0 (Dist. : Load Only)
":'2.0 . \'"
o 0.5 1.0 o 0.5 1.0
DEFLECTI.ON DEFLECTION
Fig.3.8 The First Buckling Modes of Columns with
Distributed Axial Loads
I-'
I.>J

RELATIVE
STIFFNESS
FACTOR
(3
\
\1 7] = 0.1 tOW]
--- Til11oshenko
-+- 10 Element
-0- 2 Element
2.0
PeR= B ~ ~ B )
IT
1.0'-' ___ ___ ......L-___ """"" ___ """"" ___ ""
o 0.2 0.4 0.6 0.8
RATIO OF MOMENT OF INERTIA, 7]=IT/1a
Fig. 3.9 Comparison of Finite Element and
Analytical Solutions of Tapered Columns
1.0
~ .
I
....
C..Q
--<J
6.0
RELATIVE
STIFFNESS 4.0
FACTOR
f3
1 'T}= 1
-- 10 Elements
4 Elements
-0- 2 Elements
2.0' "
o 2.0 4.0 6.0 8.0
J PeR =f3
IT ' L2
Is
10.0
RATIO OF MOMENT OF INERTIA, 'T}=IT/Ia
Fig. 3.10 Finite Element Solution of Critical Loads of Tapered Columns
t-"
W
00
20
10
0/0 ERROR
OF CRITICAL
LOAD
- 1 l ~
0 ~ 2 ~ ~ ~ ~ ~ ~ ~ ~ ~
I 8 10
NUMBER OF ELEMENTS
177< 1.0\
- 10 18
Fig. 3.11 Convergence of Finite Element Solution
of.Tapered Columns
200
. "

..
a::

u

(f)
(f)
w
Z
IJ...
IJ...
'-'100
(f)


w
co
I
(f)
::::>
..J
:::>
o
o

z
o
.-

o
z
:::>
o
IJ...
o
--Timoshenko, Solution
10 Elements }
o 4 Elements Finite Element Solutions

n=2
n=1
n=4 n=5
20 30-
BUCKLING LOAD FACTOR, /3 = pCr/P
E
tan-I rr2/2
p
Spring
ConStant
k (x) = k
40
Fig. 3.12 Comparison of Finite Element and Analytical Solutions of Columns on Elastic Foundations
,

c;
L
L
\
y'
\
\
X
Y
Z
-
z
--
Note: Column is
rigidly connected.
Fig. 3.13 The Generalized Forces and Displacements
of a Two-Component Column
--l4.lj
r''''''

. a.
u
II
>-.
..
a::
0
....
U
<t
LL
en
en
w
z
LL
LL
....
en
w
>
....
<t
...J
W

-142
1.3
L
t
"7 = II 1I2 :: 1000
100
1.2
f
A-A
----
Analytical (H su)
o '
Finite Element
1.1 (N = 4)
ANGLE BETWEEN COLUMN SEGMENTS,a
Fig. 3. 14 Buckling Loads of TwoComponent
Columns Pin Ends
Q
0...0
II
)0...
..
0::
0
r
u
<t
LL.
C/)
C/)
w
z
LL.
LL.
t-
C/)
W
>
t-

-.J
w
0::
L
L
A-A
2.0
'T} =II 1I2 = 1000
1.5
1.0
0 30
---- Analytical (Hsu)
--(0)--- Finite Element
(N =4)
....I)


. 60
-143
0


90
ANGLE BETWEEN COLUMN SEGMENTS, a (degrees)
Fig. 3.15 Buckling Loads of Two-Component
Columns with Fixed Ends
o.lLI
........
~
0.0
II
Ql.
..
Q::
0
I-
u

LL
en
en
w
z
LL
LL
I-
en
, w
I
> I
I-

...J
w
Q::
3
2
M = 2 Segments M = 4 Segments
L
End Conditions: Both Ends Fixed
M=CO
Analytical (Fischer)
~ Finite Element
(N = 4 Elements)
-144
10_ 10 20 30
RATIO OF MOMENT OF INERTIA I 7J = II/I2
Fig. 3.16 Buckling Loads of Piecewise Prismatic
Columns with Pin Ends
II
a::
o
.....
u

L1..
CI)
CI)
W
Z
L1..
L1..
.....
CI)
W
>
.....

-1
:W
a::
L/[
2.0
M=2 Segments M=4 Segments
jL/4
II I2
L/4
L
I2
,
II
End Conditions: Spherical Pins at Both Ends
M=CO
M=4
-hS'
M=2
~ ~ ~ ~
1.5
Ana Iytical (Fischer)
---0-- Finite Element
(N = 4 Elements)
1.0 L...L. __ -'-__ ---L ___ -L-__ -J... ___ I--__ --.L.
o 10 20 30
RATIO OF MOMENT OF INERTIA, 7J = II.lI2
Fig. 3.17 .Buckling Loads of Piecewise Prismatic
Columns with Fixed Ends
Fig. 3.18 Preparation of a Pretwisted Column
Fig. 3.19 Buckling Test of a Pretwisted Column
with Knife Edge Conditions
-146
5.0
4.0
CRITICAL LOAD
FACTOR
{3= PcR/P
E
3.0
2.0
Shape: 25/S
1
x I V2
11
'IF
End Condition: Knife Edge Along Web
'fJ = 11/12 = 6.06
Height = 75 n
A
180 360
ANGLE OF PRETWIST, a (degrees)
A Experimenta I
-- Theoretical Solution
(Finite Element Method)
\
540
Fig. 3.20 Comparison of Finite Element Solution and Experimental Results of Pretwisted Columns
I
1-"

-....J'
End Condition:
Knife Edge about Web
Length: 75"
Shape: 2
5
/8 X 1'/2 WF
Angle of Pretwist, a O ~ 90
180,270
o Theory
Experimental
I Col. 02\
a =90
0
1.0
I Col. 0 I I
a=Oo
1.0
o .
Y/8
m
1.0
Fig. 3.21 Comparison of Theoretical and Experimental Buckling
Modes of Pretwisted Columns
. .
'-. ,
t
Y
12 1.0
x/L
z'-1
i.
y
y
,t
~
o 1.0
8
z
/8
m
. 1.0
Icol. 031
a = 180
0
x/L
o 1.0
8
y
/8
m
ICol.041
x/L
a =360
0
1.0 8
y
/8
m
Fig. 3.21 (cont'd) Comparison of Theoretical and Experimental
Buckling Modes of Pretwisted Columns .
-149
End Condition: Knife Edge
Along Minor
5.0r Axis
CRI1JCAL LOAD 4.0
FACTOR
{3 = PeR / P
E
3.0
2.0

Finite Element Solution
90
0
180
0
270
0
ANGLE OF PRETWIST, a (degrees)
L
360
0
Fig. 3.22 Finite Element Solution of Pretwisted Columns with n i f ~ Edge Conditions
I
......"
In
o
a..1JJ
........
a:::
0..0
II
co.
a::
0
.--
u
<t
lL.
Q
<t
0
--'
--'
<t
U
.--
a::
u
1",-
2----
2.0
1.5
1.0
0
__ 2
"I
.",'
,,""
d;'
/
20
Principal Axes
Minor Axis: I-I
Major Axis : 2 - 2
End Conditions
Top: Pinned about 0-0
Bottom: Pinned about I-I
/
-/
/
I
I
I
I
I
/
,
,
f
,
)'
/
/
;I
/
-151
/
)1'
- - - Theory (Ashwell) for "1=1
--<>- Finite Element Solution
(N = 4 Elements)
40 60 80
ANGLE BETWEEN PINS, a (degrees) _
Fig. 3.23 Buckling Loads of Prismatic Columns
with Crossed Pins (Pin Ends)
L
..
0::
o
t-
u
it
a

o
--1
--1

u
t-
o::
u
8
I" CR 2
...... "...... ./:::..-y---a
a -----::../ "" ........ ../. a _
2 ........ ,
Principal Axes
Major Axis : I-I
Minor Axis: 2--2
End Conditions
-152'
Top: Pinned about 0-0
Bottom: Fixed
-0- Finite Element Solution
4 (N = 4 Elements)
3
20 40 60 80
ANGLE BETWEE N PINS a (degrees)
Fig. 3.24 Buckling Loads of Prismatic Columns
with Crossed Pins Ends)
1.0
0.8
%
ERROR
0.6
OF
CRITICAL
MOMENT
0.4
0.2
-0.8
21.6%
MCR MCR
0
;VI-l'
L
1 V
1-
PcR PcR

;p;
;))..
4 5 6 7
NUMBER OF ELEMENTS
o Cubic Polynomial Formulation (Iw:l= 0)
D Cubic Polynomial Formulation (I
w
=0)
Strain Field Formulation (I
w
::. 0)

Fig. 3.25 Convergence of. Finite Element Solution of
Torsional and Lateral-Torsional Buckling
-153
8
100
80
60
CRITICAL
LOAD
PARAMETER
Y2 40
20
n
w PcR
-fR - :IT\'
Jj--- - =i I 4> ;
L _I
I ..
i1:R= X2 /EIy GK
T
/L
2
Timoshenko
o Finite Element
- I
(N = 6 Elements)
0' ""I "" "' '"'
0.1 10 100 1000
BEAM PARAMETER, GK
T
L2/EIw
Fig. 3.26 Comparison of Finite Element and Analytical solutions of Lateral Torsional Buckling
I
1-'.,
U1-

-155
1.242"

EI
y
=372
I}742"
EIw= 764
400
GK
r
= 7.98
'0.1225
-
0
0
.ci
-
a.'"
300
0
..
Z
<t
a.
C/)
I-
200
::I:

IP2
(!)
0::
J3 1

;;>;
1-
60"
.1 ..
60"
-I
Q
100
<t
Trahair Test 0
0
-1

Finite Element
0 100 200 300 400
LOAD AT LEFT SPAN t P (lb.)
Fig. Buckling Loads for Two-Span Continuous Beams
-156
yT
aDr
- -
-r- ....
r-
D a
- -
- -


}

I 1----__ L __
ELEVATION

f3 8 :.:- .=- -8 ---=- =
u y
PLAN
D
END SECTION
1.-._-
8
---t-IJT
1"
w
MID SECTION
Fig. 3.28 Simply Supported Tapered I-Beam
..
200
CRITICAL
LOAD
PCR (lb)
100
0
I.
o
o
Taper
Depth, a
Width, {3
Thic kness, y
0.2
TAPER
60" _I
o
o
Trahair Finite
Experiment Element (N=6)
1:1
_ _ _
0 -0--
A --6-
0.4 0.6 0.8
CONSTANT (a,{3,y)
Fig. 3.29 Comparison of Finite Element Solution and
Experimental.Results of Lateral Torsional
Buckling of Tapered Beams
-157
1.0
! i
200
CRITICAL
LOAD
P
CR
(lb)
100
o
,PeR
F ~ : : ' < : ' ~ ~
I. 60
11
--f
--- --------
---
I
~
Trahair Finite
Taper Theory
Element (N= 6)
. Depth. a
----
-----
Width. f3 ---
-eo-
Thickness. y ---
-A--
0.2 0.4 0.6 0.8
TAPER CONSTANT (a, /3. y)
Fig. 3.30 Comparison of Finite Element Solution and
Finite Integral Solution of Tapered Beams
-158
I
1.0
-
z
-
o
<t
o
-l
-l
<t
()
-
.....
0::
()
200
(890)
100
(445)
o
PcR
-! -'=Jf-
Ie
77>-
a=0.8
0.6
0.4
0.2
60
11
(1.52 m)
0.2 0.4
-I
Variabley a:;:: 1.0
Variable a , y = 1.0
0.6 0.8
FLANGE WIDTH TAPER PARAMETER. P
Fig. 3.31 Finite Element Solutions of Beams of Combined Tapers
-159
1.0
6
N Elements
L

4 1-
L
-I
Fixed At
0/0 ERROR
OF CRITICAL
LOAD
---<>- Pinned At
2
NUMBER OF ELEMENTS, N
Fig. 3.32 Convergence Characteristics of
Knee-Bent Frames
-160
0.1
-161
P
CR
P
CR
5
4
(E III
6
(ED
1
h
h
7
A = (E1)2 L
(E1)2 8
9
Analytical
r
L
-I
Finite Element
PeR /P
E
4.0
-------
H
2.0
~ l 1
I.0L..__-o------0----------<;>- -
H
H
0.2 0.5 1.0 2 5 10
FRAME PARAMETER,A
Fig. 3: 33 Buckling of Single-Story Portal Frames

p
6
J @ @
I ~ r ~
5
9

10

L
r-------- --:
H-f I
t
---- ---/-
Sway Mode

,@
I
e
12
I
I
"
I
,... .... --"
~ ........
WIOx49
L=60 ry
h=20 ry
.....
I
..,
I
,
I ,
c.. ... ,
.... .... .... I
.... ,
.... ....,
Twist Mode
Fig. 3.34 Space Frame Buckling by Finite Element Method
-162
7 h
II
-163
(MX>CR
7i" p
.. .
(0)




----0-
(MX)CR tJ
[j
(MX)CR
-I
i

0
U
By
By
(b)
,-.......
,
SZI
-+-
, ,.
B
z


' ... _ ....

:;?"'"" <)
-0
Fig. 3.35 Buckling Modes of Shafts under pure Torsion
F
Incremental
Solution,
. F i ~ ~ ~
L Exact
Solution
F
1
Fj _I 1-----'-----------,,,
(Error) j
"Solution
Fig. 4.1 Basic Incremental. Procedure
-164
8
-165
F
Newton - Raphson
Fj ~ ~ ~ = ~ \
'-Exact
. Solution
F i _I J--------------.r
c
F.
Modified Newton - Raphson
2 34
F i ~ ~ ~ ~ ~ ~ ~ ~ ~
Fi-I f
8
L-__________ ________________ _____ ___
BotFo
8
Fig. 4.2 Newton-Raphson Methods
(a)
(b)
1.0
AXIAL
lOAD
PIPE
0.5
-166
~ o
__ ----.---e>
... --- ...... -
--
-.-
~ 1
~ P
~
I I
, I
l
I "
, I
1/
. .1.
D 5 Increments
o 10 Increments
20 Increments
Analytical
8
0
= L/500
!\1=10
... ~ ... ___ .. L- ___ --1-.,_. __ . _--,-_______ .. L
1.0 2.0
lATERAL DEFLECTION, 8/l x 10-
2
Flg 4.3 Effect of Number of Load Increments
when Us ing tJ:le Incremental Method
!."
AXIAL
LOAD
PIPE'
0.5
L
Euler Load,
a
~
J
I I
I I
I ,
I I
I '
-167
N = 10 Elements
No. of Load Inc. = 20
o 1.0
LATERAL DEFLECTION I all x 10-
2
Fig. 4.4 Use of the Incremental Method to Determine
Load-Deflection Curves of Columns
2.0
_ 1.6
o
<t
g
o
W
...J
o.S
8: 0.6
<t
a
L
- Timoshenko
-0- Finite Element
(8
0
= L/500)
0.2 0.4 0.6 0.8
lATERAL DEFLECTION, 8/l
Fig. 4.5 Use of the Incremental Method for the
Analysis of Large Deflections '
, i
1.0
-169
9,. REFERENCES
1. Todhunter, I. and Pearson, K.
A HISTORY FROM GALlLEI TO LORD KELVIN, OF THE THEORY OF
EIASTICITY AND THE STRENGTH OF MATERIALS, Dover Pub., Inc.,
New York, 1960.
2. Love, A. E. H.
A TREATISE ON THE MATHEMATICAL THEORY OF EIASTICITY, Dover
Publications, 4th Ed., 1944.
3. Truesdell, C. and Toupin, R.
THE CIASSICAL FIELD THEORIES, Encyclopedia of Physics,
Edited by S. F1Ugge, Springer-Verlag, 1960.
4. Medick, M. A.
ONE-DIMENSIONAL THEORIES OF WAVE PROPAGATION AND VIBRATIONS
INELASTIC BARS OF RECTANGULAR CROSS SECTION, Journal of
Applied Mechanics, Vol. 33, Sept. 1966.
5. Soler, A. I.
HIGHER-ORDER THEORIES FOR STRUCTURAL ANALYSIS USING LEGENDRE
POLYNOMIAL EXPANSIONS, Journal of Applied Mechanics, Dec.
1969.
6. Hertelendy, P.
AN APPROXIMATE THEORY GOVERNING SYMMETRIC MOTIONS OF ELASTIC
RODS OF RECTANGUIAR OR SQUARE CROSS SECTIONS, Journal of
Applied Mechanics, June 1968.
7. Timoshenko, S. P.
HISTORY OF STRENGTH OF MATERIALS, McGraw-Hill Book Company,
Inc., New York, 1953.
8. Wagner, H.
TORSION AND BUCKLING OF OPEN SECTIONS, 25th Anniversary
Publication, Technische Hochschu1e Danzig, 19-04-1929,
Translated in NACA TM. 807, 1936.
9. Goodier, J. N.
TORSIONAL AND FLEXURAL BUCKLING OF BARS OF THIN WALLED OPEN
SECTIONS UNDER COMPRESSIVE AND BENDING LOADS, Journal of
Applied Mechanics, ASME, Vol. 64, Sept. 1942.
10. Timoshenko, S.
THEORY OF BENDING, TORSION AND BUCKLING OF THIN-WALLED MEMBERS
OF OPEN CROSS SECTION, Journal, Franklin Institute, Vol. 239,
No.3, March 1945.
11. Vlasov, V. Z.
THIN WALLED BEAMS, English translation published for U.S.
Science Foundation by Israel Program for Scientific Trans-
lations, 1961.
-170
12. Chen, W. F. and Santathadaporn, S.
REVIEW OF COLUMN BEHAVIOR UNDER BIAXIAL LOADING, Journal of
the Structural Division, ASCE, Vol. 94, ST12, Dec. 1968.
13. Hrennikoff, A.
SOLUTION OF PROBLEMS IN ELASTICITY BY THE FRAMEWORK METHOD,
Journal of Applied Mechanics, Vol. 8, 1941.
14. Ang, A. H. S. and Newmark, N. M.
A NUMERICAL PROCEDURE FOR THE ANALYSIS OF CONTINUOUS PLATES,
Second Conference on Electronic Computation, ASCE, 1960.
15. Zienkiewicz, O. C.
THE FINITE ELEMENT METHOD IN ENGINEERING SCIENCE, McGraw-Hill
Book Company, Inc., London, 1971.
16. Oden, ~ T.
FINITE ELEMENTS OF NONLINEAR CONT INUA, McGraw-Hill Book
Company, Inc., 1972.
17 . Oden, J. T.
A GENERAL THEORY OF FINITE ELEMENTS, Part 1 and Part 2, Int.
Journal Numerical Methods Eng., Vol. 1, No.2 and No.3, 1969.
18. Archer, J. S.
CONSISTENT MATRIX FORMULATIONS FOR STRUCTURAL ANALYSIS USING
INFLUENCE COEFFICIENT TECHNIQUES, AIAA Journal, Vol. 3, Oct.
1965.
19. Prezemieniecki, J. S.
THEORY OF MATRIX STRUCTURAL ANALYS IS, McGraw-Hill, New York,
1968.
20. Kabaila, A. P. and Fraeijs de Veubeke, B.
STABILITY ANAIXSIS BY FINITE ELEMENTS, Air Force Flight
Dynamics Laboratory, Tech. Report AFFDL-TR-70-35, March
1970.
21. Thomas, J. M.
A FINITE ELEMENT APPROACH TO THE STRUcTURAL INSTABILITY OF
BEAM COLUMNS, FRAMES AND ARCHES, NASA TN D-5782, May 1970.
22. Barsoum, R.
A FINITE ELEMENT FORMULATION FOR THE GENERAL STABILITY
ANALYSIS OF THIN-WALLED MEMBERS, doctoral dissertation,
Cornell University, 1970, University Microfilms, Ann Arbor,
Michigan.
23. Gallagher, R.
A CORRELATION STUDY OF METHODS OF MATRIX STRUCTURAL ANALYSIS,
Pergamon Press, Oxford, 1964.
24. Washizu, K.
SOME REMARKS ON BASIC THEORY FOR FINITE ELEMENT METHOD,
Proc. Japan-U.S. Seminar on Matrix Methods in Structural
Analysis and Design, Tokyo, 1969.
25. pian, T. H. H.
-171
FORMUlATION OF FINITE ELEMENTS FOR SOLID CONTINUA, Proe.
Japan-U.S. Seminar on Matrix Methods in Struetuetura1 Analy-
sis and Design, Tokyo, 1969.
26. Washizu, K.
MARIATIONAL METHODS IN ELASTICITY AND PLASTICITY, Pergamon
Press, 1st Ed., 1968.
27. Hu, H. C.
ON SOME VARIATIONAL PRINCIPLES IN THE THEORY OF EIASTICITY
AND PLASTICITY, Scintia Sinica, Vol. 4, No.1, March 1955.
28. Prager, W.
VARIATIONAL PRINCIPLES OF LINEAR ELASTOSTATICS FOR DIS-
CONTINUOUS DISPLACEMENTS, STRAINS, AND STRESSES, The
Folke-Odquist Volume, edited by B. Brobert, J. Hult, and
F. Niordson, Almquist & Wikse11, Stockholm, 1967.
29. Reissner, E.
NOTE ON THE METHOD OF COMPLEMENTARY ENERGY, Journal of
Mathematics and Physics, Vol. 27, 1948.
30. Pian, T ~ H. H.
ELEMENT STIFFNESS MATRICES FOR BOUNDARY COMPATIBILITY AND
FOR PRESCRIBED BOUNDARY STRESSES, Proc. Conf. on Matrix
Methods in Structural Mechanics, AFFDL-TR-66-80, 9 6 5 ~
31. Timoshenko, S. P.
ON THE CORRECTION FOR SHEAR OF THE DIFFERENTIAL EQUATION FOR
TRANSVERSE VIBRATIONS OF PRISMATIC BARS, Philosophical
Magazine, Vol. 41, 1921.
32. Timoshenko, S. P. and Goodier, J. N.
THEORY OF ELASTICITY, McGraw-Hill Book Company, Inc., New
York, 3rd Edition, 1970.
33. Soko1nikoff, I. S.
MATHEMATICAL THEORY OF ELASTICITY, McGraw-Hill Book Company,
Inc., New York, 2nd Edition, 1956.
34. Cowper, G. R.
THE SHEAR COEFFICIENT IN TIMOSHENKO'S, BEAM THEORY, Journal
of Applied Mechanics, Vol. 33, No.2, June 1966.
35. Mason, W. E. and Hermann, L. R.
EIASTIC SHEAR ANALYSIS OF GENERAL PRISMATIC BEAMS, Journal
of Engineering Mechanics Division, ASCE, EM4, August 1968.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
-172
Lauterbach, G. F., Peppin, R. J., Schapiro, S. M., and Krahu1a, J. L.
A FINITE ELEMENT SOLUTION FOR SAIN-VENANT BENDING,
Journal, Vol. 9, No.3, March 1971.
Popov, E. P.
KELVIN'S SOLUTION OF TORSION PROBLEM, Journal of Engineering
Mechanics Division, ASCE, Vol. 96, EM6, Dec. 1970.
Krahu1a, J. L.
ANALYSIS OF BENT AND TWISTED BARS USING THE FINITE ELEMENT
METHOD, AIAA Journal, Vol. 5, No.6, June 1967.
Krajcinovic,' D.
A CONSISTENT DISCRETE ELEMENTS TECHNIQUE FOR THIN WALLED
ASSEMBLAGES, Int. J. Solids Structures, Vol. 5, 1969.
Janssen, J. D.
DE MOGELIJKHEDEN VAN DE METHODE DER EINDIGE ELEMENTEN BIJ
DE BEREKENING VAN DE STERKTE EN STIJFHEID VAN DUNWANDIGE
BALKEN (THE APPLICATION OF THE FINITE ELEMENT METHOD FOR
THE EVALUATION OF THE STRENGTH AND STIFFNESS OF THIN-WALLED
BEAMS), De Ingenieur, Vol. 82, Nr. 43, October 1970.
Heins, 'C. P., and Seaburg, P. A.
TORSION ANALYSIS. OF ROLLED STEEL SECTIONS, AIA File. No.
13-A-1, Bethlehem Steel, Bethlehem, Pa., 1963.
Basler, K. and Ko1lbrunner, C. F.
TORSION IN STRUCTURES, Springer-Verlag, New York, 1969.
Martin, H. C.
ON THE DERIVATION OF STIFFNESS MATRICES FOR THE ANALYSIS OF
IARGE DEFLECTION AND STABILITY PROBLEMS, Proc. AFIT Conf.
on Matrix Methods in Struct. Mech., AFFDL TR 66-80, 1965.
Mallett, R. and Marcal, P. V.
FINITE ELEMENT ANALYSIS ON NONLINEAR STRUCTURES, Journal of
Str'uctura1 Division, ASCE, Vol. 94, ST9, Sept. 1968.
Oden, J. T.
FINITE ELEMENT APPLICATIONS IN NONLINEAR STRUCTURAL ANALYSIS,
Proc. Symposium -Application of Finite Element Methods in
Civil Engineering, Vanderbilt University, 1969.
Martin, H. C.
FINITE ELEMENTS AND THE ANALYSIS OF GEOMETRICALLY NONLINEAR
Proc. Japan-U.S. Seminar on Matrix Methods in
Struc,tural Analysis and Design, Tokyo, 1969.
47. Gallagher, R.
GEOMETRICALLY NONLINEAR FINITE ELEMENT ANALYSIS, Proc.
Specialty Conf. - Finite Element Method in Civil Engineering,
McGill University, 1972.
48.
49.
-173
Novoshilov V. V.
. FOUNDATIONS OF THE NONLINEAR THEORY OF ELASTICITY, Graylock
Press, Rochester, N. Y., 1953.
Prezemieniecki, J. S.
DISCRETE-ELEMENT METHODS FOR STABILITY ANALYSIS OF COMPLEX
STRUCTURES, The Aeronautical Journal of the Royal Aeronautical
Society, Vol. 72, Dec. 1968.
50. Dzhanelidze, G. Y.
ON THE THEORY OF THIN AND THIN-WALLED RODS, NACA TM 1309,
Oct. 1951.
51. Golenveizer, A. L.
THEORY OF THIN-WALLED RODS, NACA TM 1322, Oct. 1951.
52. Gallagher, R. and padlog, J.
DISCRETE ELElNT APPROACH TO STRUCTURAL INSTABILITY ANALYSIS,
AIAA Journal, Vol. 1, No.6, June 1963.
53 .. Timoshenko, S. P. and Gere, J. M.
THEORY OF ELASTIC STABILITY, 2nd Edition, McGraw-Hill Book
Company, Inc., New York, 1961.
54. Bleich, F.
BUCKLING STRENGTH OF METAL STRUCTURES, McGraw-Hill Book
Company, Inc., New York, 1952.
55. Ziegler, H. .
PRINCIPLES OF STRUCTURAL Blaisdell Publishing
Company, London, 1968.
56. Wang, C. T.
APPLIED ELASTICITY, McGraw-Hill Book Company, Inc., New York,
1953.
57. Dalal, T. S.
COLUMNS WITH DISTRIBUTED AXIAL LOADS, The Structural Engineer,
Vol. 48, No. 11, Nov. 1970.
58. Gere, J. M. and Carter, W. O.
CRITICAL BUCKLING LOADS FOR TAPERED COLUMNS, Journal of
Structural Division, ASCE, Vol. 88, STl, February 1962.
59. Wang, L. R.
PARAMETRIC MATRICES OF SOME STRUCTURAL MEMBERS, Journal of
Structural Division, ASCE, Vol. 96, ST8, Aug. 1970.
60. Girijavallabhan, C. V.
BUCKLING LOADS OF NONUNIFORM COLUMNS, Journal of Structural
Division, ASCE, Vol. 95, STll, 1969.
61. Hetenyi, H.
BEAMS ON ELASTIC FOUNDATION, The University of Michigan
Press, Ann Arbor, Michigan, 1964.
62. Toak1ey, A. R.
-174
BUCKLING LOADS FOR ELASTICALLY SUPPORTED STRUTS, Journal of
Engineering Mechanics Division, ASCE, Vol. 91, EM3, June
1965.
63. Wieghardt, K.
UBER DEN BALKEN AUF NACHGIE-BIGER UNTERLAGE (BEAMS ON FLEX-
IBLE FOUNDATIONS), Zeitschrift FUr Angewandte Mathematic
und Physic, Vol. 2, 1922.
64. Fi1enko-Borodich, M. M.
SOME APPROXIMATE THEORY OF ELASTIC FOUNDATIONS, Uchenye
Zapiski, MGU, No. 46, 1940.
65. V1asov, V. Z. and Leontev, N. N.
BEAMS, PLATES AND SHELLS ON ELASTIC FOUNDATIONS (Trans 1ated
from Russian), NASA and NSF, Washington, D. C., 1966.
66. Biot, M. A.
BENDING OF AN INFINITE BEAM ON ELASTIC FOUNDATION, Journal
of Applied Mechanics, ASME, 36, 1937.
67. Hsu, Y.
ELASTIC BUCKLING OF TWO-COMPONENT COLUMNS, Journal of Engin-
eering Mechanics, ASCE, Vol. 95, EM3, June 1969 .
. 68. Fischer, W.
VERGLEICH DER IDEALEN KNICKLASTEN VERWUNDENER UND VERSETZTER
STABE (Comparison of Ideal Buckling Loads of Discontinuuosly
Twisted Columns), Stalbau, Vol. 39, No.2, February 1970.
69. Ziegler, H.
DIE KNICKUNG DES VERWINDEN STABES (Buckling of Twisted Bars),
Schweizerische Bauzertung, Vol. 66, 1948.
70. Zickel, J.
PRETWISTED BEAMS AND COLUMNS, Journal of Applied Mechanics,
Vol. 23, 1956.
71. Fischer, W.
DIE KNICKLAST DES SCHEIDENGELAGERTEN VERWUNDENEN STABES
(Buckling Load of Pretwisted Rods on Knife Edge Supports),
Ingenieur-Archiv, Vol. 39, 1970.
72. Tebedge, N. and Tall, L.
PROCEDURE FOR TESTING CENTRALLY-LOADED COLUMNS, Fritz
Engineering Lab Report No. 351.6, Dec. 1971.
73. Ashwell, D.
THE EULER BUCKLING LOAD OF A STRUT WITH CROSSED PINS,
Engineering, 1950, pp. 709.
74. Trahair, N. S.
-175
ELASTIC STABILITY OF CONTINUOUS BEAMS, Journal of Structural
Division, ASCE, Vol. 95, ST6, June 1969.
75. Kitipornchai, S. and Trahair, N. S.
ELASTIC STABILITY OF TAPERED I-BEAMS, Journal of Structural
Division, ASCE, Vol. 98, ST3, March 1972.
76. Lu, L. W.
A SURVEY OF LITERATURE ON THE INSTABILITY OF FRAMES, Welding
Research Council Bulletin Series No. 81, 1962.
77. Morino, S. and Lu, L. ,W.
ANALYSIS OF SPACE FRAMES, Fritz Engineering Lab Report
No. 331.13, March 1971.
78. Marca1, P. V. and King, I. P.
ELASTIC-PLASTIC ANALYSIS OF TWO-DIMENSIONAL STRESS SYSTEMS
BY THE FINITE ELEMENT METHOD, Int. J. Mech. Science, Vol.
9, 1967.
79. Yamada, Y., Yoshimura, N. and Sakurai, T.
PLASTIC STRESS-STRAIN MATRIX AND ITS APPLICATION FOR THE
SOLUTION OF ELASTIC-PLASTIC PROBLEMS BY THE FINITE ELEMENT
METHOD, Int. J. Mech. Science, Vol. 10, 1968.
80. Marca1, P. V.
FINITE ELEMENT ANALYSIS WITH MATERIAL NONLlNEARITIES--
THEORY AND PRACTICE, Proc. Specialty Conf.--Finite Element
Method in Civil Engineering, McGill University, 1972.
81. Turner, M. J., Dill, E. H., Martin, H. C. and Melosh, R. J.
LARGE DEFLECTION OF STRUCTURES SUBJECTED TO HEATING AND EXTERNAL
LOADS, Journal of Aerospace Sciences, Vol. 27, 1960.
82. Haisler, W. E., Stricklin, J. A. and Stibbins, F. J.
DEVELOPMENT AND EVALUATION OF SOLUTION PROCEDURES FOR GEO-
METRICALLY NONLINEAR STRUCTURAL ANALYSIS, AIAA Journal,
Vol. la, No.3, March 1972.
83. Carnegie, W.
STATIC BENDING OF PRE-TWISTED CANTILEVER BLADING, Proceedings
of the Institution of Mechanical Engineers, Vol. 171, 1957.
84. Houbolt, J. C. and Brooks, G. W.
DIFFERENTIAL EQUATIONS OF MOTION FOR COMBINED FLAPWISE BENDING,
CHORDWISE BENDING, AND TORSION OF TWISTED NON-UNIFORM ROTOR
BLADES, TN 1346, 1958, NACA.
85. Lawrence, K. L.
TWISTED BEAM ELEMENT MATRICES FOR BENDING, AIAA Journal, Vol.
8, No.6, June 1970.
la. VITA
The author was born the first son of Debre Boga1e and
Tebedge Guddaye in Debre Tabor, Ethiopia on October 16, 1944.
-176
He graduated from Medhane A1em High School in Addis Ababa
in 1960. From September 1960 to June 1965 he studied at the College
of Engineering in Haile S"e1assie I University. There he received the
degree of BacHelor of Science in Civil Engineering with Distinction
and was awarded the Haile Se1assie I University Chancellor's Medal.
In June 1965 he joined the College of Engineering in H.S.I.U.
as a Graduate Assistant and continued working in the same university
as an Assistant Lecturer until August 1967.
In September 1967 the author received a scholarship, sponsored
by the African-American Institute, to study in the Department of Civil
Engineering at "Lehigh University. In June 1968 he joined the staff
at Fritz Engineering Laboratory as a Research Assistant. He has been
associated with the research projects on Residual Stresses in Thick
Welded plates and also on European Column Studies. He was awarded
the degree of Master of Science in Civil Engineering in June 1969.