Static N Dynamic Analysis of Structures Computer Matrix
Static N Dynamic Analysis of Structures Computer Matrix
by
JAMES F. DOYLE
Aeronautics and Astronautics Department,
Purdue University, West Lafayette,
Indiana, U.S.A.
"
~
Preface xi
2 Rod Structures 14
2.1 Rod Theory . . . . . . . . . . 14
2.2 Rod Element Stiffness Matrix 20
2.3 Structural Stiffness Matrix .. 22
2.4 Boundary Conditions . . . . . 26
2.5 Member Distributions and Reaction 29
2.6 Distributed Loads. 32
Problems 35
Exercises . 36
VII
Vlll Contents
5 Structural Stability 97
5.1 Elastic Stability . 97
5.2 Stability of Truss Structures . . . . . . 99
5.3 Matrix Formulation for Truss Stability 105
5.4 Beams with Axial Forces . . . . . . .. 109
5.5 Beam Buckling . . . . . . . . . . . . . 115
5.6 Matrix Analysis of Stability of Beams 120
5.7 Stability of Space Frames 128
Problems 133
Exercises . 135
Problems 381
Exercises 381
Appendices
References 437
Index 439
Preface
This book is concerned with the static and dynamic analysis of structures. Specifi-
cally, it uses the stiffness formulated matrix methods for use on computers to tackle
some of the fundamental problems facing engineers in structural mechanics. This is
done by covering the Mechanics of Structures, its rephrasing in terms of the Matrix
Methods, and then their Computational implementation, all within a cohesive setting.
Although this book is designed primarily as a text for use at the upper-undergraduate
and beginning graduate level, many practicing structural engineers will find it useful
as a reference and self-study guide.
Several dozen books on structural mechanics and as many on matrix methods are
currently available. A natural question to ask is why another text? An odd devel-
opment has occurred in engineering in recent years that can serve as a backdrop to
why this book was written. With the widespread availability and use of comput-
ers, today's engineers have on their desk tops an analysis capability undreamt of by
previous generations. However, the ever increasing quality and range of capabilities
of commercially available software packages has divided the engineering profession
into two groups: a small group of specialist program writers that know the ins and
outs of the coding, algorithms, and solution strategies; and a much larger group of
practicing engineers who use the programs. It is possible for this latter group to use
this enormous power without really knowing anything of its source. Therein, in my
opinion, lies the potential danger - the engineer is seduced by the power, the litany
of capabilities, the seeming ease of use, and forgets the old computer adage: garbage
in, garbage out, or its more recent sinister form: garbage in, gospel out. We use,
and we should use, commercial packages when they are available. But to make safe,
efficient, and intelligent use of them we need to have some idea of their inner workings
as well as the mechanics foundations on which they are built.
It would seem reasonable, therefore, that the following abilities must be considered
a minimum for any structural engineer:
Know how to idealize structures in a sensible manner, and know when (and
why) a structure is beyond the capabilities of a particular computer program.
Know how to use consistency and cross checks to validate the computer output.
XI
Xli Preface
Although these are the input/output brackets to the program, the form they take is
dictated by the internal modeling. Hence, it is impossible to avoid some consideration
of the inner workings of the programs. That is, to be an intelligent user of these pro-
grams requires some appreciation of the full range of assumptions and procedures on
which they are based. Without doubt, an understanding of the mechanics principles
is essential, but it is not sufficient, because these principles are transformed in subtle
ways when converted into algorithms and code.
With the foregoing in mind, this book sets as its goal the treatment of struc-
tural dynamics starting with the basic mechanics principles and going all the way to
their implementation on digital computers. I believe that only by studying this in
its complete extent do the unique difficulties of computational mechanics manifest
themselves. I have made an effort to ensure that this book is not just a collection of
disparate topics - I sought a unity that would give meaning to the pieces as part
of a coherent whole. This is achieved by covering completely the analysis of 3-D
frame structures and using it as a paradigm for treating other structural systems. In
the same vein, rather than discuss particular commercial packages, I have developed
STADYN: a complete (but lean) program to perform each of the standard procedures
used in commercial programs. Each module in this program is reasonably complete
in itself, and all were written with the sole aim of clarity plus a modicum of effi-
ciency, compactness and elegance. (I have included complete source code listings to
these modules in an appendix; their electronic form can be obtained through ikayex
SOFTWARE TOOLS, 615 ELSTON ROAD, LAFAYETTE, INDIANA 47905, USA.)
This book takes a bootstrapping approach to developing the material. That is,
the elemental blocks are developed on first principles; these are then combined to
model more complicated problems. It is only then that the general principles are
established. For example, Chapter 2 looks at the simplest of structures, the rod. The
analysis is developed fully from the governing differential equations all the way to the
matrix formulation. Chapter 3 discusses beam structures and develops the analysis
in nearly identical fashion. The following chapter then shows how these elementary
structures can be combined to form complicated 3-D structures. The twin concepts
of compatibility and equilibrium are emphasized throughout these three chapters.
Chapter 5 refines the concept of equilibrium in the process of discussing structural
stability. Chapter 6 introduces the energy methods, and shows how the concept of the
stationary potential energy (coupled with the Ritz method) can be used to recover
the results of the previous chapters in a unified and consistent fashion. The first
part of the book concludes with an introduction to the computational aspects of the
matrix methods in Chapter 7.
The second part of the book begins with a summary of the vibration of a simple
spring-mass system. Chapters 9,10&11 are the dynamic analogs of Chapters 2,3&4;
they achieve the matrix formulation for the dynamic analysis of 3-D continuous struc-
tures. Chapter 11 introduces modal analysis as a means of understanding the com-
plicated dynamics on a structural or global level. These results are put into a unified
Preface Xlii
form in Chapter 12 using the energy concepts. Again, the Ritz method emerges as a
powerful technique for converting continuous systems into discrete matrix form. The
computational methods for direct integration and modal analysis are developed in
Chapter 13.
Appendices dealing with matrix algebra, spectral analysis, and computer source
code round out the coverage of the material. Admittedly, many problems (such as
those associated with thermal loading, non-linear material and non-linear geometric
effects) have been left out, even though most of them are treatable by matrix meth-
ods. Consequently, an effort is made to supplement each chapter with a collection of
pertinent problems that indicate extensions of the theory and the applications. These
problems, combined with selected references to relevant literature, can form the basis
for further study.
A book like this is impossible to complete without the help of many people, but
it is equally impossible to properly acknowledge them all individually. However, I
would like to single out: Graham Gladwell for his very many helpful suggestions and
his understanding of what I was trying to r-chieve; Albert Danial, Shiv Joshi, Tim
Norman, and Steve Rizzi for their construc~ive criticisms of the almost final version
of the manuscript; and all my former students who suffered through the early drafts,
their feedback was invaluable. The errors and inaccuracies remaining in this book are
purely my own doing; I would deem it a kind service to be informed of them.
I used a combination of UTEX and PostScript to typeset this manuscript; I thank
all those people who made these wonderful systems available for the desktop com-
puter.
"I am quite convinced that the central question is how to make computers
more usable, how to make their software more comprehensible, and how to
avoid the dangers imposed by the complexities of standard software in the
current generation."
C. A. R. HOARE
Chapter 1
Structural mechanics is concerned with the behavior of solid objects (or assemblies
of them) under the action of applied loads. The behavior is usually described by
idealized models from which the internal forces and displacements are found. We
present, in the following chapters, the stiffness formulated matrix methods as models
to tackle some of the fundamental problems facing engineers in structural mechanics.
It is a happy coincidence that matrices are also a convenient means for carrying
out calculations on a computer and therefore it is only natural that an intimate
relation has evolved between structural analysis and computers. The two are linked
and modern courses must reflect that. We will attempt to cover the Mechanics of
Structures, its rephrasing in terms of the Matrix Methods, and then the Computational
implementation, all within a cohesive setting. This chapter sets out to describe the
context within which we will do this. References [11,23, 33, 35] are very useful sources
for additional background material.
1
2 Chapter 1. Background and Scope
lions). It is necessary to determine the first in order to know whether the structure is
capable of withstanding the applied loads. The second must be determined to assure
that excessive displacements do not occur. To understand the construction of an an-
alytical structural model, therefore, requires an understanding of the function of the
particular structure as well as the requirements of the analysis. In other words, if the
structure does not experience dynamic loads then a dynamic model need not be devel-
oped. Similarly, if a stability analysis is of interest then this capability must be built
into the model from the beginning. It is not surprising that structural models abound;
we will attempt to keep them to a minimum by stressing the underlining principles
of model building. We will, however, set for ourselves the tasks of developing models
to describe the static, stability, and dynamic responses of structures.
The modeling of the dynamic response of structures introduces many additional
considerations probably not anticipated from a static analysis. It is therefore worth
our while at this juncture to say a few words about structural dynamics. The subject
of rigid-body dynamics treats the physical objects as bodies that undergo motion
without any change of shape. This has many applications: the movement of machine
parts; the flight of an aircraft or space vehicle; the motion of the earth and the plan-
ets. In many instances, however, the primary concern is dynamic response involving
changes of shape. This is particularly so in the design of structures and structural
frames as encountered in automobiles, ships, aircraft, space vehicles, offshore plat-
forms, buildings, and bridges.
Dynamic response involving deformations is usually oscillatory in nature; the
structure vibrates about a configuration of stable equilibrium. For example, sup-
pose a building structure is in a state of static equilibrium under the gravity loads
acting on it. When subjected to wind, the structure will oscillate about this position
of static equilibrium. An airplane provides an example of oscillatory motion about
an equilibrium configuration that involves rigid-body motion. When in flight, the
whole system moves as a rigid body but is also subjected to oscillatory motion due
to engine and aerodynamic loads.
The analysis of vibration response is of considerable importance in the design of
structures that may be subjected to dynamic disturbances. Under certain situations
vibrations may cause large displacements and severe stresses in the structure. This
may happen when the frequency of the exciting force coincides with a natural fre-
quency of the structure and resonance ensues. A related problem is that fluctuating
stresses, even of moderate intensity, may cause material failure through fatigue and
wear. Also, the transmission of vibrations to connected structures may lead to un-
desirable consequences: delicate instruments may malfunction or human occupants
may suffer discomfort.
With the increasing use being made of lightweight, high-strength materials, struc-
tures today are more susceptible than ever before to critical vibrations. Modern build-
ings and bridges are lighter, more flexible, and are made of materials that provide
much lower energy dissipation; all of these may contribute to more intense vibra-
1.2 Types of Structures Considered 3
The individual members of frames are modeled in terms of the following three
elemental forms: Rod, Shaft, and Beam. The difference between these is not so much
their geometric shape but rather the types of loading they support and the type of
resulting deformation. In this book
Rods support only axial loads (tensile or compressive) and the resulting defor-
mation is only along the length.
Shafts support only a torque acting along its length resulting in only axial twist.
Beams are designed to carry both moments and transverse loads.
Fundamental to the description of 3-D frames is the assumption that the response
of a general member is a simple superposition only of the above three actions. As a
result, these actions can be developed separately and only as a final step need they
be combined.
A truss consists of a collection of arbitrarily oriented rod members that are inter-
connected at pinned joints. They are usually loaded only at their joints and (because
the joints cannot transmit bending moment) must be triangulated to avoid collapse.
Plane trusses are loaded in their own plane, whereas the joints of a space truss can
be loaded from any direction.
A frame structure, on the other hand, is one that consists of beam members which
are connected rigidly or by pins at the joints. The members of a frame can support
bending (in any direction) as well as axial loads, and at the rigid joints the relative
positions of the members remain unchanged after deformation. Rigidly jointed frames
4 Chapter 1. Background and Scope
Space frame
are often loaded along their members as well as at their joints. Plane frames, like
plane trusses, are loaded only in their own plane. In contrast, grids are always loaded
normal to the plane of the structure. Space frames can be loaded in any plane. The
space frame is the most complicated type of jointed framework. Each member can
undergo axial deformation, torsional deformation, and flexural deformation (in two
planes). Its supports may be fixed, pinned, elastic, or there may be roller supports.
We will concentrate on the space frame because all other types of jointed frame-
works are special cases that can be obtained by reduction from it. Further, it is
assumed throughout most of the subsequent discussions that the structures have pris-
matic members; that is, each member has a straight axis and uniform cross-section
throughout its length. The non-prismatic case is developed as part of the exercises.
Equilibrium
Loads can be regarded as being either internal or external. External loads consist
of applied forces and moments and the corresponding reactions. The applied loads
usually have preset values, whereas the reactions assume values that will maintain
the equilibrium of the structure. Internal forces are the forces generated within the
structure in response to the applied loads. We will sometimes refer to them also as
member loads. The applied loads may be concentrated forces, distributed loads, or
couples.
There are two types of equilibrium encountered in this book: static and dynamic.
To illustrate these two concepts, consider a single mass point under the action of a
system of forces. Newton's second law gives
This can be read as sayin&. that the resulting forces ('L, F) cause the mass (m) to
accelerate by the amount (17). If there are no resulting accelerations, then
Constitutive Behavior
We shall restrict ourselves to materials that are linearly elastic. By elastic we mean
that the stress-strain curve is the same for both loading and unloading. The restric-
tion to linear behavior allows us to use the very important concept of superposition.
Fortunately, most materials of interest in structural applications are modeled ade-
quately by this and so we will always assume that the materials behave according to
Hooke's law. For example, the one-dimensional stress-strain behavior is given by
or F= EA u
L
The second form is that of the force-deflection relation.
There is one exception to this. The dynamic response of structures usually exhibit
some type of dissipation of energy. (This is referred to as damping.) Since the amount
of damping is usually not significant, and since (for structures) it is difficult to get
precise experimental values, we will generally describe the dissipation using a linear
viscoelastic model. In so doing, we are allowing the material behavior to be time
dependent.
The principle of superposition is one of the most important concepts in structural
analysis. In general terms the principle states that the effects produced by several
causes can be obtained by combining the effects due the individual causes. It may be
used whenever linear relations exist. This occurs whenever the following requirements
are satisfied: (1) the displacements are small; (2) the material obeys Hooke's law; and
(3) there is no interaction between flexural and axial effects in the members. This
principle is fundamental to the stiffness method of analysis and therefore it will always
be assumed that the structure being analyzed meets the stated requirements.
caused by the applied loading. Obtaining the responses at points other than nodes
is then properly viewed as part of the post-processing operation. Thus an essential
part of the description is determining the minimum number of unknowns necessary
for an adequate characterization of the structural response.
The number of possible displacement components at each node is known as the
nodal degree of freedom (DoF); the nodal degree of freedom for different framework
types is shown in the following table:
where
Px , Py , Pz are external force components,
Tx , T y , T z are external moment (torque) components,
u, v, w a r e global linear displacement components,
<Px, <Py, <Pz are global rotational displacement components.
The degree of freedom of a structure is the number of displacement components to
be found during the analysis. Finding these displacements by the stiffness method
involves the solution of a system of equations relating the unknown nodal displace-
ments to the known applied loads. We will distinguish between unconstrained (i.e.,
unknown) and constrained (i.e., known in some way) degrees offreedom. The number
of equations involved is equal to the total number of unconstrained degrees of freedom
of the structural system.
The loads and displacements are conveniently stored in one dimensional arrays
known as the nodal load vector {p} and nodal displacement vector {u}, respectively.
The number of entries in each vector is equal to the nodal degree of freedom of the
system as shown above. When these vectors have the entries associated with the
constrained degrees of freedom removed they are often referred to as generalized force
8 Chapter 1. Background and Scope
and displacement vectors, respectively. In our description the loads will be thought
of as being "applied," while the displacements are thought of as "resulting." This
distinction may seem pointless since we can obviously apply a known displacement to
a structure and seek the resulting force. But the distinction is actually very useful in
the stiffness method because we will show that all structural problems can be phrased
as
[K]{u} + [C ]{u} + [M]{il} = {p}
where [K 1 is the stiffness matrix, [M 1the mass matrix, and [c 1is the damping
matrix. In fact, we will even rephrase the imposed displacement problem as an applied
force situation.
Load
..........,.,................ ; .........., .._...............
: "
... ...
Time
........._ .
pressur~~
/
.
, ..
..........
Time
A dynamic load that varies in magnitude with time, and that repeats itself at
regular intervals, is called a periodic load. The harmonic load caused by an unbalanced
rotating machine is an example of a periodic load. Loads that do not seem to repeat
are called non-periodic loads. Such a load may be of a comparatively long duration,
in which case the response is usually modeled as being static or quasi-static. Non-
periodic loads may also be of short duration and are then called transient or impulsive
loads; an explosive blast striking a building is an example.
We will investigate the response of structures to periodic and transient loads.
Actually, by use of spectml analysis we will show that a transient signal can be
represented as a collection of periodic signals of different period. (In other words, it
can be represented by a spectrum of frequency components.)
instruments called seismometers (or seismographs), and then analyzing the structure
for the recorded motion. Such analytical studies play an important role in the design
of structures expected to undergo seismic vibrations.
.4 Acceleration [in/ s2]
.2
.0 ......~-A---.,I-+-J+l.f....II-\l-tWJ,I ......+4lf+-!-#-.,l"tMl>4,-IM"......,---......
-.2
-.4 L.-&. '-""-'- ~ .....
O. 2. 4. 6. 8. 10.
[s] Time
Figure 1.4: Surface accelerations recorded during an earthquake.
Imperial units. We use both sets of units in the examples and exercises in the following
chapters.
In the SI system, the basic unit of length is the meter (m), of mass the kilogram
(kg), of time the second (s). The unit of force is the Newton (N), defined as the force
that will produce an acceleration of 1 m/s 2 on a mass of 1 kg. An alternative set of
units for the Newton is therefore
kgm
N=--
S2
Multiples and submultiples in the metric system go in powers of 10 3 . Thus for stress,
there is the Pascal (Pa), the kiloPascal (1 kPa = 1000 Pa), the MegaPascal (1 M Pa =
106 Pa), and the GigaPascal (1 CPa = 109 Pa). Lengths are in multiples of the
meter as, for example, the millimeter (1 mm = 10- 3 m), and the micrometer (l/im =
10- 6 m).
In the Imperial system, the basic unit of length is the foot (It), of force the pound
(lb), of time the second (s). The unit of mass (called a slug) is defined as the mass
that will be accelerated at a rate of 1 ft / S2 when acted on by a force of lIb. Hence,
an alternate set of units for the slug is
lb.
It
S2
slug =
Inches (in) are often taken as the standard unit of length. Also, multiples may go in
thousands. For example, one thousand pounds is referred to as a kilopound or as a
kip, while a thousand pounds per square inch (1000 psi) is called a ksi.
The weight of a body is the gravitational force exerted on the mass. The value of
the gravitational constant is usually taken as
Thus, a mass of 1 kg weighs 9.81 N or 2.2Ib, whereas a mass of 1 slug weighs about
32lb or 143 N. The SI units are an absolute system of units because the fundamental
quantity of mass is independent of where it is measured. On the other hand, the
Imperial system is referred to as a gravitational system because the fundamental unit
of force (defined as the weight of a certain mass) varies with location on Earth.
The algorithms and source code presented do not have a built-in system of units
and do not utilize any dimensional conversion constants. Therefore, any consistent
system of units may be used for input and the corresponding calculated results will be
in the same units. For example, if the Young's modulus is specified as 10,000 ksi, the
reported values of force will be in kips. Since no particular system of units is preferred
in this book, we generally present the results in non-dimensional form. When needed
two useful conversion factors are 1 in = 25.4 mm and lIb = 4.448 N.
Some typical material values that can be used in the examples and exercises may
be taken as
Exercises 13
Exercises
1.1 Model a rocket as a simple structure with two nodes. What degrees of freedom
would you give each node? Why?
1.2 How many degrees of freedom are needed to describe the truss in Figure 1.1?
1.3 A two-story building has relatively rigid floors. If the wind is assumed as a
static load on one face only, how would you model the building? What degrees
of freedom would you give each node?
1.4 Suppose the building of the previous problem is modeled as a space frame with
nodes at each corner. What is the total number of degrees of freedom? How
does it compare to your previous model?
1.5 If the grid of Figure 1.1 is stiffened by adding diagonal members, what is the
increase in the number of degrees of freedom?
1.6 A rod (length = 100 in, area = 1 in 2 ) has a force of 1000 lb applied to it at one
end only. What is the variation of stress along the rod? How does this compare
to when the force is applied equally at both ends?
where L is the length. Show that the same frequency (in radians per second)
is obtained using either SI or Imperial units.
1.8 The fundamental frequency for a vibrating simply supported beam is
_ 11"2 [iii
w- L2Vp;.
where L is the length, A the area, and I the second moment of area. Show that
the same frequency is obtained using either SI or Imperial units.
Chapter 2
Rod Structures
A rod is a slender member that supports only axial loads. It is one of the simplest
components and therefore a very suitable vehicle for introducing the basic concepts
of the matrix analysis of structures. In this chapter, we consider only structures
composed of rods arranged longitudinally; the more general arrangement is left to
the chapter dealing with frames and trusses. We first review the basics of rod theory
and derive the governing differential equations. These are then used to obtain the
stiffness of a single rod element. The scheme for forming the structural stiffness
matrix is established by using equilibrium and compatibility conditions at each node.
Finally, the effects of boundary conditions are incorporated by simple row and column
reduction of the matrices.
The developments in this chapter act as a blueprint for the corresponding devel-
opments in most of the other chapters; we introduce the major themes that will recur
in the subsequent chapters.
14
2.1 Rod Theory 15
segment by
(u + ~u) - (u) ~u
t- --
- ~x - ~x
t = lim ~u = du
~X'-'O ~x dx
This differential relation tells us that, for example, a linear distribution of displace-
ment along the axis results in a constant distribution of strain.
Equilibrium
Let the distribution of the externally applied axial load per unit length be q(x). The
equilibrium of the small element shown in Figure 2.1 therefore leads (in the limit as
~x becomes very small) to
dF
F + ~F - F + q ~x = 0 or -
dx
=-q
where F is the resultant axial force. This says that the distribution of internal force
is constant if there is no applied loading.
q(x)
__11_
IL-E....;,.,A
I---> X u
~'-
u+~u
_
W F+~F
F
du
u = Et = E-
dx
where E is the Young's modulus. The axial resultants on the cross-section of area A
are
F =
J du '
udA = EA dx M = - uydA = -EA dx J
du ydA = 0 J
This last relation is true since the x-axis is taken to go through the centroid of the
section. The combination constant EA is called the axial stiffness.
16 Chapter 2. Rod Structures
Summary
The governing relationships for the structural quantities in the rod may now be sum-
marized as
Displacement: u = u(x)
Force: F = EA du
dx
(2.1)
d?u
Loading: q=-EA-2 (2.2)
dx
It is seen from these that the displacement (via its derivatives) is the quantity that
connects them together. That is, the displacement function u(x) can be viewed as the
fundamental unknown from which all the other quantities can be determined. The
significance of this point will become apparent later when we develop the displacement
method for solving structural problems.
Example 2.1: Find the displacement shape for a rod with a uniformly applied
load and fixed at both ends as shown in Figure 2.2.
u(x) ~_~_-=_
F(X)~
q( x) = qo = constant F====-----=:::::::::::-==---.--
Illilllll-----------llllli:illl
EA,L q(X)f ~
Figure 2.2: A fixed-fixed rod with uniformly distributed load.
Since the loading is constant, we will take this as the starting point. Thus
d2 u
EA-2 = -q = -qo = constant
dx
This is easily integrated to give
du
EA
dx
EAu =
where CI and C2 are constants of integration. To determine these constants, it is
necessary to impose some additional conditions. For the present problem, we know
that the boundary conditions are
at x = 0: u=O C2
This gives two equations for two unknowns, allowing us to solve for the coefficients.
They are
Cl = ~qoL,
EAu(x) = ~qox[L - x]
These functions are shown plotted in Figure 2.2. Notice that while the internal
(member) force goes from positive to negative, the displacement is always positive.
The maximum displacement occurs at x = L/2 and is
Example 2.2: Reconsider the last example, but this time using an applied load-
ing of a single concentrated force P at a distance a from the first boundary.
u( x) f-------------== .
F(x) t__--,
P f .....1 '
11!11!111---a---...J.IIIL.----L---a-4I!I!I!!~!l q( x) f------L.t _00 __
EA,L
It is possible to treat the rod as one complete piece of length L, but that will
require us to write the loading q as a function of x. However, we encounter a
singularity at the point of the applied force because the loading (force per unit
length) is infinite. There are ways of mathematically handling such behavior but
it is conceptually simpler to divide the rod into separate sections not including
the loading point. This approach is somewhat cumbersome because later we must
impose additional conditions when we recombine the sections back together. Thus
we can generate very many simultaneous equations. However, it will be seen that it
has the significant advantage of being adaptable to computer implementation.
18 Chapter 2. Rod Structures
Since the loading is zero for points between the boundaries and the applied load,
we will take the loading relation as the starting point, thus
0< x <a
d2u
EA 2
dx
-q =0
du
c1
dx
u(x) = c1x + C2
Similarly, for the second section
a< x <L
d2u
EA 2 -q = 0
dx
du
cj
dx
u(x) cjx + c2
To determine the four unknown coefficients C1, C2, ci, C2' we start by imposing the
boundary conditions
at x =0 : u = 0 = C2
at x =L : U = 0 = cjL + c2
This gives us two equations; hence two more conditions are required. We can obtain
one of these from the compatibility requirement at x = a: obviously, the displace-
ments are the same immediately to the left and right of the applied force, requiring
that
or
We obtain the final condition by considering the forces at the connection; consider
a small free body diagram extracted from around x = a, then
O<x<a a<x<L
PL a x PL x a
u(x) = -[1
EA
- -J(- J
L L
u(x) = -[1
EA
- -J(-J
L L
With these in hand, we can calculate all other quantities of interest (such as strains,
forces, and so on). For example, the force distributions are
O<x<a a<x<L
a a
F(x) = P[I -"L J F(x) = -P[-J
L
2.1 Rod Theory 19
Notice that the force distribution is constant in each section with a jump precisely
equal to the amount of the applied force..
To draw a connection between this example and the last, suppose the distributed
load is lumped into a single concentrated force of P = qoL at the center of the rod.
The maximum displacement then occurs at x = a = L/2 and is
PL qo L2
= 4EA = 4EA
U
max
F max = ~P = ~qoL
In comparison with the previous example, you will notice that the maximum forces
are the same, but the lumping procedure has increased the maximum displacement
by 100%.
_p
Example 2.3: Find the displacement shape for the rod structure with multiple
I
sections and applied loads, shown in Figure 2.4.
BA"L, BA"L,
Figure 2.4: Rod structure with multiple sections and applied loads.
Since the loading is zero for points between the boundaries and the applied load,
we will take the displacement solution as the starting point; thus for the first section
U(x) = CIX + C2
Similarly, for the second section
U(x) = cix + c2
Note that in this case we will use the variable x to mean the local distance along
the rod segment. To determine the four unknown coefficients CI, C2, cj, c2, we start
by imposing the boundary conditions. For the first segment we have
at x = 0: U = 0 = C2
whereas, at the end of the second segment we have
at x = L2 : F(L 2) = P = EA 2ci
This gives us two equations; hence two more conditions are required. We can obtain
one of these from the compatibility requirement at the junction of the two segments;
obviously, the displacements are the same immediately to the left and right of the
junction, requiring that
or
20 Chapter 2. Rod Structures
We obtain the final condition by considering the forces at the connection; consider
a small free body diagram extracted from around this region, then
or - EA 1 Cl + EA 2 ci = 0
These four equations can now be solved to give the undetermined coefficients as
P * PL 1
Cl = EA 1 ' C2 = 0; C2 = EA 1
The displacement distributions are now written for the two segments as
Px Px PL 1
1: u(x) = EA 1 2: u(x) = EA 2 + EA 1
Again, note that the variable x is the local distance along the rod segment.
The approach we used in the previous examples is applicable to general rod struc-
tures made of multiple segments. That is, we divide the rod structure into sections
by placing nodes (joints) where the loads are applied and where there are section
discontinuities. Each section is then integrated under the condition of zero applied
load. Finally, we determine the constants of integration by imposing the boundary
conditions, compatibility between each section, and equilibrium at each joint.
While this is straightforward enough, it is seen that once the number of sections
exceeds three or four, the simultaneous equations become unwieldy. To make the
approach workable, it is necessary to marshall the equations in such a way that the
treatment of each section is highly repetitious and many patterns recur. We will do
this by considering a typical section under the action of typical applied loads. Once
the relations for this case are established, then the relations for other sections will
follow immediately as special cases. We will incorporate a significant refinement in
this by arranging the equations so that compatibility between sections is automatically
assured.
F(x) ......
3.....-
...... UI U2
---.... ....JI-- ~"L
FI ...... X F2
CD Q)
taken at the first end, say, then the forces acting on it are FI and F(O). Equilibrium
(and the sign convention for the internal forces) requires that
FI + F(O) = 0 or FI = -F(O)
This procedure is set up so that we will emphasize the difference between the external
and internal forces. It is important to realize that FI and F2 (which we will call
element nodal forces) are to be viewed as external applied forces different from the
internal force distribution F(x).
If there are no loads applied in between the two nodes (FI and F2 are applied at
the ends), then we may start by integrating the rod equation
<flu
EA dx 2 = -q = 0
We find by integrating twice that the displacement in the element is a linear function
of x:
u(x) = ao + alX
where ao and al are constants. The relation between the displacements at each node
and the coefficients is obtained by imposing
u(O) = UI = ao
u(L) = U2 = ao + alL
This gives
ao = UI
al = -(ul-u2)/L
allowing us to write the displacement distribution in terms of the nodal values as
(2.3)
The known functions fl(x) and h(x) are called the rod shape functions.
As already said, if the displacements are known, then everything else about the
problem is also known. This is now modified to read that if the nodal displacements
22 Chapter 2. Rod Structures
are known then everything else is known. For example, we can easily obtain the
member force as
du EA
F(x) = EA dx = T( -Ut + U2)
By equilibrium of a free body diagram near each node, the member and element nodal
forces are related by
F t = -F(O), F2 = F(L)
Consequently, the nodal forces are related to the nodal displacements by
(2.4)
or symbolically,
{F} = [ k ){ u}
The symmetric matrix
is called the stiffness matrix for the rod element, {F} is called the vector of element
nodal forces, and {u} the vector of nodal displacements. The latter is often referred
to as the nodal degrees of freedom.
It is noted that [ k 1is singular, i.e., det[ k 1= 0; and consequently, its inverse does
not exist. This implies that given an arbitrary force vector {F}, it is not possible to
find a unique solution for { u}. A singular stiffness matrix indicates that the structure
is not stable. In the present case, it means that the element will translate indefinitely
and boundary conditions must be specified before there is a solution.
where the superscript designates the element by reference to its nodes. It is neces-
sary for us to devise a scheme for attaching these elements together subject to the
requirement that the resulting assembled structure not violate any compatibility or
equilibrium requirements.
<D Q)
+
...-
+
...1---..
p(12)
....: ...-1-+
F.(12) p(23) F.(23)
-PI -
P2
~
P3
Figure 2.6: Simple rod structure with two members and three nodes.
Nodal Equilibrium
Given that the structure is in a state of equilibrium, it follows that each node (or
joint) must also be in a state of equilibrium. Take a free body diagram containing
Node 1 as shown in the figure. Notice that we are referring to two types of nodes: one
represents the joint and one the end of the element. This device is very useful because
it allows us to draw the distinction between the loads applied to the structure and
the element nodal loads. The applied force is denoted by PI, while the nodal force
acting on Member 1-2 is denoted by FP2). For equilibrium at this node, we must
have
or P _ F(12)I - I
P2 = FP2) + FP3)
P3 = FP3)
This is the statement of equilibrium of each node. (It is also, consequently, a statement
of the equilibrium of the structure as a whole.) There are as many equations as there
are nodes, and as many vectors on the right hand side as there are members.
24 Chapter 2. Rod Structures
Assembly
We can relate the two vectors on the right-hand side of the equilibrium relation
to the nodal displacements {UI, .. ,U3} through augmented [3 x 3) element stiffness
matrices as follows
0 } (23) [ 0 0 ]{ UI }
{ {F} = 0 [k(23) {U}(23)
Substituting these into the nodal equilibrium conditions leads us to
{
;: } =
P3
[:il:; k~~2~~2~g3)
0 k(23)
21
ki~3)] ~:
k(23)
22
{
U3
}
Example 2.4: A three element (four noded) rod is loaded as shown in Figure 2.7.
Set up the structural stiffness matrix and the system of equations to be solved.
To construct the structural stiffness matrix, it is necessary for us to find the
contributions to [ K ] from each member. Thus, with reference to the figure, we first
state each element stiffness matrix, that is
11I11t +~---"L=---t!I!I!!1
~ ~ @ @
Figure 2.7: Four noded example.
We now augment each matrix to the size of the number of unknowns. In this case
the size is [4 x 4J corresponding to the degrees of freedom {Ul, U2, U3, U4}' The three
augmented stiffness matrices are, respectively,
3 -3 00] 00 0]
EA -3 3 0 0 EA 00 01 0
-1 0]
0 EA 0 0 oo 0
L [o 0 0
0
0 0
0 0
' L [ 0 -1
o 0
1
0
0 '
0
L
[ 0 0 2
o 0 -2
-2
2
We form the structural stiffness matrix by simply collecting the entries in the aug-
mented element stiffness matrices according to their position. The result is
3 -3 0 0] [3 -3 0 0]
[K J = EA -3 3 + 1 -1 0 = EA -3 4 -1 0
L [ 0 -1 1 + 2 -2 L 0 -1 3 -2
o 0 -2 2 0 0 -2 2
Note that this matrix must be symmetric as indeed it turns out to be. Also note
that it is banded, that is, the non-zero entries are located close to the main diagonal.
This is a characteristic of structural systems, and we will say more about it later
when considering the computer aspects.
The system of equations for us to solve can now be established as
EA 3 -3 0 0] {
-3 4 -1 0
Ul
[
U2
L 0 -1 3 -2 U3
o 0 -2 2 U4
This system of equations cannot be solved directly because the right hand side con-
tains some unknowns. That this occurred is associated with the boundary conditions
and how that information is implemented. In order to proceed, it is necessary for
us to rearrange these equations, and we will treat this in a general way in the next
section.
It is worth noting in this last example that the knowns and unknowns in the
final set of equations form mutually exclusive groups. That is, if a displacement is
known at a node, then the corresponding applied load is an unknown (we refer to
this as an unknown reaction). Conversely, if the applied load is known at a node,
then the corresponding nodal displacement is unknown (we refer to this as a degree
of freedom). Understanding this exclusivity is one of the keys to understanding the
matrix approach to structural mechanics.
26 Chapter 2. Rod Structures
where
are the unknown and known nodal displacement and force vectors. The equations for
the whole structural system are expressed in the partitioned form as
We can therefore obtain the unknown displacements from the first equation by solving
(2.5)
since this has only known quantities on the right hand side. The unknown external
loads are obtained from the subsequent computation using the second equation above
(2.6)
where now, we know everything on the right hand side since {uu} was obtained from
the previous calculation.
2.4 Boundary Conditions 27
Special Case
If the known nodal displacements are zero (such as for fixed boundary conditions)
then the above take on the particularly simple form after substituting {Uk} = 0
(2.7)
For convenience, we define [K*] == [Kuu ]; the system of equations for us to solve
becomes simply
where [K*] is called the reduced structural stiffness matrix relating the unknown dis-
placements to the given external loads.
This method of rearranging the equations reduces the size of the system to be
solved. That is, [K*] is smaller than [K] because {uu} is smaller than {u}. What
we have done is posed the problem in terms of the minimum number of unknowns;
this involves the unknown displacements (degrees of freedom) but not the unknown
applied forces (reactions). In the general structural principles to be developed in
Chapters 6&12, we refer to {uu} as the generalized coordinates and the equations
{ud = 0 as constraint equations.
Example 2.5: Complete the solution of the three element rod structure of the
example in the last section.
The nodal displacements at Nodes 1 and 4 are zero, leaving the unknown nodal
displacements as
{uu} = {:~}
The corresponding known applied loads are
We obtain the rearranged structural system by interchanging both rows and columns
EA =
L -3 0 3 0 UI = 0 PI = ?
o -2 0 2 U4 = 0 P4 = ?
The reduced stiffness matrix is just the first [2 x 2J entries giving us the reduced
system of equations as
EA
L
[4-1 -1]
3
{'U ~
2 } ={ 0}
P
Once the nodal displacements are obtained, the loads in the members can be
calculated.
This scheme can be easily adapted to a computer. To show the approach laid out
somewhat in the pedantic manner followed by the computer, we will redo the previous
problem.
Example 2.6: Assemble the reduced stiffness matrix for the rod structure of
Figure 2.7.
Step 1: There are 4 nodes, hence the system size is
Step 2: There are fixed boundaries at Nodes 1 and 4, hence the reduced system is
Step 4: The reduced element stiffness is (after retaining rows and columns asso-
ciated only with Node 2)
(k*(12)] = EA[ 3 ]
L
Step 5: Add to the currently unpopulated reduced structural stiffness
[K*] = EA
L
[30 0]0
Step 6: For this step just repeat Steps 3,4 and 5 for each member. This will be
done here in one swoop each.
Steps 3,4,5:
L -1 ~1 ]
= EA
[k(23)] [ 1
L -1 ~1 ]
= EA
[k*(23)] [ 1
L -1 ~1 ]
= EA
[K*] [ 4
Steps 3,4,5:
[k(34)] 2EA [ 1
L -1 ~1 ]
[k*(34)] = EA[2]
L
[K*] = EA [ 4
L -1 ~1 ]
This last array is the desired reduced structural stiffness matrix.
It is worth mentioning that the way this is actually programmed is not by 'scratch-
ing' rows and columns. Rather, a pointer array is maintained whose content is an
identifier for each node (degree of freedom) and this identifier is a zero when the node
is fixed. A refinement is to also use the array to keep track of the equation numbers;
thus for the above problem the array is
advance the precise output desired. That is, the results of the solution are usually
investigated interactively. The post-processing arrangement lends itself to this ap-
proach, and so we will assume that the computation of member quantities needs to
be determined in this spirit.
Once we know the nodal displacements, the distribution of displacement is ob-
tained simply from
This shows that the distribution is linear and therefore the plot of displacement is a
straight line connecting the nodal values of Ut and U2 for each element.
Determining the member load distribution is only slightly more involved. Recall
that the following equation relates the member force to the global displacement
{F} = [ k ]{ u}
Thus the application of this to each member gives the nodal forces for each member.
The internal (member) load in the rod member is related to the nodal forces by
Example 2.7: Plot the displacement and member force distributions for the three
element rod structure of Figure 2.7.
U(Xl~
F(Xlt a I r
Figure 2.8: Displacement and force distributions.
The displacement distribution is obtained by plotting these four points and connect-
ing them by straight lines.
The member force distribution is obtained by first determining the nodal forces.
Hence, using the stiffness relation for each member we get
3EA
L
[1-1 -1]
1
{O} ~ = {-3} P
1 11EA 3 11
These values of nodal force are related to the internal member forces by the previ-
ously derived equilibrium conditions. Thus for the first member
We know that the distribution of force is constant, and these relations show that we
have a choice as to which node we use in order to determine its value. Proceeding
in like manner for the other elements, we get
and
F(O)(34) = -Ff34) = -8P, F(2L)(34) = +FJ34) = -8P
11 11
Note that the discontinuity in the plot at Node 3 corresponds precisely to the amount
of the applied load of P.
Boundary Reactions
As pointed out before, the unknown forces (which are usually the boundary reactions)
can be obtained from the calculated displacements as
(2.8)
32 Chapter 2. Rod Structures
'-- --3
..' ....
FJl)
-P
F1 2 )
--e'O:. -'
For fixed boundaries, we will generally take that the applied load P is zero, but the
above allows for other cases also.
Example 2.8: Find the boundary reactions for the three element rod structure
of Figure 2.7.
Following on from the last example, we can use the stiffness relation for each
member attached to a boundary of interest to get
3EA [ 1
L -1
-1]{O}
1
PL {-3}P
1 IlEA = 3 11
2EA [ 1
L -1
-11 ]{4}
IlEA
PL {8} P
= -8 11
These values of nodal force are related to the reactions at Nodes 1 and 4 by
The plus sign in both cases indicates that the action is tending to move the bound-
aries from left to right. Thus, the reaction at Node 1 is tensile, while that at Node
2 is compressive.
<flu
EA dx 2 = -q(x)
We find by integrating twice that the displacement in the element is no longer a linear
function of x but given by:
<flw
u(x) = ao + alX + w(x), EA dx 2 == -q(x)
2.6 Distributed Loads 33
where ao and at are constants. The relation between the displacements at each node
and the coefficien ts is obtained by imposing
u(O) = Ut= ao + Wt
u(L) = uz=ao+atL+wz
where we have defined Wt == w(O) and Wz == w(L). This gives
ao = Wt
Ut -
at = -(Ut-Wt-uz+wz)/L
allowing us to write the displacement distribution in terms of the nodal values as
x x x x
u(x) = (1- -)Ut + -Uz - (1- -)Wt - -Wz + w(x) (2.9)
L L L L
As already stated, if the displacements are known, then everything else about the
problem is also known. We can now obtain the member forces as
du
F(x) = EA-
dx
= -EA
( - U t + uz) -
L
EA
-(-Wt
L
dw
+ wz) + EA-
dx
By equilibrium of a free body diagram near each node, the member and nodal forces
are related by
F t = -F(O), Fz = F(L)
Consequently, the nodal forces are related to the nodal displacements by
-l]{Ut}_EA[ 1 (2.10)
1 Uz L-1
The remaining step is to substitute the applied loading q( x) into the expressions
for w( x). The results are simplified considerably if we note that
(2.11)
where the known functions It (x) and fz( x) are the rod shape functions of Equation(2.3).
From this, we see that the distributed load is replaced by concentrated loads associ-
ated with Nodes 1 and 2. Therefore, the assembled equations for the structure are
given as
{P+Q} = [K]{u}, {Q} == Lm J
{fq}m dx
34 Chapter 2. Rod Structures
where {p} are the applied concentrated loads and {Q} are the assembled distributed
loads treated as concentrated applied loads.
These results look deceptively similar to those of the point loads only case; all
we have done is replaced the distributed load with concentrated loads at the nodes.
Keep in mind, however, that the displacement functions are quite different in the two
cases; thus for a given set of nodal displacements both would give different computed
values for the member loads.
Example 2.9: Use matrix methods to find the axial force and displacements in
the uniformly loaded rod of Figure 2.2.
IIll1it'-------lI~-tllllll
<D ~
We will model the rod with two elements (three nodes). The global degrees of
freedom are then
{tI} = {til, tl2, tl3}
The stiffness relation for element 1-2 is
r'
EA[
L/2 1 tl2 Jof {(1-X2/L)} /
(x2/L)
d
qo x
EA2 [1 -1]
=
L -1 1 {til}
tl2 _qoL4 {11}
Similarly, the stiffness relation for element 2-3 is
{FFl }(23)
2
= EA2
L
[1-1 -1]
1 {til}
tl2 _qoL4 { 11}
{ ~ } = :2 [~1
The last vector in the above equation is the vector {Q}.
The boundary conditions impose that the degrees of freedom are zero at Nodes
1 and 3; further, there are no applied concentrated loads, hence the reduced system
becomes
Problems 35
EA2 [ 1
L -1
-11 ]{O} q L2 _ 'loL {I} = qoL {-I}
1 SEA
O
4 1 2 0
EA2 [ 1 -1 ] { 1 } qoL2 _ qoL { 1 } = qoL { 0 }
L -1 1 0 SEA 4 1 2-1
These values of nodal force are related to the internal member forces by the previ-
ously derived equilibrium conditions. Thus it is easy to show that
Problems
2.1 A uniform rod of length L is hung vertically under the action of gravity. Show
that the loading per unit length is q(x) = pAg, where p is the mass density
and 9 is the gravitational constant. Consequently, show that the displacement
distribution is
pAL x2
u(x)=-gx(l--)
EA 2L
2.2 Model the previous self weight problem using 3 elements. Account for the
distributed mass by 'lumping' it at the nodes. Compare the distribution of
displacement with the exact values.
2.3 Show that when the cross-sectional area or the elastic modulus changes along
the length of the rod that the relevant equations are
du d du
F(x) = EA dx ' q(x) = --[EA-J
dx dx
Consider a rod of length L that has a varying area of the form
x
A(x) = At + (A 2 - Ad!;
If this is fixed at one end and a load of P applied at the other, show that the
displacement function is
2.4 Model the previous problem using three elements of different area but each
of the same length. (Assume that At = 2A 2 .) Compare the displacement
distributions.
36 Chapter 2. Rod Structures
2.5 When a rod is heated uniformly it expands by the amount t1L = ext1TL, where
ex is the coefficient of thermal expansion, and {).T is the change in temperature.
Hooke's law for the material is modified as
u = E(f - ext1T)
Show that when temperature is taken into account, that the element stiffness
relation becomes
2.6 What modifications (if any) are needed during the assemblage stage to account
for temperature effects? If the second term on the right hand side of the
previous stiffness relation is viewed as a force vector { q}, show how the global
equilibrium equations are affected.
[Reference [35], pp. 131]
Exercises
2.1 Rework Example 2.9 but using only one element. Show that the same values
of boundary reaction are obtained.
2.2 The pendulum of a clock has a 3/b weight suspended by three parallel rods of
30 in length. Two of the rods are brass (E = 15 msi diam = .10 in) and the
third steel (diam = .05 in). What is the force in each rod?
[F. = 0.61b H = 1.21b)
2.3 A square reinforced pier 1 It X 1 It in cross-section and 4 It high is loaded
axially by 150 kip. The concrete is stiffened with eight 1 in 2 steel reinforcing
bars placed symmetrically about the vertical axis. What percentage of the force
is supported by the steel? [33%]
2.4 Show that the governing equations for the twisting of a shaft are
d d d
T(x)=GJ
dx
' q(x) = --[GJ-)
dx dx
where T is the torque and is the angle of twist.
2.5 Following the procedure for a rod, show that the shaft element shape functions
are the same as for the rod. Consequently, show that the stiffness relation for
a shaft element can be derived as
Tl }
{ T2
= GJ
L
[ 1-1] {
-1 1
1 }
2
Chapter 3
Beam Structures
A beam is a slender structural member designed to carry transverse loads and applied
couples. In response to these loads, it develops internal bending moments and shear
forces. We shall refer to a beam structure as a collection of beams arranged in a
collinear manner. These are sometimes referred to as continuous beams. Additional
aspects of the mechanics of beams can be found in strength of materials books such
as Reference [131.
This chapter develops the governing differential equations of beam theory so as to
set the basis for the introduction of the beam element. This is done by emphasizing
the application of the twin principles of compatibility and equilibrium at the junction
of each beam connection. In this way the matrix assemblage procedure is seen as just
an application of these principles. Consequently, the matrix description of these beam
structures follows very closely that of the rod structures of the previous chapter. We
do, however, introduce new procedures for handling elastic boundaries.
4>(x) = ~v ~ dv
~x dx
37
38 Chapter 3. Beam Structures
q(x)
C(JI S~M
V ~x V+~V
Since the beam deforms (locally) into the arc of a circle, then the deformed length of
a small line segment a distance y from the axis is (R - y)6.</>, where R is the radius
of curvature. The strain of the segment is therefore
(R - y)6.</> - R6.</> -y
i = R6.</> = if
This shows that the strain is linearly distributed on the section. From geometry, it
is known that radius and curvature are related by 1/R = tflv/dx 2 , hence we can also
obtain for the strain
tflv d</>
i= - y -2 = - y -
dx dx
This gives the relation between the strain (i), slope (</, and curvature (tflv/dx 2 ).
Equilibrium
The equilibrium of the small element shown in Figure 3.1 leads to the following two
equations (in the limit as 6.x becomes very small)
dV
dx +q=O, dM + V =0
dx
where V is the resultant shear force, M the resultant moment, and q the distributed
load per unit length. In contrast to the rod, the beam has two equilibrium equations.
One consequence of this is that the variety of beam behavior will be much greater
than for the rod.
Stress Resultants
Applying Hooke's law to the axial stress state of the beam gives
d</> tflv
a = Ei = -y E- = -y E -2
dx dx
3.1 Beam Theory 39
Note that the stress, too, is distributed linearly on the section. Knowing the explicit
form for this distribution allows us to determine the following stress resultants on the
cross-section
For a rectangular section (b x h) and a circular section (diameter D), this gives
I circ -- 641 D4
respectively. Similar expressions exist for other sections.
Summary
All the relationships for the structural quantities may now be collected as
Displacement v = v(x)
Slope 4> = dv (3.1)
dx
Jlv
Moment M = +EI dx 2 (3.2)
It is seen from these that the deflected shape v(x) can be viewed as the fundamental
unknown of interest; all other quantities are obtained by differentiation. This is pre-
cisely the same conclusion already drawn from the rod analysis. It is also interesting
to note that the section properties are reduced to the single combination term EI.
This is called the flexural stiffness.
When solving beam problems, we may be given information at any of the five
levels above, and have to carry out integrations (or differentiations) to obtain the
other quantities. Integration gives rise to constants of integration which must be
found from the boundary and compatibility conditions. In the general case there
are four constants of integration, twice as many as for the rod. Thus twice as many
conditions must be imposed at each section.
Example 3.1: Find the deflected shape of the fixed-pinned beam shown in Fig-
ure 3.2. The applied load per unit length, w o , is uniformly distributed.
40 Chapter 3. Beam Structures
v(x) .......
J(x)
M( x) ~-"....:.::-
. . _---':::=.......
V(x)p~'
Figure 3.2: Uniformly loaded beam.
The loading is constant and given as q(x) = -w(x) = -Wo; therefore we will
take this as the starting point. That is,
d4 v
E I -4 =q= -Wo
dx
Integrate to obtain
d3v
El =
dx 3
d2v
El
dx 2
=
dv
El
dx
=
Elv
We will now impose the boundary conditions. At the fixed end, both the deflection
and rotation (slope) are constrained to be zero. At the pinned end, only the de-
flection is constrained to be zero; in other words, it is free to rotate, which in turn
means there is no restraining moment. The four conditions are expressed as
at x =0 : v = 0, dv =0
dx
d2v
at x = L: v = 0, M = El dx 2 = 0
These give, respectively,
o = C4 o = C3
o = -i4woL4 + ~cIL3 + !c2L2 + C3L + C4 o = -!woL2 + clL + c2
After solving for the coefficients, we find the deflected shape to be
3.1 Beam Theory 41
We can obtain the other quantities of the solution by differentiation. This gives the
slope, moment, and shear distributions as
Both the fixed-fixed and pinned-pinned cases have symmetric distributions. The
maximum deflection of the fixed-pinned case lies between the maximum for these
two.
Example 3.2: Find the deflected shape of the fixed-fixed beam shown in Fig-
ure 3.3. The concentrated applied load is located a distance a from the end.
t
M(1) P M(2)
/1!/!I!If---:_tP _ ( tot)
L
- a_II!!!111
EI,L V(l) V(2)
Because of the change in loading due to the concentrated force, we will consider
the beam to be made of two sections; compatibility and equilibrium conditions will
then be imposed to match them at the connection. In the present problem, the
loading on each section is zero and so it is convenient for us to start there. That is,
for the first section
We will first impose the boundary conditions of zero deflection and slope at the near
end. That is,
at x = 0: v = 0, dv = 0
dx
42 Chapter 3. Beam Structures
In analyzing the second section, we will use the variable x as a local variable,
that is, it ranges over 0::; x ::; (L - a). With that in mind, we have
d4V
EI = q(x) = 0
dx 4
EIv ~cix3 + ~C2X2 + cjx + c~
The boundary conditions for this section are zero deflection and slope at the end.
That is
at x = L - a: v = 0, dv = 0
dx
Hence, we obtain that
lc*(L
c*4-- 3 1 - a)3 + lc*(L
22 - a)2
The four boundary conditions were insufficient to determine all eight of the
unknown coefficients of integration. It is necessary for us to impose compatibility
across the junction of the two sections. That is, at the junction
dV{I) dv(2)
dx = dx
We still need more equations and to obtain them we will consider the equilibrium
conditions in the vicinity of the junction. To this end, isolate a small segment near
the applied load as shown in Figure 3.3, the equilibrium conditions give
y{I) = y(2) + P
The two compatibility and two equilibrium conditions become (when we rewrite
them in terms of the unknown coefficients)
C2 = ~Pa, ci = ~P,
3.2 Beam Element Stiffness Matrix 43
The beam element stiffness matrix will be derived by directly integrating the
governing differential equations for the beam. Since the loading between nodes is
zero, we have
d4v
El 4 = 0
dx
giving the general solution for the deflection curve as
v(x) = ao + atX + a2x2 + a3x3
where ao, at, a2, and a3 are constants. By using the following end conditions
v(O) = Vt ,
dv(O) = 4>t
dx
v(L) = V2,
dv(L) = 4>2
dx
44 Chapter 3. Beam Structures
we can rewrite the constants in terms of the nodal displacements VI and V2, and the
nodal rotations <PI and <P2. Specifically, they are
ao = VI
al = <PI
3 2 3 1
a2 = --VI -
2
-<PI
L
+ -V2
L2
- -<P2
L
2 1 2 1
a3 = L3 VI + 2 <PI - L3 V2 + 2 <P2
Substitution of these into the expression for the deflection leads us to
The functions gn (x) are called the beam shape functions. The complete description
of the element is captured in the four nodal degrees of freedom VI, <PI, V2, and <P2
(since the shape functions are known explicitly). If, in any problem, these can be
determined, then the solution has been obtained.
The slope, moment and shear force are obtained (in terms of the nodal degrees of
freedom) simply by differentiation. For example, the moment distribution is
By considering a free body diagram of each end of the element, it is easy for us to
establish the following relations among the member loads and the nodal values
VI = -V(O), M I = -M(O)
V; = +V(L), M2 = +M(L) (3.7)
Therefore, the nodal moments are easily related to the nodal degrees of freedom. By
carrying out the indicated differentiations we can also relate the nodal forces to the
nodal degrees of freedom. To summarize the above results in matrix form, we define
3.2 Beam Element Stiffness Matrix 45
the element nodal loads vector and the corresponding nodal displacements vector as
the following column matrices
6L -12
V1I
M } EI [ 12
6L 4L 2 -6L 2L6L2 ] { </>
VI }
where [ k 1is the beam element stiffness matrix. Note that this stiffness matrix is
symmetric. Also note that the nodal loads satisfy the equilibrium conditions for the
free body diagram of the element. This is as it should be, since the relation was
obtained by integrating the differential form of equilibrium.
The beam element stiffness (unlike that for the rod) contains dimensional quan-
tities inside the brackets. This comes about because the vector terms have mixed
quantities; that is, rotation is non-dimensional but deflection has the units of length.
An interesting alternative form of the above relation is
-12
V1}
MtiL EI [126 64 -6
{ V2 = 3 -12-6 12
MdL 6 2 -6
This gives dimension of force for all load terms and dimension of length for all degrees
of freedom, and makes the terms inside the matrix dimensionless. As a result the
relation is dimensionally similar to that for rods. While this form has a certain appeal,
it is unsuited for our purpose - we wish to assemble many elements of (possibly)
different lengths to form a structure and this task is made easier by having the vector
of degrees of freedom common from element to element. In the above form, these
vectors are element dependent since they contain the element length.
When the loads or displacement vectors contain mixed terms (e.g., force and
moment) they are often referred to as generalized vectors. In Chapter 6, we will
show that it is possible to establish quite general results in terms of generalized loads
and displacements without having to specify them explicitly. Therefore, using mixed
load vectors at this stage does not lead to any difficulty when we move onto more
complicated structures.
46 Chapter 3. Beam Structures
CD
.OL-
~
+------,.
Pz = Z
+ \/,(Z3)
\1.(12)
1
Tz = MJ12) + M1 Z3 )
P3 VP3)
T3 MJZ3)
Note that the total number of equilibrium equations (6) is equal to the number of
nodes (3) times the number of equilibrium equations at each node (2), and that this
coincides with the total number of degrees of freedom of the system (6). We can
express the above equations in matrix form as:
(IZ) (Z3)
PI VI 0
T1 M1 0
Pz l'2 VI + {F}(Z3)
Tz = M z + M1
== {F}(lZ)
P3 0 l'2
T3 0 M z
3.3 Structural Stiffness Matrix 47
This is arranged so that there are as many column vectors on the right hand side
as there are members. We now replace each element nodal load vector by use of the
element stiffness relation augmented to the system size. That is,
{p} = [K]{u}
where
{p} == {u} ==
where [K ] is the structural stiffness matrix of the system, and is obviously the sim-
ple superposition of the respective element stiffnesses augmented to the full structural
size. This is precisely the same result that occurred when we analyzed the rod struc-
ture. In fact, the pattern that is emerging when forming the structural stiffness matrix
is one that is true for all linear elastic structures.
Vl=cPl=V3=cP3=O
48 Chapter 3. Beam Structures
lillllr--E-l-h-L-l----f tl l l!
CD @
and therefore the reduced system has only the following two free degrees of freedom
{IDbc} = {0,0;1,2jO,0}
The two unknown nodal displacements and the given loads are, respectively,
{Pd = { ~ } = { ~}
In order to form the stiffness matrix, we must first form the element stiffness ma-
trix and then 'scratch' the rows and columns corresponding to the zero degrees of
freedom. The contributing reduced element stiffness matrices for our problem are
determined to be
-6L 1 ]
4Lr '
The assembled reduced structural stiffness matrix is therefore
As a special case, consider when each segment is the same, that is, E 1 = E 2 = E,
h = h = I, L 1 = L 2 = L. Then
= EI
[K*]
3
[24 0]
82
Note how the off-diagonal terms are zero, indicating uncoupling of the degrees of
freedom. We can now obtain the inverse of this matrix by forming the reciprocal of
the diagonal terms. This gives
The rotation at Node 2 is zero showing that there is symmetry of the deflection.
Both results are identical to that obtained in Example 3.2. This is a reminder that
the matrix method can give the exact result.
Example 3.4: A uniform beam is pinned at two points and attached to a vertical
roller as shown in Figure 3.7. Set up the reduced structural stiffness matrix for this
problem.
~lil ilil
Figure 3.7: Three node problem.
We start by dividing the beam into two elements with three nodes, and number
them as shown. The total system degrees of freedom are
giving a total system size of [6 x 6]. The boundary conditions (or geometric con-
straints) require that
VI = 0,
{IDbc} = {0,1;0,2;3,0}
- 6L 2]
12
~~
4L2 2L2 + + ] = ~~ [
+ 2L2
4L2 2L2 ]
+
[K*] = 2L2 4L2 0- 12L/8 2L2 6L2 -3L/2
[
0- 12L/8 0+ 12/8 -3L/2 3/2
50 Chapter 3. Beam Structures
These values are used in Examples 3.6 and 3.7 to find the load distributions and
reactions.
1 l
Pl p2
q(x)
l
t U I
lumped
CD @
1 PlTl
consistent
t P2T2
where ao, all a2, and a3 are constants. By using the following end conditions
dv(O) = 1>.
v(O) = v. ,
dx
dv(L)
dx
we can rewrite the constants in terms of the nodal displacements. Substitution of
these into the expression for the deflection leads us to
where the functions 9n(X) are the beam shape functions of Equation(3.6) and w. =
w(O), w~ = dw(O)/dx, and so on.
The moment and shear force distributions are obtained (in terms of the nodal
degrees of freedom) by differentiation. After we evaluate the end values, we can
arrange them in matrix form as
Vi}
M. EI [12
-12
6L -6L
6L 4L2
{ MV2 2
= L3 -12 -6L 12
6L 22 -6L
The final step is to replace the function w( x) by the distributed load function q(x).
We can show by integration by parts that the following relation holds true
3w
fL
Jo 9. (x )q( x) dx = kn w. + k. w~ + k. 3W2 + k 14 w; -
2
E I dd 3
x
Similar results can be obtained for the products of q( x) with the other beam shape
functions. They allow us to simplify the stiffness relation to
-12
V. }
M. EI [12 6L -6L
6L 4L2
(3.10)
{ V2 = D -12 -6L 12
M2 6L 2L2 -6L
We have thus succeeded in replacing the distributed load by equivalent loads associ-
ated with the nodes and given by
p. = 1 q(X)9.(x)dx,
L
P2 = 1 q(X)93(X) dx
L
T. = L 1 q(X)92(X) dx,
L
T2 = L 1L
q(X)94(X) dx (3.11)
52 Chapter 3. Beam Structures
These loads are referred to as initial stress terms; that is, they are loads in addition
to the nodal loads {F}. The initial stress loads are also called the consistent load
representations because they involve displacement shape functions consistent with
the stiffness matrix formulation. The assembled system of equations take the form
{P+Q} = [J{]{u}
where {p} is the collection of applied nodal loads and {Q} is the assembled form of
the initial stress loads. We see from this that the distributed load can be accounted
for by using a set of equivalent nodal loads.
The consistent loads are a statically equivalent load system. This can be shown
in general, but we demonstrate it with the following special cases. For a uniform
distributed load qo, for example,
PI = ~qoL, P2 = ho L
T I = +hqoL2, T2 = --hqo2
The resultant vertical force is PI + P2 = qoL which is in agreement with the total
distributed load. The resultant moment about the first node is TI + T2+ P2L = ~qo2
again in agreement with the value from the distributed load. The corresponding
values for a linear distribution of load with a maximum of q(L) = qm are
PI = foqmL,
T I = ioqm L2 ,
The resultant force and moment are hmL and ~qmL2, respectively. Notice that the
consistent loads have moments even though the applied distributed load does not.
the sandbags to their nearest node. That is, we have lumped the total load into two
statically equivalent concentrated loads at the nodes. The mathematical expressions
for this are
Jor r q(x)dx
L 2 L
/
PI = q(x)dx, P2 = (3.12)
JL/2
This is called the lumped load approximation and is essentially a scheme for replacing
the actual distribution in terms of its average. It is evident that when q(x) is highly
irregular, many nodes are required to yield an accurate representation.
For a uniform distributed load qo, for example, the lumped load representations
are
That is, half of the total load is placed at each node. Whereas, if the distribution is
linear with a maximum of q(L) = qm, the corresponding lumped loads are
In comparison to the consistent load representation, we see that the lumped repre-
sentation does not contain any end moments.
CD (2)
./,L
q(x) = -W o
:L :::::::::::::::~::::
M(x) = -~ox (x - L)
-5woL4 -W L2
V
max
= 384/ ' M max = --8-
O
Note that the applied moments at the center cancel each other, and that the forces
at the ends do not cause deflections. Therefore, the equivalent problem is that of a
simply supported beam under the action of a concentrated force -~woL acting at
the center and two concentrated end moments fgw oL2 acting in opposite directions.
Because of symmetry we actually need model only one half of the beam. In
doing so, however, we must retain the T 2 loading for the first element even though
it cancels when assembled with the second element. Imposing the constraints that
VI = 0 and </12 = 0 we get the reduced stiffness relation as
-12
-woL/4 }
VI }
M 1 8El [12
=V
3L 3L
L2 -3L 12
3L ] { </>1O} { -W-wL2/48
O
{ V2
M2
-12
3L
-3L
~L2
12
-3L
~~~ ~2- L/4o
+wo L2/48
the nodal loads since the initial stress term is neglected in the element stiffness
relation. Therefore, we can write from above
{
~I
V2
} = woL
48
{ ~i }
-12
M2 7L
The moment at the center of the beam is in error by about 2%. A non-zero value
of shear force is also found at that location.
Lumped Load
Realizing that the length of each element is L/2, then the equivalent lumped loads
acting on the nodes of each element are
1- 2: PI = -~woL P2 = -~woL
2 - 3: PI = -~woL P2 = -~woL
Note that the loads at Nodes 1 and 3 are taken directly by the supports and do not
contribute to the deflection. Therefore the equivalent problem is that of a simply
supported beam under the action of a single concentrated force -~woL acting at
the center. The reduced system of equations are
-12
MVI } 8EI [12
3L 3L
L2 -3L
I
V
{ 2 = 13 -12 -3L 12
M2 3L ~L2 -3L
fM>i C) Ti
When a beam rests on a system of springs at some of the nodes, constraints are
exerted on the beam system. For a linear spring and a rotational spring, for example,
we have reactions of
P = -av, T = -f3>
where a is the linear spring constant and f3 is the rotational spring constant. Consider
Node i, as shown in Figure 3.10. The resultant applied linear force is (Pi - aivi),
that is, for a positive deflection the spring exerts a restoring force. Similarly, for the
torsional spring, the resultant applied moment is (Ti - f3d>i). Therefore, during the
assemblage process, the equilibrium equations are modified to become
[K]{u}=
The right hand side contains a set of unknowns (the displacements Vi and rotations >i)
and so this system cannot be solved as is. We recognize, however, that the unknowns
are simple functions of the system degrees of freedom; and their contribution to the
load may simply be written as
[Ks]{u}
o 0
o 0
[Ksl ==
o 0
o 0
3.5 Elastic Supports 57
Example 3.6: Consider the uniform beam attached to a spring at its mid-span,
as shown in Figure 3.11. Determine the resulting deflection and rotation at the point
of attachment for various values of spring stiffness.
<D
P
1111~11U"__E_I_,_L ~~ E_I_'_L_~:::I~IIII
We will model this beam with two elements (three nodes). The reduced degrees
of freedom are
{Pk}
PI = P}
= { P2 = 0
T2 = 0
The reduced element stiffness matrix for Elements 1-2 and 2-3 are, respectively,
[k*(I2)j = ~;
12 -12
12
6L]
-6L [k*(23)j = EI [12 6L]
[ sym 4L2
13 6L 4L2
The resulting reduced structural stiffness matrix is (assuming both element lengths
are the same)
[K*j = ~;
12 -12
12 + 12
6L] =EI- [
-6L + 6L
12 -12
24
[ sym 4L2 + 4L2 L3 sym
58 Chapter 3. Beam Structures
The reduced spring constant matrix [K.l is simply the diagonal term associated
with the vertical deflection at Node 2. This gives the overall system to be solved as
EI
L3
[12 24-12
+ a 6L] {
0
VI} = { P}
V2 0
sym 8L2 1/>2 0
The corresponding results in the limit of a very stiff spring (a = 00) are, respectively,
Example 3.7: Reconsider the problem of the last example, but this time place
the applied load at the point of attachment of the spring.
The beam is modeled as in the last example, except that the applied load vector
is different and given by
PI = 0 }
{Pk} ={ P2 =P
T2 =0
This gives the overall system to be solved as
EI
L3
[12 24-12
+ a 6L] {
0
VI} = {O}
V2 P 0'*
aL3
=--
- EI
sym 8L2 4>2 0
PL3 10 _PL 2 12
V2 = EI (48 + lOa.) , 4>2 = EI (48 + lOa)
Now in the limit of a very stiff spring (a = 00) these give respectively
Imposed Displacements
Not all problems of interest are posed in terms of applied loads; it often happens that
we know information in the form of imposed displacements. The exact approach to
solving these problems is to rearrange the stiffness equations explicitly. For example,
consider a [3 x 3] system where we wish to impose a known displacement at the second
node. The rearranged system of equations are
The results of the last example suggests an alternative scheme for imposing known
(non-zero) displacements as a boundary condition. The required displacements can be
imposed by adding the constraint equations that express the prescribed displacement
conditions into the structural equilibrium equations. Assume that the displacement
is to be specified at degree of freedom i, say Uj = u o , then the constraint equation
kUj = ku o
is added to the equilibrium equations as was done for the spring. This means that a
spring stiffness of k is added to the diagonal while simultaneously a load of Pj = ku o
is added to the load vector. If we insure that k ~ Kjj (the assembled stiffness at i),
then the solution of the modified equilibrium equations must now give Uj = U o ' We
note that only the diagonal element in the stiffness matrix is affected, resulting in a
numerically stable solution.
Mathematically, the procedure corresponds to an application of the penalty method.
As regards the computer implementation, it is most conveniently achieved by actually
adding a spring as an additional element.
60 Chapter 3. Beam Structures
Displacement distributions
We assume that the structural problem has already been solved and therefore we
have in hand the nodal degrees of freedom (Vi, ;) at every node. More specifically,
we know all the nodal displacements {VI, I, V2, 2} associated with the member of
interest. The distribution of deflection and slope are then obtained by using the shape
functions of Equations(3.6) separately for each element. That is,
v(x) 91(X)VI + 92(X)L1 + 93(X)V2 + 94(X)L2
(x) = 9~(X)VI + 9;(x)L1 + 9;(X)V2 + 9~(x)L2
where the prime refers to differentiation with respect to x. Because these are written
for each element, the variable x must be understood to be local to that element, and
not the global or structural coordinate; it is zero at the first node and L at the second.
Note that since v(x) is a cubic function in x, it is not accurate to connect the nodal
values with a straight line to get the distribution. This is usually done, nonetheless,
just to get a sense of the distribution pattern.
Load Distribution
Based on our particular beam element model, we have a shear force distribution that
is constant throughout the length, that is,
V(x) = V(O) = V(L)
and a bending moment which varies linearly,
x
M(x) = M(O) + L[-M(O) + M(L)]
In both cases, the distributions are plotted simply by connecting the end values with a
straight line. The task now is to obtain the end values from knowledge of the degrees
of freedom.
We have already shown that the end values of the distributions of moment and
shear force are related to the element nodal load values by
V(O) = -VI, M(O) = -MI
V(L) = +V2 , M(L) = +M2
3.6 Member Loads and Reactions 61
The nodal forces and moments of each element are, in turn, obtained from the element
stiffness matrix by
v;(1) V;(2)
2 1
~M(')2
M")~1
Reactions
Figure 3.12 shows the sign convention used to determine the boundary reactions from
the nodal loads; keep in mind that for our one-dimensional problems we have at most
two elements meeting at a support. Based on equilibrium, the total reactions are seen
to be
PR = -VP) - VP) +P
TR = -MJl) - Mi 2) +T (3.13)
The element nodal load values are obtained from the elemen t nodal degrees of freedom
as indicated above. That is, for each member meeting at the node, we use
Example 3.8: For the beam problem given in Figure 3.7 of Example 3.4, plot
the shear force and bending moment distributions.
62 Chapter 3. Beam Structures
M(x) p ..
V(x) 1----
Figure 3.13: Shear force and bending moment distributions.
VI 0
<Pt 9
V2 TL 0
</>2 = 281 -4
V3 -4L
</>3 0
The element nodal loads for the first element, Member 1-2, are determined from
l 6L -12
VI }(12 [12 4L2 30 }
M1 1 6L -6L 2L2
6L ] { 90 } TL T { 28L
{ V2 = 3 -12 -6L 12 ~~~ ~4 281 = 28L -30
M2 6L 2L2 -6L 2L
2T
V(O) = -VI = 0, M(O) = -M1 = -
28
2T
V(L)=+V2 =0, M(L) = +M2 = 28
Problems 63
Example 3.9: Find the reactions at each support for the beam of the previous
example.
We have already determined the element nodal loads in the last example, now
we substitute these values into the relations for the reactions. At the first node, we
have
(12) -30T
PRl = - VI + PI = 28L '
Notice that the moment reaction is zero, consistent with Node 1 being pinned. There
are two elements meeting at the second node, hence
P _ _ V(12) _ V(23) +P _ 30T - M(12)
T R2-- 2 -
M(23)
1
+ .'2- 0
'T' -
R2 - 2 1 2 - 28L '
Again, the zero moment reaction is in agreement with our expectation for a pinned
support. Finally, for the third node
_
PRJ - - V2
(23)
+ P3 -_ 0,
The zero vertical force is in agreement with our expectation for vertical rollers.
Problems
3.1 Show that when the cross-sectional area or the elastic modulus changes along
the length of the beam that the relevant equations are
d d2 v
<P - dv
- dx' V(x) = - dx[El dx 2 ]
[Reference (13). pp. 542)
3.2 A beam is cantilevered at x = L and has a vertical force applied at x = O. Its
moment of inertia varies as
x
l(x) = 10[1 + rI]
where r is a numerical factor. Show that the deflected shape is given by
2 2
v(x) = rPL3
3 El
0
[r x
(22 +2
rx
+l)-(I+
T
rx
)log(I+
rx ]
) T T +CI X + C2
3.3 Model the previous problem using two elements of different inertia but each of
the same length. (Take that r = 2.) Compare the displacement distributions.
3.4 An alternative way to derive an equivalent load system is to base it on the work
done by the system. Show that if we define the work done by the equivalent
force and moment system as
and on replacing the deflection with its representation in terms of the beam
shape functions that we recover the consistent load formulas.
[Reference [49]. pp. 143]
3.5 Show that the consistent load representation gives rise to a statically equivalent
system even for a general distributed applied load q( x).
3.6 If a beam is subjected to a temperature differential b..T between its top and
bottom, then it will tend curve. Show that the curvature is
d2 v
dx 2 = ab..T/h
where a is the coefficient of thermal expansion, and h is the beam thickness.
[Reference [30]. pp. 115]
3.7 Using the result of the last exercise, show that the stiffness relation for a beam
element with a temperature differential is given by
12 6L -12
6L 4L2 -6L
-12 -6L 12 -6L
2i~ ]{;: }+
V2
aE1b..T ~h ~1 }
{
0
6L 2L2 -6L 4L2 4>2 1
Exercises
3.1 An aluminum cantilever beam of uniform width 100 mm tapers from a thickness
of 150 mm to 75 mm over its length of 4000 mm. If a uniformly distributed load
of 0.5 N /mm is applied find the tip deflection using four uniform elements.
[12.42mm]
3.2 A cantilever beam of length 2L has a non-simple support at its middle. If a
load Q is applied at the end and a rotation 4> applied at the middle support,
determine the reactions at the support.
[T = QL + (4E1/L)4>, P = Q - (6E1/L2)4>]
Exercises 65
3.3 Show that for the uniformly loaded beam of Figure 3.9 that ten elements are
required before the error in maximum deflection is less that 1% when using the
lumped load approximation.
3.4 A carpenter with a power saw has a 20 ft plank of uniform weight per unit
length W o and two saw horses. He wishes to cut a 6 ft length from the plank
but in order to minimize splitting of the ends he wants to cut it at a point
where the bending moment in the plank is zero. If he places one sawhorse at
one end of the plank, where should he place the other? [~ 15.5 ft]
3.5 A bookshelf is made by placing a wooden plank on two brick supports. Where
should the bricks be placed so as to make the maximum bending moment as
small as possible. [~ .21L .79L]
3.6 A 5 m beam of cross-section 100 mm x 100 mm is built in at both end. If one end
=
slips an amount fJ 25 mm relative to the other end, determine the maximum
moment generated. [M = 6EIfJIL2]
3.7 A cantilever beam of length 2 m has a distributed load of 4800 N 1m over half of
its length beginning at the free end. It also has a concentrated force of 3000 N
at the mid point and a linear spring (stiffness 200 N 1m) attached at the tip.
The beam has an EI of 40 kN 1m2 Use two elements and the consistent load
approximation to obtain the deflections at the mid point and tip.
[11.66 mm 6.80 mm]
3.8 Timber beams 18 ft long, 12 in deep and 4 in wide, are simply supported at
the ends. How far apart, center to center, should such beams be placed when
supporting a floor loaded with 40 lbl ft 2? The beams weigh 40 lbl ft 3 and the
maximum allowable bending moment in each beam is 6800 lb ft. [3.86 ft]
3.9 A beam, 20 ft long, carries a uniform load of 1 toni ft run. It is simply sup-
ported at one end and at some other point. Find the position (from the free
end) of the other support so that the maximum bending moment may be as
small as possible. [5.87 ft]
3.10 A timber beam, 8 in wide and 6 in deep, is placed directly above another timber
beam 6 in wide and 8 in deep. The beams are held apart by three solid blocks,
one under each end and one under the center of the upper beam. The beams
are 20 ft long, and the whole rests on two supports, one under each end of the
lower beam. If the upper beam is loaded with 100 lbl ft run, find the deflection
at the center, and the maximum bending moment in each beam. (E = 1.5 msi)
[0.6 in; 1800, 4000 lb ft]
Chapter 4
A truss is a structure composed of rod members arranged to form one or more tri-
angles. The joints are pinned (do not transmit moments) so that the members must
be triangulated. A frame, on the other hand, is a structure that consists of arbi-
trarily oriented beam members which are connected rigidly or by pins at joints. The
members support bending as well as axial loads.
The essential new aspect to the study of these structures is the consideration of
the element stiffness of an arbitrarily oriented member. Since differently oriented
members are to be considered simultaneously, they must have common or global ref-
erence axes. The connection between the global and the local axes is established via
the 3-D rotation matrix. We will first develop the analysis for plane structures as an
intermediate step to introducing the general case of space frames. References [3, 45]
are excellent sources for additional details on modeling 3-D structures.
A fundamental assumption in the following developments is that the principle of
superposition holds. In this way, we can assemble the general frame by combining
the separate actions of the simpler cases developed in this and previous chapters.
66
4.1 Truss Analysis 67
LOCAL
Figure 4.1: Load and displacement components in local and global coordinates.
referred to the local coordinates x - y, is as given for the rod. That is,
{ :~ } = ELA [ _ ~ ~ 1 ] { ~~ }
The forces and displacements referred to the local coordinates are now denoted by
the barred notation. We wish to give the truss member extra degrees of freedom
that allow it to move in the plane. That is, we want to add the degrees of freedom
VI and V2. The corresponding stiffness relation is obtained by augmenting the above
equation by adding zero forces in the y-direction, and the nodal displacements VI
and V2 in the y-direction at Node 1 and Node 2, respectively. The resulting element
stiffness relation is
or {F} = [ Tc Hu}
The axial forces FI and F2 can be decomposed into global components in the x
and y-direction. At Node 1, for example,
[ k ] = [ T ]T[ k )[ T ]
That is, the global stiffness matrix is obtained from a transformation of the stiff-
ness matrix in local coordinates. The explicit expression for [ k ], after multiplying
through by the transformation matrices, is
(4.1)
where C == cos () and S == sin (). Note that while the special case of () = 0 corresponds
to a horizontal rod, the resulting stiffness is not that of the rod - the presence of the
VI, V2 degrees of freedom means that the member can move transversely. To recover
the one-dimensional rod behavior, it is necessary to suppress these degrees of freedom
by imposing some transverse constraints.
This element stiffness relation allows a simple expression for the axial force to be
determined from
P = F X2 cos () + F Y2 sin ()
By substituting for the force components we get
- EA
F = L [(U2 - udcos() + (V2 - VI) sin ()] (4.2)
This is useful as a quick means to obtain the axial force once the global degrees of
freedom are known.
4.1 Truss Analysis 69
tPyl
L....Px1
F(12) -:-e _ F(13)
xl ~ ~ xl
F(12) F(13)
yl yl
This free body diagram of the node is typical of the other nodes also because all force
components are written in terms of the global reference frame. To amplify on this
point, the specific orientation of a member does not enter these free body diagrams;
that information is included in the stiffness matrix. Therefore, we can write similar
equilibrium equations at Node 2 and Node 3. The six equations in all can be expressed
in the form
(12) (23) (31)
Pxl F X1 0 F X2
Pyl FyI 0 F Y2
Px2 F X2 F Xl 0
Py2
=
F y2 + FyI
+ 0
Px3 0 F X2 F Xl
Py3 0 F Y2 FyI
There is a force vector on the right hand side for each member in the structure. These
vectors are related to the nodal displacements through the [4 x 4] element stiffness
70 Chapter 4. Truss and Frame Analysis
{F
o
}(12) = [[k(t2)]
0
0]
0 { },
U
{F0}(23) = [00 [k(23)]
0 ]{ U }
Combining all the stiffnesses together leads to the structural equilibrium equations
The stiffness property of the truss as a whole is characterized by the square symmetric
matrix [/{] which relates the external loads to the nodal displacements. The entry
/{ij which relates Pi with Uj, is the sum of all the entries in the element stiffness
matrices that relate Pi and Uj.
Connedivities
In order to use the general form of the element stiffness matrix, we must be consistent
in our designation of the orientation. For example, a member that is upright can be
said to have an orientation of either 0 = 90 0 or 0 = -900 depending on which end is
considered the 'first' node. We will state this information in the form of a connectivity.
That is, we will state (as a property of the member) the two nodes that it connects.
Note that the actual numbering sequence of the nodes is not important, only which
of the two is to be considered the first node.
Consider the numbering of the frame shown in Figure 4.3. The numbering of the
members is quite arbitrary but, as will be seen later, the numbering of the nodes can
be very important to the computational efficiency. The connectivities can be specified
4.1 Truss Analysis 71
Example 4.1: A three-bar truss is loaded as shown in Figure 4.4. Find the
displacements at the load application point. Each member has the same EA.
The power of the matrix methods is that similar solution procedures can be
used irrespective of the particular structural shape. Thus we can use the approach
introduced for the rod and beam structures. Labelling the nodes as shown gives the
total degrees of freedom as
The displacements at Nodes 2, 3, and 4 are zero leaving the unknown nodal dis-
placements and the corresponding known loads as
EA
() = 90 0
<D 1000#
2000#
{IDbc} = {I ,2;0,0;0,0;0,0}
since only the first two degrees of freedom are non-zero.
To construct the reduced stiffness matrix [K*), it is necessary for us to find
the contributions from each member. We will state the connectivities of the three
members as, 1-2, 1-3, 1-4, respectively. In each case, then, the reduced element
stiffness is of the form
[ k* ) -_ EA
L
[C2 C 5 ]
C5 52
With reference to the figure, the orientation of the three mem bers are 120 0 , 90 0 , and
60 0 , respectively. Substituting these values of orientation gives
The reduced structural stiffness matrix is formed by collecting the elements in the
reduced element stiffness matrices according to their position vis-a-vis the global
degrees of freedom. In this case each matrix shares the same global degrees of
freedom, hence the result is
Once the nodal displacements are obtained, the loads in the truss members can be
calculated.
Example 4.2: Determine the deflections of the equilateral truss shown in Fig-
ure 4.5. One joint is on a horizontal roller.
EA,L
Both displacements at Node 1 are zero, while the vertical displacement at Node 2
is zero. This gives the unknown nodal displacements and the corresponding known
loads as
{IDbc} = {0,0;I,Oj2,3}
We will state the connectivities of the three members as, 1-2, 2-3, 3-1, respec-
tively. For Member 1-2, the orientation is 0, giving the reduced stiffness matrix
as
[k(*12)] = EA[ 1 ]
L
For Member 2-3, the orientation is 120, giving the reduced stiffness matrix as
-1
1
V3]
-V3
-V3 3
74 Chapter 4. Truss and Frame Analysis
The associated degrees of freedom are {U2' U3, V3}' Finally, for Member 3-1 the
orientation is -120", giving the reduced stiffness matrix as
for the associated degrees of freedom {U3' V3}' The reduced structural stiffness
matrix is formed by collecting the elements in the reduced element stiffness matrices
according to their position vis-a-vis the global degrees of freedom. The result is
IX] : ~: [j -~ 1]
The system of equations to be solved is
EA [ _
4L V3
~ - ;
0
~]
6
{ :~
V3
} ={ ; }
0
Once the nodal displacements are obtained, the loads in the truss members can be
calculated.
The element stiffness relation for the beam element with axial loading can be ex-
pressed as
{P} = [ k Hii}
where the following augmented matrices have been introduced
PI iii
iII
M
VI
I ~l
{P} == {u} ==
P2 U2
V2 V2
M2 ~2
and
1 0 0 -1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 12 6L 0 -12 6L
[k] = EA 0 0 0 0 0 0 EI 0 6L 4L2 0 -6L 2L 2
- L -1 0 0 1 0 0 +3 0 0 0 0 0 0
0 0 0 0 0 0 0 -12 -6L 0 12 -6L
0 0 0 0 0 0 0 6L 22 0 -6L 4L 2
The beam under combined loading is seen to be represented by a [6 X 6] stiffness
matrix. We will generally leave it in separated form like this so as to emphasize that
it arose as a simple superposition of the axial and flexural actions.
Fy2
M2
FX2
system x - fi and the global coordinate system x - y of the member shown in the
Figure 4.6. The nodal forces referred to these two systems are denoted by
PI FXI
lit FYI
{P} =
MI and {F} =
MI
P2 FX2
V2 Fy2
M2 M2
respectively. The relation between {P} and {F} is obtained by using the transfor-
mation for vector components. Thus, for example, at the first node
There are similar expressions at the second node. We can write all these equations
in the symbolic form
{P} = [ T ]{F}
where
cosO sinO 0 0 00
-sinO cosO 0 0 00
0 0 1 0 00
[TJ= 0 0 0 cosO sinO 0
0 0 0 -sinO cosO 0
0 0 0 0 0 1
It is easy to confirm that this matrix is orthogonal. Proceeding in like manner, we
obtain for the nodal displacement vectors {u} and {u}
UI UI
VI VI
<PI <PI
{u} == {U} == {u}=[T]{U}
U2 U2
V2 V2
<P2 <P2
Substitution for the barred displacements and forces into the element stiffness relation
gives
[ T ]{F} = [ k )[ T ]{ u}
from which the following is obtained
{F} = [T JT[ k )[ T ]{ u} = [k ]{ u}
4.2 Plane Frame Analysis 77
since [ T 1is orthogonal. The explicit expression for the frame element stiffness rela-
tion is given as
FXI C2 sym
FYI CS S2
MI EA 0 0 0
-C2 -CS 0 C2 (4.3)
FX2 L
Fy2 -CS -S2 0 CS S2
M2 0 0 0 0 0 0
12S2 sym UI
-12CS 12C2 VI
E1 -6LS 6LC 4L2
+v
<PI
-12S 2 12CS 6LS 12S2 U2
12CS -12C2 -6LC -12CS 12C2 V2
-6LS 6LC 2L 2 6LS -6LC 4L2 <P2
where, as before, the abbreviations C == cos 0, S == sin 0 are used. Note that both of
these matrices reduce to the respective matrices when 0 = O. But also note that the
first matrix is the augmented global stiffness for the truss element. Therefore, in a
sense, the truss behavior is embedded in the frame.
This relationship shows that even for the arbitrarily oriented frame member that
the axial and flexural actions are uncoupled. This decomposition of the frame stiffness
is due to our initial assumption.
The boundary conditions at the first and third nodes require that
UI = VI = <PI = 0 , V3 =0
78 Chapter 4. Truss and Frame Analysis
~----~,.;z:-~p
EA,EI,L
{Pd = {O,O,O,P,O}
The corresponding I Dbc matrix
{IDbc} = {O,O,Oj1,2,3;4,O,5}
We will take the connectivities for the members as 1-2 and 2-3. The reduced
element stiffnesses are obtained by substituting 8 = 90 and 8 = 0 for Members 1-2
and 2-3, respectively. This gives for Member 1-2
The unknown nodal displacements can now be solved for and are
4.2 Plane Frame Analysis 79
-6P [EA + 3 El ) P
L4>2 = D T L3 EI
p[ EA
U3 = El)P
D 7T +30];3 EI
3P [EA 6 El ) P
L4>3 = D L
_
];3 EI
D 12[4 + 3ft)
This form emphasizes that the two stiffness terms interact in a complicated fashion
to give the final deflections. It also shows that the sense of the rotation at Node 3
depends on the relative values of the axial and flexural stiffnesses.
Example 4.4: A two-member frame is loaded as shown in Figure 4.8. One end
is fixed while the end with the applied load is fixed to a horizontal rollers. Both
members have the same material and section properties, and the joint connecting
them is at 90 0
Find the deflections.
EA,EI,L
The unknown nodal displacements and known applied loads are easily identified as,
respectively,
and
{IDbc} = {0,0,0;1,2,3;4,0,0}
80 Chapter 4. Truss and Frame Analysis
We will take the connectivities for the two members as 1-2 and 2-3, respectively.
For member 1-2 we have (J = 45, hence the reduced stiffness is
[k*(12l] = EA
2L [110] + [6
1 1 0
0 0 0
E;
L
-6
-6
6
3LV2 -3LV2
3LV2]
-3LV2
4L2
6LV2
EA 02 02 00 -1]
1 E/ [120 0
12 o
[
2L [ 0 0 0 0 + L3 6LV2 0 8L2
-1 1 0 1 -6 -6 -3LV2
rP2 = 0
The unknown nodal displacements are
U2 }
V2
= !:.. { 2 } P L3 {~1 }
4EA 1 + 48E/
{
U3 4 4
Note that if the axial flexibility is small enough, this result predicts that the center
can actually move upwards.
EA 0 0 0 0 0
0 12EIz/ L2 0 0 0 -6EIz/ L
~33 ~34 +1 0 0 12EIy/ L2 0 6EIy/L 0
[ ] =
k 43 k 44 L 0 0 0 GIx 0 0
0 0 6EIy/L 0 4EIy 0
0 -6EIz/L 0 0 0 4EIz
(4.4)
For regularity of notation, the torsional stiffness is written as GIx = GJ. The special
cases of plane frame and grid, for example, can be obtained by setting the appropriate
degrees of freedom to zero. This will be done in a later section.
A point to note is the regularity in the repetition of terms. For example, the
first quadrant of [ k ] is almost the same as the fourth except for the signs of the
off-diagonal terms.
ITI~[r ~ ; l1
is a [12 x 12] matrix. Substituting for the barred vectors into the element stiffness
relation allows us to obtain the global stiffness as
This is formally the same relation as obtained for the truss and plane frame. In
practice, however, we take advantage of the special nature of [ T] to reduce this
further to
[kll ] = [R f[kllH R], [k12 ] = [R f[kn][ R] or
4.4 Determining the Rotation Matrix 83
where the summation is over each member. In practice, there is no need to augment
the member stiffness since we assemble the reduced global stiffness directly. It is
important to realize that the transformation occurs before the assembly.
where lx, m x and n x are the cosines of the angles that the x axis makes with the
global x, y, Z axes, respectively, as shown in Figure 4.9. The other terms give the
orientations of the iJ and z axes, respectively. Since the member axis of a frame or
truss coincides with x, then the direction cosines of the first row can also be written
as
Ix = (Xj - xi)/L ij , m x = (Yj - y;)/L ij , n x = (Zj - zi)/Lij (4.5)
where (Xi, Yi, Zi) and (xj, Yj, Zj) are the coordinates of the first and second nodes,
respectively, and L ij is the length of the member. The problem here is to find the
remaining elements of [ R ], as lx, m x , and n x define only the orientation of the member
x axis.
In order to determine the rotation matrix we will assume that the member arrived
at its current position by successive rotations of axes. Consider the member to be
initially oriented along the X axis. The first rotation is through an angle Q' about the
Z axis. The second is a rotation through an angle f3 about the iJ axis. (This sequence
leaves the member iJ-axis always oriented so as to lie in the global X - Y plane.) The
resulting rotation matrix is therefore
member
Equating the first row to the direction cosines of the member gives
Therefore, the functions cos a, sin a, cos {3 and sin {3 may be expressed in terms of the
direction cosines of the member by
lx . mx
cos a = D' slna = D' cos {3 = D, sin {3 = n x , D =- .)12x + m x
2
~]
This is the rotation matrix [ R 1for a space truss member and is also valid for space
frames when the member has a symmetric cross-section.
Non-symmetric Sections
A space frame member may have its principal axes in general directions. There
are various ways in which the orientation of these cases can be defined and the one
chosen involves specifying the orientation by means of an angle of rotation about the
member i; axis. In order to visualize clearly how such an angle is measured, consider
three successive rotations from the structure axes to the member axes. The first
two rotations through the angles a and {3 are exactly the same as before. The third
transformation consists of a final rotation through the angle, about the member
i; axis, resulting in the fj and z axes coinciding with the principal axes of the cross
section.
4.4 Determining the Rotation Matrix 85
The rotation of axes through the angle I about the member axis requires the
introduction of a rotation matrix [Ry 1 given by
Multiplication of the three successive rotations gives [R] = [Ry][R,a][Rc.] and when
multiplied out in terms of the direction cosines this becomes
mx
Ix cos I - mxn x sin I
Dsin l (4.6)
. D
-Ix Sill 1- mxn x cos I
Dcos l
D
This rotation matrix is expressed in terms of the direction cosines of the member
(which are readily computed from the coordinates of the joints, Equation(4.5)) and
the angle I' which must be given as part of the description of the structure itself.
Note that if this angle is equal to zero, the matrix [ R] reduces to the form given
previously for a space truss member.
Special Case
When the member axes are specified in the manner just described, there is no am-
biguity about their orientations except in the special case of a member oriented
along the global z-axis. There is no unique rotation to get to that orientation, i.e.,
0' = 0, f3 = 90 or 0' = 90, f3 = 90. To overcome this difficulty, the additional
specification will be made that the member y-axis is always taken to be along the
global y-axis for these cases. That is 0' = 0, f3 = 90.
The complete set of direction cosines is therefore
10 0][00 o1 n
x o
[o -
]
[Rzl = 0 CO~I sin l 0 cos I
Sill I cos I -n x o 0 -Sllli
This is for the general case of a member of non-symmetric section. All that is necessary
is to substitute for the direction cosine n x its appropriate value, which is either 1 or
-1.
This gives D = V5l6. The structure is specified as being a truss, hence the angle I is
zero, therefore the rotation matrix can be obtained by substituting into Equation 4.6.
The result is
.4082 -.8165 -.4082]
[ R J= .8944 .4472 .0000
[ .1826 -.3652 .3727
The length of Member 3-4 is
This is the special case with n x = -1, hence the rotation matrix is
00-1]
[RJ=
[0 1
1 0
0
0
Global Reductions
There are six degrees of freedom at each node in the space frame. Many problems,
however, do not need this many; for example, the plane frame only requires three,
while the plane truss uses two. Obviously to analyze a 2-D structure as a 3-D frame
is a waste of computer resources.
The key to understanding the reduction of the general case is the idea of imposing
constraints. We saw in the case of fixed boundary conditions that we specify the
degree of freedom as zero, and consequently 'scratch' the associated rows and columns
in the stiffness relation. In essence we do the same here; we specify constraints on the
4.5 Special Considerations 87
u = v = 0, <Pz = 0
at each node. The non-zero degrees of freedom are the out of plane displacement
w, and the two in plane rotations <Px, <py. The member nodal forces and degrees of
freedom are
FZ1 WI
M X1 <PxI
M yl <Pyl
{F} = {u} =
FZ2 W2
M X2 <Px2
M y2 <Py2
0 0 0 0 0 0 12 0 -6 -12 0 -6
0 1 0 0 -1 0 0 0 0 0 0 0
[k] = GJ 0 0 0 0 0 0 EI -6 0 42 6 0 22
0 0 0 0 0 0 +3 -12 0 6 12 0 6
0 -1 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 -6 0 22 6 0 42
where GJ is the torsional rigidity of the member and EI is the flexural rigidity for
bending out of the plane. The rotation matrix
[R] =
Ix
-mx
m x
Ix
0]
0
[
001
is used to obtain the stiffness matrix for an arbitrarily oriented grid member.
Similar global reductions can be done to recover the cases already considered as
well some new ones such as the shaft in torsion.
For large scale problems, however, it is much simpler to find the boundary reac-
tions by adding a 'boundary element' to the structure. The basic idea is to connect
the structure to a rigid support using a small frame element. The nodal loads of this
element are the required boundary reactions. This is feasible since many elements
are already being used to model the structure.
Actually this idea of the boundary element adds a lot of flexibility to the matrix
analysis method. We already saw how it allowed the specification of displacements
as boundary conditions, we next show how it can also be used to implement oblique
supports.
We have formulated the stiffness approach in terms of a global coordinate system.
Therefore the allowable constraints must also be in terms of the global coordinates.
Consider the case of a frame with oblique supports, that is, the frame is attached
to rollers on an inclined surface. The boundary condition is that the displacement
normal to the surface is zero. This is a constraint condition written as
(a)
c:::::=:::::=~11111111
(b)
~
ttmm~mm
(c)
411111!111
Figure 4.10: Equivalent boundary conditions.
infer this information from the symmetry (or anti-symmetry) of the geometry and
loading conditions. We can then implement these constraints by use of equivalent
boundary conditions and thereby simplify the problem.
Figure 4.10 shows a few examples for a beam, but the idea is quite general. The
first structure shown has a loading symmetric with respect to the center of the beam.
Thus, at the center, 1> = 0 but v # O. An equivalent boundary support is as shown
on the right hand side. It should be noted that only half of the load acts on the
equivalent model. In fact, the rigidity of members lying in a plane of symmetry must
also be halved in order to divide the structure into two equal parts. The second case
corresponds to an anti-symmetric loading condition achieved with either a couple or
moment. This corresponds to a pinned support condition. The symmetric version of
the problem corresponds to a fixed support.
The use of symmetry and anti-symmetry does not involve any approximation and
therefore when the opportunity arises, advantage should be taken of it. It is worth
keeping in mind that this can be done as long as the structure is symmetric, because
any unsymmetrical loading can be decomposed into the sum of a symmetric and
anti-symmetric load.
ty
Lx
! I !I~
p p p
CD
Due to the symmetry, only half of the structure need be analyzed. This is done
by modeling it as shown where the vertical displacement at the apex is left free while
the rotational and horizontal degrees of freedom are suppressed. Number the nodes
as shown, then the total degrees of freedom are
VI = 0, U2 = 0,
The unknown nodal displacements and known applied loads are easily identified as,
90 Chapter 4. Truss and Frame Analysis
respectively,
and
[K*] = EA
2L [ 1]1 + ~
1
1
[ 6)2
)213 -6)2
-6L
-6)2 ]
6L
6)2
The unknown nodal displacements are
1 { P } P {L} P L2 { 2L }
:~ = [K*t ~ = 2EA ~ + 12El ~fl
Ut }
{
Inextensible Structures
Many frame structures are more flexible in bending than in axial deformation, and,
consequently, the axial deformation can often be neglected. The structure is then
referred to as being in extensible. Sometimes, the complexity of a structures problem
may be reduced by making use of this approximation.
it is usual to assume that the structure is inextensible. For example, reconsider the
frame problem in Figure 4.7. In the limit as the axial stiffness is very large then we
have
7PL3
48E1
o
_PL3
8E1
7PU
U3 = 48E1 = U2
PL3
16E1
In that case, V2 ~ 0 and the horizontal displacements at Nodes 2 and 3 are the same.
The inextensibility assumption is implemented in the form of a set of new bound-
ary constraints. Figure 4.12 shows the modeling of the above problem. The vertical
deflection at Node 2 (being assumed small) is forced to be zero by attaching the
node to a horizontal roller by pins. The inextensibility of the horizontal member is
implemented by the constraint relation
It should be noted that when axial deformation is neglected, the axial force cannot
be calculated from the nodal displacement solution. Instead, equilibrium must be
used.
.-----~---.----------- .
./
,
,
J----i-----J-
, -- -- --~
,
1--+-----1---1
I---'-----~---I
/ EI 2EI
There are three degrees of freedom {u, v, 4>} at each joint. The inextensibility of
the vertical members results in
v=o
92 Chapter 4. Truss and Frame Analysis
at each joint. The inextensibility of the floor, however, gives the constraint
{VVI} = EI
2
[ 12
13 -12
-12] {
12
UI }
U2
= 12EI [1
13-1
This resembles the rod stiffness relation, but actually the displacements are trans-
verse to the axis. This element is called a shear element.
In this way we reduced the system size from 18 down to 3. This is a very useful
approximation in the dynamic analysis of buildings.
4.6 Substructuring
We conclude this chapter with an introduction to the approach used to analyze com-
plex structures. The division of a structure into components for separate analysis is
called substructuring. The main advantage of the approach is the replacement of a
large problem by many smaller ones each of which may be more manageable for a
given set of resources. Also, design changes on a component do not affect the work
already done on the other components. As compared with a single pass analysis, the
computational cost is usually much greater, the coding is usually much more compli-
cated and needs to be customized to specific problems. When the structure is made
of repetitive substructures, this approach is very efficient. The basic idea, as will
be shown, is essentially that of a super-element - the substructure is used in the
same way as an individual finite element with internal degrees of freedom that are
eliminated prior to the element assemblage process.
Consider the case of a structure divided into two. Now look at one of those
substructures. We can partition the arrays as
where the subscript i refers to interior nodes, and c refers to common or connection
nodes. The equations for the substructural system are expressed in the partitioned
form as
J(ii J(iC] { Ui } = { Pi }
[
J(ci J(cc uc Fc
4.6 Substructuring 93
For concreteness, we will assume that we have already reduced the zero degrees of
freedom from the system. Multiplying this out gives
(4.7)
Substitute this into the second equation to give
[IC]{ u} = {P;}
The load term {Qe} = [Kei]{u} contains the reactions at the fixed connection points.
Example 4.8: Analyze the truss structure shown in Figure 4.14 by substructur-
ing. Each member has the same length and section properties.
The original structure will be divided at Node 2; thus there is only one connection
node and two connection degrees of freedom. Both substructures are similar and
therefore have the same reduced structural stiffnesses. These are formed in the usual
way and can easily be shown to be
-y'3
[J(*] = EA [ -~ 3 -4o 0]
0
4L -4 o sy'3
o o y'33
94 Chapter 4. Truss and Frame Analysis
The [6 X 6] stiffness matrix has been reduced by scratching the rows and columns
associated with the degrees of freedom at Node 3 or Node 5.
For Substructure 1, the partitioned arrays are
, EA [
[Ri;] = 4L
5
-v'3 -v'3
3
] .
[Ric] 0 ~]
EA [ -4
= 4L
EA
[K ] = 4L [5v'3 v'3]
3 -
EA
4L [40 00]
The stiffness relation for the substructure is (the load term {Qe} is zero)
{ Fxe } = EA [ 1 v'3] { U e }
Fye 4L v'3 3 Ve
,
[Ai;]
EA
= 4L 3
[5
v'3 v'3] ,
[Aie] =
EA
4L [-40 0]
0
[K .]
Ct
= EA
4L
[-4 0] 0 0 [Kee] = 4L -v'3 -v'3
EA [ 5
3
]
Problems 95
{Fe} = {~::}
The inverse matrix is
[ K-.)-l _ ~ [ 3 -V3 ]
It - 3EA -V3 5
5 -V3]
EA [ -J3
[K) = 4L 3 EA [40 00]
- 4L
The loading term is non-zero and determined as
EA[-4
{Qe} = 4L 0
0] L [ 3
0 3EA -J3
-J3]{Px}
5 Py = 121[-120 4J303]{~y}
r.
Ve =0
Once the nodal displacements are obtained, the loads in the members can be calcu-
lated.
Problems
4.1 Maxwell's reciprocal theorem states that for an elastic structure the deflection
produced at some point i by a unit load at point j is equal to the deflection at j
due to a unit load at i. Show that this statement is true for systems described
by the stiffness relation
{P}=[K){u}
[Reference (45). pp. 43)
4.2 Show that a frame attached rigidly to rollers on an incline can be modeled
using the idea of two boundary elements separated a distance and connected
by a cross- bar.
[Reference (45). pp. 375]
96 Chapter 4. Truss and Frame Analysis
Exercises
4.1 Give an example that proves it is possible that the diagonal term ki i of the
element stiffness matrix can be zero. Show, by a simple example, that it cannot
be less than zero.
V2 = VI + </>IL,
Show that for the plane frame, the forces required to produce such a motion
are zero.
4.3 For the truss of Figure 4.3, renumber the nodes in an attempt to reduce the
bandwidth. State the connectivities and assemble the connectivity matrix.
4.4 Solve the problem of Figure 4.4, but choose the connectivities as 2-1, 3-1,4-1.
4.6 Solve the problem of Figure 4.4, but make use of symmetry and anti-symmetry.
4.7 For the truss problem of Figure 4.14, solve it making use of symmetry and
anti-symmetry.
4.9 A steel portal frame ABCD, fixed at A and D, is 8 it high and 12 it wide. It has
a distributed load of 500 fbi it along BC and a 3000lb horizontal force applied
at B. Determine the deflections. [{ UB} = {0.092 in, -0.001 in, - .0014 T} J
Structural Stability
[k ]u= Q, (5.1)
where k b is the axial stiffness of the bar. (Identical results are obtained if the matrix
methods of Chapter 4 are used.) For the purpose of the later discussion, let the bar
97
98 Chapter 5. Structural Stability
be very stiff, from which we conclude that the displacement is only horizontal, i.e.,
v = O. Keep in mind that all displacements are considered to be very small.
1!1/!I!!
We will now look at equilibrium again, but this time based on the deformed con-
figuration. Consider the situation when the bar has already displaced by an amount
u as shown in Figure 5.1. Since it is in equilibrium, summing the moments about the
base gives
-QL-Pu+kuL=O
This can be put in the form of a stiffness relation
[k - P/L}u =Q (5.2)
In general, we can solve this for the unknown displacement u if both loads P and Q
are specified. Contrast this stiffness relation with the one obtained in Equation(5.1);
using the deformed position as the reference state for application of the equilibrium
equations has coupled both applied loads in the stiffness relation. The consequences
of this, as we shall see, are what dictates the stability behavior of a structure. For
example, in the present case a critical situation arises when the applied load P ap-
proaches a value near kL since then an inordinately large displacement is indicated.
Indeed, even if Q is almost zero very large deflections are calculated. Obviously, the
structure is experiencing some critical behavior; in fact, it is becoming unstable. The
critical condition is when P = kL and this value of load is called the critical load Per.
The critical load is the borderline (or transition) between the structure being stable
and unstable, and is the primary concern in stability analysis. Additional insight
can be gained by inquiring as to what force (Q) is necessary to a achieve a given
displacement u. At or near the critical load, it takes virtually no load Q to move the
structure any amount. We call the situation neutral equilibrium.
Some important points can be drawn from this example. First, the question of
stability arose only after we used the deformed configuration to establish the equi-
librium conditions. Second, the critical behavior occurred when the stiffness tended
5.2 Stability of Truss Structures 99
to zero. That is, instability is associated with a loss of stiffness; there is no longer
a unique displacement solution for a given set of loads. Third, while we call the
term in square brackets in Equation(5.2) a stiffness, it is unlike any of the stiffnesses
encountered in the previous chapters because it includes one of the applied loads.
Actually, this is the total stiffness and is comprised of an elastic stiffness k and a
geometric stiffness P / L. The latter is due to the geometry change of the equilibrium
configuration, hence its name and why it is load dependent. Finally, the critical value
of P does not depend on the value of Q. But it is important for us to realize that Q is
necessary to actually cause instability. That is why our concept of stability involves
imagining the consequences as the structure is displaced slightly, in other words, we
will assume there is always an upsetting force present.
The concept of stability we will develop is rooted in the linear, small deflection
analysis of structures. After instability has occurred, it may well be that the structure
achieves a new equilibrium position corresponding to the large deflection position.
What we are concerned with is occurrence of a neutral equilibrium state under small
deformation conditions.
Fy F
q( x) Lc:::: Fx
~t- F+b.F
EA
.... 1.., .
........ A.. _
~ 'I' -
dv
dx
F-!V---
V b.x V + b.V
b.x
Consider a typical member shown in Figure 5.2; for convenience the reference axis
is taken along the member and the deflection is assumed only in the plane. This
figure is similar to that of Figure 2.1 except that the infinitesimal element is shown
slightly displaced in the y direction and slightly rotated about the z axis. The axial
force F acts along the member, but since we want a set of equations written with
100 Chapter 5. Structural Stability
- - - dv
Fx = F cos ~ F , Fy = Fsin ~ F ~ F dx
This last relation shows that the member has a shear force that depends on the
amount of rotation. Thus the equilibrium relations must also take this shear into
account. To make the relations look similar to those of the previous chapters, we
introduce the notations
F == Fr , V == Fy
Summing the forces on the element gives (in the limit of ~x becoming very small)
dF dV dv
dx = -q, -=-q-
dx dx
A slightly more convenient form of the second equation is obtained by summing the
moments, this results in
dv
V =F dx
These equilibrium equations are actually coupled (through F) and this makes them
very difficult to use. In fact, they are a system of nonlinear equations. Since our
whole development is for small displacements, we take advantage of this and solve
the equations approximately in two steps. First we solve them with respect to the
original configuration and compute an estimate of the axial force. Call this force the
initial force Fo and for reasons to become apparent later call this analysis the pre-
buckle analysis. We then re-solve the equations, this time using Fo as a parameter in
the shear expression. This effectively uncouples the above equations.
All the relationships for the structural quantities may now be collected as
Note that the initial force term Fo plays a role similar to that of the axial stiffness EA.
As we have seen so many times before, the displacements (u and v) can be viewed as
the fundamental unknowns of interest. In the special case when the loading is zero,
we have for these displacements
5.2 Stability of Truss Structures 101
Since the displacements are linear, we conclude that the axial force and shear force
are constant along the section.
The manner in which we have posed the above problem makes it a continuation
of the methods already developed in the previous chapters. That is, the solution
to a particular problem is obtained by determining the constants of integration for
each section by imposing equilibrium and compatibility at the junctions. We give an
example to demonstrate the approach.
~-..:..(1~)- - & r - P
III!I:[::+:
tQ
F(I)--Ie-p
V(I)~n(2)
+F(2)
(2)
EA,L
Normally, this stage of the solution can be quite extensive, requiring a full analysis
of the static problem. We will see this in some of the later examples.
For the horizontal member, the first boundary is pinned, hence the deflections
are
u(x) = CIX,
Choosing a local x for the vertical member, we have
C1 = -bi,
102 Chapter 5. Structural Stability
Example 5.2: Determine the sensitivity of the joint displacements of the last
example to small variations in the applied loading.
We will first rewrite the above results in the form of stiffness relations for the
joint (x = L)
EA pJ2) EA pJ1)
[L+T]UL=P, [L+T]VL=Q
It is important to remember that the initial loads pJ1) and pJ2) are obtained from
the undeformed configuration and therefore they will always have the relation
PJ1) = P, PP) =Q
irrespective of the actual displacements UL and VL. That is, we view the pre-buckle
analysis as establishing the stiffness properties of the structure which can then be
treated as constant for the subsequent analysis.
It should be pointed out that the contribution to the stiffness from the initial
load is quite small in the above example. For example, for the contributions to be
equal would require the axial stress to be of magnitude E. This is much too large
since the yield stress of a material is usually of the order E /100.
Example 5.3: Determine the relation between the displacements at the apex of
the truss shown in Figure 5.4 and the applied loads. Use equilibrium based on the
deformed configuration.
5.2 Stability of Truss Structures 103
Proceeding as in the last example, we first solve the equilibrium equations in the
undisturbed configuration to determine the initial loads. This analysis gives that
F(l) cos (J - V(1) sin (J - V(2) = P, F(1) sin (J + V(1) cos (J + F(2) = Q
where we have taken (J = 45 and I = V2L. Since we actually want the stiffness
relation, we can avoid having to solve this system by writing the coefficients in terms
of the displacements. That is,
While the stiffness matrix has a linear dependence on the initial loads, the resultant
effect these loads have on the displacements is not clear.
We can solve for the displacements only if the determinant associated with this
system of equations is non-zero. Thus the question arises: Is it possible for the
determinant of this stiffness matrix to be zero? As a particular case, suppose that
P = 0, then the initial loads are PjIl = 0 and PJ2l = Q. The determinant of the
stiffness matrix is
It is apparent that Q can be chosen such that the determinant is zero. This value is
given by
-EA
Q = 1 +2V2
Note that Q must act so as to put one of the members in compression. Even in
the general case, we can find values of P and Q so that the determinant is zero. In
contrast to the earlier example where individual stiffness terms went to zero, we see
that for a general system the indication of neutral equilibrium is that the determinant
of the stiffness matrix goes to zero.
To see the consequences of this, recall that we always set up a system of simulta-
neous relations as
[J<]{U}={p}
where {p} represents the known applied loads. (Generally, we will have already incor-
porated the known displacement boundary conditions.) A condition used implicitly
in solving all the simultaneous equations of the last few chapters is that det[ J< ] =I o.
As shown above, however, it is possible to choose initial loads to make det[ J<] = o.
At precisely those values there is not a unique connection between the displacements
and the boundary conditions. This means that two or more of the equations are
linearly dependent and therefore the solution is not unique. In other words, there is
a second solution {u} that also satisfies the equations
Let the difference of the two solutions be {4>} = {u} - {u}; it therefore satisfies the
homogeneous equations
[J<]{4>} = {OJ
In general, such an homogeneous equation admits only the trivial solution {4>} =
0, but when det [ J<] is zero there are other solutions. These are the ones we are
interested in.
From the above it is clear that we need only consider the homogeneous problem in
order to determine the critical points at which the structure attains neutral stability.
5.3 Matrix Formulation for Truss Stability 105
[K(A)]{</>} = {oJ
where Ais some parameter (in our case associated with the critical load), special values
of which cause the determinant of [ K I to be zero. This is known as an eigenvalue
problem. We obtain these special values of A (called eigenvalues) by setting the
determinant to zero. Corresponding to each eigenvalue, we can find a solution for
{</>}; this is called an eigenvector or mode shape. For all truss problems, the critical
loads appear in the equations in a linear fashion and therefore we will deal with a
simpler version of the eigenvalue problem in the form
Since there are standard solution procedures available for this equation, we will even
approximate the frame analysis so as to cast it in this form.
Similar expressions can be determined for F 2 and V2 This is put into matrix form as
The first term is recognized as the augmented element stiffness for a rod and is referred
to as the elastic stiffness. The second term is associated with the change in geometry
of the structure and is referred to as the geometric matrix. For future reference, we
will write this as
106 Chapter 5. Structural Stability
In obtaining this relation we have considered that the transformation is applied only
to the vector displacement {u} and vector load {F}, and that the axial force Po is
left intact. We do this because the axial force is a member quantity and therefore
independent of our choice of global coordinates. The expression for [ k ] is a general
one and if the appropriate form for [T] is used it can be valid even for the 3-D
truss. As a special case, we get the total stiffness matrix for the plane truss in global
coordinates as
sc
C2 ]
-SC
C2
(5.6)
where, as before, the abbreviations C == cos (), S == sin () are used.
Assembling the global stiffness matrix is now a matter of adding all the element
stiffnesses. The only significant point is that during assemblage there is an initial force
Po for each member. Thus the assemblage occurs in two stages. First, we assemble
only the elastic stiffness and use it to find the initial force in each member by solving
[I<E]{U} = {p}
where [I<E] is the assembled elastic global stiffness. This corresponds to the pre-buckle
analysis. In the analysis of stability problems, we do not know the loading {p}; we
do, however, know the distribution of loads, therefore, let these be normalized as
P {p }. That is, all loads are assumed proportional to a single parameter P. Thus the
pre-buckle problem is solved using {p} as the load vector. From this, the axial force
in each member is obtained. Let the axial forces obtained in this way be called 10,
then the axial forces due to {p} are
In this manner the eigenvalue problem associated with truss stability can be estab-
lished as
[[I<El + P[loI<GJ] {u} = {oJ
where [lJ<G] is the assembled global geometric stiffness. The notation [JoI<G] is a
reminder that each element may have a different initial axial load. Indeed, some
members may even have zero axial forces. These points will be made clearer by the
following example.
5.3 Matrix Formulation for Truss Stability 107
Example 5.4: Find the critical load for the truss shown in Figure 5.5. All joints
are assumed pinned and both members have the same properties.
~ P
CD
We first solve the pre-buckle problem in order to get the initial axial force in
each member. Since the only non-zero degrees of freedom are {UI' VI}' the reduced
element elastic stiffnesses are
= EA
L
[00 0]1 ' [k*(13)] =2V2L
EA [1 11]
1
where we have ratioed the applied loads to P. The structural system of equations
are therefore
L~~D 1+~V2]{::}={~I}
Note that the scaling factor P does not appear. Solving this system gives the
displacements at Node 1 as
L
UI = EA'
The axial loads in each member can now be obtained from the stiffness relation or
by using the general formula
-("")
F'] = (EA)
L ij[coSOij(Uj-ui)+sinOij(Vj-Vi)]
and gives
Both of these values could have been obtained by inspection, but it is instructive to
go through the general proced ure.
108 Chapter 5. Structural Stability
[k*(12)j = EA
L 01
[0 0] _~LOO'
[1 0]
giving the assembled total stiffness matrix as
[K*] = EA
2.../2L
[1- P2.../2/EA
1
1]
1 + 2.../2
For the stability problem, we want to know what value of P causes buckling. This is
obtained by finding the value of P that makes the determinant of the total stiffness
matrix zero. That is,
P P
det =1- - 2V2 + 2V2
EA
- -EA 8- 1 =0
or
p _ EA
cr - (1 + 2.../2)
Note that only one value of the buckling load is obtained even though there are two
degrees of freedom in the system. Actually the second critical load is infinite.
The value of critical load we obtained is very large, and in a real structure it is
more than likely that either plastic yielding or a column type buckling would have
occurred. However, there are two circumstances when this type of buckling could
be significant. First, if the members forming the triangulated system form very
shallow angles the structure becomes susceptible to this type of failure. Second, as
the structure becomes more complex, made of a greater number of cells of about the
same size, then the load to cause buckling decreases.
Example 5.5: Find the critical load for the truss shown in Figure 5.5 but apply
the force P horizontally in the x direction.
We first solve the pre-buckle problem in order to get the initial axial force in
each member. It is straightforward to show that
fJ 12 ) = -P,
In this case both members have an axial component of force.
The total element stiffness is almost identical to the last example except for the
addition of the geometric stiffness of Member 1-2.
[k*(12)j = EA
L
[0010] _~LOO'
[1 0] 1]1 + 2LP[-11
giving the assembled stiffness matrix as
(1- V2A)(1 + A) = 0
This gives the two critical values of
Per = -EA/V2 or EA
The negative sign in the first value means that the horizontal force to cause buckling
acts in the negative x direction. Note that in this problem two finite values of the
buckling load are obtained.
Lx
.... 'Ir-=.~.:..-~.- trd_. )_1
... -....... .". .--.. V _(
x
~x
It is apparent that the axial force causes a bending action and therefore will enter
the beam equilibrium equations. Thus, balance of force and moment on the small
segment of Figure 5.6 gives (for small slopes and deflections), respectively,
dV
-=-q,
dx
110 Chapter 5. Structural Stability
The governing differential equation for the deflection shape is given by combining
these and using the moment-deflection relation to get
(5.7)
We shall refer to this equation as the coupled beam equation. When the initial axial
force Po is compressive only, then it is referred to as the beam-column equation.
All the relationships for the structural quantities may now be collected as
Displacement: v = v(x)
Slope: </> = dv (5.8)
dx
d2v
Moment: M = +EI dx 2 (5.9)
cFv - dv
Shear: V = - E I dx 3 + Fo dx (5.10)
d4 v _ d2 v
Loading: q = +EI dx 4 - F o dx 2 (5.11)
The only difference in comparison with the beam equations of Chapter 3 is the addi-
tion of the Po related terms in the expressions for the loading and shear. Thus with
Po treated as a parameter, it is seen from these that the transverse displacement v(x)
(via its derivatives) can be viewed as the fundamental unknown of interest.
When solving the coupled beam equations, information may be given at any of the
five levels above, thus requiring integrations (or differentiations) to obtain the other
quantities. Integration gives rise to constants of integration which must be found
from the boundary and compatibility conditions. In the general case, there are four
constants of integration (since the highest derivative is four).
We reiterate that for the purpose of integrating these governing equations, Po
is considered known from the pre-buckle analysis. As an example, to integrate the
loading relation, we will assume that the material properties EI, the axial force Fo ,
and the loading q, are all constant. Integrating twice gives
Example 5.6: Consider the cantilevered beam shown in Figure 5.7 with the
combination of axial and transverse forces. Determine the deflection shape.
t Q
I!I!II!II....--E-I..;..'-L--------,I- p
Let us first solve this problem in the manner of the previous chapters (this will
constitute the pre-buckle analysis). The axial behavior gives a displacement and
force distribution of
Px
u( x) = - E A ' F( x) = - P
That is, the axial analysis simply gives that the axial force is the same everywhere
along the beam and is P compressive.
Now consider the flexural behavior. Since the loading is zero, the general de-
flected shape is
v(x) = CIX 3 + C2x2 + C3X + C4
Imposing the boundary conditions of zero displacements and slope at one end, zero
moment and an applied force of Q at the other, allows the constants of integration
to be determined. Specifically, at x = 0, we have
o=
v=o
dv = 0
* C4
dx * 0= C3
At x = L, we have
d2v
M = EI = 0 + 6CIX]
dx 2 * 0= EI[2c2
d3v
V = _EI 3 = Q Q = -EIct
dx *
112 Chapter 5. Structural Stability
In essence, we have solved for the transverse displacements independent of the axial
behavior.
We will now solve the coupled beam equations. The applied loading is zero and
the axial force is compressive, therefore, the general deflected shape is
v =0 => 0 = Ct + C4
dv =0 => 0 = C2k + C3
dx
At x = L, we have
d2v
M = El = 0 => 0 = EI[ -Ctk2 cos kL - c2k2 sin kL)
dx 2
d3v
V = _El _ pdv = Q => Q = -Elk2c3
dx 3 dx
Solving directly gives (noting that P = k 2 EI)
Q sinH Q Q sinH
Ct =-----,
k 3EI coskL c3=-k2EI' C4= ----
k 3EI coskL
As a result, the deflected shape is determined to be
Example 5.7: Determine the effect of the initial axial load on the stiffness rela-
tion for the cantilever beam.
Since we are eventually interested in stiffness type relations, consider the above
solution for a particular point. For convenience, choose the end point x = L, then
Figure 5.8 shows a plot of the stiffness as a function of the initial axial load. The
axial load is seen to have a profound effect on the stiffness. Indeed, there are multiple
points when the stiffness is zero and other points where it is infinite (both positive
and negative).
To show that these results are consistent with the straight beam, choose P small
so that the parameter ~ == kL is small also, then
=3Q- [ e/3 .. ]
Q [~-e/6"'-~(1-e/2"')] QL3
v(L):::::-
k 3 EI 1- e /6 .. k EI 1 - /6 .. e 3EI
2.
o.
-2.
-4. L..L....L-..........L..~-'-..&......ol--'-.........................................L-&....1.-'-~--I-J
O. 2. 4. 6. 8. 10.
~ = kL = JPP/EI
Figure 5.8: Stiffness as a function of axial load.
The zero crossings in this example correspond to critical loads because at these
points the structure has complete loss of stiffness and is in a state of neutral equilib-
rium. That is, even for a nominal Q the deflections are indicated to be very large.
(An alternative view point is that any specified position can be obtained with very
little Q.) These critical points occur at
There are many critical loads. The minimum load to cause buckling is at n = 0,
hence
114 Chapter 5. Structural Stability
When we analyzed the truss, we saw that the critical loads were of the order EA,
but in the present case we see that
p
cr
= EA.!.-(~)2
A 2L
~ EA(!!:')2
L
The ratio (hi L) is called the slenderness ratio and relates the thickness h of the
beam h to its length L. When the length of a beam is 10 times it thickness, we see
that the critical load has already decreased to 1% of EA.
Example 5.8: Determine the deflected shapes corresponding to the critical loads
for the beam loaded as in Figure 5.7.
v(x)
v(x) I l"mOde~x
Figure 5.9: Deflected shapes for the first and second buckling modes.
The deflection at any load is obtained by substituting for kL into the expressions
for v(x) as
1
v( x) = k3~I cos kL [- sin kL cos kx + cos kL sin kx - kx cos kL + sin kL]
However, if we do this for the critical values kL = !(n + 1)11" we see that the cosine
term goes to zero and the deflection will become infinite. To overcome this, we will
normalize the deflection with respect to the tip deflection
v(x) = [sinkL :(~l coskL] [- sinkL coskx + coskLsin kx - kx coskL + sin kL]
This expression is valid for any value of kL, in particular, if we substitute the critical
values we obtain the deflection as
These are shown plotted in Figure 5.9 for the first two modes. It is important to
realize that the absolute position of any point on the beam is unknown because the
structure is in neutral equilibrium; however, the relative position of points (mode
shapes) take on very definite forms.
[ A ]{c} = {B}
Let the difference of the two solutions be {} = {c} - {c*}; it satisfies the homoge-
neous equations
[A]{}={O}
In general, this admits only the trivial solution, but when det [ A I is zero there are
other solutions.
From the above it is clear that we need only consider the homogeneous problem
in order to determine the critical loads. That is, we consider the problem
[A(,\)]{} = {oJ
where A is some parameter (in our case the critical load), special values of which
cause the determinant of [ A I to be zero. This is known as an eigenvalue problem;
actually, since [ A I contains transcendental functions it is usually referred to as a
transcendental eigenvalue problem. We obtain the special values of A (called eigen-
values) by setting the determinant to zero. Corresponding to each eigenvalue, we can
find a solution for {}, this is called an eigenvector.
116 Chapter 5. Structural Stability
v(x) ~:~_..:~.~~...
...-.....
!"'- - - - - -.-.,."'-,- : : - - - - - - , . . . . . . X
................ _ __ ..
1
Lx
==========~~ I/~------
lililllil=!
E I, L
P
JJ~~J~t~
v(x) .
x
Example 5.9: Determine the buckling loads for the pinned-pinned beam shown
in Figure 5.10.
Except for P the problem does not state any other applied loads, therefore, by
inspection get that Fo = -Po The boundary conditions are
at x =0 : v =0 ~ 0 = Ct + C4
d2v
M = El 2
dx
=0 ~ 0 = -Ctk2
at x =L : v =0 ~ 0 = Ct cos kL + C2 sin kL + C3L + C4
d2v
M = El 2
dx
=0 ~ 0 = -Ctk2 cos kL - c2k2 sin kL
We now inquire if the determinant of this matrix can be zero. On multiplying out
get
det = k 4 LsinkL = 0
There are many values of k when this equation is satisfied. The obvious one of k =0
corresponds to the trivial case of zero axial load. The other possibilities are
There are many critical loads, and corresponding to each there is a different deflected
shape. To determine these, let us reconsider the relation among the coefficients. At
5.5 Beam Buckling 117
0 C1
!]{
[ _k
12 0 0 C2
1 0 L C3 } =0
_k 2 0 0 C4
M(1) M(2)
I!I!II,L...)_ _E_I,_L
I ..-
'iffiiiM"
E_I,_L_ _ Jft!
mmm
( ~Dt)
V(1) l V(2)
We will assume that the spring does not exert any horizontal restraint on the
beam, and therefore conclude that the axial force is the same over the complete
length and given by Fo = -Po The spring, however, acts as a concentrated force (of
magnitude avo) hence the beam must be divided into two regions of integration.
For the first section we have at x = 0
v =0 ~ 0 = C1 + C4
d2v
M = El dx 2 =0 ~ 0 = -c1k2
For the second section, we have at x =L (we use x to mean the local distance
for the section)
where k 2 = PIEl is the same as for the first section. From these we get for the
deflected shape
It does not matter which set of coefficients we use for Vo ' From the first and third
of these equations we get
c'}=c2sinkL,
The remaining equations for both beam sections can be put into a reduced matrix
form as
kcoskL 2 -k]
-a*sinkL 2k 2 - a*L 0
[
sin kL cos kL 0 sin kL
where a* == al El. The determinant of this is
det = 4k3 cos kL sin kL + 2a*[sin kL - kL cos kLI sin kL =0
This equation cannot be solved explicitly for k. There are multiple modes but
computing them requires solving a transcendental equation.
The complete results for an arbitrary spring are shown in the Figure 5.12. These
were computed numerically. We see that increasing the spring stiffness causes the
buckling loads to increase. That is, the structure becomes less sensitive to buckling.
Example 5.11: For the previous example, consider the special limits of a very
flexible spring, and a very stiff spring, respectively.
5.5 Beam Buckling 119
~ = ";Pcr L2/E1
6. 0 0 0 0 0 0
5.
2nd sym. mode 0
4. 0
0
3.
l,t anti-sym. mode
2.
1. o two element results
O.
O. 10. 20. 30. 40. 50.
spring stiffness, aL2 / E1
Figure 5.12: Effect of spring stiffness on the first few critical loads.
sin kL = 0 or kL=1r,21r,
tan~ =~
Now make a plot of the functions h (0 == tan ~ and 12(0 == ~ for different values of
~ as shown in Figure 5.13. The intersections of these plots are at the roots of the
characteristic equation. The first few roots are
kL = ~ = 4.493,7.725,10.904
120 Chapter 5. Structural Stability
8.
6.
4.
2.
O.
-2.
-4. L.....&--'-.........~:..L- ........10..00.1..........L--'-.................l......JL.......l.---'-..........-......J
O. 2. 4. 6. 8.
~ = kL
Figure 5.13: Graphical scheme for solving transcendental equation tan ~ = ~.
For future reference, Table 5.1 gives results for some other boundary conditions.
The critical load is given by
z zEI
Per = 0'11" (kL) Li
The numerical parameter 0' is given in the table.
Table 5.1: Critical load factor for the first three buckling modes.
integration in terms of the nodal degrees of freedom. Using the displacement function,
we are then in a position to determine the member moment and shear distributions.
These, in turn, can be related to the nodal loads.
With reference to Figures 5.6 and 3.4, we assume that the loading is zero over
the element length and that the axial force is constant and compressive. The general
deflected shape is
The system of equations to determine the remaining coefficients can now be arranged
as
[ (1 - ~ (g ~ g~ ]{ ~: } = { VI4>~l~L4>~Lv2 }
where we have used the notations ( == kL, and C == cos (, S == sin (. We will use
Cramer's rule to solve this. First we get the determinant and rearrange it as
~ = ([2 - 2C - (S]
At this stage we can rewrite the deflection function just in terms of the nodal
degrees of freedom. From this, the moment and shear force distributions can be
obtained. For example, the moment distribution is
l
eS ((1 - C) -eS ((1 - C) ]
[ k ] = EI ((1 -~) -2((C - S) -((1 -~) 2(( - S) e (5.14)
3 -( S -((1 - C) (S -((1 - C) ~
((1 - C) 2(( - S) -((1 - C) 2(S - (C)
This is symmetric and exhibits many of the same repetitions as the beam element
stiffness.
Example 5.12: Solve for the tip deflection of the cantilever problem in Figure 5.7
using the exact stiffness matrix.
The procedure we follow is the same as for other matrix methods. Number the
nodes with 1 at the fixed end and 2 at the other, then, the unknown degrees of
freedom are
{uu} = {V2' 4>2}
There is only one element, hence the reduced structural stiffness matrix is just the
fourth quadrant of the element stiffness. The total stiffness relation is therefore
EI [ e5 -(L(1 - C) ] e { V2 } { Q }
3 -(L(1 - C) L2(5 - (C) ~ 4>2 = 0
The second of these equations gives
~ (1 - C)
4>2 = L (5 _ ~C) V2
Using this in the first equation and canceling the determinant term gives
QL 3 (5 - ~C)
V2 = EI f,3C
This is the same as previously obtained.
5.6 Matrix Analysis of Stability of Beams 123
~ - 2C - ~S]
~[2
~ ~[2 - 2(1 - e/2 + e/24 - ~6/720 + ...) - ~(~ - e/6 + e/120 _ ...)]
~ e[l - e
/15 + .. ]/12
It is important to realize that the series has been truncated precisely as indicated
above because that is the first non-trivial point at which the determinant can exhibit
a zero. We will need the reciprocal of the determinant and this is approximated as
1 12 ( 2
~ ~ ~5 1 + ~ / 15 +... )
We now do the expansion on the stiffness terms out to the same order. For example
EI [~ 4( ~
kn ~ D - ~3/6 + ...)] 12
~5 (1 + ~ 2/15 + ...) = EI [ -2~ / 10 + ... ]
3121
EI Po 12
kn = -[12] + -[-]
3 L 10
The first term is recognized as the kn of the beam element stiffness; the second term
is very similar to that of truss buckling. In like manner, we can expand for all the
stiffness terms to finally get the complete approximate element stiffness matrix as
12 6L -12 -36 3L ]
EI 6L 4L2 -6L 226L ] P. [36 3L -3L
3L 42 _L 2
[ k ] = 3
[
-12 -6L 12 -6L + 30~ -36 -3L 36 -3L (5.15)
6L 2L 2 -6L 4L 2 3L _L 2 -3L 4L2
124 Chapter 5. Structural Stability
This is symmetric, and the axial load appears in a linear fashion. We can make this
stiffness matrix resemble the results for the truss by writing
where [kE] is the element elastic stiffness matrix, and Fo [ kG] is the element geomet-
ric stiffness. Note that Fo may be positive or negative, the former means that the
structure gets stiffer with tensile loading.
If we compare the diagonal terms of the approximate total stiffness matrix with the
exact values, we see they go through a zero only once as Fo is made more compressive.
Therefore, at most, only four buckling loads are obtained for this element. Contrast
this with the infinite number obtainable from the exact stiffness matrix. However,
we can improve the approximate solution simply by using many elements for a given
member length.
Assemblage
Assembling the global stiffness matrix in stability problems is done in exactly the
same way as in the other chapters. The only significant difference is that during
assemblage there is an initial force Fo for each member. These axial forces are not
known; indeed, these are precisely what we are trying to determine in a buckling
analysis.
Following the procedure used for the truss, we assume there is a load or a group
of loads {p} applied and they are normalized as P {p }. We solve the pre-buckle
problem so that the initial axial force can be obtained as
for each member. In this manner the assembled matrices of the eigenvalue problem
can be established as
where [fJ<G] is the assembled global geometric stiffness. That is, we have extracted a
common P from each member, and now the problem resembles an eigenvalue problem.
The notation [lJ<G] is a reminder that each element may have a different axial load.
Example 5.13: Obtain an approximation for the buckling load for the clamped-
clamped beam shown in Figure 5.14. Use two elements.
The first step is to solve the pre-buckled equations in order to determine the
distribution of initial axial load. In this problem, Fo is the same throughout the
beam and simply equal to the applied axial load as Fo = -Po Therefore, we can go
directly to the eigenvalue problem.
Number the nodes as shown, then, the unknown degrees of freedom are
5.6 Matrix Analysis of Stability of Beams 125
II!Illt E/,L
Ef,JI-p
Figure 5.14: Clamped-clamped beam modeled with two elements.
The reduced element stiffness matrices are (using the fact that Fo =-P for both
elements)
[k*(12)] = E/
13
[12 -6L] P [36 -3L]
-6L 4L2 - 30L -3L 4L2
[k*(23)] E/ [ 12 6L] P [36 3L]
13 6L 4L2 - 30L 3L 4L2
In this case, there are two values of the critical load. Solving these gives the critical
loads as
P = 24 E/30 = lOE/ p = E/ 30 = 30E/
cr L272 L2' cr L2 L2
Compare these with the exact solutions (keeping in mind that the beam is of length
2L)
p _ 47r 2 E/ _ 9.87 E/ p _ 167r 2 E/ = 39.48E/
cr-~--v cr - 4L2 L2
The difference in the lowest critical load is less than 1%. The second corresponds to
the antisymmetric mode.
It must be remembered that even though there is no distributed load, the matrix
solution is an approximate one. Hence, to obtain a better solution it is necessary
to use more elements. Recall that for a given length of beam, this makes L of the
element smaller and therefore extends the range of validity of the approximation.
Further, the number of eigenvalues obtained is directly proportional to the number
of elements used. Again, this is equivalent to extending the range of the Taylor
series expansion. A rule of thumb is that if N eigenvalues are computed, then
approximately the first N/2 of them are reasonably accurate. Table 5.2 indicates the
convergence of the computed critical loads as a function of the number of elements
used; more or less the rule of thumb seems to be substantiated.
126 Chapter 5. Structural Stability
2 2 4.05 12.16
4 6 4.03 8.39 16.21 30.44 52.17 80.05
8 14 4.00 8.20 16.12 24.57 37.20 50.61
16 30 4.00 8.18 16.00 24.21 36.09 48.40
Exact 4.00 8.18 16.00 24.20 36.00 48.20
Vi = 0, V3 =0
giving the reduced system as
{llJbc} = {0,1;2,3;0,4}
A simple static analysis shows that the axial force is the same in each element
(and of compressive value P), hence we can proceed directly to assembling the elastic
and geometric stiffnesses. The reduced element stiffnesses are
[K*] = El
4L2 -6L
-6L 24 + a*
2L2
o6L
0] P
[ -3L
4L2
-3L
72 o
3L ]
13 [ 2L2 o 8L2 2L2 - 30L - L2 o _L2
o 6L 2L2 4L2 0 3L 4L2
5.6 Matrix Analysis of Stability of Beams 127
11111===-) -===--==--I~--==-----==---~:Ht-
EI, L EI, L flttiJ
P
~
Figure 5.15: Beam with an elastic support.
where a* == aLJ / EI. Note that the elastic stiffness [KE) changes by the presence of
the spring, but [KG) remains the same as in the previous example.
Generally, we solve for the critical loads by expanding the determinant. However,
the equation to be solved is then quartic in the load, so we will take a different tack.
Because of the geometric symmetry of the problem, we know that the results will
contain separate symmetric and anti-symmetric modes. We will simplify the problem
by considering each separately.
For the symmetric modes, we have that
The second critical value is hardly affected by the spring. The lowest critical load
increases linearly by the presence of the spring. In the case of a stiff spring, we have
EI
Per = L2 [12.7] , Per = ~~ [0.450:*)
The first critical value is not affected by the spring, whereas, the second critical
load increases to infinity. Indeed, if the spring stiffness is large enough this load
may exceed that of the second mode. What this means is that buckling is being
128 Chapter 5. Structural Stability
prevented from occurring in that mode. The full results are shown in Figure 5.12
where they are compared to the exact results.
For the anti-symmetric modes we have that
[[ 42 2]
4 - 30
oX [4-1 -1]] {</>1</>2-
4
} _0
where oX == P L2 / El. The spring does not enter these equations since the displace-
ment at the attachment is zero. The determinant of this system is
oX
2
- oX 72 + 720 = 0
This is solved to give
oX = 12 or 60
The critical loads are therefore
El El
PeT =V [12] , PeT =V [60]
What these results show is that for a flexible spring the first antisymmetric mode
lies between the first and second symmetric modes. However, as the spring stiffness
increases, it is possible for the anti-symmetric mode to become the lowest buckling
load.
where [KEl is the assembled elastic global stiffness and {p} is the normalized load
distribution related to the actual loading by {p} = P{p}. The initial problem is
then solved to give the axial force in each member (represented by 10)' The axial
forces due to {p} are therefore
Po = pIa
In this manner the eigenvalue problem can be set up as
Plane Frames
The total global stiffness matrix for an arbitrary frame element in two-dimensions
reduces to
C2
GS S2 sym
[k 1=
EA o 0 0
(5.16)
L -C2 -cs 0
-cs -S2 0
o 00 o
12S 2
-12CS 12G 2 sym
E1 -6LS 6LC 4L2
+ L3 -12S 2 12CS 6LS 12S 2
12GS -12G2 -6LG -12CS 12C2
-6LS 6LC 2L2 6LS -6LG 4L2
36S 2
-36CS 36C2 sym
Po -3LS 3LC 4L2
+ 30L -36S 2 36CS 3LS 36S 2
36CS -36G2 -3LC -36CS 36C2
-3LS 3LC -L2 3LS -3LG 42
where the abbreviations G == cos 0, S == sin 0 are used. The first two matrices are
recognized as the elastic stiffness for a plane frame member.
Example 5.15: Find the buckling loads for the plane frame shown in Figure 5.16.
Each member has the same material and section properties. The load of V2,P is
applied at an angle of 45 0 to the horizontal.
5.7 Stability of Space Frames 131
V2P~ Q)
p;;:::::::::::;.:!!II/I!
El, EA, L
CD fttJt
We will use two elements to model the problem. Numbering the nodes as shown,
the total degrees of freedom are
{u}={~,~,~;~,~,~;~,~,~}
UI = VI = </>1 = 0 and
giving the reduced system as
[k(*23)j = EA 1 0 0
0 00] + E: [00 12
0 6L
0]
L [ 0 0 0 L 0 6L 4L2
The reduced structural stiffness relation is therefore obtained as
[k~12)]=Plor(12) [36 0
0 0
3L]
0 [k*(23)] = PIar(23) [00 36
0 0] 3L
30L 3L 0 4L2 30L 0 3L 4L2
The reduced structural stiffness matrix can therefore be assembled. Introducing the
notations
_ JJ12) _ 11 23 ) a
,== 30T = 30L = (a + 12(3)L
we get the eigenvalue problem as
and
This is the symmetric mode. In the inextensible case (0' very large) we get
Per = 0.8330'L
which shows that the frame will not buckle in this mode.
The other two critical loads occur at
Problems
5.1 In reference to Figure 5.5, let the inclined member be at an angle of 8. Show
that the dependence of the buckling load on 8 is given by
p _ EAsin8cos 2 8
cr - (1 + sin38)
From this it is clear that as the inclined member becomes more vertical that
the buckling load decreases dramatically.
[Reference (41), pp. 145)
5.2 In reference to the previous exercise, at what angle 8 is the critical load 1% of
EA?
5.3 Show that for a cantilevered beam loaded by an axial force and concentrated
moment applied at the free end
v(X) = MoL2
EI
[1-k coskL
COSkX]
2
5.4 Recover the uncoupled solution from the previous exercise, and show that it is
in agreement with the straight-beam solution.
5.5 If the cantilever beam of Figure 5.7 has a tensile axial force applied, show that
5.7 Show that the lowest critical load of the second set of solutions of the previous
exercise is given by
5.8 Consider a cantilever beam with a rigid 'handle' of length a welded to the free
end. The handle is oriented along the length of the beam with its tip closest
to the fixed end of the beam. Show that if a force is applied to the handle and
pointing toward the fixed end of the beam that the characteristic equation for
the buckling load is
kL tan kL = L/a
[Reference [51). pp. 56)
5.9 If the force is reversed in the previous exercise so that the beam is in tension,
show that buckling occurs when
kL tanh kL = L/a
5.10 In reference to the spring supported beam of Figure 5.11, an interesting special
case arises when the spring has the special values a* = 2k 2 / L. Show that the
critical loads and the corresponding mode shapes are given by
Per
E1
= n 2 1r 2 V' v(x) = C2 [Sin n; x n1r 2~]
5.11 Consider a pinned-pinned beam of length 2L with a second axial load applied
at its middle. Show that the displacement shape in each beam segment can be
taken as
sin kL L -1 o L
[
k cos kL
Elk 2 sin kL
o
1
0
-Elk 3
o
-Elk*2
o
-k*
o
o
1
o
Elk*3 { c2
~~ } = 0
o 0 cosk* L sin k* L -1 c3
5.12 Explore some of the special cases of the last exercise (such as k* = k or P2 = 0)
and show that the reduce to the known results. An interesting situation arises
when P2 = - P3
5.13 A long homogeneous timber of rectangular cross-section [a X b] floats in water
with its top face horizontal. Prove that a necessary and sufficient condition for
the timber to be in stable equilibrium is 8 2 - 8 + b2 /6a 2 > 0, where 8 is the
specific gravity of the wood.
[Reference [27]. pp. 38]
Exercises
5.1 Plot the shear force and bending moment diagrams for a beam loaded as in
Figure 5.7.
5.2 Compare the kl l term from the approximate beam stiffness with its exact val-
ues. At what value of kL do they differ by 5%?
5.3 A simple beam carries a concentrated load at the center. If the beam is sub-
jected to a tensile force equal to three times the buckling load, by what per-
centage is the deflection reduced? [74.26%]
5.4 A steel tube oflength 4 ft has a 1 in outer diameter and a thickness of 0.036 in.
Find the load which may be applied axially so as not to exceed one-third of the
Euler buckling load. [544Ib]
5.7 A steel truss similar to Figure 4.5 has a length of 5 m and cross-sectional area
of 250 mm 2 Use one element per member to obtain the buckled mode shapes.
What would happen if more elements were use? [30.34 M N]
The previous chapters developed the analysis of particular structural systems by the
consistent use of the twin concepts of compatibility and equilibrium. In each case, we
derived a set of governing differential equations, integrated them, and then determined
the constants of integration by satisfying the appropriate boundary and compatibility
conditions. This approach will work in every case, but as a practical tool it suffers
from a number of drawbacks. Primary among these is that the solution can become
very cumbersome when more than one integration region is involved. It may also
happen that the governing differential equations cannot be integrated in a closed
form. When approximations are used the differential equations are not in a suitable
form for direct manipulation. For example, we obtained the approximate geometric
matrix in Chapter 5 by first obtaining its exact form, and then using a Taylor series
expansion to obtain the linearized version. But if the cross-sectional area varies, for
example, this approach is not feasible at all; it would be much more convenient if we
could go directly to the approximate solution.
The purpose of this chapter is to introduce alternative methods for establish-
ing equilibrium and compatibility requirements for structural behavior. Since the
equilibrium conditions involve forces, and the compatibility conditions involve dis-
placements, then a process through which both conditions are satisfied involves a
quantity dependent on both force and displacement. This quantity is called work,
and we will use it to show that, when a structure is in equilibrium, certain structural
quantities achieve a stationary value. This powerful idea will then form the basis for
our approximate analysis of structures. References [20, 27] are excellent introductions
to these energy approaches.
136
6.1 Work and Strain Energy 137
We will restrict the following analysis to structures for which the strains and
displacements are small and there is no dissipation of energy during loading. Although
the material behavior is elastic, it is not necessarily linear. A typical force-deflection
or stress-strain curve is shown in Figure 6.1. The elasticity requirement is that both
the loading and unloading paths coincide.
P,u
B
----_ ...
AW
L....... AU ~ u, f
Work
The work done by a force at a point is the scalar product of the force and the
displacement at the point. For example, in terms of our global coordinate system,
dW = F. du = Fxdu + Fydv + Fzdw
where the arrow indicates a vector. The force is understood to be constant during
the infinitesimal displacement duo When the force moves along a path from State A
to State B, the work done is
[8 [8 _ [8
~W = JA dW = JA F du = JA (Fxdu + Fydv + Fzdw)
The loading curve of Figure 6.1 can be interpreted as a sequence of possible equilib-
rium states as the load (or deflection) is changed. Thus, A and B are two possible
equilibrium states with different load conditions. We use the ~ symbol to signify
that an increment of work is performed in moving from one state to the other - the
change in configuration may be small or large. In the special case when the initial
configuration is the unstressed, unstrained, virgin state we will simply use W without
the ~ for the work done in reaching a certain state.
The systems we are interested in have multiple forces and moments, so we will
generalize the above expression for work to
(6.1)
E
yy
U yx
U yz
j.x
z
Strain Energy
If external forces {p} act on the structure, stresses are set up inside the body. Con-
sider a typical small volume of dimensions D.x, D.y, D.z, as shown in Figure 6.2. This
volume is under the action of the following system of stresses
Note that there are only six independent components of stress and strain, since the
shear components are related through relations such as
T xy = T yx , IXY = IYX
{t} = [V ]{u}
where [V 1is the matrix of partial differential operators and {U} = {u, v, w}.
A line that was originally of length D.x in the undeformed state has a new length
D.x(l + t xx ) in the deformed state. Now consider an additional infinitesimal deforma-
tion added to the body; a typical change in length is D.x(dt xx ). The work done by
the U xx component of stress during this change is
6.2 Linear Elastic Structures 139
where the symbol V represents the volume. Hence the work increment for the whole
body is obtained by allowing the typical volume to become infinitesimal and then
integrating over the total volume. That is
AW = lB i axxdtxxd:V
Because of perfect elasticity, the work expended can be regained if the loads are
gradually decreased. That is, the work is stored in the elastically distorted body in
the form of energy known as strain energy or elastic energy. We will use the definition
AU = lB i{aV{dt}d:V (6.4)
to mean the increment of strain energy. Although related, work and energy are
distinctly different concepts; we can say that forces perform work, but the system
possesses energy. Therefore, when work is performed on a system, a change of energy
occurs.
The formal equivalence of the two measures of work, Equation 6.1 and Equa-
tion 6.3, can be shown; we will take it as self evident. The principle of conservation
of energy can now be stated as
AU - AWe = 0 or (6.5)
where A We is the external work done by the forces. This says that for elastic systems,
the external work done is equal to the strain energy stored. It is to be kept in mind
that the changes in going from State A to State B could be large.
Strain Energy
Let the material obey Hooke's law in the form
1
t xx = E[(l + v)axx - v(axx + ayy + azz )]
140 Chapter 6. General Structural Principles I
1
t yy E[(1 + v)O'yy - v(O'xx + O'yy + O'zz)]
1
t zz = E[(1 + v)O'zz - v(O'xx + O'yy + O'zz)]
1 1 1
IXY = G t xy , IYz = G t yx , IZX = G 7zx (6.6)
where E is the Young's modulus, v the Poisson's ratio, and G = E/2(1 + v) the shear
modulus. The reciprocal of these for the normal components are
v
= 2G[f xx + -1--(t
+v
xx + f yy + fzz)]
v
= 2G[f yy + -1--(f xx + f yy + fzz)]
+v
v
= 2G[f zz + -1--(f xx + f yy + fzz)]
+v
These stress-strain relations can be summarized in the matrix forms
{t}=[c]{O'}, {O'}=[]{f},
(6.7)
The above relations will now be particularized to the structural systems of interest
by writing the distributions of stress and strain in terms of resultants. We have
already developed these structures in the previous chapters, so the assumptions and
restrictions as stated there are still assumed to apply here.
For the rod member, there is only an axial stress present and it is uniformly
distributed on the cross-section. Let F be the resultant force; then 0' = F / A = Ef
and
axial:
6.2 Linear Elastic Structures 141
For the beam member in bending, there is only an axial stress but it is distributed
linearly on the cross-section in such a way that there is no resultant axial force. Let
M be the resultant moment; then a = -My/ I = Et and
bending:
The shear forces in a beam can also do some work. Let the shear stress be assumed
to be uniformly distributed on the cross-section, and the resultant shear force be V,
then T = V/A = G, and
shear:
For a circular shaft in torsion, the shear stress is linearly distributed on the radius
giving T = Tr/J = G, and
torsion:
where </> is the twist per unit length, T is the resultant torque, G the shear modulus,
and J is the polar moment of inertia.
There are, of course, other types of structures, and an energy expression can be
set up for these also. However, they will all have a similar form. For instance, for the
cases considered above, there are the resultant loads
F, M, V, T
du d?v dv d</>
, ,
dx dx 2' dx dx
L(. lL
l
[load)2
energy =~ o stzf fness
) dx =~ 0
(stif fness) [deformation]2 dx
Note that even the general expression, Equation(6.7), follows this form.
142 Chapter 6. General Structural Principles I
Thus, C ll is the displacement at Point 1 caused by application of a unit load (in the
PI direction) at Point 1 (with all other loads vanishing); C 21 is the displacement at
Point 2 caused by the application of a unit load at Point 1. In general, Cij can be
interpreted as the displacement at Point i due to a unit load applied at Point j.
If the flexibility matrix [C] is not singular, i.e., det[Cij ] i- 0, then its inverse
[ C ]-1 exists. Let
then we obtain
N
{p}=[]{]{u} or Pi = L ]{ijUj (i = I,N)
j=I
The matrix [ ]{] is called the stiffness matrix. This is the same stiffness matrix that
has already been encountered in the earlier chapters.
P'l
Sequence 1
P2! Sequence 2
...:;~::;...... ....u2
For load sequence 2 on the other hand, P2 is applied first and then followed by Pl.
These are also shown in Figure 6.3. The work done by both forces during the entire
loading is
Wpl = [!C22Pi + (C 2I PdP2 ] + [!Cll PI2 ]
which is different from that given by load sequence 1. The total work done by PI
and P2 with load sequence 1 and with load sequence 2 should be the same, however.
That is,
That is, the flexibility matrix [C] is symmetric (and by inference, so also is the
stiffness matrix). This property is referred to as Maxwell's reciprocal theorem and
states that, in a linear elastic structure, the displacement at point i in direction nl
due to a unit force at point j in direction n2, is equal to the displacement at point j
in direction n2 due to a unit force at point i in direction nl. This is a special case of
some more general principles obeyed by linear elastic structures as described by the
fundamental matrices [ C ] and [1<]. We look at two such principles next.
144 Chapter 6. General Structural Principles I
Castigliano's Principles
Using the relations involving the stiffness and flexibility matrices allows the work and
energy expressions to be written alternatively as
U = ~L{LKikUk}Ui=~LLUiKikUk=Hu}T[K]{U}
i k i k
Thus, the work and energy are quadratic in the forces and displacements. These
expressions also show that the strain energy can be considered either as a function of
displacements, or of the loads. Use will be made of this in the following.
Consider the partial derivative of the strain energy with respect to the displace-
ments. For example,
where use has been made of the symmetry property of K ij . In general, it is possible
to write
p _ aU (6.8)
- aUi
This is Castigliano's first theorem and is essentially an alternative statement of equi-
librium. That is, it determines those forces that are consistent with a given compatible
displacement field. In a similar manner, consider the derivatives of the strain energy
with respect to the forces. That is,
where, again, use is made of the symmetry of Cij . In general, it is possible to write
aU
Ui = aP (6.9)
i
Castigliano's principles are our first indication that elastic structures obey certain
minimum principles. For example, consider a point on the structure where the applied
load is zero. Now imagine that this point is displaced in a compatible fashion by the
application of a force. As this is done, the strain energy of the structure changes (and,
of course, we expend work in moving the point), but it is clear that as we traverse
back through the original equilibrium point the force exerted is zero and we have
8U
Po = - =0
8u o
That is, the strain energy achieves a minimum. We will develop this concept more
fully in the next section.
It is also possible to show by further differentiation that
8 2U 8u 8 2U 8u
8 8
UI UI
= Lj I<Ij ~ )
UUI
= I<11 , 8 8
UI U2
= Lj I<Ij ~ )
UU2
= I<I2
Thus, in general, the elements of the stiffness matrix are related to the second deriva-
tives of the strain energy as
82U
I<;j = - - -
8u;8uj
This indicates that if the strain energy can be established in terms of nodal displace-
ments (degrees of freedom) associated with a compatible displacement field, then
the stiffness matrix can be obtained simply by differentiation. Indeed, this will be
demonstrated later in this chapter.
Example 6.1: Determine the total strain energy for the cantilever beam that is
loaded as shown in Figure 6.4. Then determine the vertical deflection and rotation
at the tip of the beam.
T.{\--;x----~!I!II!I ~t)
I
T V(x) M(x)
I
EI,L .....- x
This force system is in equilibrium since we obtained it from use of a free body
diagram. The strain energy is therefore
Comparing the various contributions of P l l it is seen that when the length of the
beam is large then the shear contribution is small. More specifically, the shear effect
can be neglected when
GA EI
L~ 13
Indeed, this assumption is inherent in the Bernoulli-Euler beam model. In most of
the work in this book, the shear effect is assumed to be small and usually neglected.
The deflections and rotations are obtained by use of Castigliano's second prin-
ciple as
au
VI = aPI'
That is, the deflections are determined by differentiating the energy expression.
Alternatively, it is usually more convenient to do the differentiating before the inte-
gration over the length. This will be illustrated here giving
{ 4>1VI} = 6EI
L [2L2
-3L
-3L] {
-6
PI }
TI
Example 6.2: Determine the deflection at a point a distance a from the tip of a
cantilevered beam. Neglect the strain energy due to shear deformation.
A force is not present at the location where we desire to know the deflection,
consequently we cannot use Castigliano's principle directly. A way around this is
to introduce a fictitious force at the point of interest (as shown in Figure 6.5) and
then use Castigliano's principle. Subsequently putting the force to zero will give the
solution. This is sometimes referred to as the dummy load method.
6.2 Linear Elastic Structures 147
The moment distribution is piece-wise linear, hence the integral for the strain
energy can be divided up as
Again, this force system is in equilibrium because we obtained the moment distri-
bution from application of the free body diagram. Note that if a = L then only PI
contributes to the energy; forces at fixed supports do not contribute to the strain
energy. On the other hand, if a = 0 then the energy is
It is apparent from this that the total strain energy of a combined loading system
is not the sum of the energy of the individual components.
We now use the general expression for the strain energy to obtain the deflection
at 2 as
V2 = -PI
1 (1-L 3- 1 2 + -a
-aL 1 3) = --PI
1 ( L - a)2( 2L + a )
EI 3 2 6 6EI
If P2 is retained then by differentiation of U with respect to both PI and P2 , we get
Now the relationship between deflection and force is quite clear; we have introduced
the concept of a flexibility at a point even if there is no concentrated force present.
148 Chapter 6. General Structural Principles I
The displacement distribution is piece-wise linear, hence the integral for the strain
energy can be divided up as
The force necessary to cause the deflection Uo is now obtained from Castigliano's
first theorem:
aU L
Po = auo = 2EAuo(a(L _ a))
This result is in agreement with what would be obtained using two rod elements as
done in Chapter 2.
applied loads. This subtle idea provides additional techniques for establishing equi-
librium and compatibility conditions; and will lead to some very powerful minimum
theorems.
J-.x
z
Now imagine that the particle is displaced a small amount ba and the forces move
with it, but with no change in the magnitude or the direction of any of the forces. It
is important to reiterate that we are not speaking of the actual displacement of the
particle caused by the action of the applied forces; we are speaking of an imagined
displacement in which the forces are imagined to behave as stated. For this reason,
the displacement ba is referred to as a virtual displacement to distinguish it from the
actual displacement. Throughout this chapter, the symbol b before any quantity will
indicate that it is a virtual or imaginary quantity.
The work done on the particle during the virtual displacement is the sum of all the
components of force in the direction of ba, times ba. This is called the virtual work.
Let the virtual displacement have components bu, bv, bw in the coordinate directions;
we have for the virtual work
By virtue of the fact that the particle is already in equilibrium (that is, the terms in
parenthesis are zero) we conclude that the virtual work is zero,
bW =0 (6.10)
150 Chapter 6. General Structural Principles I
Thus, if a particle is in equilibrium under the action of a set of forces, the total virtual
work done by the forces during a virtual displacement is zero. Note that the converse
statement is not true; just because the virtual work is zero we cannot conclude that
the body is in equilibrium. (This is easily demonstrated by a particle acted on by
a horizontal force but given a virtual displacement in the vertical direction.) It is
correct to say that if the virtual work is zero for every possible virtual displacement
then the body is in equilibrium, but this would be an impractical approach to testing
for equilibrium. We can simplify matters considerably by introducing the idea of
independent virtual displacements; an arbitrary displacement of a point can always
be written as the sum of three independent displacements. This leads to the much
more useful statement that if the virtual work is zero for every independent virtual
displacement then the body is in equilibrium. With this seemingly simple result we
have managed to interpret Equation(6.10) as a condition that defines equilibrium
rather than a result of equilibrium. Furthermore, some unknown forces that would
appear in a free body diagram (reactions, for example) will not appear in the virtual
work statement if they do no work.
This simple example identifies some important features for the general statement
of virtual work. First note that there are three states of the particle:
(A) an original state where there are no loads and the spring is unstrained,
(B) a final equilibrium state after the loads are applied and the spring has been
stretched an amount u, and
(C) an imaginary state where a virtual displacement ba is present.
We emphasize that the true displacement of the particle is the difference in displace-
ment between States A and B, and that the virtual displacement is the difference
between States Band C. Equation(6.10) can be viewed as either a consequence of
equilibrium or a definition of equilibrium. In the former case, the virtual displace-
ments bU, etc., need not actually be 'displacements' - they could be velocities or
pressure or just a purely mathematical function. In the latter case of defining equi-
librium, we must restrict the statement to independent virtual displacements - es-
sentially we refer to the generalized displacements or the degrees of freedom of the
body. A final point to note is that the condition that the forces do not change during
the virtual displacement also applies to the internal spring force ku, even though it
is clearly dependent on the real displacement.
P (1
~W* ~U*
8P
~W ~u
of the body. In State A, constraints are in place, but the body is unloaded. To keep
things simple, we consider the case where the constraints all impose zero displace-
ments and so the body is undeformed and unstressed. State B is the true state of
loading and deformation of the body; that is, the true loads are applied to the body
defined in State A and it is allowed to deform to its true state of deformation and
achieve a state of static equilibrium. Finally, State C is the virtual state of the body
obtained by assigning virtual displacements. It is in this state that the principle of
virtual work provides us with a suitable statement of equilibrium.
We take as given that the geometric constraints imposed in State A can never be
violated by the true displacements in State B or the virtual displacements in State C.
A deformable body has an infinite number of material points and hence has an infinite
number of degrees of freedom. We therefore must assume that all real and virtual
displacements are continuous (differentiable) functions of the coordinates x,y, and z.
Functions that have these properties are called kinematically admissible functions.
In the following development we will concentrate on a single component of force
P and a single component of displacement u; later we will generalize the results to a
system of arbitrary size.
With reference to Figure 6.8, the curved path can be viewed as the sequence of
possible equilibrium states as the load is increased. (Although we deal only with
linear elastic systems, it is more insightful to consider the non-linear elastic case at
this stage.) The states A and B can therefore be any two points along this curve.
Suppose B is relatively close to A then a Taylor series expansion about A gives
~W ~ bW + 1b2W + ... = Pbu + 1bPbu + ...
where we used the trapezoidal rule to compute the area ~W. A similar expansion
for the strain energy gives
~U ~ bU + 1b
2
2U +... = Jv[[(1bt + 1b(1bt
2
+ ... J dV'
The virtual strains are related to the virtual displacements by
au a
Dt xx = b ax = ax (bu)
152 Chapter 6. General Structural Principles I
and so on. Consider only the contributions of the first terms in the work and energy,
then
~W~<SW=P<Su,
<SP = 0, <Sf = 0
during the deformation. This is not physically realizable since the load-deformation
curves show a monotonic increase. For this reason, the small displacements <Su are
called virtual displacements, and the corresponding small strains <Sf. virtual strains.
For a deformable body the principle of virtual work may be stated as follows;
A deformable solid body is in equilibrium if and only if the total virtual work
is zero for every independent kinematically admissible virtual displacement.
This says that for every virtual displacement, a necessary and sufficient condition for
equilibrium of an elastic body is that the virtual work done by the external loads be
identically equal to the virtual elastic strain energy stored. The principle of virtual
work is as valid for establishing equilibrium as Newton's laws.
It is important to realize the difference between this principle and that of the
conservation of energy since both resemble each other closely. In essence, the latter is
concerned with total energies whereas the former is concerned with small changes or
what are called variations. As a result, virtual strain energy and work are calculated
on the assumption that the forces remain unchanged during the variation of the state
of strain.
The virtual work of a deformable solid body is divided into two parts as follows:
where <SWe is the virtual work of external forces and can be broken further into a
contribution <Sws of the surface forces and a contribution <SWb of the body forces.
These terms are given explicitly as
where the integrations are carried out over the surface and volume, respectively, of
the body, and the vector {U} contains the three displacement components u, v, w.
6.3 Virtual Work 153
Concentrated and distributed line forces on the surface can be taken as special forms
of these equations.
The other term bWi is the virtual work of the internal forces, in this case the work
due to the elastic straining of the body. That is, the internal virtual work is given by
In summary, the statement for virtual work can be put in matrix form as
bW = bWS + bW b - bU
= is {bUV {fS}dS + J) bU}T {l}eN -1 {bE}T {a}eN = 0 (6.13)
where we call bW* complementary work and bU* complementary strain energy cor-
responding to the virtual force bP and virtual stress ba, respectively. Virtual forces
are referred to as statically admissible forces because they must satisfy the conditions
of static equilibrium. We can now state the principle of complementary virtual work:
For an elastic structure in equilibrium, a necessary and sufficient condition
for compatibility of the deformation is that the complementary virtual work be
zero for any virtual force.
It is worth recalling that the principle of virtual work establishes the equilibrium
conditions when the deformation is already compatible, the complementary virtual
work principle establishes the conditions for compatibility when the loads are already
equilibrated. This duality can also be seen in Castigliano's principles. To show this,
apply the virtual work principles to a system with a collection of applied loads
If we now substitute these into the virtual work relations, and realize that the varia-
tions bUi and bPi are arbitrary, then we conclude that
Pi
au
=-,
au
aUi Ui = api
For linear elastic structures we have that U = U; the only difference is that we write
the energy as a function of displacement or load.
The matrix methods we have developed in the previous chapters are called 'dis-
placement formulated' because we always take the displacements as the fundamental
set of unknowns. In so doing we always enforce compatibility and use the stiffness
relation as the statement of equilibrium. Even when the displacements are approxi-
mate (as in beam buckling, for example) compatibility is assured. Hence, in the later
sections we will not make use of complementary work since it assumes equilibrated
load systems.
Example 6.4: Use the principle of virtual work to determine the forces in mem-
bers AB and AC of the truss shown in Figure 6.9.
2L
~---------."B
Consider a virtual displacement fJv at the point of application of the load. Since
this is a virtual displacement, we will specify that all other nodes remain in the same
position. Consequently, the new lengths of members AB and AC are (assuming small
deflections)
AB ~ v0.L + fJv/v0., AC~L
Consequently, the new lengths of members AB and AC are (again assuming small
deflections)
AB:::::: .../2L - cuj.../2, AC:::::: L - CU
The virtual strains in the two members are
The virtual work done by P in the virtual displacement cu is zero, hence we have
for all four members
or
We have already obtained a value for FAB, hence we conclude that FAG = P.
hW = hWe - hU = 0, hW =hW;-hU =0
The new concept introduced here is to see these relations (and consequently, equilib-
rium and compatibility) as the achievement of a minimization of certain quantities in
the structure. While this very powerful idea can be used to recover all of the previous
results, it turns out that its most fruitful application is in the approximate analysis
of structures.
6W = L B
Fdu = LB
-dV = -(VeE) - V(A)) or VeE) = V(A) - 6W
If State A is considered to be a zero state, then we see that we can recover work of
amount mgh.
Assume that we can obtain the external forces from a potential function; specifi-
cally, let it be a function of the displacements such that
aV
Pi = - - or {p} = - { - }
aV
aUi au
We will take as the potential for the external forces V = -{p)T{u}. The external
work term now becomes
aV
hWe = {p}T {hu} = -{ au}T{hu} = -hV
v= -Pu
n = ~Ku2 - Pu
These terms are shown plotted in Figure 6.10 for different values of displacement
u. It is apparent that n can achieve a minimum; the principle indicates that this
minimum occurs at the equilibrium position.
. . . . n
.......... u Ilil!l~
Ku= P
Note that the second derivative of the potential gives
This is positive thus confirming that the stationary value is actually a minimum in
this case.
The significant point of this example is that the minimizing of a potential can
be equivalent to deriving the system equations of equilibrium.
Example 6.6: Use the principle of stationary potential energy to determine the
equilibrium equations for a system of springs and applied loads.
The strain energy and potential energy of the applied forces are, respectively,
158 Chapter 6. General Structural Principles I
[K]{U}={P}
m m m
where [K(m)] is the augmented element stiffness. This relation shows that an al-
ternative view of the assemblage process is one of summing the strain energies of
the individual components. In the previous chapters we achieved assemblage by im-
posing equilibrium at the connections of the components, here equilibrium is taken
care of by our stationarity principle. We reiterate that this is possible only because
we have imposed compatibility between the components (by use of a set of common
nodal degrees of freedom {u}).
Boundary Conditions
When we apply the principle of stationary potential energy, we need to identify two
classes of boundary conditions, called essential and natural boundary conditions. The
essential boundary conditions are also called geometric boundary conditions because
they correspond to prescribed displacements and rotations. The natural boundary
conditions are sometimes called the force boundary conditions because they corre-
spond to prescribed boundary forces and moments.
To see the different roles played by the boundary conditions, we will reconsider
the simple problem of a rod with body forces and an applied load. We will show that
by invoking the stationarity condition on II, that the governing differential equation
of the problem and the corresponding natural boundary condition can be derived.
Example 6.7: The total potential energy and essential boundary conditions for
the rod shown in Figure 6.11 are
11111111L: :A,L
Figure 6.11: Rod with axial load.
where fb is the body force per unit length of the rod. Note that the highest derivative
in the potential is one.
The stationarity condition gives
8 (dU) d (8u)
8a dx = dx 8a
assuming that EA is constant, and using integration by parts, yields
i
L 2
-
o
[ EA-
d u
dx 2
+ fb ] 8u
-dx
8a
+ [dU
EA-lx=L -
dx
P ] -8uL - [ EA-Ix=o
8a
du
dx
] -8uo
8a
=0
To obtain the governing differential equation and natural boundary conditions, we
use the argument that 8u/8a is arbitrary at all points interior to the region, and
therefore the integrand must be zero, giving
82u
EA
8x 2 + fb = 0
It is also required that the other terms be zero, separately. Thus we must have
When an essential boundary condition is specified then its variation (with respect to
the state variable a) is zero, hence we see that this is the case in the first boundary
condition but not in the second, hence the boundary condition at x = L is
This is a natural boundary condition. We notice that this relation could also be
derived by considering a small free body diagram of the vicinity of the applied load.
160 Chapter 6. General Structural Principles I
Thus the natural boundary conditions are connected with the equilibrium of the
body.
Example 6.8: Use the stationarity of the total potential energy functional to
recover the differential equation governing the static buckling of the beam shown in
Figure 6.12.
::::;:::L El, L
111~~111t=:-x----------I--' oil p
~
Figure 6.12: Beam with buckling load.
There are three contributions to the strain energy, giving the total potential
energy as
L1 2V - I V 12
n = 10 2 El (ddx2 )2 dx + Fo 1L
0 2
(d)2
dx dx + 2 avL
Note that the axial force Po = -p does no work, but instead appears in the form
of a strain energy. The essential boundary conditions at x 0 are =
v = 0, dv =0
dx
Again, we will invoke the stationarity of n to derive the governing differential equa-
tions and also the natural boundary conditions at x L. =
The stationarity condition yields
10o Elv"'!.!!-dx
oa oa
El
0 oa
If we continue to use integration by parts, we eventually obtain
Since the variations on v and v' must satisfy the essential boundary conditions, we
have that at x = 0
vv= 0, VV ' = 0
va va
It follows that the third and fifth terms are therefore zero. The variations on v and
v' are arbitrary at all other points (including x = L), hence to satisfy bIT = 0 we
conclude that the following equations must be satisfied:
where the 9i are linearly independent trial functions, and the ai are multipliers to
be determined in the solution. The trial functions satisfy the essential (geometric)
boundary conditions but not necessarily the natural boundary conditions. Hence,
in the Ritz solution, there is error in the satisfaction of the differential equations of
equilibrium and the natural boundary conditions but we attempt to minimize this
error. By substituting the trial functions into n we can generate N simultaneous
equations for the parameters ai using the stationarity condition of n,
Example 6.9: Use the Ritz method to formulate the equations from which can
be obtained an approximate buckling load for the beam of Figure 6.12.
6.5 Ritz Approximate Analysis 163
The functional and associated boundary conditions governing the problem were
already given as
v(O) = 0, dv(O) =0
dx
oa2
all = 1 EI( 2a2 + 6a3x)2dx -1 P ( 2a2X + 3a3x2) 2xdx + a(a2 L2 + a3 L3 )L 2
oa3
all = 1 EI(2a2 + 6a3x)6dx -1 P ( 2a2X + 3a3x 2) 3x 2dx + a(a2 L2 + a3 L3 )L 3
From these, we obtain the system of equations
Not surprisingly, we end up with an eigenvalue problem. The solution of this eigen-
problem gives two values of P for which v( x) is nonzero. The smaller value of P
represents an approximation to the lowest buckling load of the structure.
Suppose we take as the displacement shape the following reduced function
[
E1(4L) + aL 4 (1)- 30
PL3
(40) ] a2 = 0
Notice that the reduced eigenvalue problem corresponds to considering only the first
entry in the [2 X 2) matrices. That is, as we increase the expansion of v( x), then each
164 Chapter 6. General Structural Principles I
additional term adds a row and column to the matrices but otherwise the existing
matrices are unaffected.
Example 6.10: Consider a non-uniform rod fixed at one end and subjected to
an axial concentrated force at the other end, as shown in Figure 6.13. The variation
of axial stiffness is EA(x) = EA o (1 + x/L)2. Obtain a Ritz approximate solution,
and compare with the exact solution.
u(x)
F(z) 1-_____ _ x
The exact solution is easily calculated using the methods of Chapter 2, and gives
(X IX P PL x/L
u(x) = Jo f(x)dx= Jo EA o (1+x/LF dx = EA o (1+x/L)
The exact force distribution in the bar is
du
F(x) = EA dx =P
This last result is in agreement with what would be obtained from the equilibrium of
simple free body diagrams. Both the displacement distribution and force distribution
are shown plotted in Figure 6.13. We will use these results to evaluate the quality
of the Ritz approximate solutions. Specifically, we wish to investigate the use of
different trial functions.
The total potential energy of the structure is
II =1 {L EA (du)2 dx - PUL
2 Jo dx
We will calculate the displacement and force distributions using the following as-
sumed form for the displacement:
This must satisfy the essential boundary condition, hence ao = O. Note that the
remaining polynomial does not necessarily satisfy the natural boundary condition.
Substituting the assumed displacements into the total potential energy expression,
we obtain
Invoking the stationarity of II with respect to the coefficients an, we obtain the
following equations for al and a2 after differentiating
E A o [70L 852] { al } {P L }
30 85L2 1243 a2 = P L2
Note that this is symmetric. Solving this system gives for the two coefficients
78P 10 2
u(x) = 97EAo1x - 26Lx ]
Substituting this into the potential energy expression and minimizing, we obtain
166 Chapter 6. General Structural Principles I
displ force
x/L exact 1 term 2 term hi-linear exact 1 term 2 term hi-linear
0.0 0.0 0.0 0.0 0.0 1.0 0.428 0.804 0.6316
0.5 0.3333 .2143 .3247 .3158 1.0 0.96 1.113 1.421
0.5 0.3333 .2143 .3247 .3158 1.0 0.96 1.113 0.729
1.0 0.5 .4285 .4948 .4779 1.0 1.714 0.740 1.297
Tahle 6.1: Displacement and force results for the non-uniform rod.
EA o [70LJaJ = PL
30
We recognize this as precisely the first term in the above [2 X 2J matrix form. That
is, as we increase the expansion of u( x), then each additional term adds a row and
column to the matrices but otherwise the existing matrices are unaffected.
Solving for al gives
3 P
al = ---
7 EA o
The approximate solution for the displacement and force are, respectively,
3P
u(x) = 7EA o [x], F(x) = 3: [IJ(1 + xl L)2
These results are also shown in Table 6.1 as the I-term columns. Note that this
force distribution does not satisfy the differential equation of equilibrium.
Example 6.11: As a second Ritz solution to the rod problem of Figure 6.13,
assume that the displacements are given in a piece-wise linear form as
2x
u(x) = yU 2 0$ x $ LI2
where U2 and U3 are the displacements at points mid-way and the end of the rod.
These displacements satisfy the essential boundary condition at x 0, and also the =
continuity of displacement condition at x = L12. There is no continuity of the first
derivative duldx at x = L12, but that is permissible since the highest derivative in
the potential is duldx.
Using these trial functions in the potential energy gives
IT = ! 1
L 2
/ EA o (1 + xl L)2 C~2) 2
dx +! 1:2
EA o (1 + xl L? ( _ + 2~2 2~3) 2
dx - PU3
In this case, the displacements U2 and U3 are the state variables (degrees of freedom).
Invoking that IT is stationary with respect to these, we obtain
EA o [56 -37] { U2 } ={ 0 }
6L -37 37 U3 P
6.5 Ritz Approximate Analysis 167
These results are shown in Table 6.1 as the bi-linear columns. Note again that the
displacements are quite accurate, but the forces are significantly off. Indeed, there
is not even equilibrium at the joint.
This last example is very important in demonstrating the relationship between a
Ritz analysis and a finite element analysis. Indeed, a Ritz analysis can be regarded as
a finite element analysis and vice versa, and we make use of that in the next section
to derive some elements.
At this stage, it is worthwhile to summarize some of the characteristics of the Ritz
method. The most important ones are:
Usually, the accuracy of the assessed displacement is increased with an increase
in the number of trial functions.
While fairly accurate expressions for the displacements are obtained, the corre-
sponding forces may differ significantly from the exact values.
Equilibrium is satisfied in an average sense through minimization of the total
potential energy. Therefore, forces (computed on the basis of the displacements)
do not, in general, satisfy the equilibrium equations.
The approximate system is stiffer than the actual system and therefore buckling
loads and vibration resonances are overestimated.
Since the Ritz solution is approximate, sometimes it can be confusing to know
what check conditions the solution should satisfy. It often helps to realize that the
approximate solution is actually the exact solution to some other problem. For ex-
ample, in the bi-linear approximation of the rod, we could think of the solution as
replacing the original problem with a piece-wise constant rod, each section of length
L/2 and of stiffnesses
19 EA 37 EA
6 0, 6 0
respectively. This problem now resembles a typical one from Chapter 2 and thus
we expect conditions such as equilibrium at the joints to be satisfied. Indeed, when
viewed in this manner, we find that the axial force is constant and given by F(x) = P.
Thus a safe interpretation of the Ritz approximation is that we replace the original
system with a new physical system comprised of simpler elements; the effective prop-
erties of these elements are obtained by use of the principle of stationary potential
energy.
168 Chapter 6. General Structural Principles I
We can view the finite element method as an application of the Ritz method, where
instead of the trial functions spanning the complete domain, the individual functions
span only subdomains (the finite elements) of the complete region. Figure 6.14 shows
an example of a bar with a hole modeled as a collection of many triangular regions.
The use of relatively many functions in regions of high strain gradients is made pos-
sible simply by using many elements as shown around the hole in the figure. The
combination of domains with different kinds of strain distributions may be achieved
by using different kinds of elements to idealize the domains.
In order that a finite element solution be a Ritz analysis, it must satisfy the es-
sential boundary conditions. However, in the selection of the displacement functions,
no special attention need be given to the natural boundary conditions, because these
conditions are imposed with the load vector and are satisfied approximately in the
Ritz solution. The accuracy with which these natural boundary conditions are satis-
fied depends on the specific trial functions employed, and on the number of elements
used to model the problem. This idea was demonstrated in the last example where
it was shown that improvement in the approximate solution could be obtained either
by increasing the number of functions, or by increasing the number of sub-domains.
6.6 The Finite Element Method 169
Element Stiffnesses
In order to develop the element stiffness relation we view the element itself as a
structure with essential boundary conditions at the nodal points (these are the degrees
of freedom). An application of the stationarity principle then gives the approximate
equilibrium conditions that the nodal loads satisfy.
Let the distributed displacement fields {U(x, y, z)} be represented by an expression
of the form
u(x,y,z) = Lg(x,y,z)J{u}
with similar expressions for the other displacements v(x, y, z) and w(x, y, z). Here
Lg(x, y, z)J is a set of known admissible functions of the coordinates, and {u} is a
set of constants. Let us insert this into the expression for the potential energy for a
general elastic body
where [ B ] == [ 1) ][ 9 ] and
[k){u}={F}
This is exactly the same set of equations we derived in the early chapters using the
direct approach. What is different is that we can talk about the accuracy of a given
choice of functions and the convergence of the solutions when more admissible func-
tions are added. Generally, adding a term will reduce the value of the potential energy,
thus improving the answer. At worst, an admissible but otherwise inappropriate term
(for example, an antisymmetrical function in a solution which is clearly symmetrical)
will have little or no effect because the Ui for that term will prove to be very small or
even zero.
170 Chapter 6. General Structural Principles I
We are reminded that in the Ritz method, we are replacing an infinite degree of
freedom system by one which has finite degrees of freedom. This tends to overestimate
the stiffness of the system so stress, strain, and deflection err on the side of a stiffer
structure.
This displacement function satisfies the essential boundary conditions in that it gives
the displacements at both nodes. We therefore identify 11 (x) and h( x) as the Ritz
functions, and Ul and U2 as the state variables. The normal strain in the element is
EA(
U = 21 JfL
o
(dU)2
EA dx dx = 2 U2 - ut}
2
which is the stiffness matrix for the rod element with respect to the local coordinates.
We will now use the general expression to derive the general form of the rod
stiffness. From the shape functions we get
d
[B] = [V ][g(x)] = [dxllf1(x) , h(x)j = If;(x) , f~(x)j
k jj = i EAfJ(x)fj(x) dx
The corresponding element nodal forces are
Note that if q(x) is replaced by concentrated forces at the end and a distributed load
along the length, then we recover the usual nodal forces, but in addition, we get the
rule for treating distributed loads. This example shows the natural manner in which
the approximate stiffness can be obtained for complicated problems.
This satisfies the essential boundary conditions. The total strain energy stored in the
beam element is (neglecting shear)
II = UB + Uc + V
The entities in the total stiffness matrix can be obtained by extremizing the total
potential to give
r r
L L
kij = !:J82~ = EIg:'(x)gj'(x)dx+ Fog:(x)gj(x)dx
UUiUUj Jo Jo
Carrying out these integrations under the condition of constant EI and constant Fo
gives the beam element stiffness matrix and the beam geometric stiffness matrix as
already obtained. In the present form they are applicable to non-uniform sections.
Note that it is the symmetry of the terms gi'( x )g'/(x) and gi( x )g; (x) that insures the
symmetry of the stiffness matrices.
We have already shown that the condition 6II = 0 governs equilibrium; a study of
the second order terms 62 rr therefore governs the nature of the equilibrium, that is,
Actually, if 62 II = 0 then we must check the higher order terms also. We conclude
that for stable equilibrium the potential energy is a minimum.
To clarify the meaning of the second variation, let us reconsider the example of
Figure 5.1. The total potential for the problem is
II = !ku 2 - Pv - Qu
6.7 Stability Reconsidered 173
Since we assume the vertical bar to be rigid, then there is a geometric constraint that
relates the displacements u and v. Specifically, we have
u2
V = L - "/2 - u 2 :=::::-
2L
This last approximation is consistent with our assumption that all deflections are
small. Consequently, the potential for the problem becomes
Note that while P originally appeared in the work term Pv, it now appears in a strain
energy like term. The different orders of variation of the potential are
bil =
all
au bu = [ku - PulL - Q]bu
2
1a Il bu 2 = l[k _ PI L]bu 2
2 au 2 2
= 0
The first equation gives the equilibrium condition; when the system is in equilibrium
this term is zero, therefore it is the second equation that determines the sign of ~Il.
Depending on the value of P this second term can be either positive or negative. What
is most important to note is that the relation does not contain the other applied load
Q. In fact, the potential of the applied loads is always linear in displacement, hence
they will never contribute to the second variation of the potential and thus affect
stability. (However, they do affect the initial load term Fo .)
Note that the second variation in potential energy is related to the total stiffness
as
Hence the conditions for stability of equilibrium can be applied to the total stiff-
ness matrix; that is, we inspect whether this matrix is positive definite, positive
semi-definite, or negative definite. A necessary and sufficient condition that b2 Il be
positive definite is that the determinant of the matrix [I<] and all its principal mi-
nors be positive. Usually, it is unnecessary to examine a sequence of minors since the
determinant itself vanishes before any of its minors. Thus our criterion
det[ I<] = 0
Problems
6.1 A uniform rod is held between two rigid walls. Determine the strain energy
and complementary strain energy stored in the rod as it is heated.
[Reference [10]. pp. 87]
d2 u
dx 2 + 100 = 0 o~ x ~ 10
6.3 Obtain a Ritz solution to the previous exercise using the trial function f( x) =
x(lO - x). Compare this approximate solution with the exact solution.
6.4 In general, it is not possible to obtain a functional for problems whose governing
differential equation contains odd-power derivatives. A special exception is the
following case
d2 u du
dx 2 + a dx + bu = 0
where a and b are constants. Show that the functional
6.6 Suppose we wish to derive a "higher-order" rod element and to that end we
take the deflection function in the cubic form
where the functions 9n( x) are identical to those for the beam shape functions
of Equation(3.6). Show that the higher-order rod element stiffness matrix is
36 3L -36 3L ]
EA 3L 4L2 -3L -L2
[ k ] = 30L -36 -3L 36 -3L
[
3L _L2 -3L 4L2
Note that this is almost identical to the geometric stiffness matrix for the beam;
Why?
[Reference (49). pp. 207)
6.7 Suppose that only the axial forces F 1 , F2 are taken as the nodal loads in the
previous exercise, show that the elementary rod stiffness relation is recovered.
Propose some nodal loads that do not give the trivial result.
6.8 If in the derivation of the beam geometric stiffness matrix, we use the rod shape
functions instead of the beam shape functions, i.e.,
[ kG ] = F
10-1 0]
[
o 0 0 0 0
L -1 0 1 0
o 0 0 0
Where might such an element be more useful than its consistent counterpart?
[Reference [46]. pp. 319]
6.9 With reference to Figure 6.6, use the Ritz method to find an approximate
solution using
Exercises
6.1 A particle in the (x, y) plane moves in the force field Fx =
-ky, Fy = kx (k
is a constant). Prove that when the particle describes any closed path in the
counterclockwise sense, the work performed on the particle is 2kA where A is
the area enclosed by the path.
176 Chapter 6. General Structural Principles I
6.2 With reference to Figure 4.14, verify Maxwell's reciprocal theorem for a hori-
zontalload at Node 1 and a vertical load at Node 4.
6.3 Use Castigliano's second principle to find the reactions for the rod in Figure 6.6.
The procedure is to replace the fixity constraints by a system of unknown forces
and then use the principle to determine the displacements at the end points.
By setting these displacements to zero a sufficient number of equations are
obtained to determine the reactions.
6.4 Consider a cantilever beam, fixed at the end x = 0 and subjected to a given
displacement VL = c at the other. Show that the following is a set of admissible
displacements and obtain the corresponding Ritz solution.
6.5 Consider a cantilever beam, fixed at the end x = 0 and subjected to a concen-
trated lateral applied force at the other. Using the Ritz method, show that the
following displacements
leads to the exact solution. Show that the addition of extra terms have zero
contributions.
6.6 A uniformly loaded beam is simply supported at both ends. Find the deflection
and bending moment at the center using the Ritz method. First use a quadratic
function in x and then use a sine function in x. Compare the results with the
exact solution and say why the second solution is better than the first.
6.7 Suppose the rotation distribution in a circular shaft can be written in terms of
the nodal rotations as
6.8 For a uniform column with clamped ends, assume v = ax 2( L - x)2 and deter-
mine the critical load. [Pcr = 42EI/ L2]
6.9 For a uniform column with clamped at one end and free at the other, assume
v = ax 2 and determine the critical load. [Pcr = 2.5EI/ L2]
6.10 For a uniform column with clamped at one end and free at the other, assume
the mode shape is the static deflection shape due to a point load at the tip and
determine the critical load. [Pcr = 42EI/17L2]
Chapter 7
Computer Methods I
The use of the computer is essential for modern structural analysis; therefore, it is
important that the engineer have some understanding of how the computer is actually
used to accomplish this analysis. This chapter introduces some of the basic computer
methods and algorithms used in structural analysis; these include schemes for solving
simultaneous equations, solving eigenvalue problems, and for efficient data storage.
In addition, it attempts to describe the computer environment in which the analysis
takes place. References [4, 32, 45] consider many of the general aspects of using
computers for structural analysis.
It is difficult to describe computer algorithms without describing the complete
programming context. Therefore, most of the discussion will refer to the algorithms
implemented in the program STADYN; the source code is listed in Appendix B. The
availability of this program also gives the opportunity to discuss some of the design
aspects of a large structural analysis program.
177
178 Chapter 7. Computer Methods I
RAM available affects the program in two ways; first, in the accuracy of the numbers
computed, and second in the size of problems that can be solved. Much of the art
in programming revolves around minimizing the number of computations and the
memory requirements.
Data Representation
To get a clearer idea of how RAM is used to store data, consider a floating point num-
ber 123456.78 represented in the form .12345678E6. We call the portion .12345678
the mantissa and the portion E6 the exponent. Numbers such as these are represented
in a computer in terms of bits; a bit is the smallest unit of memory and has a value
of either 0 or 1. According to the IEEE standard [32], the standard real variable
(float) is stored as a 32 bit (4 byte) word. In FORTRAN, this is referred to as single
precision. According to the standard, 23 bits are reserved for the mantissa resulting
in approximately 2310g 10 (2) = 6.9 decimal digit accuracy. Eight bits are reserved for
the exponent and this gives a dynamic range of about 2255 = 1O38. The remaining
bit is the sign bit.
It is apparent that because of the manner in which numbers are stored in the
computer, only a finite set of real numbers can be represented. This gives rise to
side effects that would not be present with exact arithmetic. One such effect is that
if a computation obtains a number outside the dynamic range, for example, 10+ 40
or 10- 40 then an overflow or underflow error occurs, respectively. As an illustration
of how easy it is to cause an overflow, consider computing the determinant of the
stiffness matrix. If the matrix is diagonal the determinant is the product of all the
diagonal terms. Since each of these terms is of order 106 in imperial units (power
11 in metric units), then a system size of only 4 to 8 is sufficient to cause overflow.
Obviously computations like this must be done with some circumspection; usually
this means scaling the numbers before multiplication. Another consequence of the
finite representation of numbers in the computer is that every operation (addition,
multiplication, and so on) results in loss of accuracy due to round-off error. It must
be kept in mind that these errors are quite separate from errors due to modeling.
To help alleviate both of these problems, programs may make use of higher precision
representation of numbers. For example, the double precision model represents float-
ing point numbers as 64 bit (8 byte) words. This results in approximately 16 decimal
digit accuracy and a dynamic range of 10308. However, each number requires double
the amount of computer storage.
It appears that all these difficulties may be overcome by making more RAM avail-
able. Sometimes the RAM limit is due to the operating system; this is the case with
the DOS limit of 640kB on the PC microcomputer. More often the limit is due to
the cost of memory. Thus, for example, workstations running Unix often have from
4 to 8 MB of RAM even though the operating system can handle much more. Some
programs (and operating systems as in the case of Unix) can give the effect of having
7.1 Computers and Data Storage 179
large memory by using the disk drive as a virtual memory source. However, the speed
with which data can be accessed is the critical performance factor of an external stor-
age device; these devices are always much slower than internal memory, and thus their
use is usually restricted. For example, this approach is adequate for static problems
but can be very slow for dynamic problems because of the very frequent writing and
reading of the disk.
Matrix Storage
A matrix is banded if all nonzero coefficients cluster about the diagonal. The stiffness,
mass, and geometric matrices encountered in structural analyses are often of this type.
Banded storage is a simple way to exploit matrix sparsity because zeros outside the
band need not be stored or processed.
# + + + # + 0 0 + +
+ # + + + # 0 + + + 0
# + + + + # + + + + 0
+ + # + + # + + 0 0 0
+ + + + # + + + # + 0 0 + +
+ + + + + # + + + => # 0 + + + 0
+ # + + +
+ # + + + + 0
+ + # + + # + + 0 0 0
+ + + + # + # + 0 0 0 0
+ + + + + # + # 0 + 0 0 0
+ # + # + 0 0 0 0
+ + # # 0 0 0 0 0
Full Storage N X N Band Storage N X B
A simple storage format for a banded symmetric matrix is shown above. The N
principal diagonal elements are stored in the first column, the (N - 1) elements of the
next diagonal are stored in the first (N - 1) positions of the second column with zero
in the last position, and so on for the diagonals which have non-zero elements. The
minimum number of columns needed is called the semi-bandwidth, B. It is six in the
above case. The total bandwidth includes coefficients to the left of the diagonal and
is therefore 2B - 1. The entire information content of a symmetric matrix resides in
the [N x B) matrix. In practice, the bandwidth B may be significantly smaller, by
a factor of ten or more, than the matrix order N, so there is merit in storing just
N B coefficients rather than all N2 matrix coefficients. Also, as will be seen in the
next section, the equation-solving expense decreases by a factor of about B2/ N2 if
we operate on only N B coefficients instead of on all coefficients in the upper triangle.
Because of the fewer operations, the solutions are more accurate.
Other schemes for efficient data storage are available and a clear discussion of
them can be found in Reference [46].
180 Chapter 7. Computer Methods I
IP"-Pme",,o,
Structure
DataFile
X,Y,Z
I connectivity
material
loads
. Bes
I'
t oo ......... O' . . _ _ . . . . . _ .. __ .....
Processor
.
.. : {p}
[ J( ]
{u}
F(x), M(x), ...
Post- Processor
u(x), v(x), ...
In many large scale problems, the volume of possible output can be gigantic.
The ability to interpret the results quickly and use them effectively is of critical
importance. The post-processor manages and presents the results of the computations.
Typical tasks are the plotting of deformed shapes, or the contours of constant values
of stress. Sometimes error estimates are provided with recommendations for mesh
refinement. The post-processor may also have self-adaptive features incorporated
with appropriate feed-back to the pre-processor and processor units.
Design Considerations
In any design project compromises must be made between competing goals; in this
respect, the design of a large scale computer program is no exception. Ideally, we
would like a program that is blindingly fast, can have access to as much memory as
needed, utilizes every hardware enhancement available, and gives output in just the
right form. This is not to be achieved (yet) in either the mainframe or micro computing
world. The following specific guidelines were used to coordinate the design decisions
when writing STADYN:
STADYN is to be a special purpose program for the static and dynamic analysis
of 3-D frame structures.
Dynamic problems are emphasized over static ones.
The essential data must be kept in-core during the computationally intensive
stages.
A disk drive may be used as a virtual memory source, but only for secondary
and seldom used variables.
STADYN is to be developed as a processing engine; therefore it should have
hooks to pre- and post-processors.
The program must run reasonably well on a variety of platforms from
micro to mainframe computers.
All modules must be written in a high level language, in a style that is
portable and easily maintained. The modules are to be re-usable for other
special purpose programs.
These guidelines, to a greater or lesser extent, shaped the choice of the algorithms to
be described later. The resulting design has proven to be robust and quite suitable
for microcomputers.
It is emphasized that all programs (intentionally or not) are developed with par-
ticular guidelines in mind and this can profoundly affect their suitability for certain
applications. For example, a program designed to handle very large static problems
may make copious use of a disk drive as a secondary memory source; this program
would probably perform very poorly if applied to a dynamic problem. On a more
global level, some general purpose programs are designed to make available to the user
the complete library of elements and routines; these programs demand vast computer
182 Chapter 7. Computer Methods I
resources and usually are restricted to mainframe environments. On the other hand,
specialized programs can target specific applications (structural analysis of buildings,
say) and thus provide performance, problem capacity, and accuracy exceeding that of
the general purpose programs. Invariably, these programs can also run in a smaller
computing environment.
The point is that no single program is the 'best' for all situations, and the user
must be aware of the inherent limitations in each program and not attempt to use it
inappropriately.
...........
Dataln Mat:
{Be}
......................
DataCnv
BC I 1
..~:.!.. j
ElmStif I [K] :
MAIN
ASSEMble
o
{u}
................o
UDU I
SOLVEr 1
o
Bak
u(x)
................o
0
Disp 0
:
OutPut F(x) ~
Member ~ ...... _...... ~
freedom, and the boundary conditions imposed, in module DataCnv. The global
system of equations must be rearranged and renumbered to account for this. It is
at this stage that STADYN does any special global reductions such as restricting the
frame to be 2-D only. The book-keeping for the reduced degrees of freedom is taken
care of in the array {IDbc}. This is a vector array whose elements identify each node
(in sequence) and the six degrees of freedom at each node. The constrained degrees of
freedom contain a 0 while the non-zero ones contain the equation number associated
with the global system of equations.
During the input of the structure data file, every element is scanned and the
largest nodal number difference (for nodes with non-zero degrees of freedom) is taken
as the bandwidth. Thus on exit from the input routine, the storage requirements
are known even though the stiffness matrix has not yet been assembled. The module
ASSEMble forms the stiffness matrix. For each element, the local element stiffness
is established and then transformed to global coordinates. This is entered into the
structural stiffness matrix by associating the appropriate nodal numbers. Those nodes
which have a zero degree of freedom are not entered at all. After assemblage, the
stiffness matrix is stored to disk in compact form ready for later re-use.
The solution of the system of equations is accomplished in module SOLVEr. This
first decomposes the stiffness matrix and then uses forward and back substitution to
obtain the solution. The obtained nodal degrees of freedom are stored to disk for
use by other modules. The routines in module OutPut convert the nodal degrees of
freedom into other informational forms (such as member loads). Members can be
184 Chapter 7. Computer Methods I
Step 1: Pick a reference point in space, typically away from the structure.
Step 2: Calculate the distance of each node from the reference point.
Step 3: Order the nodes according to their proximity to the reference point.
This is a very effective scheme when the mesh is relatively uniform as is the case in
2-D and 3-D continuum problems.
7.3 Node Renumbering 185
Cuthill-McKee Algorithm
The Cuthill-McKee method [121 is a significantly more sophisticated renumbering
scheme than the wave front method. It explicitly considers the connectivities at each
node and is therefore good for both uniform and highly non-uniform meshes.
The algorithm is as follows:
Step 1: Scan all the nodes and order them according to their degree. The degree of
a node is the number of other nodes connected to it.
Step 2: Pick a starting node, usually one of low degree, and number it 1. Node 1 is
said to be in Level 1.
Step 3: Nodes connected to Node 1 are said to be in Level 2. These are then num-
bered 2, 3, and so on, in order of increasing degree.
Step 4: Un-numbered nodes connected to nodes in Level 2 are then numbered, again
in order of increasing degree.
Step 7: Repeat from Step 1 starting with a different node until a satisfactory mini-
mum bandwidth is achieved.
A variation on this algorithm, called the reverse Cuthill-McKee method, proceeds as
above but then reverses the numbers. That is, the last numbered node becomes Node
1, the second last becomes Node 2, and so on. This does not affect the bandwidth,
but it does affect the distribution of zeros inside the band.
Example 7.1: Use the Cuthill-McKee method to renumber the nodes in Fig-
ure 7.3 so as to decrease the bandwidth.
ill()
CD ill
Figure 7.3: Renumbering of the Nodes.
The initial numbering sequence follows the perimeter and then includes the in-
ternal nodes. This may seem an unnatural or artificial numbering scheme, but is
186 Chapter 7. Computer Methods I
quite likely to occur if the interior nodes were generated automatically using a mesh
generator. The purpose of this example is to show how, in a fairly automatic way,
the algorithm can give significant improvement in the bandwidth.
Recall that the renumbering does not take the degrees of freedom at each node
into account. Therefore, the connectivity matrix (introduced in Chapter 4) is very
similar to a stiffness matrix if there is only one degree of freedom at each node. For
the initial numbering, this matrix is
1 2 3 4 5 6 7 8 9 10 11 12
1 # + +
2 + # + + +
3 + # + + + +
4 + # +
5 + + # + +
6 + # + +
7 + # + +
8 + # + + +
9 + # + +
10 + + + # +
11 + + + + + # +
12 + + + + + + #
This gives a semi-band B = 10. It must be emphasized at this stage that the above
numbering is not necessarily bad; the 'goodness' of a numbering depends on the type
of algorithms used to operate on the arrays. The above array is considered sparse,
and there are routines designed explicitly to handle this type of array. The solver
we present next works best with tightly banded arrays.
The first step in the renumbering is to determine the degree of each node. That
is, we determine for each node the number of other nodes connected to it; this is its
degree.
degree 1: None
degree 2: 1,4
degree 3: 6, 7, 9
degree 4: 2,5,8, 10
degree 5: 3
degree 6: 11,12
Pick Node 1 as the starting node (since it is of lowest degree) and run through each
level renumbering the nodes in sequence.
1 2 3 4 5 6 7 8 9 10 11 12
1 #-+ +
2 + #- + + +
3 + + #- + +
4 + # + + +
5 + #- + + + + +
6 + + + + #- + + +
7 + # + +
8 + + + + #- + +
9 + + + #- + + +
10 + + + + + #- + +
11 + + + #- +
12 + + + #-
It is possible that another renumbering could reduce the bandwidth even more,
but the reduction will be minimal compared to what has already been achieved. Thus
the main benefit of this renumbering is in the context of automatic mesh generators.
That is, it makes it possible to piece the mesh (frame) together using any convenient
numbering scheme; then a pass through this algorithm will give consecutive numbers
with a small bandwidth. This is especially useful if the mesh is made of separate
pieces with separate numbering schemes. For example, each segment could have a
numbering beginning with a different multiple of 1000, say.
The primary mathematical task associated with matrix analysis of frame structures
consists of solving a set of N simultaneous linear algebraic equations for N unknowns.
For small scale problems, a wide variety of schemes can be used; but for the large
systems usually associated with structural analyses, great care must be exercised in
choosing the correct method.
The following introduces an effective scheme which is called the factorization
method (also referred to as the method of decomposition). This approach is particularly
well suited for matrix analysis of structures because it provides the efficiency of the
well-known Gaussian elimination process with a matrix format. (See Reference [34]
for more discussion of the Gaussian elimination scheme.) Since the stiffness matrices
of linearly elastic structures are always symmetric, a specialized type of factorization
known as the Cholesky method will be developed for symmetric matrices.
188 Chapter 7. Computer Methods I
[I{J = [U fr D J[ U]
This, in expanded form, is
Kl l K 12 K 13 KIN 1 0 0 0
K 21 K n K 23 K 2N U12 1 0 0
K 31 K 32 K 33 K 3N U13 U23 1 0 X
K Nl KN2 K N3 K NN U1N U2 N U3 N 1
Dn 0 0 0 1 U12 U13 U1 N
0 D 22 0 0 0 1 U23 U2 N
0 0 D23 0 0 0 1 U3 N
o 0 0 D NN 0 0 0 1
This technique is known as the modified Cholesky method; the 'usual' Cholesky
method uses the decomposition [K] = [U V[ Uland requires all diagonal terms
r
to be positive. We will see that the diagonal matrix D J has some very interesting
and useful properties.
The elements of the matrix [K ] essentially are obtained as inner products of the
rows and the columns of [ U ]. Thus, the inner products of the first column with itself
and subsequent columns produce
... ,
Similarly, the inner products of the second column with itself and subsequent columns
are
... ,
and so on. From this, the diagonal terms in [ K 1 can be written as
or, compactly,
i-I
i-I
I<ij = DiiUij +L DkkUkiUkj (1 ~ i ~ N) (i +1 ~ j ~ N)
k=1
Since we actually want the elements of r D J and [ U ], then these can be found
by rearranging the above as follows:
i-I
Dii = I<ii - L DkkUkiUki (1 ~ i ~ N)
k=1
These three equations constitute the recurrence formulas in an algorithm for factoring
the matrix [I<] into rUT DU] form. In these equations, the rows are traversed as
i = 1, Nj and for each row the column is traversed as j = i, N.
A significant attribute of this decomposition algorithm is that it can be done
in place. That is, the decomposed matrix is stored in the same memory locations
as the original matrix. This gives the minimum storage requirements. An additional
consequence is that the decomposed array inherits the same bandedness as the original
array.
The diagonal coefficients D ii have a very useful structural interpretation which
can be seen by re-arranging the system
Example 7.2: The following array is typical (except for size) of those usually
encountered in structural mechanics. It is desired to decompose it into its triangular
form.
1 -1 0 0]
[o
-1 2 -1 0 T
[K] = 0 -1 3 -2 = [U H D J[ U]
0 -2 3
D 33 K 33 - (DIlU13U13 + D22U23U23) = 3 - 1 = 2
U34 [K 34 - (DlIU12U14 + D22U23U24)]/D33 = [-2]/2 =-1
Finally, increment i to 4 and obtain
1 -1 0 0]
[ oo
[UTDU] = -1 1 -1 0
-1 2 -1
0 -1 1
It is significant that this array has the same bandedness as the original.
7.4 Solving Simultaneous Equations 191
Column-wise Sequence
The above recurrence equations imply that the diagonal term D jj is computed first,
followed by the calculation of the other terms in the i-th row of [ U ). This row-wise
generation of terms requires twice as many multiplications as actually needed. This
increase in the number of operations can be avoided by changing to the following
column-wise sequence:
Note that the product DkkUkj appears in both equations, so let this product be
and perform the calculations for Uij and D jj in the following manner:
i-I
U i} = J(jj - L UkiUkj (1 '5: j '5: N) (1 '5: i '5: j)
k=1
j-I
D jj = J(jj - L UkjUkj (1 '5: j '5: N) (7.2)
k=1
where U kj = Uk) D kk Thus, for column j the intermediate product Uij is obtained
for each off-diagonal term after the first. Then, the diagonal term D jj is computed,
during which calculation the final value of each off-diagonal term is also found. By
this sequence of operations, the number of multiplications is reduced to that for the
usual Cholesky method; in addition, the calculation of square roots is avoided. In
fact, it is slightly more efficient for that reason.
Example 7.3: Decompose the array of the last example, but this time use the
column-wise sequence.
Start with j = 1, obtain
Dn = Kn =1
Now increment j to 2, and take i = 1 to obtain
K12 =-1
(Ui2)/D n = -1
(/(22 - U12 Ui2) = 1
192 Chapter 7. Computer Methods I
[J(]{U}={p}
in which {u} is a column vector of N unknowns and {p} is a column vector of N
known constant terms. As a preliminary step, replace [ J( l with its decomposition to
obtain
[UfrDJ[U]{U}={p}
Now define the vector {y} to be
[U]{U}={Y}
In expanded form this expression is
1
o
o
o o 0 1
rDJ{Y}={z}
7.4 Solving Simultaneous Equations 193
) I I
Dl l 0
~
0 0 Yl ZI
0 D22 0 0 Y2 Z2
0 0 D33 0 Z3
1
Y3
0 0 0 DNN YN ZN
[uf{z}={p}
I I
1
)~I
0 0 0 ZI PI
U 12 1 0 0 Z2 P2
U 13 U 23 1 0 Z3 P3
U 1N U 2N U 3N 1 ZN PN
The original vector of unknowns {u} may now be obtained in three steps using
the above three intermediate products. First solve for the vector {z}. Since [ U jT
is a lower triangular matrix, each element can be calculated in a series of forward
substitutions. For example, the first few terms are
The second step consists of solving for the vector {y}. Since D J is a diagonal ma- r
trix, each element can be found simply by dividing terms in {z} by the corresponding
r
diagonal term of D J, as follows:
(1 ~ i ~ N)
and so on. In general, the elements of {u} may be calculated from the recurrence
formula:
N
UN = YN, Ui = Yi - L UikUk (N ~ i ~ 1)
k=i+l
This step completes the solution of the original equations for the unknown quantities.
The solution procedure involves both a forward and a backward substitution.
It is apparent that once a decomposition has been obtained, the solution to many
different load problems requires only the forward and back substitution. We will make
use of this feature when solving dynamic problems; for example, in Chapter 13 we
will pose transient problems as a sequence of pseudo-static problems with the load
vector changing but the stiffness remaining constant.
~1o -1
~1 ~13 -2~] { :~ } = { ~3 } U3
[
o 0 -2 3 U4 4
This is the same array as in the previous examples, so using its rUT DUj form of
decomposition, gives
10
-1 1
o -1
~ ~] [~ ~ ~ ~] [~ ~1 ~1
100020001
[
o 0 -1 1 0 0 0 1 0 0 0
This is broken into a series of substitution problems corresponding to the three
square arrays. The first of these requires solving
Zl = 1, Z3 = 6,
These form the right hand side of the second problem
Yl = 1, Y4 = 10
7.4 Solving Simultaneous Equations 195
1
-1 -1
2 0
-1 0]
0 { 17}
16 ={ 17 - 16
-17+32-13 } = { I2 }
[ o -1 3 -2 13 -16+39-20 3
o 0 -2 3 10 -26 + 30 4
This process of verifying the solution is often used to assess if significant round-
off errors have accumulated. Note, however, that the original matrix (and not its
decomposed form) must be used, thus requiring additional storage.
K{i,j) = KUb{i,j - B + 1)
where the superscript ub means upper banded.
It happens that multiplication and division are far more computationally costly
than addition or subtraction, hence the computational cost is usually based only on
the former. The cost for the modified Cholesky method is
cost -- !2 N B 2 - !3 B 3 + !2 B 2 - '31 B
For large systems this shows that the computational cost of the decomposition is
about !N B2. It is seen that if advantage is not taken of the bandedness, then the
method is essentially on the order of iN3; since B is typically less than 10% of N
this translates in a computational difference of at least 100. The cost for the solution
stage is about
cost = 2NB + N
from which it is apparent that the decomposition is the significant part of the total
solution cost.
To put this computational cost into perspective, consider solving a medium sized
problem (I 000 degrees of freedom) on our benchmark machine of 1 MFlops. Assuming
a 10% bandwidth, we have
On the other hand, a large problem with 10,000 degrees of freedom would take about
one hour.
More details on the properties of eigenvalues and eigenvectors can be found in Chap-
ter 11; here we use only those properties necessary to help us estimate the first few
lowest eigenvalues of the system.
We will concentrate on the following form of the eigenvalue problem
[/{]{u} -'x[M]{U} =0
2
In vibration problems ,X = w (where w is the vibration frequency) and therefore ,X is
always positive. On the other hand, in stability problems ,X = -P and therefore it
can be positive or negative (as we have already seen in Example 5.5). To avoid loss
of generality in the following, we will therefore allow ,X to be positive or negative.
This scheme requires pre-processing of the data which adds overhead to the com-
putational cost. This cost is on the order of N B2 for the decomposition and inversions;
it is quite significant when only a few eigenvalues are to be computed. The scheme
also has difficulties if there are zero masses since then the decomposition will fail.
The method to be presented next works directly on the general eigenvalue problem
and is efficient when only a few of the eigenvalues and eigenvectors are required.
Chapter 13 will present some other schemes also.
198 Chapter 7. Computer Methods I
Let us assume a trial vector for {}, say {u} 1, and introduce on the right-hand side
an equivalent load vector of
{Rh = [M]{uh
The original eigenvalue problem will now resemble an equilibrium equation as en-
countered in a static analysis, which we may write as
Since the roots are ordered so that Al < A2 < A3'" the ratios AI/An are less than
or equal to unity. Hence, with each iteration each guess will become more and more
dominated by the vector {4> h. That is, there will be convergence to the fundamental
mode. It must be realized that for each iteration, effectively new coefficients {c} are
calculated.
Convergence depends on the ratio Ad A2' the smaller this is, the faster the process
converges. It also depends on the relative strengths of the modes in the initial guess,
thus if {u h has a large component of the first eigenvector {4> h, i.e., CI ~ C2, C3, . ,
relatively few iteration steps are required for convergence. On the other hand, if the
first guess is totally deficient of the first mode, i.e., CI = 0, then theoretically the
process will not find {4> h. In practice, however, rounding errors arising from finite
precision in the computer will eventually produce a component of the first mode in one
of the subsequent guesses, and this will eventually become dominant. It is expected,
therefore, that convergence will always be to the fundamental mode, although the
rate depends on the initial guess. For this reason, it is usual to take the initial guess
to be fully populated with ones. Repeated roots (A2 = AI, say) is obviously a problem
but we postpone its solution until Chapter 13.
Computer Algorithm
The iteration is begun by assuming a starting iteration vector { u h. The basic step in
the iteration scheme is the solution of the pseudo-static problem in which we evaluate
a vector {it} k+I from
This normalization calculation merely assures that the length of the new iteration
vector is unity with respect to the mass. That is,
If this scaling is not included in the iteration, the elements of the iteration vectors
grow (or decrease) in each step and do not converge to {4> h but to a multiple of it.
n
Provided that {u h is not M -orthogonal to {4> h, meaning that {u} M]{ 4> h i= 0,
we then have the following limits
as k-+oo
[K]=[ufrDj[ul
Step 6: Use the Rayleigh quotient to get an estimate of the eigenvalue from
{Qh+l
{}
II k+l = II Z Il I / 2
Step 8: Check for convergence. Convergence can be said to have been achieved
when successive values of P satisfy
A= Pk+I ,
Example 7.5: Determine the lowest eigenvalue and the corresponding eigenvec-
tor for the following system.
Step 1 of the algorithm requires the decomposition of the stiffness matrix. This
is the same array as in the previous example, so using its rUTnUl form of decom-
position gives
~, ]
0 0 0 0 -1 0
~JU Wi
[ -1
1 1 0 1 0 1 -1
[K)= ~ -1 1 0 2 0 1
0 -1 0 0 0 0
202 Chapter 7. Computer Methods I
Step 2 chooses an initial guess for the effective load vector. Take
{R}f = {I, 1, 1, I}
Step 3 solves the static problem
[ K ]{ u} = {R}
using the initial guess for the load vector. The solution is
{uV{R}
P = IIZII = 0.0751104
In anticipation of the need for another iteration, we compute in Step 7 another load
vector
{R}
{Q} = {.4613, .8140, .8954, .8683}
=~
IIZII
We now check for convergence according to
This criterion is based on successive estimates of P, and we only have one, we must
repeat the process.
Just the major numbers will now be stated.
Iteration 2 :
{u} {5.8614, 5.4001 , 4.1247, 3.0392}
{Q} = {5.8614 , 10.8002, 12.3741, 12.1570}
IIZII 13.441248
P = 0.0743521
{R} = {.4360, .8035, .9206, .9044}
CONV 1.019899 X 10- 2
Iteration 3 :
{u} {5.8204, 5.3843,4.1447, 3.0646}
{Q} {5.8204 , 10.7687 , 12.4342, 12.2586}
7.5 Solving Eigenvalue Problems 203
IIZII = 13.452346
P = 0.0743352
{R} = {.4326, .8005, .9243, .9112}
CONV = 2.268693 x 10- 4
Iteration 4 :
{u} = {5.8133 , 5.3806 , 4.1475 , 3.0687}
{Q} = {5.8133, 10.7613, 12.4425, 12.2750}
IIZII = 13.452642
P = 0.0743348
{R} = {.4321 , .7999, .9249, .9124}
CONV 6.365913 x 10- 6