Plan Stress PDF

6
Plane Stress
Transformations
6–1
Lecture 6: PLANE STRESS TRANSFORMATIONS
TABLE OF CONTENTS
Page
§6.1. Introduction 6–3
§6.2. Thin Plate in Plate Stress 6–3
§6.3. 2D Stress Transformations 6–5
§6.3.1. Why Are Stress Transformations Important? . . . . . . . 6–5
§6.3.2. Method of Equations . . . . . . . . . . . . . . 6–5
§6.3.3. Double Angle Version . . . . . . . . . . . . . . 6–6
§6.3.4. Principal Stresses, Planes, Directions, Angles . . . . . . 6–6
§6.3.5. Maximum Shear Stresses . . . . . . . . . . . . . 6–7
§6.3.6. Principal Stress Element . . . . . . . . . . . . . 6–8
§6.3.7. Mohr’s Circle . . . . . . . . . . . . . . . . . 6–8
§6.4. What Happens in 3D? 6–10
§6.4.1. Including the Plane Stress Thickness Dimension . . . . . 6–10
§6.4.2. 3D Mohr Circles . . . . . . . . . . . . . . . . 6–10
§6.4.3. Overall Maximum Shear . . . . . . . . . . . . . 6–11
§6.4.4. Plane Stress Revisited . . . . . . . . . . . . . . 6–12
§6.4.5. The Sphere Paradox . . . . . . . . . . . . . . . 6–13
6–2
§6.2 THIN PLATE IN PLATE STRESS
§6.1. Introduction
This Lecture deals with the plate stress problem. This is a two-dimensional stress state, briefly
introduced in §1.5 of Lecture 1. It occurs frequently in two kinds of aerospace structural components:
1. Thin wall plates and shells; e.g., aircraft and rocket skins, and the pressure vessels of Lecture 3.
2. Shaft members that transmit torque. These will be studied in Lectures 7–9.
The material below focuses on thin flat plates, and works out the associated problem of plane stress
transformations.
§6.2. Thin Plate in Plate Stress
In structural mechanics, a flat thin sheet of material is called a plate. The distance between the plate
faces is the thickness, which is denoted by h. The midplane lies halfway between the two faces.
The direction normal to the midplane is the transverse z
direction. Directions parallel to the midplane are called Top surface
in-plane directions. The global axis z is oriented along
the transverse direction. Axes x and y are placed in the
midplane, forming a right-handed Rectangular Carte-
sian Coordinate (RCC) system. Thus the equation of y
the midplane is z = 0. The +z axis conventionally
defines the top surface of the plate as the one that it x
intersects, whereas the opposite surface is called the Figure 6.1. A plate structure in plane stress.
bottom surface. See Figure 6.1.
A plate loaded in its midplane is said to be in a state of plane stress, or a membrane state, if the
following assumptions hold:
1. All loads applied to the plate act in the midplane direction, as pictured in Figure 6.1, and are
symmetric with respect to the midplane.
2. All support conditions are symmetric about the midplane.
3. In-plane displacements, strains and stresses are taken to be uniform through the thickness.
4. The normal and shear stress components in the z direction are zero or negligible.
The last two assumptions are not necessarily consequences of the first two. For those to hold, the
thickness h should be small, typically 10% or less, than the shortest in-plane dimension. If the plate
thickness varies it should do so gradually. Finally, the plate fabrication must exhibit symmetry with
respect to the midplane.
To these four assumptions we add an adscititious restriction:
5. The plate is fabricated of the same material through the thickness. Such plates are called
transversely homogeneous or (in aerospace) monocoque plates.
The last assumption excludes wall constructions of importance in aerospace, in particular composite
and honeycomb sandwich plates. The development of mathematical models for such configurations
requires a more complicated integration over the thickness as well as the ability to handle coupled
bending and stretching effects. Those topics fall outside the scope of the course.
6–3
y
Midplane
Mathematical
idealization
Plate
Γ
x
Figure 6.2. Mathematical model of plate in plane stress. (Symbols and , used to
denote the plate interior and the boundary, respectively, are used in advanced courses.)
In-plane internal forces y + sign conventions

for internal forces,
dx dy stresses and strains
h pyy
Thin plate in plane stress y x
pxx pxy dy
z x dx
dx dy In-plane stresses In-plane body forces

y dx dy
dx dy
h
σyy y h
x y
x σxx τxy = τ yx by
bx
x
In-plane strains In-plane displacements

dx dy dx dy
h h
εyy y uy y
x ε xx γ xy = γ yx x u x
Figure 6.3. Notational conventions for in-plane stresses, strains, displacements and
internal forces of a thin plate in plane stress.
Remark 6.1. Selective relaxation from assumption 4 leads to the so-called generalized plane stress state, in
which nonzero σzz stresses are accepted; but these stresses do not vary with z. The plane strain state described
in §5.4.2 of Lecture 5 is obtained if strains in the z direction are precluded: zz = γx z = γ yz = 0.
Remark 6.2. Transverse loading on a plate produces plate bending, which is associated with a more complex
configuration of internal forces and deformations. This topic is studied in graduate-level courses.
Remark 6.3. The requirement of flatness can be relaxed to allow for a curved configuration, as long as
the structure, or structure component, resists primarily in-plane loads. In that case the midplane becomes a
midsurface. Examples are rocket and aircraft skins, ship and submarine hulls, open parachutes, boat sails
and balloon walls. Such configurations are said to be in a membrane state. Another example are thin-wall
members under torsion, which are covered in Lectures 8–9.
The plate in plane stress idealized as a two-dimensional problem is illustrated in Figure 6.2.
6–4
§6.3 2D STRESS TRANSFORMATIONS
σyy σtt
(a) (b) τtn
τyx τnt σnn
τxy
P σxx P
y t y
n
Global axes Local axes
x,y stay fixed x θ x n,t rotate by θ
z z
Figure 6.4. Plane stress system referred to global axes x, y (data) and to local rotated axes n, t.
(Locations of point P in (a,b) coincide; they are drawn offset for visualization convenience.)
In this idealization the third dimension is represented as functions of x and y that are integrated
through the plate thickness. Engineers often work with internal plate forces, which result from
integrating the in-plane stresses through the thickness. See Figure 6.3.
In this Lecture we focus on the in-plane stresses σx x , σ yy and τx y and their expressions with respect
to an arbitrary system of axes
§6.3. 2D Stress Transformations

The stress transformation problem studied in this Lecture is illustrated in Figure 6.4. Stress com-
ponents σx x , σ yy and τx y at a midplane point P are given with respect to the global axes x and y,
as shown in Figure 6.4(a). The material element about P is rotated by an angle θ that aligns it with
axes n, t, as shown in Figure 6.4(b). (Note that location of point P in (a,b) coincide — they are
drawn offset for visualization convenience.)
The transformation problem consists of expressing σnn , σtt , and τnt = τtn in terms of the stress
data σx x , σ yy and τx y , and of the angle θ. Two methods, one analytical and one graphical, will
be described here. Before that is done, it is useful to motivate what are the main uses of these
transformations.
§6.3.1. Why Are Stress Transformations Important?

The transformation problem has two major uses in structural analysis and design.
• Find stresses along a given skew direction. Here θ is given as data. This has several
applications. Two examples:
1. Analysis of fiber reinforced composites if the direction of the fibers is not aligned with the
{x, y} axes. A tensile normal stress perpendicular to the fibers may cause delamination,
and a compressive one may trigger buckling.
2. Oblique joints that may fail by shear parallel to the joint. For example, welded joints.
• Find max/min normal stresses, max in-plane shear and overall max shear. This may be
important for strength and safety assessment. Here finding the angle θ is part of the problem.
Both cases are covered in the following subsections.
6–5
§6.3.2. Method of Equations

The derivation given on pages 524-525 of Vable or in §7.2 of Beer-Johnston-DeWolf is based on
the wedge method. This will be omitted here for brevity. The final result is
σnn = σx x cos2 θ + σ yy sin2 θ + 2 τx y sin θ cos θ,

σtt = σx x sin2 θ + σ yy cos2 θ − 2 τx y sin θ cos θ, (6.1)
τnt = −(σx x − σ yy ) sin θ cos θ + τx y (cos2 θ − sin2 θ).
Couple of checks are useful to verify these equations. If θ = 0◦ ,
σnn = σx x , σtt = σ yy , τnt = τx y , (6.2)
as expected. If θ = 90◦ ,
σnn = σ yy , σtt = σx x , τnt = −τx y . (6.3)
which is also OK (can you guess why τnt at θ = 90◦ is − τx y ?). Note that
σx x + σ yy = σnn + σtt . (6.4)
This sum is independent of θ , and is called a stress invariant. (Mathematically, it is the trace of the
stress tensor.) Consequently, if σnn is computed, the fastest way to get σtt is as σx x + σ yy − σnn ,
which does not require trig functions.
§6.3.3. Double Angle Version

For many developments it is convenient to express the transformation equations (6.1) in terms of
the “double angle” 2θ by using the well known trigonometric relations cos 2θ = cos2 θ − sin2 θ
and sin 2θ = 2 sin θ cos θ, in addition to sin2 θ + cos2 θ = 1. The result is
σx x + σ yy σx x − σ yy
σnn = + cos 2θ + τx y sin 2θ,
2 2 (6.5)
σx x − σ yy
τnt = − sin 2θ + τx y cos 2θ.
2
Here σtt is omitted since, as previously noted, it can be quickly computed as σtt = σx x + σ yy − σnn .
§6.3.4. Principal Stresses, Planes, Directions, Angles

The maximum and minimum values attained by the normal stress σnn , considered as a function of
θ, are called principal stresses. (A more precise name is in-plane principal normal stresses, but
the qualifiers “in-plane” and “normal” are often dropped for brevity.) Those values occur if the
derivative dσnn /dθ vanishes. Differentiating the first of (6.1) and passing to double angles gives
dσnn
= 2(σ yy − σx x ) sin θ cos θ + 2τx y (cos2 θ − sin2 θ)
dθ (6.6)
= (σ yy − σx x ) sin 2θ + 2τx y cos 2θ = 0.
6–6
This is satisfied for θ = θ p if

2 τx y
tan 2θ p = (6.7)
σx x − σ yy
There are two double-angle solutions 2θ1 and 2θ2 given by (6.7) in the range [0 ≤ θ ≤ 360◦ ] or
[−180◦ ≤ θ ≤ 180◦ ] that are 180◦ apart. (Which θ range is used depends on the textbook; here
we use the first one.) Upon dividing those values by 2, the principal angles θ1 and θ2 define the
principal planes, which are 90◦ apart. The normals to the principal planes define the principal
stress directions. Since they differ by a 90◦ angle of rotation about z, it follows that the principal
stress directions are orthogonal.
As previously noted, the normal stresses that act on the principal planes are called the in-plane
principal normal stresses, or simply principal stresses. They are denoted by σ1 and σ2 , respectively.
Using (6.7) and trigonometric relations it can be shown that their values are given by

2
σ1,2 = ± + τx2y . (6.8)
2 2
To evaluate (6.8) it is convenient to go through the following staged sequence:

1. Compute
2
σav = and R=+ + τx2y . (6.9)
2 2
Meaning of these values: σav is the average normal stress at P (recall that σx x + σ yy is an
invariant and so is σav ), whereas R is the radius of the Mohr’s circle, as described in §6.3.7
below, thus the symbol. Furthermore, R represents the maximum in-plane shear stress value,
as discussed in §6.3.5 below.
2. The principal stress values are
σ1 = σav + R, σ2 = σav − R. (6.10)
3. Note that the a priori computation of the principal angles is not needed to get the principal
stresses if one follows the foregoing steps. If finding those angles is of interest, use (6.7).
Comparing (6.6) with the second of (6.5) shows that dσnn /dθ = 2τnt . Since dσnn /dθ vanishes for
a principal angle, so does τnt . Hence the principal planes are shear stress free.
§6.3.5. Maximum Shear Stresses

Planes on which the maximum shear stresses act can be found by setting dτnt /dθ = 0. A study
of this equation shows that the maximum shear planes are located at ±45◦ from the principal
planes, and that the maximum and minimum values of τnt are ±R. See for example, §7.3 of the
Beer-Johnston-DeWolf textbook.
This result can be obtained graphically on inspection of Mohr’s circle, covered later.
6–7
principal
directions
σ2 =10 psi |τmax |= R=50 psi
σyy = 20 psi θ2= 108.44
(a) (b)
τ xy = τyx =30 psi σ1 =110 psi
(c) 18.44 +45
= 63.44
P σxx =100 psi P θ1=18.44 P
x x
principal principal
planes planes
y t y
principal
n planes
x θ (d) 45 principal
stress
x element
plane of max 45
inplane shear P
Figure 6.5. Plane stress example: (a) given data: stress components σx x , σ yy and τx y ; (b) principal stresses and
angles; (c) maximum shear planes; (d) a principal stress element (actually four PSE can be drawn, that shown is one
of them). Note: locations of point P in (a) through (d) coincide; they are drawn offset for visualization convenience.
Example 6.1. This example is pictured in Figure 6.5. Given: σx x = 100 psi, σ yy = 20 psi and τx y = 30 psi, as
shown in Figure 6.5(a), find the principal stresses and their directions. Following the recommended sequence
(6.9)–(6.10), we compute first
2
100 + 20 100 − 20
σav = = 60 psi, R=+ + 302 = 50 psi, (6.11)
2 2
from which the principal stresses are obtained as
σ1 = 60 + 50 = 110 psi , σ2 = 60 − 50 = 10 psi. (6.12)
To find the angles formed by the principal directions, use (6.7):
2 × 30 3
tan 2θ p = = = 0.75, with solutions 2θ1 = 36.87◦ , 2θ2 = 2θ1 + 180◦ = 216.87◦ ,
100 − 20 4 (6.13)
whence θ1 = 18.44 , θ2 = θ1 + 90◦ = 108.44◦ .
◦
These values are shown in Figure 6.5(b). As regards maximum shear stresses, we have |τmax | = R = 50 psi.
The planes on which these act are located at ±45◦ from the principal planes, as illustrated in Figure 6.5(c).
§6.3.6. Principal Stress Element

Some authors, such as Vable, introduce here the so-called principal stress element or PSE. This is
a wedge formed by the two principal planes and the plane of maximum in-plane shear stress. Its
projection on the {x, y} plane is an isosceles triangle with one right angle and two 45◦ angles. For
the foregoing example, Figure 6.5(d) shows a PSE.
There are actually 4 ways to draw a PSE, since one can join point P to the opposite corners of
the square aligned with the principal planes in two ways along diagonals, and each diagonal splits
the square into two triangles. The 4 images may be sequentially produced by applying sucessive
rotations of 90◦ . Figure 6.5(d) shows one of the 4 possible PSEs for the example displayed in that
Figure. The PSE is primarily used for the visualization of material failure surfaces in fracture and
yield, a topic only covered superficially in this course.
6–8
2θ2 = 36.88 +180 = 216.88

τ = shear
(a) Point in plane stress stress τmax = 50
50
σyy = 20 psi 40
(a) H
τ xy = τyx =30 psi 30
Radius R = 50
20
σ = normal
P σxx =100 psi 10
0 20 40 60 80 100 stress
0
σ2 = 10 −10
C
y
−20 σ1 = 110
−30
−40 V
2θ1 = 36.88
−50
x τmin = −50
(b) Mohr's circle
coordinates of blue points are
H: (20,30), V:(100,-30), C:(60,0)
Figure 6.6. Mohr’s circle for plane stress example of Figure 6.5.
§6.3.7. Mohr’s Circle

Mohr’s circle is a graphical representation of the plane stress state at a point. Instead of using the
methods of equations, a circle is drawn on the {σ, τ } plane. The normal stress σ (θ) and the share
stress τ (θ ) are plotted along the horizontal and vertical axes, respectively, with θ as a parameter.
All stress states obtained as the angle θ is varied fall on a circle called Mohr’s circle.1
This representation was more important for engineers before computers and calculators appeared.
But it still retains some appealing features, notably the clear visualization of principal stresses and
maximum shear. It also remains important in theories of damage, fracture and plasticity that have
a “failure surface”.
To explain the method we will construct the circle corresponding to Example 6.1, which is re-
produced in Figure 6.6(a) for convenience. Draw horizontal axis σ = σnn (θ) to record normal
stresses, and vertical axis τ = τnt (θ) to record shear stresses. Mark two points: V (for ”vertical
cut”, meaning a plane with exterior normal parallel to x) at (σx x , −τx y ) = {100, −30} and H (for
“horizontal cut”, meaning a plane with exterior normal parallel to y) at (σ yy , τx y ) = (20, 30).
The midpoint between H and V is C, the circle center, which is located at ( 12 (σx x + σ yy ), 0) =
(σav , 0) = (60, 0). Now draw the circle. It may be verified that its radius is R as given by (6.9);
which for Example 6.1 is R = 50. See Figure 6.6(b).
The circle intersects the σ axis at two points with normal stresses σavg + R = 110 = σ1 and
σav − R = 10 = σ2 . Those are the principal stresses. Why? At those two points the shear stress τnt
vanishes, which as we have seen characterizes the principal planes. The maximum in-plane shear
occurs when τ = τnt is maximum or minimum, which happens at the highest and lowest points of
the circle. Obviously τmax = R = 50 and τmin = −R = −50. What is the normal stress when the
shear is maximum or minimum? By inspection it is σav = 60 because those points lie on a vertical
line that passes through the circle center C.
1 Introduced by Christian Otto Mohr (a civil engineer and professor at Dresden) in 1882. Other important personal
contributions were the concept of statically indeterminate structures and theories of material failure by shear.
6–9
Other features, such as the correlation between the angle θ on the physical plane and the rotation
angle 2θ traversed around the circle will be explained in class if there is time. If not, one can find
those details in Chapter 7 of Beer-Johnston-DeWolf textbook, which covers Mohr’s circle well.
§6.4. What Happens in 3D?
Despite the common use of simplified 1D and 2D structural models, the world is actually three-
dimensional. Stresses and strains actually “live” in 3D. When the extra(s) space dimensions are
accounted for, some paradoxes are resolved. In this final section we take a quick look at principal
stresses in 3D, stating the major properties as recipes.
§6.4.1. Including the Plane Stress Thickness Dimension

To fix the ideas, 3D stress results will be linked to the plane stress case studied in Example 6.1. and
pictured in Figure 6.5. Its 3D state of stress in {x, y, z} coordinates is defined by the 3 × 3 stress
matrix
100 30 0
S = 30 20 0 (6.14)
0 0 0
It may be verified that the eigenvalues of this matrix, arranged in descending order, are
σ1 = 110, σ2 = 10, σ3 = 0. (6.15)
Now in 3D there are three principal stresses, which act on three mutually orthogonal principal
planes that are shear stress free. For the stress matrix (6.14) the principal stress values are 110, 10
and 0. But those are precisely the eigenvalues (6.15). This result is general:
The principal stresses in 3D are the eigenvalues of the 3D stress matrix
The stress matrix is symmetric. A linear algebra theorem says that a real symmetric matrix has a full
set of eigenvalues and eigenvectors, and that the eigenvalues are guaranteed to be real. The normals
to the principal planes (the so-called principal directions) are defined by the three orthonormalized
eigenvectors, but this topic will not be pursued further.
In plane stress, one of the eigenvalues is always zero because the last row and column of S are null;
thus one principal stress is zero. The associated principal plane is normal to the transverse axis z,
as can be physically expected. Consequently σzz = 0 is a principal stress. The other two principal
stresses are the in-plane principal stresses, which were those studied in §6.3.4. It may be verified
that their values are given by (6.8), and that their principal directions lie in the {x, y} plane.
Some terminology is needed. Suppose that the three principal stresses are ordered, as usually done,
by decreasing algebraic value:
σ1 ≥ σ2 ≥ σ3 (6.16)
Then σ1 is called the maximum principal stress, σ3 the minimum principal stress, and σ2 the inter-
mediate principal stress. (Note that this ordering is by algebraic, rather than by absolute, value.)
In the plane stress example (6.14) the maximum, intermediate and minimum normal stresses are
110, 10 and 0, respectively. The first two are the in-plane principal stresses.
6–10
§6.4 WHAT HAPPENS IN 3D?
Overall All possible

τ = shear max shear stress states lie
in the shaded τ = shear
(a)
stress area between stress τmax
overall
= 55
(b)
the outer and τmax
inplane
= 50
Inner Mohr's inner circles 50
circles 40
30
20
σ = normal σ3 = 0 10 σ = normal
stress 0 20 40 60 80 100 stress
0
σ3 σ2 −10
σ1 σ2 = 10
−20 σ1 = 110
−30
Outer Mohr's
circle −40
−50
Yellow-filled circle is the inplane one
Figure 6.7. Mohr’s circles for a 3D stress state: (a) general case; (b) plane stress example of Figure 6.5. In (b)
the Mohr circle of Figure 6.6 is the rightmost inner circle. Actual stress states lie on the grey shaded areas.
§6.4.2. 3D Mohr Circles

More surprises (pun intended): in 3D there are actually three Mohr circles, not just one. To draw
the circles, start by getting the eigenvalues σ1 , σ2 and σ3 of the stress matrix, for example through
the eig function of Matlab. Suppose they are ordered as per the convention (6.16). Mark their
values on the σ axis of the σ vs. τ plane. Draw the 3 circles with diameters {σ1 , σ2 }, {σ2 , σ3 } and
{σ1 , σ3 }, as sketched in Figure 6.7(a).2 For the principal values (6.15) the three Mohr circles are
drawn in Figure 6.7(b), using a scale similar to that of Figure 6.6.
Of the three circles the wider one goes from σ1 (maximum normal stress) to σ3 (minimum normal
stress). This is known as the outer or “big” circle, while the two others are called the inner circles.
It can be shown (this is proven in advanced courses in continuum mechanics) that all actual stress
states at the material point lie between the outer circle and the two inner ones. Those are marked
as grey shaded areas in Figure 6.7.
§6.4.3. Overall Maximum Shear

For ductile materials such as metal alloys, which yield under shear, the 3D Mohr-circles diagram
is quite useful for visualizing the overall maximum shear stress at a point, and hence establish the
factor of safety against that failure condition. Looking at the diagram of Figure 6.7(a), clearly the
overall maximum shear, called τmax overall
, is given by the highest and lowest point of the outer Mohr’s
circle, marked as a blue dot in that figure. (Note that only its absolute value is of importance for
safety checks; the sign has no importance.) But that is also equal to the radius Router of that circle.
If the three principal stresses are algebraically ordered as in (6.16), then
σ3 − σ1
τmax
overall
= Router = (6.17)
2
Note that the intermediate principal stress σ2 does not appear in this formula.
2 In the general 3D case there is no simple geometric construction of the circles starting from the six x, y, z independent
stress components, as done in §6.3.7 for plane stress.
6–11
If the principal stresses are not ordered, it is necessary to use the max function in a more complicated
formula that picks up the largest of the 3 radii:

σ1 − σ2 σ2 − σ3 σ3 − σ1
τmax = max
overall , , (6.18)
2 2 2
Taking absolute values in (6.18) is important because the max function picks up the largest algebraic
value. For example, writing τmax
overall
= max(30, −50, 20) = 30 picks up the wrong value. On the
other hand τmax max(|30|, | − 50|, |20|) = 50 is correct.
overall
§6.4.4. Plane Stress Revisited

Going back to plane stress, how do overall maximum shear and in-plane maximum shear compare?
Recall that one of the principal stresses is zero. The ordering (6.16) of the principal stresses is
assumed, and one of those three is zero. The following cases may be considered:
(A) The zero stress is the intermediate one: σ2 = 0. If so, the in-plane Mohr circle is the outer
one and the two shear maxima coincide:
σ1 − σ3
τmax
overall
= τmax
inplane
= (6.19)
2
(B) The zero stress is either the largest one or the smallest one. Two subcases:
σ1
(B1) If σ1 ≥ σ2 ≥ 0 and σ3 = 0 : τmax
overall
= ,
2 (6.20)
σ3
(B2) If σ3 ≤ σ2 ≤ 0 and σ1 = 0 : τmax
overall
=− .
2
In the plane stress example of Figure 6.5, the principal stresses are given by (6.12). Since both
in-plane principal stresses (110 psi and 10 psi) are positive, the zero principal stress is the smallest
inplane
one. We are in case (B), subcase (B1). Consequently τmax overall
= 12 σ1 = 55 psi, whereas τmax =
1
(σ − σ2 ) = 50 psi.
2 1
6–12
§6.4 WHAT HAPPENS IN 3D?
τ = shear τ = shear
stress stress
(a) (b) τmax
overall
= 40
50 50
40 40
30 τinplane
max = 0
30 τinplane
max = 0
20 σ3 = 0 20
10 σ = normal 10 σ = normal
0 20 40 60 80 100 stress 0 20 40 60 80 100 stress
0 0
−10 −10
−20 σ1 = σ 2=80 −20 σ1 = σ 2=80
−30 −30
−40 −40
−50 −50
inplane
Figure 6.8. The sphere paradox: (a) Mohr’s in-plane circle reduced to a point, whence τmax = 0;
(b) drawing the 3D circles shows that τmax
overall = 40.
§6.4.5. The Sphere Paradox

Taking the third dimension into account clarifies some puzzles such as the “sphere paradox.”
Consider a thin spherical pressure vessel with p = 160 psi, R = 10 in and t = 0.1 in. In Lecture 3,
the wall stress in spherical coordinates was found to be σ = p R/(2 t) = 80 ksi, same in all
directions. The stress matrix at any point in the wall, taking z as the normal to the sphere at that
point, is
80 0 0
0 80 0 (6.21)
0 0 0
The Mohr circle reduces to a point, as illustrated in Figure 6.8(a), and the maximum in-plane shear
stress is zero. And remains zero for any p. If the sphere is fabricated of a ductile material, it should
never break regardless of pressure. Plainly we have a paradox.
The paradox is resolved by considering the other two Mohr circles. These additional circles are
pictured in Figure 6.8(b). In this case the outer circle and the other inner circle coalesce. The
overall maximum shear stress is 40 ksi, which is nonzero. Therefore, increasing the pressure will
eventually produce yield, and the paradox is solved.
6–13
17
Free
Single-DOF
Oscillator
17–1
Lecture 17: FREE SINGLE-DOF OSCILLATOR
TABLE OF CONTENTS
Page
§17.2. Notation 17–3
§17.3. Free Vibrations of Undamped SDOF Oscillator 17–3
§17.3.1. Equation of Motion Of Undamped Oscillator . . . . . . . 17–4
§17.3.2. Undamped Response in Terms of Complex Exponentials . . 17–4
§17.3.3. Undamped Response in Terms of Trigonometric Functions . . 17–5
§17.3.4. Effect of Initial Conditions . . . . . . . . . . . . 17–6
§17.4. Free Vibrations of Viscous-Damped SDOF Oscillator 17–6
§17.4.1. Equation of Motion of Damped Oscillator . . . . . . . . 17–6
§17.4.2. Damped Response in Terms of Trigonometric Functions . . . 17–7
§17.4.3. Underdamped Case: ξ < 1 . . . . . . . . . . . . . 17–8
§17.4.4. Critically Damped Case: ξ = 1 . . . . . . . . . . . 17–9
§17.4.5. Overdamped Case: ξ > 1 . . . . . . . . . . . . . 17–10
17–2
§17.3 FREE VIBRATIONS OF UNDAMPED SDOF OSCILLATOR
§17.1. Introduction
This lecture begins Part VI of the course. This Part provides a quick introduction to structural
dynamics and in particular free and forced vibrations. Nowadays the term dynamics has acquired
several interpretations, especially in the social sciences. It is used here in the traditional sense:
“The study of the relationship between motion and the forces affecting motion”
This is in fact meaning (1a) given in the American Heritage Dictionary of the English Language.
The corresponding adjectives are dynamic1 and its equivalent dynamical.2
But even the traditional definition is far too broad in two respects. First, Mechanics embodies
a wide range of scales that span from cosmological through atomic and sub-atomic. Second, the
study can focus on three aspects: theoretical, applied and computational. Our focus is restricted to a
particular flavor: Classical Mechanics, which is the subset of Mechanics that obeys Newton’s laws.
This restriction allows the use of continuum (field) models as well as certain “lumped” idealizations
(for example, point masses) that can be derived directly from such laws.
An important application area of Classical Mechanics is to structures, a subject naturally called
Structural Dynamics. This is the application considered here. More advanced treatments, which
fall under the purview of Solid and Continuum Mechanics, are studied in dedicated courses at the
senior and graduate level.
§17.2. Notation
As just observed, Structural Dynamics studies the motion of structures under time-dependent forces.
The time will be always denoted by t. Derivatives with respect to t will be abbreviated by superposed
dots. For example, if u(t) is a scalar motion, the associated velocity and acceleration are compactly
written
du(t) d 2 u(t)
u̇ ≡ , ü ≡ . (17.1)
dt dt 2
Sometimes we will denote velocity by v and acceleration by a where appropriate. Note that partial
differentiation with respect to time is not needed: the ordinary differential symbol d suffices.
§17.3. Free Vibrations of Undamped SDOF Oscillator

We start up this study with the simplest dynamic system that obey the laws of Classical Dynamics:
the unforced spring-mass oscillator undergoing free vibrations.
A weight of mass m > 0 hangs from an extensional spring of stiffness k ≥ 0 under gravity
acceleration g directed along the spring direction. The weight will be treated as a point mass. The
weight force is W = m g. We will assume small displacements throughout. The static deflection
of the mass is
δs = W/k = m g/k. (17.2)
This configuration is called a Single Degree of Freedom oscillator, or SDOF oscillator for short.
See Figure 17.1(a,b).
1 Etymology: French dynamique, from Greek dunamikos: powerful, from dunamis, power, from dunashai, to be able.
2 The variant dynamical tends to be often used in a more abstract sense. For example, MathWorld defines dynamical
system as “a means of describing how a state evolves into another over the course of time.”
17–3
Apply now an initial displacement u 0 and initial velocity v0 to the mass at t = 0, so that
u(0) = u 0 , u̇(0) = v0 . (17.3)
These are called the initial conditions or ICs. Then release the mass and do not apply any force
for t > 0. The oscillator is set in motion, and a free vibrations response ensues. Let u(t) be the
deviation from the static equilibrium position so that the total motion is δs + u(t), as shown in
Figure 17.1(c). The function u(t) is called the dynamic response. The only DOF is u = u(t). We
next study the equation of motion (EOM) and its solution.
§17.3.1. Equation of Motion Of Undamped Oscillator

Remove the spring and replace it by the force Fs = k(δs + u(t)), as drawn in Figure 17.1(d). The
other two forces acting on the mass are its weight W = m g, which is independent of time, and the
inertia force m ü(t), which acts in the opposite direction to the acceleration ü = ü(t) (because it
resists it). The resulting picture is called a Dynamic Free Body Diagram, or DFBD.
Equilibrium of forces in the x direction requires m ü = W − k(δs + u) = m g − kδs − ku. On
canceling m g − kδs = 0 we get
m ü + k u = 0. (17.4)
This is the equation of motion (EOM) of the free, undamped, single DOF, linear oscillator. It is
a second order, linear ODE in u(t). Notice that the static deflection δs has disappeared from this
EOM; in fact the mass will oscillate about that position.
To convert to canonical form, divide (17.4) through by m and introduce ωn2 = k/m:

k
ü + ωn2 u = 0, in which ωn = + . (17.5)
m
This parameter ωn is called the undamped circular natural frequency, or natural frequency for short.
Its units are radians per second (rad/s).
Key features of the equation of motion (17.5) are:
1. EOM is linear: superposition applies
2. EOM has constant coefficients: solutions are elementary circular functions
3. EOM is second-order in time: two constants of integration ⇒ two initial conditions required
§17.3.2. Undamped Response in Terms of Complex Exponentials

Let u(t) = C eλt , where C
= 0 and λ are generally complex. Since ü = λ2 C eλt , substitution into
(17.5) requires
(λ2 + ωn2 ) C eλt = 0 (17.6)
But since C
= 0 (else the solution is null) and the exponential never vanishes, the term in parenthesis
must be zero. This gives the characteristic equation
λ2 + ωn2 = 0, ⇒ λ1,2 = ±i ωn , (17.7)
17–4
§17.3 FREE VIBRATIONS OF UNDAMPED SDOF OSCILLATOR
;;; ;; ;;
(a) (b) (c) (d)
k k g g
k
Fs = k (δs+u)
δs ..
u=u(t) mu
δs = W/k
W ..
x W W mu
Figure 17.1. Undamped, unforced spring-mass SDOF oscillator undergoing free vibrations.
√
in which i = −1 denotes the imaginary unit. Because both eiωn t and e−iωn t satisfy the linear ODE
(17.4), so does any linear combination of them. We thus arrive at the general solution in terms of
complex exponentials:
u(t) = C1 eiωn t + C2 e−iωn t , (17.8)
in which C1 and C2 are generally complex numbers.
§17.3.3. Undamped Response in Terms of Trigonometric Functions

We know that the motion u(t) is supposed to be real. Therefore it is convenient to recast (17.8) as
a real expression. Using Euler’s relation for the complex exponential
e±iθ = cos θ ± i sin θ, (17.9)
the solution (17.8) becomes
u(t) = (C1 + C2 ) cos ωn t + i(C1 − C2 ) sin ωn t. (17.10)
To simplify the linkage of this expression to the initial conditions, introduce A1 = C1 + C2 and
A2 = i(C1 − C2 ). Then we compactly express the response in terms of trigonometric functions as
u(t) = A1 cos ωn t + A2 sin ωn t. (17.11)
Here A1 and A2 are real constants that can be directly determined from the initial conditions (17.3):
u(0) = u 0 = A1 , u̇(0) = v0 = A2 ωn , (17.12)
Replacing into (17.11) yields

v0
u(t) = u 0 cos ωn t + sin ωn t. (17.13)
ωn
17–5
This is the free vibration response of an undamped SOF system expressed in terms of trigonometric
functions and initial conditions. Although (17.13) holds only for free vibrations, this solution will
be later reused as the homogeneous portion of the response of a forced SDOF oscillator. Another
common version of (17.13) is the phased response (also called phase-shifted response) form
u(t) = U cos(ωn t − α). (17.14)
Here the amplitude U and the phase angle α are linked to the initial conditions (17.3) by

2
v0 v0
U= u 20 + , tan α = . (17.15)
ωn ωn u 0
§17.3.4. Effect of Initial Conditions
To study the effect of initial conditions, consider first the case where the mass is displaced from its
static equilibrium position by u 0 and released. If so v0 = 0, and (17.13) reduces to
u(t) = u 0 cos ωn t. (17.16)
This is plotted in Figure 17.2(a). The resulting response is a simple harmonic motion with amplitude
u 0 , undamped natural frequency f n and undamped natural period Tn , given by
ωn 1 2π
fn = , Tn = = , (17.17)
2π fn ωn
Frequency f n is expressed in cycles per second or Hertz, abbreviated to Hz (1 Hz = cycle/s). The

period is given in seconds per cycle, or simply seconds (s). It s easily verified that in this case the
phased response (17.14) coalesces with (17.16) since U = u 0 and α = 0.
Figure 17.2(b) shows a response plot when u 0 = 0 and v0 is nonzero. In this case the phased
response (17.14) gives a similar harmonic curve but shifted in time. If both u 0 and v0 are nonzero,
the response has amplitude U related to the initial conditions by the first of (17.15). The maxima
and minima occur at t = α/ωn + 12 nTn , n = 1, 2, . . . .
§17.4. Free Vibrations of Viscous-Damped SDOF Oscillator
The response plots for the previous case convey the idea that, once an undamped SDOF oscillator
is set in motion, that motion will continue for all times t > 0 without any change in amplitude.
In reality all mechanical systems exhibit some damping, which dissipates energy and causes the
motion eventually to die out. In the present section we consider the case of a SDOF oscillator under
the simplest damping model: viscous damping, which is linearly proportional to the velocity. It is
modeled by a dashpot element of viscosity c ≥ 0, as illustrated in Figure 17.3(a).
17–6
§17.4 FREE VIBRATIONS OF VISCOUS-DAMPED SDOF OSCILLATOR
(a) (b)
2 2
1.5 u(t) u0 = 1, v 0 = 0 1.5 u(t) u0 = 0, v0 = 1
1 1
0.5 0.5
0 0
−0.5 −0.5
−1 −1
−1.5 −1.5
−2 −2
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
Figure 17.2. Response of undamped spring-mass oscillator with k = m = 1: (a) unit initial
displacement and zero initial velocity; (b) unit initial displacement and unit initial velocity.
§17.4.1. Equation of Motion of Damped Oscillator

From the Dynamic Free Body Diagram (DFBD) shown in Figure ?(d) we obtain the physical form
of the EOM:
m ü + c u̇ + k u = 0. (17.18)
Divide this equation through by m and denote

k k c c
ωn2 = , ωn = + , c = 2ξ ωn m, ξ= = √ . (17.19)
m m 2ωn m 2 km
Here ωn is the undamped√ circular natural frequency introduced in the previous section, while
ξ = 2 c/(ωn m) = 2 c/ k m is the viscous damping factor, also called damping ratio and damping
1 1
coefficient. We thus convert (17.18) to the canonical form of the EOM:
ü + 2ξ ωn u̇ + ωn2 u = 0. (17.20)
This is a linear, second-order ODE with constant coefficients. The only difference with respect to
the undamped case is the presence of the velocity term 2ξ ωn u̇.
§17.4.2. Damped Response in Terms of Trigonometric Functions

As usual in solving ODEs with the canonical form (17.20), assume a solution
u(t) = A eλt . (17.21)
where A and λ are generally complex values to be determined. Inserting into (17.20), we obtain
the characteristic equation
λ2 + 2ξ ωn λ + ωn2 = 0. (17.22)
17–7
;; ;; ;;
(a) (b) (c) (d)
k c k c g c g
k .
Fs = k (δs +u) Fd = cu
δs ..
u=u(t) mu
δs = W/k
W
x W ..
W mu
Figure 17.3. Damped spring-dashpot-mass oscillator undergoing free vibrations: (a) initial unloaded position;
(b) static equilibrium position; (c) dynamic equilibrium position; (d) Dynamic Free Body Diagram (DFBD).
which is quadratic in λ. Its roots are given by

λ1,2 = −ξ ωn ± ωn ξ 2 − 1. (17.23)
Note that for the undamped case: ξ = 0, the roots reduce to ±iωn , as found previously. The
magnitude of the damping factor ξ compared to unity can be used to distinguish three cases:
ξ <1 Underdamped case. Damping is called subcritical. The roots λ1,2 in (17.23) are complex
conjugate. The motion is oscillatory with decreasing amplitude. This is the most common
case in typical structures, and thus the most practically important one.
ξ >1 Overdamped case. Damping is called overcritical. The roots λ1,2 in (17.23) are negative
real and distict. The motion is non-oscillatory. Its amplitude decays monotonicallyexcept
possibly for one zero crossing.
ξ =1 Critically damped case. The roots λ1,2 are negative real and coalesce. The motion
is non-oscillatory. Its amplitude decays monotonically except possibly for one zero
crossing. Since the response exhibits the most rapid decay toward the static position,
critical damping is of interest in the design of devices such as instruments and suspensions.
We now proceed to study those three cases in more detail. Emphsis is placed on the underdamped
case, which is the most important one for structural dynamics.
§17.4.3. Underdamped Case: ξ < 1
For convenience the roots (17.23) of the characterist equation may be expressed as follows
λ1,2 = −ξ ωn ± iωd , (17.24)
in which ωd denotes the damped circular natural frequency given by

ωd = ω n 1 − ξ 2 . (17.25)
17–8
(a) (b)
1.5 1.5
u(t) u0 = 1, v 0 = 0 u(t) u0 = 0, v0 = 1
1 1
ξ=0
0.5 0.5
ξ=1
ξ=1 ξ=0.5
0 0
ξ=0.5 ξ=0.2
−0.5 −0.5 ξ=0.2
ξ=0
−1 −1
−1.5 −1.5
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
Figure 17.4. Response of underdamped (and critically damped) spring-dashpot-mass oscillator

with k = m = 1, for various values of the damping factor ξ that range from 0 to 1: (a) unit initial
displacement and zero initial velocity; (b) unit initial displacement and unit initial velocity.
Like ωn , this is expressed in radians per second. The corresponding damped period is
2π
Td = . (17.26)
ωd
With the help of this symbols and Euler’s formula the general solution can be expressed as
u(t) = e−ξ ωn t (A1 cos ωd t + A2 sin ωd t) . (17.27)
Again using the initial conditions u(0) = u 0 and u̇(0) = v0 , we obtain A1 = u 0 and A2 =
(v0 + ξ ωn u 0 )/ωd , which substituted into (17.27) gives the response

−ξ ωn t v0 + ξ ωn u 0
u(t) = e u 0 cos ωd t + sin ωd t . (17.28)
ωd
This equation may in turn be rewritten in the phased form
u(t) = U e−ξ ωn t cos(ωd t − α) (17.29)
in which

2
v 0 + ξ ωn u 0 v0 + ξ ωn u 0
U= u 20 + , tan α = . (17.30)
ωd ωd u 0
This kind of response is illustrated in Figure 17.4.
17–9
(a) (b)
0.5
1.4 u(t) u(t)
u0 = 1, v 0 = 0 u0 = 0, v0 = 1
1.2 0.4 ξ=1
1
0.3
0.8 ξ=10 ξ=2
0.6 ξ=5 0.2
0.4 ξ=2 ξ=5
0.1
0.2 ξ=1
ξ=10
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
(c)
1.5
u(t) u 0 = 1, v0 = −5
1
ξ=10
0.5
ξ=5
0
ξ=2
−0.5
ξ=1
−1
0 2 4 6 8 10 12
Time t
Figure 17.5. Response of overdamped and critically spring-dashpot-mass oscillator with k = m = 1, for various
values of the damping factor ξ that range from 1 to 10: (a) unit initial displacement and zero initial velocity; (b)
unit initial displacement and unit initial velocity; (c) unit initial displacement and high negative initial velocity.
§17.4.4. Critically Damped Case: ξ = 1

In this case the characteristic roots coalesce: λ1 = λ2 = −ωn . The second order ODE theory says
that the solution is u(t) = (A1 + A2 t)e−ωn t . When the initial conditions are taken into account we
get
u(t) = [u 0 + (v0 + ωn u 0 ) t] e−ωn t . (17.31)
Plots of critically damped system responses are illustrated in Figures 17.4 and 17.5 by taking ξ = 1
(red curve response).
§17.4.5. Overdamped Case: ξ > 1

If ξ > 1 the characteristic equation has two distict negative real roots. For convenience define

ω∗ = ωn ξ 2 − 1. (17.32)
Then the solution may be written in terms of the hyperbolic sine and cosine as
u(t) = e−ξ ωn t (A1 cosh ω∗ t + A2 sinh ω∗ t), (17.33)
17–10
in which A1 and A2 depend on initial conditions. Introducing these we arrive at the response

−ξ ωn t ∗ v0 + ξ ωn u 0 ∗
u(t) = e u 0 cosh ω t + sinh ω t . (17.34)
ω∗
For a pictorial effect of damping factor on the response of an overdamped system (including the
transition critically damped case ξ = 1) see Figure 17.5.
17–11
18
Harmonically Forced
SDOF Oscillator
18–1
Lecture 18: HARMONICALLY FORCED SDOF OSCILLATOR
TABLE OF CONTENTS
Page
§18.2. Harmonically Excited Viscous-Damped SDOF Oscillator 18–3
§18.2.1. Equation of Motion and Steady-State Solution . . . . . . 18–4
§18.2.2. Magnification Factor and Phase Lag . . . . . . . . . 18–4
§18.2.3. Frequency Response Features . . . . . . . . . . . . 18–5
§18.2.4. Total Response . . . . . . . . . . . . . . . . 18–6
§18.3. Response to Harmonic Base Excitation 18–6
18–2
§18.2 HARMONICALLY EXCITED VISCOUS-DAMPED SDOF OSCILLATOR
This Lecture generalizes the solutions presented in Lecture 17 for the SDOF oscillator by considering
an applied harmonic force for t ≥ 0. The frequency of this harmonic force is called the exciting
frequency. This forced SDOF configuration is called the forced harmonic oscillator.
When an external force is applied, the EOM, while still a second-order ODE, becomes nonhomogenous.
The total response is the superposition of the homogeneous solution, which is defined by the initial
conditions, and the particular solution, which is defined by the prescribed harmonic force. Structural
engineers call these two components the starting transient or simply transient, and the steady-state
components of the dynamic response, respectively.
Those names suggest that the particular (= steady-state) solution will eventually dominate. And indeed
if even a slight amount of damping is present, the homogenous component becomes negligible after
sufficient time has elapsed, and we can focus our attention on the particular or steady-state solution.
This focus will allow us to exhibit the important phenomenon of resonance. This term quantifies the
response amplification that occurs when the exciting frequency approaches the natural frequency of
the oscillator.
We will consider a damped oscillator from the start since the presence of damping is crucial in estab-
lishing the peak magnitude of resonance amplification. The undamped case can be easily recovered
by setting the damping factor to zero.
;;
(a) (b)
k c
Static equilibrium
position .
Fs =k u Fd = c u
Dynamic u=u(t)
motion Mass m Mass m
p(t) =p0 cos Ω t ..

p0 cos Ω t FI = m u
Figure 18.1. Damped spring-dashpot-mass oscillator under prescribed harmonic force

excitation: (a) dynamic equilibrium position; (b) dynamic FBD.
§18.2. Harmonically Excited Viscous-Damped SDOF Oscillator
We consider again the mass-dashpot-spring system studied in §17.3, but now subjected to a harmonic
excitation force
p(t) = p0 cos t, (18.1)
18–3
(a) (b)
100 200
50 100
Phase lag α in degrees

Magnification factor Ds
ξ=0.01 ξ=5
50
10 ξ=1
ξ=0.1 ξ=0.5
5 20
ξ=0.2 ξ=0.2
10
1 ξ=0.5 ξ=0.1
5
0.5 ξ=1 ξ=0.01
ξ=5
2
0.1 0.2 0.5 1 2 5 10 0 0.5 1 1.5 2

Frequency ratio r = Ω/ωn Frequency ratio r = Ω/ωn
Figure 18.2. Axially loaded ductile bar: failure by yield due to crystal slip at ±45◦ .
in which is the excitation frequency. See Figure 18.1(a). Note that all static actions, such as gravity,
are omitted from the drawing to improve clarity, being tacitly understood that u = u(t) is the dynamic
motion measured from the static equilibrium position.
§18.2.1. Equation of Motion and Steady-State Solution

The DFD shown in Figure 18.1(b) gives us the dynamic equations
m ü + c u̇ + k u = p0 cos t. (18.2)
As remarked in the Introduction, the complete response will be the sum of the transient (= homogenous)
and steady-state (= particular)components. If damping is present, after a while the transient part, which
is that dependent on initial conditions, becomes unimportant. This solution, denoted by u p (t), with
subscript p for “particular,” is conveniently expressed in the phased form
u p (t) = U cos(t − α). (18.3)
The associated velocity and acceleration are
u̇ p (t) = − U sin(t − α), ü p (t) = −2 U cos(t − α). (18.4)
§18.2.2. Magnification Factor and Phase Lag

Inserting this into the EOM (18.2), and representing its components as a force vector polygon in the
phase space (see Craig-Kurdila, pages 87–88) we obtain
c
(k U − m2 U )2 + (cU )2 = p02 , tan α = , (18.5)
k − m2
For convenience introduce the static displacement U0 and the frequency ratio as
p0
U0 = , r= . (18.6)
k ωn
18–4
§18.2 HARMONICALLY EXCITED VISCOUS-DAMPED SDOF OSCILLATOR
(a) (b)
10 175

8 ξ=0.01 150
125
6
100
4 ξ=0.1 75
ξ=5
ξ=0.2 50 ξ=0.1
ξ=1
2 ξ=0.2
ξ=0.5 25 ξ=0.5
ξ=1 ξ=0.01
ξ=5
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 18.3. Frequency response log plots for harmonically forced viscous-damped oscillator: (a) plot of
magnification factor versus frequency ratio; (b) plot of phase lag angle versus frequency ratio.
(a) (b)
100 200
50 100
ξ=0.01 ξ=5
50
10 ξ=1
ξ=0.1 ξ=0.5
5 20
ξ=0.2 ξ=0.2
10
1 ξ=0.5 ξ=0.1
5
0.5 ξ=1 ξ=0.01
ξ=5
2
0.1 0.2 0.5 1 2 5 10 0 0.5 1 1.5 2

Figure 18.4. Frequency response log plots for harmonically forced viscous-damped oscillator: (a) log-log
plot of magnification factor versus frequency ratio; (b) log plot of phase lag angle versus frequency ratio.
Observe that U0 is the displacement that the mass would undergo if a force of magnitude p0 were to
be applied statically.
With the help of (18.6) the the expressions in (18.5) can be compactly expressed as
def U (r ) 1 2ξ r
Ds (r ) = = , tan α(r ) = . (18.7)
U0 (1 − r 2 )2 + (2ξr )2 1 − r2
Here Ds is the steady-state magnification factor, also called gain. These are plotted in Figures 18.2(a,b).
The peaks observed as r is near unity and ξ << 1 identify resonance.
Since Ds may be quite large near resonance whereas r may cover a broad range of frequencies, the
plots are often shown in log-log format, as shown in Figure 18.3(a,b). The log-log frequency plot is
known in control applications as the Bode plot.
18–5
Prescribed
base motion
ub =Ub cos ( Ω t − αb )
;;;
;;;
(a) (b)
Static equilibrium k c
position for zero
base motion . .
Fs =k (u − ub) Fd =c (u − ub)
u=u(t)
Mass m
FI = m u..
Figure 18.5. Damped spring-dashpot-mass oscillator under prescribed harmonic base
motion: (a) dynamic equilibrium position; (d) dynamic FBD.
§18.2.3. Frequency Response Features

The combination of amplitude versus frequency and phhase versus frequency information is called the
frequency response of the system. Graphs such as shown in Figures 18.2 and 18.3 are called frequency
response plots. The following significant features can be observed:
1. The steady state motion (18.3) is also harmonic and has the same frequency as the excitation.
2. The amplitude of the steady-state response is a function of both the amplitude and frequency
of the excitation, as well as of the natural frequency and damping factor of the system. The
magnification factor Ds may be considerably greater than unity (near resonance) or may be less
than unity.
3. The excitation p0 cos t and the steady-state response u p = U cos(t − α) are not in phase;
meaning that they do not attain their maximum values at the same times. The steady-state
§18.2.4. Total Response
Transcribing results of Lecture 17 for the free damped oscillator, the total response may be expressed
as
u(t) = U0 Ds cos (t − α) + u h (t), in whichu h (t) = e−ξ ωn t (A1 cos ωn t + A2 sinωn t).
(18.8)
where Ds and α are given by (18.7) whereas A1 and A2 in the homogeneous solution u h (t) are
determined by the initial conditions. Since u h (t) dies out with time if ξ > 0, that justifies the name
starting transient in commonly uses by structural engineers.
§18.3. Response to Harmonic Base Excitation

Suppose now that the SDOF oscillator is not subject to an applied force, but instead the base displaces
by a prescribed harmonic motion specified in the form
u b (t) = Ub cos(t − α), (18.9)
18–6
§18.3 RESPONSE TO HARMONIC BASE EXCITATION
See Figure 18.4(a). This kind of excitation is commonly called ground motion or ground excitation
in Civil structures subject to seismic loads. It also appears in vibration isolation design problems in
which equipment, instruments or payloads must be protected from base vibratory motions caused by
launches, maneuvers, etc. It is also the configuration pertinent to vibration resonance tests such as the
“Airplane Shaker” Lab 3 demos on November 9.
Assume that the only excitation is base motion. From the FBD of Figure 18.4(b) we get FI +Fs +Fd = 0,
or
m ü + k(u − u b ) + c(u̇ − u̇ b ) = 0, (18.10)
Passing the prescribed motion terms to the RHS we arrive at the EOM:
m ü + cu̇ + ku = ku b + cu̇ b . (18.11)
If the base excitation is harmonic as in (18.9),

m ü + cu̇ + ku = Ub k cos(t − αb ) − c) sin(t − αb ) . (18.12)
Comparing to (18.2) shows that this is similar to the equation of motion of a harmonically forced
oscillator, except that both cosine and sine terms are present in the forcing function. If the damping
vanishes: c = 0, (18.12) and (18.2) are identical if U0 is replaced by k Ub .
18–7
19
.
MDOF
Dynamic
Systems
19–1
Lecture 19: MDOF DYNAMIC SYSTEMS 19–2
TABLE OF CONTENTS
Page
§19.2. A Two-DOF Lumped-Parameter System 19–3
§19.2.1. Equations of Motion . . . . . . . . . . . . . . . 19–3
§19.2.2. Matrix Form of EOM . . . . . . . . . . . . . . 19–4
§19.2.3. Undamped Natural Frequencies and Modes . . . . . . . 19–5
§19.3. Generalization to Arbitrary Number of DOF 19–6
§19.3.1. Linkage to the Algebraic Eigenproblem . . . . . . . . 19–6
§19.3.2. Generalized Mass and Stiffness . . . . . . . . . . . 19–7
§19.3.3. Eigenvector Normalization Criteria . . . . . . . . . . 19–9
19–2
19–3 §19.2 A TWO-DOF LUMPED-PARAMETER SYSTEM
We move now to dynamic structural systems with Multiple Degrees of Freedom (abbreviation:
MDOF). For simplicity we focus on linear, viscous-damped, lumped parameter models. All explicit
examples have only two DOF, but that is sufficient to show how such models are treated using matrix
algebra. Once the matrix equations are exhibited in compact form, the extension to an arbitrary,
but finite, number of DOF is readily done.
;; (a) (b)
k1 c1
Static equilibrium
position .
Fs1 = k1 u1 Fd1 = c1 u 1
u 1 = u1(t)
Mass m1
p1 (t) p1 (t) ..
FI1 = m 1 u1
Fs2 Fd 2
Static equilibrium k2 c2
position
Fs2 = k2 (u2−u1) . .
Fd 2 = c2 (u 2 −u1)
u 2 = u2(t)
Mass m2
p2 (t) ..
p2 (t) FI2 = m 2 u 2
x
Figure 19.1. Two-DOF, lumped parameter, spring-mass-dashpot dynamic system.
§19.2. A Two-DOF Lumped-Parameter System

Consider the lumped-parameter, mass-spring-dashpot dynamic system shown in Figure 19.1(a). It
has two point masses m 1 and m 2 , which are connected by a spring-dashpot pair with constants k2
and c2 , respectively. Mass m 1 is linked to ground by another spring-dashpot pair with constants k1
and c1 , respectively. The system has two degrees of freedom (DOF). These are the displacements
u 1 (t) and u 2 (t) of the two masses from their static equilibrium positions. These displacements are
assumed small. Known dynamic forces p1 (t) and p2 (t) act on the masses. Our goal is to formulate
the equations of motion (EOM), and to express them in matrix form.
19–3
§19.2.1. Equations of Motion

Because we have two DOF, we need two dynamic equilibrium equations. These can be formed
using two Dynamic Free Body Diagrams (DFBD). To get equations with a minimal number of
terms, is convenient to isolate the two masses as pictured in Figure 19.1(b). Spring and dashpot
forces are denoted by Fsi and Fdi for the i th spring and i th dashpot, respectively, with i = 1, 2. The
inertial force acting on the i th mass is FI i .
The positive force conventions are those shown in Figure 19.1(b). Summing forces along x we get
the equilibrium conditions

DFBD #1: Fx at mass 1 : − FI 1 − Fs1 − Fd1 + Fs2 + Fd2 + p1 = 0,
(19.1)
DFBD #2: Fx at mass 2 : − FI 2 − Fs2 − Fd2 + p2 = 0.
Next express the spring, dashpot and inertial forces in terms of displacement, velocities and ac-
celerations. Note that for the spring and dashpot forces we must use relative displacements and
velocities, respectively, whereas for inertia forces total accelerations must be used:
Fs1 = k1 u 1 , Fs2 = k2 (u 2 − u 1 ), Fd1 = c1 u̇ 1 , Fd2 = c2 (u̇ 2 −u̇ 1 ), FI 1 = m 1 ü 1 , FI 2 = m 2 ü 2 .

(19.2)
Replacing into (19.1) we get
−m 1 ü 1 − k1 u 1 − c1 u̇ 1 + k2 (u 2 − u 1 ) + c2 (u̇ 2 − u̇ 1 ) + p1 = 0,
(19.3)
−m 2 ü 2 − k2 (u 2 − u 1 ) − c2 (u̇ 2 − u̇ 1 ) + p2 = 0.
Finally, collect all terms depending on the DOF u 1 and u 2 and their temporal derivatives in the
LHS while moving everything else to the RHS. Sign convention used in mechanics: mass-times-
acceleration terms in LHS should be positive. Following that convention we get
m 1 ü 1 + c1 u̇ 1 + c2 u̇ 1 − c2 u̇ 1 + k1 u 1 + k2 u 1 − k2 u 2 = p1 ,
(19.4)
m 2 ü 2 − c2 u̇ 1 + c2 u̇ 2 − k2 u 1 + k2 u 2 = p2 .
§19.2.2. Matrix Form of EOM

Rewrite (19.4) in matrix notation as

m1 0 ü 1 c1 + c2 −c2 u̇ 1 k + k2 −k2 u1 p1
+ + 1 = . (19.5)
0 m2 ü 2 −c2 c2 u̇ 2 −k2 k2 u2 p2
Passing to compact notation,
M ü + C u̇ + K u = p. (19.6)
Here M, C and K denote the mass, damping and stiffness matrices, respectively, while p, u, u̇ and
ü denote the force, displacement, velocity and acceleration vectors, respectively. The latter four
are function of time: u = u(t), etc., but the time argument will be often omitted for brevity.
Note that matrices M, C and K are symmetric, and that M is diagonal. The latter property holds
when masses are “point lumped” at the DOF locations, as in this model.
19–4
19–5 §19.2 A TWO-DOF LUMPED-PARAMETER SYSTEM
Example 19.1. Suppose m 1 = 2, m 2 = 1, c1 = 0.1, c2 = 0.3, k1 = 6, k2 = 3, p1 = 2 sin 3tk, p2 = 5 cos 2t.

Then

2 0 ü 1 0.4 −0.3 u̇ 1 9 −3 u1 2 sin 3t
+ + = . (19.7)
0 1 ü 2 −0.3 0.3 u̇ 2 −3 3 u2 5 cos 2t
The data matrices and vectors are

2 0 0.4 −0.3 9 −3 2 sin 3t
M= , C= , K= , p= . (19.8)
0 1 −0.3 0.3 −3 3 5 cos 2t
§19.2.3. Undamped Natural Frequencies and Modes

Suppose that the foregoing system is undamped: c1 = c2 = 0, whence C = 0, the null 2 × 2 matrix.
It is also unforced: p1 = p2 = 0, whence p = 0, the 2 × 1 null vector. The matrix system (19.6)
reduces to
M ü + K u = 0. (19.9)
This is the MDOF generalization of the free, undamped SDOF oscillator covered in Lecture 17.
Next assume that this undamped, unforced system is undergoing in-phase harmonic motions of
circular frequency ω. This assumption can be mathematically stated using either trigonometric
functions, or complex exponentials:
u(t) = U cos (ωt − α), or u(t) = U eiωt , (19.10)
in which U is a nonzero 2-vector that has the amplitudes of the motions of the point masses as
entries, and (in the first form) α is a phase shift angle. Both expressions lead to identical results.
For the ensuing derivations we select the first form: u(t) = U cos(ωt − α). The corresponding
velocity and accelerations are u̇(t) = −ω U sin(ωt − α) and ü = −ω2 U cos(ωt − α). Substitute
these in (19.9), and extract U cos(ωt − α) as common post-multiplier factor:

M ü + K u = K − ω2 M U cos(ωt − α) = 0. (19.11)
If this vector expression is to be identically zero for all time t, the product of the first two terms
(matrix times vector) must vanish. Thus

K − ω2 M U = D(ω) U = 0. (19.12)
def
Here D(ω) = K − ω2 M is called the dynamic matrix.
Equation (19.12) states the free vibrations eigenproblem for an undamped MDOF system. For
nontrivial solutions U = 0 the determinant of the dynamic matrix must vanish, whence

C(ω2 ) = det D(ω) = det K − ω2 M = 0. (19.13)
This is the characteristic equation for free undamped vibrations. For a two-DOF system, C(ω2 ) is
a quadratic polynomial in ω2 , which will yield two roots: ω12 and ω22 .
19–5
The next section shows that under certain conditions on K and M, which are satisfied here, those
roots are real and nonnegative. Consequently the positive square roots

ω1 = + ω12 , ω2 = + ω22 , (19.14)
are also real and nonnegative. These are called the undamped natural circular frequencies. (Qual-
ifiers “undamped” and “circular” are often omitted unless a clear distinction needs to be made.) As
usual they are measured in radians per second, or rad/sec. They can be converted to frequencies in
cycles per second, or Hertz (Hz), by scaling through 1/(2π ), and to natural periods in seconds by
taking their reciprocals:
ωi 1 2π
fi = , Ti = = , i = 1, 2. (19.15)
2π fi ωi
How about the U vector in (19.12)? Insert ω2 = ω12 there and solve for U = 0. (This operation
is possible because K − ω12 M is singular by definition.) Call that solution vector U1 . Repeat with
ω2 = ω22 and call that vector U2 .
Vectors U1 and U2 are mathematically known as eigenvectors. Since U1 and U2 may be scaled by
arbitrary nonzero factors, we must normalize them to make them unique. Normalization criteria
often used in dynamic analysis are discussed in §19.3.3.
For this particular problem the eigenvectors are called undamped free-vibrations natural modes, a
mouthful often abreviated to just modes, on account of the physical interpretation discussed later.
The visualization of an eigenvector as a motion pattern is called a mode shape.
§19.3. Generalization to Arbitrary Number of DOF
Thanks to matrix notation, the extension to an undamped, unforced MDOF system with n degrees
of freedom is immediate.
In this general case, matrices K and M are both n × n. The characteristic polynomial is of order
n in ω2 . This polynomial provides n roots labeled ωi2 , i = 1, 2 . . . n, arranged in ascending order.
Their square roots are the natural frequencies of the system.
The associated vibration modes are the eigenvectors of the vibration eigenproblem. They are called
Ui before normalization. Normalization schemes are studied in §19.3.3.
§19.3.1. Linkage to the Algebraic Eigenproblem

To facilitate linkage to the subject of of algebraic eigenproblems covered in the Linear Algebra
sophomore courses, as well as to capabilities of computer tools such as Matlab, we effect some
notational changes. First, rewrite (19.12) as
K Ui = ωi2 M Ui , (19.16)
where i is a natural frequency index that ranges from 1 through n. Second, rename the symbols
that appear in (19.16) as follows: K → A, M → B, ωi2 − > λi , Ui → xi . This allows us to rewrite
(19.16) as
A xi = λi B xi , i = 1, 2, . . . n. (19.17)
19–6
19–7 §19.3 GENERALIZATION TO ARBITRARY NUMBER OF DOF
This is the canonical form of the generalized algebraic eigenproblem that appears in Linear Algebra
textbooks. If both A and B are real symmetric, (19.17) is called the generalized symmetric algebraic
eigenproblem. If, in addition, B is positive definite (PD), it is shown in those textbooks that
(i) All eigenvalues λi are real
(ii) All eigenvectors xi have real entries
Property (i) can be further strengthened if A enjoys additional properties:
(iii) If A is nonnegative definite (NND) and B is PD, all eigenvalues are nonnegative real
(iv) If both A and B are PD, all eigenvalues are positive real
The nonnegativity property is obviously important in vibration analysis since natural frequencies
are obtained by taking square roots of the eigenvalues, which appear as squared frequencies in
(19.16). If a squared frequency is negative, its square roots are imaginary. The possibility of
negative squared frequencies arise when studying dynamic stability and active control systems, but
those topics are beyond the scope of this course.
If matrix B in (19.17) is the identity I, (19.17) reduces to the standard algebraic eigenproblem
A xi = λi xi , i = 1, 2, . . . n. (19.18)
If A is real symmetric, (19.18) is called the standard symmetric algebraic eigenproblem, or symmet-
ric eigenproblem for short. Since the identity matrix is obviously symmetric and positive definite,
the aforementioned properties (i) through (iv) hold for (19.18).
Remark 19.1. Some higher level programming languages such as Matlab, Mathematica, and Maple can solve
eigenproblems directly through their built-in libraries. Of these 3 commercial products, Matlab has the most
extensive capabilities, being able to process the generalized eigenproblem form (19.17) directly. On the other
hand Mathematica only accepts the standard eigenproblem (19.18) as built-in function. Similar constraints
may apply to hand calculators with matrix processing abilities. To overcome those limitations, there are
transformation methods that convert (19.17) into (19.18). Only a couple of such schemes will be mentioned:
(I) If B is nonsingular, premultiply both sides of (19.17) by its inverse: B−1 to get B−1 A xi = B−1 λi B xi =
λi xi . This has the standard form Ã xi = λi xi , in which Ã = B−1 A. The eigenvalues and eigenvectors
of Ã are those of (19.17).
(II) If A is nonsingular, premultiply both sides of (19.17) by its inverse: A−1 to get A−1 A xi = A−1 λi B xi =
λi xi . This has the standard form B̃ xi = (1/λi ) xi , in which B̃ = A−1 B. The eigenvectors of B̃ are those
of (19.17) but its eigenvalues are reciprocated.
Note that matrices Ã or B̃ produced by these schemes will not generally be symmetric, even if A and B are.
However the eigenvalue-eigenvector properties listed previously will be preserved. There are fancier reduction
methods that preserve symmetry (for example, those based on the Cholesky factorization of either A or B),
but those noted above will be sufficient for this course.
Remark 19.2. For small number of DOF, say 4 or less, it is often expedient to solve for the eigenvalues directly
by finding the roots of the characteristic polynomial introduced in (19.13):
C(λ) = det(D) = 0. (19.19)
in which λ = ω2 and D = K − ω2 M = K − λM is the dynamic stiffness matrix. This procedure will be
followed in subsequent Lectures for 2-DOF examples. It is particularly useful with hand calculators that can
compute polynomial roots but cannot handle eigenproblems directly. Eigenvectors my be then obtained as
illustrated in the next Lecture.
19–7
§19.3.2. Generalized Mass and Stiffness

We go back to the stiffness-mass notation of (19.17), which is reproduced for convenience:
K Ui = ωi2 M Ui , (19.20)
For specificity we assume that:

(I) The mass matrix M is positive definite (PD) whereas the stiffness matrix K is nonegative
definite (NND). As a consequence all eigenvalues λi = ωi2 are nonnegative real, and so are
their positive square roots: ωi ≥ 0.
(II) All frequencies are distinct: ωi = ω j if i = j.
Under assumptions (I) and (II), the theory behind the algebraic eigenproblem says that for each ωi
there is one and only one eigenvector Ui , which is unique except for an arbitrary nonzero scale
factor.
An eigenvector Ui scaled as per some normalization condition (for example, unit Euclidean length)
is called a normal mode. We will write
Ui = ci φi , i = 1, 2, . . . n. (19.21)
Here φi is the notation for the i th normal mode. This is obtained from Ui on dividing through by a
scale factor ci , chosen according to some normalization criteria.
Normal modes enjoy the following orthogonality properties with respect to the mass and stiffness
matrices†
φiT M φ j = 0, φiT K φ j = 0, i = j. (19.22)
If i = j, the foregoing quadratic forms provide two important quantities called the generalized
mass Mi and the generalized stiffness K i :
Mi = φiT M φi > 0, K i = φiT K φi ≥ 0, i = 1, 2, . . . n, (19.23)
These are also called the modal mass and the modal stiffness, respectively, in the structural dynamics
literature. The squared natural frequencies appear as the ratio of generalized stiffnesses to the
corresponding generalized masses:
Ki
ωi2 = , i = 1, 2, . . . n. (19.24)
Mi
This formula represents the generalization of the SDOF formula ωn2 = k/m to MDOF.
The definition of generalized mass Mi in (19.22) provides one criterion for eigenvector normaliza-
tion, which is of particular importance in modal response analysis. If the scale factors in (19.21)
are chosen so that
Mi = φiT M φi = 1, i = 1, 2, . . . n, (19.25)
† These properties are proven in Section 10.1.5 of the Craig-Kurdila textbook.
19–8
19–9 §19.3 GENERALIZATION TO ARBITRARY NUMBER OF DOF
the normal modes φi are said to be mass-orthonormal. If this is done, (19.24) reduces to
ωi2 = K i = φi K φi , i = 1, 2, . . . n. (19.26)
Remark 19.3. If there are coincident frequencies (i.e., multiple eigenvalues appear), assumption (II) no longer
holds. Eigenvector extraction becomes more complicated. For example if ω12 = ω22 , any linear combination
of U1 and U2 is also an eigenvector. Instead of an a single vector we have to deal with a subspace spanned by
two or more vectors. This case is treated in detail in the Linear Algebra literature, and will not be reviewed
further here. Suffice to say that, if assumption (I) holds, it is always possible* to construct a distinct set of
eigenvectors that satisfy the orthogonality properties (19.22) as well as the definitions (19.23) of generalized
mass and stiffness.
§19.3.3. Eigenvector Normalization Criteria

Three normalization criteria are commonly used in structural dynamics. They are summarized
here for convenience. The i th unnormalized eigenvector is generically denoted by Ui , whereas a
normalized one is called φi . We will assume that eigenvectors are real and that each entry has the
same physical unit. Furthermore, for method (III) we assume that the mass matrix M is PD.
(I) Unit Largest Entry. Search for the largest entry of Ui in absolute value, and divide Ui by it.
(This value must be nonzero, else Ui = 0, contradicting the definition of eigenvector.)

(II) Unit Length. Compute the Euclidean length of the eigenvector l = + UiT Ui , and divide Ui
by this value. The normalized eigenvector satisfies φiT φi = 1.
√i = Ui MUi , which must be positive under
T
(III) Unit Generalized Mass. Get the generalized mass M
the assumption that M is PD, and divide Ui by + Mi . The normalized eigenvector satisfies
φiT Mφi = 1. As noted above, the resulting eigenvectors are said to be mass-orthogonal.
If M is the identity matrix I, the last two methods coalesce.
For just showing mode shapes pictorially, the simplest normalization method (I) works fine. For
the modal analysis presented in the next two Lectures, however, the Unit-Generalized-Mass nor-
malization offers computational and organization advantages, so it will be preferred one.
Example 19.2. Consider the two-DOF system developed in §19.2.1. As numerical values we take‡
m 1 = 2, m 2 = 1, c1 = c2 = 0, k1 = 6, k2 = 3, p1 = p2 = 0. (19.27)
The mass and stiffness matrices are

2 0 9 −3
M= , K= , (19.28)
0 1 −3 3
while the damping matrix C and dynamic force vector p vanish. The free vibrations eigenproblem associated
with (19.28) is

9 −3 U1 2 0 U1 9 − 2ω2 −3 U1 0
= ω2 , or = . (19.29)
−3 3 U2 0 1 U2 −3 3 − ω2 U2 0
* Using a procedure known as Gram-Schmidt orthogonalization.

‡ Note that no units are specified. It is tacitly understood that a consistent set of physical units, SI or English, is used
througout.
19–9
The characteristic polynomial equation is

9 − 2ω2 −3
det = 18 − 15ω2 + 2ω4 = (3 − 2ω2 )(6 − ω2 ) = 0. (19.30)
−3 3 − ω2
The roots of (19.30) give the two undamped squared frequencies
3
ω12 = = 1.5, ω22 = 6. (19.31)
2
√ √
The undamped circular natural frequencies (to 4 places) are ω1 = 3/2 = 1.225 and ω2 = 6 = 2.449. If
the physical system of units is consistent, these
√ are expressed in radians per second.√To convert to cycles per
−1 −1
second, or Hertz, divide by 2π: f 1 = (2π) 3/2 = 0.1949 Hz, and f 2 = (2π) 6 = 0.3898 Hz.
The calculation of eigenvectors is carried out in the next Lecture.
19–10
20
Modal Analysis
of MDOF Unforced
Undamped Systems
20–1
Lecture 20: MODAL ANALYSIS OF MDOF UNFORCED UNDAMPED SYSTEMS
TABLE OF CONTENTS
Page
§20.2. Modal Analysis of Unforced Undamped MDOF System 20–3
§20.2.1. Natural Frequencies . . . . . . . . . . . . . . . 20–4
§20.2.2. Vibration Modes . . . . . . . . . . . . . . . . 20–4
§20.2.3. Generalized Mass and Stiffness . . . . . . . . . . . 20–5
§20.2.4. Eigenvector Mass Orthonormalization . . . . . . . . . 20–5
§20.2.5. Generalized Coordinates and Modal Matrix . . . . . . . 20–6
§20.2.6. Modal Equations of Motion . . . . . . . . . . . . 20–6
§20.2.7. Modal Initial Conditions . . . . . . . . . . . . . . 20–7
§20.3. Unforced Response Example 20–7
20–2
;;
§20.2 MODAL ANALYSIS OF UNFORCED UNDAMPED MDOF SYSTEM
(a) (b)
k1 = 6
Static equilibrium
position Fs1 = k1 u1
u 1 = u1(t)
Mass m 1 = 2
..
Fs2 FI1 = m 1 u1
Static equilibrium k2 = 3
position
Fs2 = k2 (u2−u1)
u 2 = u2(t)
Mass m 2 = 1
..
x FI2 = m 2 u 2
Figure 20.1. Two-DOF, unforced, undamped spring-mass example system: (a) configuration, (b) DFBD.
This subject of this Lecture and of the next one is modal analysis. This is a technique by which the
equations of motion (EOM), which are originally expressed in physical coordinates, are transformed to
modal coordinates using the eigenvalues and eigenvectors gotten by solving the undamped frequency
eigenproblem. The transformed equations are called modal equations. In a mathematical context,
modal coordinates are also called generalized coordinates, or principal coordinates. For structural
dynamics they can be interpreted as response amplitudes of orthonormalized vibration modes.
The distinguishing feature of modal equations is that for an undamped system they uncouple. Con-
sequently, each modal equation may be solved independently of the others. Once computed, modal
responses may be transformed back to physical coordinates and superposed to produce the physical
response of the original system.
The method is a particular case of what is known in applied mathematics as orthogonal projection
methods. Instead of tackling the original equations directly, they are projected onto another space (the
‘modal space” in the case of dynamics) in which equations decouple.
§20.2. Modal Analysis of Unforced Undamped MDOF System

Here we retake the two-DOF example introduced in §19.1. As numerical values we take*
m 1 = 2, m 2 = 1, c1 = c2 = 0, k1 = 6, k2 = 3, p1 = p2 = 0. (20.1)
* Note that no units are specified; being tacitly understood that a consistent set of physical units, SI or English, is used througout.
20–3
The resulting undamped and unforced system is displayed in Figure 20.1. The known matrices and
vectors in the general EOM: Mü + Cu̇ + Ku = p, become

2 0 0 0 9 −3 0
M= , C= , K= , p= . (20.2)
0 1 0 0 −3 3 0
and the equations become

2 0 ü 1 9 −3 u1 0
+ = . (20.3)
0 1 ü 2 −3 3 u2 0
§20.2.1. Natural Frequencies

The free-vibration eigenproblem associated with (20.3) is

9 −3 U1 2 2 0 U1 9 − 2ω2 −3 U1 0
=ω , or = . (20.4)
−3 3 U2 0 1 U2 −3 3 − ω2 U2 0
The characteristic polynomial equation is

9 − 2ω2 −3
det = 18 − 15ω2 + 2ω4 = (3 − 2ω2 )(6 − ω2 ) = 0. (20.5)
−3 3 − ω2
The roots of (20.5) give the two squared frequencies
3
ω12 = = 1.5, ω22 = 6. (20.6)
2
§20.2.2. Vibration Modes

To get U1 , replace ω12 = 3/2 into the second of (20.4), set its first entry U11 to one, and solve for the
second entry U12 :

9−2 × (3/2) −3 U11 6 −3 U11 0 1
= = ⇒ U1 = . (20.7)
−3 3−3/2 U12 −3 3/2 U12 0 2
To get U2 , replace ω22 = 6 into the second of (20.5), set its first entry U21 to one, and solve for the
second entry U22 :

9−2 × 6 −3 U21 −3 −3 U21 0 1
= = ⇒ U2 = . (20.8)
−3 3−6 U22 −3 3 U22 0 −1
We now apply the first normalization method described in §19.3.3, by forcing the largest entry to be
+1:
1/2 1
φ1 = 2 U1 =
1
, φ2 = U2 . (20.9)
1 −1
The vibration modes are pictured in Figure 20.2. Note that in the first mode the masses oscillate in
phase while in the second one they are 180◦ out of phase, moving opposite to each other.
20–4
;; ;;
§20.2 MODAL ANALYSIS OF UNFORCED UNDAMPED MDOF SYSTEM
Mode 1, ω 21 = 3/2 Mode 2, ω22 = 6
1/2 Mass m 1 1
Red lines mark the

static equilibrium
positions of the masses
−1
1 Mass m 2
Figure 20.2. Vibration modes for example system.
It is convenient at this point to verify the orthogonality property: φiT Mφ j = 0 and φiT Kφ j = 0 for
i
= j. Here only one combination, namely i = 1 and j = 2, needs to be tested:

2 0 1 9 −3 1
φ1 M φ2 = [ 1/2 1 ]
T
= 0, φ1 K φ2 = [ 1/2 1 ]
T
= 0. (20.10)
0 1 −1 −3 3 −1
Note that is no need to explicitly check that φ2T M φ1 = 0 and φ2T M φ1 = 0 because (φ1T M φ2 )T =
φ2T MT φ1 = φ2T M φ1 , since M = MT if M is symmetric. Likewise for K.
The vibration modes (20.9) are orthogonal but not orthonormal with respect to the mass matrix. To
achieve that property it is necessary to rescale them so that generalized masses are unity, as done next.
§20.2.3. Generalized Mass and Stiffness
Next we compute the generalized masses and stiffnesses associated with the modes (20.10):

2 0 1/2
M1 = φ1 M φ1 = [ 1/2 1 ]
T
= 3/2,
0 1 1

2 0 1
M2 = φ2 M φ2 = [ 1 −1 ]
T
= 3,
0 1 −1
(20.11)
9 −3 1/2
K 1 = φ1 K φ1 = [ 1/2 1 ]
T
= 9/4,
−3 3 1

9 −3 1
K 2 = φ2 K φ2 = [ 1 −1 ]
T
= 18.
−3 3 −1
Quick check:
K1 9/4 K2 18
ω12 = = = 3/2, ω22 = = = 6, (20.12)
M1 3/2 M2 3
as expected. Note that should eigenvectors be normalized in a different way, M1 , M2 , K 1 , K 2 will be
generally different, but the ratios ωi2 = K i /Mi remain invariant.
20–5
§20.2.4. Eigenvector Mass Orthonormalization

Next we renormalize φ1 and φ2 to get unit generalized masses. These rescaled vectors will be initially
denoted by φ̃1 and φ̃2 , respectively. After renormalization the tildes are dropped.
T
Suppose φ̃1 = c1 φ1 . To find c1 , insert this into 1 = φ̃1 M φ̃1√= (c1 φ1T ) M (c1 φ1T ) = c12 φ1T M φ1T =
√ √
c12 × 3/2, whence c1 = 1/ 3/2 = 2/3. Likewise c2 = 1/ 3. Dropping the tildes for brevity, the
mass-orthonormalized vibration modes are

2 1/2 0.4088 1 1 0.5773
φ1 = = , φ2 = = (20.13)
3 1 0.8165 3 −1 −0.5773
After this renormalization,

3
φ1T M φ1 = 1, φ2T M φ2 = 1, φ1T K φ1 = = ω12 , φ2T K φ2 = 6 = ω22 . (20.14)
2
§20.2.5. Generalized Coordinates and Modal Matrix

We will express† the displacement vector u(t) in terms of normal coordinates η1 (t) and η2 (t), as
follows:

u 1 (t) def η1 (t)
u= = φ1 η1 (t) + φ2 η2 (t) = [ φ1 φ2 ] = Φ η. (20.15)
u 2 (t) η2 (t)
In a mathematical context, the ηi (t)s are called generalized coordinates or principal coordinates. They
represent the amplitude of the orthonormalized mode shapes, or modal amplitudes for short. Equation
(20.15) is an instance of modal superposition, which is a general feature of linear dynamical systems.
Two new matrix symbols appear in (20.15). The normal coordinate vector η collects η1 (t) and η2 (t) as
its entries. The modal matrix Φ is formed by stacking the mass-orthogonal eigenvectors as columns:
 1 
√ √1
0.4082 0.5773
Φ = [ φ1 φ2 ] =  6 3 = . (20.16)
√2 − √1 0.8165 −0.5773
6 3
Note that u̇ = Φ η̇ and ü = Φ η̈, because Φ does not depend on time. With the help of Φ, the
orthonormality conditions can be expressed compactly as‡
2
1 0 ω1 0 3/2 0
Φ M Φ = Mg = I =
T
, Φ K Φ = Kg = diag[ωi ] =
T 2
= . (20.17)
0 1 0 ω22 0 6
Here Mg and Kg denote the generalized mass matrix and generalized stiffness matrix, respectively.
If Φ is built with mass-orthogonalized eigenvectors, as in (20.16), Mg reduces to the identity matrix
while Kg becomes a diagonal matrix with squared frequencies stacked along its diagonal.
† This is a consequence of the so-called expansion theorem for MDOF systems; in the applied mathematics context it is known
as the spectral decomposition. See Section 10.1.9 of the Craig-Kurdila textbook for details.
‡ The Craig-Kurdila textbook uses bold italics symbols for generalized mass and stiffness matrices. Since those fonts are not
available to the document processor used for these Lectures, we shall use g subscripts to denote those quantities.
20–6
§20.3 UNFORCED RESPONSE EXAMPLE
§20.2.6. Modal Equations of Motion

To transform the undamped, unforced EOM: Mü + Ku = 0 to modal coordinates, replace u = Φη
and ü = Φη̈ into it, and premultiply by ΦT :
ΦT M Φ η̈ + ΦT K Φ η = ΦT 0 = 0. (20.18)
Using the generalized matrices introduced in (20.17) we get
η̈ + Kg η = 0, (20.19)
For the example problem,

1 0 η̈1 (t) 3/2 0 η1 (t) 0
+ = . (20.20)
0 1 η̈2 (t) 0 6 η2 (t) 0
Because the matrices in (20.20) are diagonal, (20.20) uncouples into two homogeneous, second-order
ODE already in canonical form:
η̈1 (t) + (3/2) η1 (t) = 0, η̈2 (t) + 6 η2 (t) = 0. (20.21)
Each ODE models an unforced, undamped SDOF oscillator. As such it may be treated using the
methods of Lecture 20 as long as initial conditions (IC) are known. Once solutions η1 (t) and η2 (t)
are available, they can be combined via the mode superposition relation (20.15) to get the physical
response u(t) = Φ η(t). But to solve (20.21) we need their IC in modal coordinates.
§20.2.7. Modal Initial Conditions
Suppose that the initial conditions for a MDOF system are
u(0) = u0 , u̇(0) = v0 . (20.22)
Here u0 and v0 denote vectors of initial displacements and velocities, respectively. Because u(t) =
Φ η(t) and u̇(t) = Φ η̇(t) for any time t, setting t = 0 we get
u0 = u(0) = Φ η(0), v0 = u̇(0) = Φ η̇(0). (20.23)
These can be solved for η0 = η(0) and η̇0 = η̇(0) by inverting the modal matrix:
η0 = η(0) = Φ−1 u0 , η̇0 = η̇(0) = Φ−1 v0 . (20.24)
or equivalently, solving the linear systems (20.23). But there is a more elegant scheme that only requires
matrix multiplication. Postmultiplying both sides of ΦT MΦ = I by Φ−1 gives Φ−1 = ΦT M, whence
η0 = ΦT M u0 , η̇0 = ΦT M v0 . (20.25)
Remark 20.1. Should the eigenvectors φi not be mass-orthonormal, we have the more complicated relations
Mg η0 = ΦT M u0 and Mg η̇0 = ΦT M v0 , in which Mg = ΦT MΦ is the generalized mass matrix. Consequently
η0 = M−1
g Φ M u0 ,
T
η̇0 = M−1
g Φ M v0 .
T
(20.26)
Since Mg is always diagonal, so is M−1

g . Consequently (20.26) merely amounts to scaling entries of η0 and η̇0 by
reciprocals of the generalized masses.
20–7
1 .
0.6 u1(t)
0.4 u1(t) 0.5
Displacements
0.2
Velocities
0
0
−0.2 −0.5
−0.4 .
u2(t) −1 u2(t)
−0.6
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
Figure 20.3. Response histories for example system under initial conditions u̇ 1 (0) = 1, others zero.
§20.3. Unforced Response Example

The following initial conditions (IC) are specified for the 2-DOF mass-spring system of Figure 20.1:
mass m 1 is given a unit initial velocity along +x, while all other initial values are zero. In matrix form:

0 1
u0 = , v0 = . (20.27)
0 0
The modal IC are obtained from (20.25) as
 1 
√ √2
η10  6 6  2 0 0 0
η0 = = Φ M u0 =
T
= ,
η20 1
√ −√ 1 0 1 0 0
3 3 (20.28)
 1   2 
√ √2 √
η̇10  6 6  2 0 1  6  0.8165
η̇0 = = Φ M v0 =
T
= = .
η̇20 √1 − √1 0 1 0 √2 1.1547
3 3 3
The solution of the modal equations (20.20) for the initial conditions (20.28) are
η̇10 2 √ 2
η1 (t) = sin ω1 t = √ sin( 6 t) = sin( 3/2 t) = 0.6667 sin(1.2247 t),
ω1 3 3 3
√ (20.29)
η̇20 2/ 3 √ 2 √
η2 (t) = sin ω2 t = √ sin( 6t) = √ sin( 6 t) = 0.4714 sin(2.4495 t).
ω2 6 3 2
These can be transformed back to physical coordinates via the modal matrix, as per u(t) = Φη(t).
Here are the calculations for the mass displacement responses, with numerical expressions listed to
4-place accuracy:
 1 
 √ 
√ √1 2 sin(√3/2 t) 1 2
sin(
√
3/2t) + sin( 6t)
u 1 (t) 3  3
= 6 √ = 
3 3
u 2 (t) 2
√ sin( 6 t) √ √
√2 − √1 1 2
2 sin( 3/2t) − sin( 6t)
3 2 3 3
6 3 (20.30)

0.2722 sin(1.2247t) + sin(2.4495t)
=
0.2722 2 sin(1.2247t) − sin(2.4495t)
20–8
§20.3 UNFORCED RESPONSE EXAMPLE
Mass velocities are obtained by direct numerical differentiation with respect to t:

1 cos(√3/2t) + 2 cos(√6t)
0.3333 cos(1.2247t) + 2 cos(2.4495t)
u̇ 1 (t) √ =
3
= √
u̇ 2 (t) 2
cos( 3/2t) − cos( 6t) 0.6667 cos(1.2247t) − cos(2.4495t)
3
(20.31)
The results (20.30) and (20.31) are plotted in Figure 20.3 over 0 ≤ t ≤ 12. Since there is no damping,
the harmonic oscillation patterns visible in the Figure will repeat with no distortions or amplitude
decay. In other words, the free vibrations will go on forever.* Observe that the initiaal conditions
(20.27) are correctly reproduced.
* The reason behind the repeating oscillation patterns in Figure 20.3 is that frequency ω2 is exactly twice ω1 in the example.
20–9
21
Modal Analysis
of Undamped Forced
MDOF Systems
21–1
Lecture 21: MODAL ANALYSIS OF UNDAMPED FORCED MDOF SYSTEMS
TABLE OF CONTENTS
Page
§21.2. Modal Analysis of Forced Undamped MDOF System 21–3
§21.2.1. Modal Forces . . . . . . . . . . . . . . . . . 21–4
§21.2.2. Modal Equations of Motion . . . . . . . . . . . . 21–4
§21.2.3. Forced Response Example . . . . . . . . . . . . . 21–5
§21.2.4. Frequency Response Functions . . . . . . . . . . . 21–7
§21.3. MDOF Under Prescribed Base Motion 21–8
§21.4. Base Motion EOM 21–8
§21.4.1. Modal Analysis . . . . . . . . . . . . . . . . 21–9
§21.4.2. Frequency Response Functions . . . . . . . . . . . 21–11
21–2
§21.2 MODAL ANALYSIS OF FORCED UNDAMPED MDOF SYSTEM
The modal analysis method presented in the previous Lecture for unforced, undamped, MDOF dy-
namical systems, is extended here to the case of applied forces as well as prescribed base motions.
It will be seen that extensions are relatively minor as long as the system is undamped. The main
new concept required is that of modal forces. This ingredient makes the modal equations non-
homogeneous, and generally requires computation of both homogeneous (transient) and particular
(steady-state) solutions. (In the unforced case, the homogeneous solution is sufficient.)
If the prescribed force or base motion is harmonic, a comparison of the response amplitude to that
of the static response will allow us to introduce the frequency response functions for MDOF. This
generalizes the SDOF results of Lecture 17. The FRF provide a quick visualization of resonances and
antiresonances. The latter is a new phenomenon that only arises in MDOF systems.
§21.2. Modal Analysis of Forced Undamped MDOF System
;;(a) (b)
k1 = 6
Static equilibrium
position Fs1 = k1 u1
u 1 = u1(t)
Mass m 1 = 2
p1 = p1(t) .. p1 = p1(t)
FI1 = m 1 u1
Static equilibrium Fs2
position k2 = 3
Fs2 = k2 (u2−u1)
u 2 = u2(t)
Mass m 2 = 1
p 2 = p2(t) .. p 2 = p2(t)
FI2= m2 u 2
x
Figure 21.1. Two-DOF, forced, undamped spring-mass example system: (a) configuration, (b) DFBD.
We again consider the two-DOF example of §22.1. As numerical values we take
m 1 = 2, m 2 = 1, c1 = c2 = 0, k1 = 6, k2 = 3, p1 = p1 (t), p2 = p2 (t). (21.1)
where p1 (t) and p2 (t) are specified force histories applied to masses m 1 and m 2 , respectively. The
resulting undamped but forced system is displayed in Figure 21.1. The known matrices and vectors
21–3
in the general MDOF equation of motion (EOM): Mü + Cu̇ + Ku = p, become

2 0 0 0 9 −3 p1 (t)
M= , C= , K= , p= . (21.2)
0 1 0 0 −3 3 p2 (t)
and the matrix EOM becomes

2 0 ü 1 9 −3 u1 p1 (t)
+ = . (21.3)
0 1 ü 2 −3 3 u2 p2 (t)
Comparing to the EOM in §23.1, it is seen that the only difference is the presence of specified forces
on the RHS. As a result the general solution is the sum of the homogeneous and particular solutions.
The modal analysis techniques described there can be reused with only minor modifications described
next.
§21.2.1. Modal Forces

It is convenient to write a system such as (21.2) in the compact matrix form
M ü(t) + K u(t) = p(t). (21.4)
Suppose that the modal eigenanalysis of the unforced version of (21.4), that is, M ü + K u = 0, has
been carried out as described in the previous two Lectures, and that the modal matrix Φ has been
constructed with the mass-orthogonalized eigenvectors φi as columns. Modal analysis relies on the
assumption introduced in §23.1.5 of Lecture 23:
u(t) = Φ η(t), (21.5)
where η(t) is the vector of modal amplitudes ηi (t). (The relation (21.5) is also called modal superpo-
sition, modal decomposition, or spectral decomposition in the literature.) Because Φ does not depend
on time, u̇ = Φ η̇(t) and ü = Φ η̈(t). Insert these into (21.4) and then premultiply through by Φt to
get
ΦT M Φ η̈(t) + ΦT K Φ η(t) = ΦT p(t). (21.6)
But ΦT M Φ = Mg and ΦT K Φ = Kg are the generalized mass and stiffness matrices, respectively,
which are diagonal. Furthermore if the eigenvectors stacked as columns of Φ have been mass-
orthonormalized as stated, Mg = I, which is the identity matrix, whereas Kg = diag(ωi2 ) has the
squared frequencies in its diagonal. To accomodate the RHS of (21.6) we introduce the definition
f(t) = ΦT p(t), (21.7)
The entries of f(t), denoted by f i (t), are called modal forces and also (in a wider context) generalized
forces. To show how those forces are obtained, consider the example system, in which the modal
matrix was displayed in §23.1.5 of Lecture 23. Applying (21.7) yields
 1   p1 (t) + 2 p2 (t) 
√ √2 √
f 1 (t)  6 6  p1 (t) =  6 .
= Φ p(t) =
T
(21.8)
f 2 (t) √1 − √1 p 2 (t) p 1 (t)√− p 2 (t)
3 3 3
21–4
§21.2.2. Modal Equations of Motion

Substitute (21.7) into (21.6), and account for the diagonal configuration of Mg = I and Kg = diag(ωi2 )
noted above. The matrix EOM in modal coordinates becomes
η̈(t) + diag(ωi2 ) η(t) = f(t). (21.9)
These are called the modal EOM. For a n-DOF dynamical system governed by the physical EOM
(21.4), (21.9) represents a set of n uncoupled, nonhomogeneous equations. To display this decoupling
clearly, consider the example system. Using results of the previous Lecture as well as (21.8), we
obtain the modal EOM in matrix form as

1 0 η̈1 (t) 3/2 0 η1 (t) f 1 (t)
+ = , (21.10)
0 1 η̈2 (t) 0 6 η2 (t) f 2 (t)
which uncouples into two nonhomogeneous, second-order ODEs:
√
η̈1 (t) + (3/2) η1 (t) = f 1 (t) = p1 (t) + 2 p2 (t) / 6,
√ (21.11)
η̈2 (t) + 6 η2 (t) = f 2 (t) = p1 (t) − p2 (t) / 3.
The solution of each ODE is the sum of its homogeneous solution, which depends on the initial
conditions (which can be obtained by the method presented in §23.1.7 of Lecture 23) and the particular
solution, which depends on the forcing terms f 1 (t) and f 2 (t).† Once the modal solutions are on hand,
they can be transformed to physical coordinates via u = Φ η.
§21.2.3. Forced Response Example

To illustrate the analysis process for the example problem, assume rest initial conditions and a harmonic
force of amplitude F2 , circular frequency , and cosine dependency, acting on mass 2:

0 0 0
p(t) = , u0 = , v0 = , (21.12)
F2 cos t 0 0
Here F2 and will be kept as variables during the analysis to help the√construction of frequency
√ function (FRF) graphs later. The modal forces are f 1 (t) = (2/ 6) F2 cos t and f 1 (t) =
response
−(1/ 3) F2 cos t. The modal equations are
√ √
η̈1 (t) + (3/2) η1 (t) = (2/ 6) F2 cos t, η̈2 (t) + 6 η2 (t) = −(1/ 3) F2 cos t. (21.13)
Because all IC are zero, the homogeneous solution of each ODE vanishes and only its particular
solution contributes. To find that quickly, note that the particular solution of η̈ + ω2 η = A cos t is‡
η p = A cos t/(ω2 − 2 ). Thus
√ √
(2/ 6) F2 (1/ 3) F2
η1 (t) = 2 cos t, η2 (t) = − 2 cos t, (21.14)
ω1 − 2 ω2 − 2
† Recall from Lecture 21 that this are called “transient” and “steady state” solutions, respectively, by structural engineers
when the system includes damping.
‡ To get that result, insert the guess η
guess = B cos t in η̈ + ω η = A cos t, and solve for B.
2
21–5
3 3
Displacements 2 u2(t) (a) Ω = 0 2 u2(t) (b) Ω = 1
Displacements
1 1
0 0
u1(t)
−1 −1
u1(t)
−2 −2
−3 −3
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
3 3
2 (c) Ω = 3/ 2 = 2.1213 ... 2 (d) Ω = 3
Displacements
Displacements
1 u2(t) = 0 1
u2(t)
0 0
−1 u1(t) −1 u1(t)
−2 −2
−3 −3
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
Figure 21.2. Response of harmonically excited example system for 4 values of the excitation frequency
: (a)
√= 0, a static force, (b) = 1, below both natural frequencies: masses oscillate in phase; (c)
= 3/ 2 ≈ 2.12132, which produces mass-2 antiresonance; and (d) = 3, higher than both natural
frequencies: masses oscillate 180◦ out of phase.
Here we have kept the natural frequencies ω1 and ω2 as generic to make the configuration of the
particular solutions more visible.
Combining the solutions (21.14) as per u(t) = Φ η(t) yields
 1   1 1 
√ √1 2 − 2
u 1 (t) 3  η1 (t) = 1 F cos t  ω1 −
2
ω2 − 
2
u(t) = = 6
u 2 (t) 2 1 η (t) 3 2 2 1
√ −√ 2 + 2
2
6 3 ω1
2
− ω2 − 2

F2 cos t ω22 − ω12 F2 cos t 6
= 2 2 = .
(ω1 − 2 )(ω22 − 2 ) ω1 + 2ω2 − 3 (6 − 2 )(3 − 2 2 ) 2(9 − 2 2 )
2 2
(21.15)
The last expression is obtained upon substituting ω1 = 3/2 and ω2 = 6.
2 2
Velocities are obtained by direct differentiation with respect to t:

u̇ 1 (t) F2 sin t 6
u̇(t) = =− . (21.16)
u̇ 2 (t) (6 − 2 )(3 − 2 2 ) 2(9 − 2 2 )
Figure 21.2 plots displacement responses (21.15) for F2 = 2, over time range 0 ≤ t ≤ 12, for the
following four values of the excitation frequency.
21–6
(a) (b)
1000
|SS magnification factor|

SS magnification factor
15 Resonance
Mass 1 100 at Ω = ω1
10 Mass 2 Mass 1
10 Resonance
5 Mass 2 at Ω = ω2
0 1
−5 0.1
Antiresonance
−10 0.01 of mass 2 at
−15 Ω = 2.12132
0.001
0 1 2 3 4 5 0.1 0.2 0.5 1 2 5 10
Excitation frequency Ω Excitation frequency Ω
Figure 21.3. Frequency Response Functions (FRF) of harmonically excited example system:
(a) natural scale plot, (b) log-log plot.
=0 Force p2 is static and equal to F2 = 2 for all time t. The static deflections are u 1 =
F2 /k1 = 1/3 and u 2 = F2 /k1 + F2 /k2 = 1.
√
=1 This excitation
√ is low frequency in the sense that < ω1 = 3/2 ≈ 1.224 and <
ω2 = 6 ≈ 2.449. The masses oscillate in phase.
√
=3/ 2 This excitation, with √numerical value ≈ 2.12132,√is intermediate frequency in the
sense that > ω1 = 3/2 ≈ 1.224 and < ω2 = 6 ≈ 2.449. Under this particular
it may be verified that u 2 (t) = 0 for all time since the second entry of the rightmost
vector in (21.15) is identically zero. Thus mass 2 never moves. This phenomenon is
called antiresonance.†
√
=3 √ is high frequency in the sense that
This excitation > ω1 = 3/2 ≈ 1.224 and
> ω2 = 6 ≈ 2.449. The masses oscillate iat 180◦ out of phase, moving in opposite
directions.
Note that vertical plot range: −3 to +3 is the same for all four plots in Figure 21.2, to facilitate visual
amplitude comparison.
§21.2.4. Frequency Response Functions
Recall from Lecture 21 that a Frequency Response Function (FRF) is the ratio Ds of the steady-state
displacement response amplitude to the static displacement, computed as a function of excitation
frequency. The FRF, also called steady-state magnification factor, pinpoints at a glance the location
of resonances. The log-log version of this plot is obtained by taking the absolute value of the
magnification factor (to avoid taking the log of negative numbers).
For a MDOF system, there are as many FRF as DOF. For our example the static displacements are
u st1 = F2 cos t/k1 = F2 cos t/6 and u st2 = F2 cos t (1/k1 + 1/k2 = F2 cos t/2. (These
values may be verified by making = 0 in (21.15).) The amplitude of the steady strate response
is obtained by setting cos t = 1 in (21.15), which gives U1 = 6F2 /(6 − 2 )(3 − 2 2 ) and
U2 = 2F2 (9 − 2 2 )/(6 − 2 )(3 − 2 2 ). Consequently the magnification factors at the point mass
† Antiresonance has important technical applications, such as vibration and noise mitigation.
21–7
;;; (a)
u b (t) : prescribed
(b)
k1 = 6
Static equilibrium
position Fs1 = k1 (u1 − u b )
u 1 = u1(t)
Mass m 1 = 2
..
FI1 = m 1 u1
Static equilibrium Fs2
position k2 = 3
Fs2 = k2 (u2−u1)
u 2 = u2(t)
Mass m 2 = 1
..
x FI2 = m2 u 2
Figure 21.4. Two-DOF spring-mass system under prescribed base excitation.
locations are
U1 18 U2 2(9 − 2 2 )
Ds1 = = , Ds2 = = . (21.17)
u st1 (6 − 2 )(3 − 2 2 ) u st2 (6 − 2 )(3 − 2 2 )
Figure 21.3(a) plots the FRF for masses 1 and 2 as the excitation frequency sweeps over the range
0 ≤ ≤ 5. Figure 21.3(b) is a log-log plot version that covers a wider excitation spectrum. This
version pinpoints resonances and the antiresonance better.
§21.3. MDOF Under Prescribed Base Motion

The two-DOF example problem is taken to be unforced, as in the previous Lecture, but under prescribed
base motion u b (t). See Figure 21.4(a).
§21.4. Base Motion EOM

Using the DFBD shown in Figure 21.4(b) we construct the matrix EOM

m1 0 ü 1 k1 + k2 −k2 u1 k1 u b (t)
+ = . (21.18)
0 m2 ü 2 −k2 k2 u2 0
As specific values we take those used in previous cases:
m 1 = 2, m 2 = 1, c1 = c2 = 0, k1 = 6, k2 = 3. (21.19)
21–8
§21.4 BASE MOTION EOM
The base motion, however, is for now left generic as well as k1 on the RHS, so (21.18) becomes

2 0 ü 1 9 −3 u1 k1 u b (t)
+ = . (21.20)
0 1 ü 2 −3 3 u2 0
§21.4.1. Modal Analysis

Since the LHS of (21.20) is the same as before, the modal analysis procedure described in S24.1 for
forced case applied with only cosmetic changes, since the prescribed base motion is equivalent to a
force p1 (t) = k1 u b (t) on mass 1 whereas p2 (t) = 0. In particular, the previously computed natural
frequencies and mass-orthogonalized mode shapes, as well as the modal matrix formed with the latter,
can be reused without change. Consequently the result (21.12) for the modal forces is applicable.
Replacing there p1 (t) = k1 u b (t) and p2 (t) = 0 gives the modal forces
 p1 (t) + 2 p2 (t)   1 
√ √
f 1 (t)
= Φ p(t) = 
T 6  = k1 u b (t)  6  (21.21)
f 2 (t) p1 (t)√− p2 (t) √1
3 3
The modal equations, leaving ω1 and ω2 symbolic for the moment, are
1 1
η̈1 (t) + ω12 η1 (t) = √ k1 u b (t), η̈2 (t) + ω22 η2 (t) = √ k1 u b (t). (21.22)
6 3
We will assume a state of rest at t = 0. Thus all IC (displacement and velocities) are zero, and
both modal solutions reduce to the particular, a.k.a. steady-state, component. To solve for this
we need a specific form for the RHS. To facilitate comparsion with the forced problem, we take a
harmonically varying base motion of excitation frequency , amplitude Ub , and cosine dependence,
whence u b (t) = Ub cos t and the modal EOM (21.22) become
1 1
η̈1 (t) + ω12 η1 (t) = √ k1 Ub cos t, η̈2 (t) + ω22 η2 (t) = √ k1 Ub cos t. (21.23)
6 3
The steady-state solutions are
√ √
(1/ 6) (1/ 3)
η1 (t) = 2 k1 Ub cos t, η2 (t) = 2 k1 Ub cos t. (21.24)
ω1 − 2 ω2 − 2
Combining these solutions as per u(t) = Φ η(t), and finally replacing k1 = 6, ω12 = 3/2 and ω22 = 6
yields
 1   2 1 
√ √1 2 + 2
u 1 (t) 3  η1 (t) = k U cos t  ω1 −
1 2
ω2 − 
2
u(t) = = 6
u 2 (t) 2 1 η (t) 1 b 1 1
√ −√ 2 − 2
2 6
6 3 ω 2
− ω − 2

1 2
(21.25)
2 + 1 
3/2 − 2 6 − 2  6Ub cos t 3 − 2
= Ub cos t  = .
1 − 1 (6 − 2 )(3 − 2 2 ) 3
3/2 − 2 6 − 2
21–9
4 4
Displacements 3
u1(t) u2(t) (a) Ω = 0 3 u2(t) (b) Ω = 1
Displacements
2 2
1 1
0 0
−1 −1 u1(t)
−2 −2
−3 −3
−4 −4
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
4 4
3 (c) Ω = 3 = 1.732 ... 3 (d) Ω = 3
Displacements
Displacements
2 2 u2(t)
1 1
0 0
−1 −1
−2 −2
u1(t)
u1(t)=0 u2(t)
−3 −3
−4 −4
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time t Time t
Figure 21.5. Response of harmonically base-excited example system for 4 values of the excitation
frequency : (a)√ = 0, a static force, (b) = 1, below both natural frequencies: masses oscillate in
phase; (c) = 3 ≈ 1.732, which produces mass-1 antiresonance; and (d) = 3, higher than both
natural frequencies: masses oscillate 180◦ out of phase.
Figure 21.5 plots displacement responses (21.24) for base amplitude motion Ub = 1, over time range
0 ≤ t ≤ 12, for the following four values of the excitation frequency.
=0 Base motion is static and equal to Ub = 1 for all time t. The static deflections are
u 1 = u 2 = Ub = 1, since the system moves as a rigid body. This can be readily verified
by setting = 0 in the last of (21.25).
√
=1 This excitation
√ is low frequency in the sense that < ω1 = 3/2 ≈ 1.224 and <
ω2 = 6 ≈ 2.449. The two masses oscillate in phase.
√
= 3 with numerical value ≈ 1.732,
This excitation, √ √ is intermediate frequency in the sense
that > ω1 = 3/2 ≈ 1.224 and < ω2 = 6 ≈ 2.449. Under this particular it
may be verified from (21.25) that u 1 (t) = 0 for all time, i.e., mass 1 never moves. As
noted in §24.1.3, this phenomenon is called antiresonance.
√
=3 √ is high frequency in the sense that ◦ > ω1 = 3/2 ≈ 1.224 and
This excitation
> ω2 = 6 ≈ 2.449. The masses oscillate at 180 out of phase, moving in opposite
directions.
Note that the vertical range: −4 to +4, is the same for all four plots in Figure 21.5 so as to facilitate
visual amplitude comparison.
21–10
§21.4 BASE MOTION EOM
|SS magnification factor|

SS magnification factor 20 Resonance
100 at Ω = ω1
Mass 1
10 Mass 2 10 Mass 1 Resonance
Mass 2 at Ω = ω2
1
0
0.1 Antiresonance
−10 of mass 1 at
0.01 Ω = 1.732
−20 0.001
0 1 2 3 4 5 0.1 0.2 0.5 1 2 5 10
Excitation frequency Ω Excitation frequency Ω
Figure 21.6. Frequency Response Functions (FRF) of harmonically base-excited example

system: (a) natural scale plot, (b) log-log plot.
§21.4.2. Frequency Response Functions

The steady-state amplification factors Ds1 and Ds2 are the ratio of the amplitudes of u 1 (t) and u 2 (t)
to their static values. As noted above, the latter are u st1 = u st2 = Ub = 1. The amplitudes can be
obtained from (21.25) by setting cos t = 1. Hence
3 − 2 3
Ds1 = , Ds2 = . (21.26)
(6 − 2 )(3 − 2 2 ) (6 − 2 )(3 − 2 2 )
The graph of Dsi versus excitation frequency are called the frequency response plots or FRF. For
this example they are shown in Figure 21.6. Figure 21.6(a) shows the two FRF for masses 1 and 2
as sweeps over the range 0 ≤ ≤ 5. Figure 21.6(b) is a log-log plot version that covers a wider
excitation spectrum; in this case the absolute value of the amplitude is taken. The log-log displays
both resonances and the antiresonance better.
21–11
22
Modal Analysis
of MDOF Forced
Damped Systems
22–1
Lecture 22: MODAL ANALYSIS OF MDOF FORCED DAMPED SYSTEMS
TABLE OF CONTENTS
Page
§22.2. What is Mechanical Damping? 22–3
§22.2.1. Internal vs. External Damping . . . . . . . . . . . . 22–3
§22.2.2. Distributed vs. Localized Damping . . . . . . . . . . 22–3
§22.3. Modeling Structural Damping 22–4
§22.4. Matrix Equations of Motion 22–5
§22.4.1. Equations of Motion Using Undamped Modes . . . . . . 22–5
§22.4.2. Three Ways Out . . . . . . . . . . . . . . . . 22–6
§22.5. Diagonalization by Modal Damping 22–7
§22.5.1. Damping Factor Guessing . . . . . . . . . . . . . 22–8
§22.5.2. Energy Equivalent Damping Factor . . . . . . . . . 22–8
§22.6. Diagonalization by Rayleigh Damping 22–8
§22.7. Remainder of Lecture 22–9
§22.8. Damped matrix EOM 22–9
§22.9. Diagonalization By Energy Balance, aka RQ Disgonalization 22–10
§22.10. Comparison of Exact Versus Diagonalized Responses 22–11
22–2
§22.2 WHAT IS MECHANICAL DAMPING?
The present lecture introduces damping within the context of dynamic modal analysis. After a brief
overview of a very complex but “fuzzy” subject, the lecture focuses on a very specific damping
model. Namely linearly viscous and light.
Linearly viscous damping is proportional to the velocity. Light damping means a damping fac-
tor (also called damping ratio and damping coefficient) that is small compared to unity. In the
terminology of Lecture 17, a lightly damped mechanical system is said to be underdamped.
Accounting for damping effects brings good and bad news. All real dynamic systems experience
damping because energy dissipation mechanisms are like death and taxes:* inevitable. Hence
inclusion makes the dynamic model more physically realistic. The bad news is that it can seriously
complicate the analysis process. Here is where the assumption of light viscous damping helps:
it allows the reuse of major parts of the modal analysis techniques introduced in the previous 3
lectures.
§22.2. What is Mechanical Damping?
Damping is the (generally irreversible) conversion of mechanical energy into heat as a result of
motion. For example, as we scratch a match against a rough surface, its motion generates heat
and ignites the sulphur content. When shivering under cold, we rub palms against each other to
warm up. Those are two classical examples of the thermodynamic effect of friction. In structural
systems, damping is more complex, appearing in several forms. These may be broadly categorized
into internal versus external, and distributed versus localized.
§22.2.1. Internal vs. External Damping

Internal damping is due to the structural material itself. Sources are varied: microstructural defects,
crystal grain relative motions, eddy currents in ferromagnetic media, dislocations in metals, and
molecular chain movements in polymers. The key macroscopic effect is the production of a
hysteresis loop in stress-strain plots. The hysteresis loop area represents the energy dissipated
per unit volume of material and per stress cycle. This kind of damping is intimately related to cyclic
motions such as vibration.
External damping come from boundary effects. An important form is structural damping, which is
produced by rubbing friction: stick-and-slip contact or impact. That may happen between structural
components such as joints, or between a structural surface and non-structural solid media such as
soil. This form is often modeled by Coulomb damping, which describes the energy dissipation of
rubbing dry-friction.
Another form of external damping is fluid damping. When a material is immersed in a fluid, such as
air or water, and there is relative motion between the structure and the fluid, a drag force appears. This
force causes an energy dissipation through internal fluid mechanisms such as viscosity, convection
or turbulence. This dissipation is collectively known as fluid damping. One well known instance
is a vehicle shock absorber: a fluid (liquid or air) is forced through a small opening by a piston.
* One exception: tax-exempt vampires.
22–3
§22.2.2. Distributed vs. Localized Damping

All damping ultimately comes from frictional effects, which may however take place at different
scales. If the effects are distributed over volumes or surfaces at macro scales, we speak of distributed
damping. But occasionally the engineer uses damping devices designed to produce beneficial
effects. For example: shock absorbers, airbags, parachutes, motion mitigators for buildings or
bridges in seismic or hurricane zones, active piezoelectric dampers for space structures.
Those devices can be often idealized as lumped objects, modeled as point forces or moments, and
said to produce localized damping. The distinction appears primarily at the modeling level, since
all motion-damper devices ultimately work as a result of some kind of internal energy conversion
at the molecular level.
Localized damping devices may be in turn classified into passive or active, with the latter responding
to motion feedback. But this would take us to far into control systems, which are not part of this
course.
§22.3. Modeling Structural Damping

The foregoing summary should make clear that damping is a ubiquitous but complicated business.
In structures containing joints, for example, Coulomb (dry friction) damping often dominates; this
model is partly nonlinear because the damping force depends on the sign of the velocity. Fluid
damping tends to be highly nonlinear if the interacting flow is turbulent, since in that case the drag
is nonlinear in the relative velocity between solid and fluid. Another modeling complication is that
friction may depend on fabrication or construction details that are not easy to predict; e.g., bolted
versus welded connections.
Balancing those complications is the fact that damping in most structures, especially metallic ones,
is light in the sense that the damping factor ξ introduced in Lecture 17 is much smaller than one. In
addition, the presence of damping is usually beneficial to safety in the sense that resonance effects
are mitigated. This gives structural engineers some leeway to simplify the dynamic analysis, by
proceeding as follows.
• A simple damping model, such as linear viscous damping, can be assumed without much
concern.
• Mode superposition is applicable because the EOM is linear. Moreover the frequencies
and mode shapes for the undamped system can be reused if additional assumptions, such
as Rayleigh damping, or modal damping, are made.
It should be stressed that such simplifications are not recommended if precise modeling of damping
effects is important to safety and performance. This occurs in the following scenarios:
• Damping is crucial to function or operation. Think, for instance, of a shock absorber. For-
tunately damper devices can be modeled more accurately than, say, dry friction, using manu-
facturer data.
• Damping may destabilize the system by feeding energy instead of removing it. This can
happen in active control systems and aeroelasticity.
The last two scenarios are beyond the scope of this course. In this lecture we will focus on linear
viscous damping. Moreover, damping levels will usually assumed to be light in the sense that the
damping factor ξ << 1.
22–4
;;
§22.4 MATRIX EQUATIONS OF MOTION
(a) (b)
k1 c1
Static equilibrium
position .
Fs1 = k1 u1 Fd1 = c1 u 1
u 1 = u1(t)
Mass m1
p1 (t) p1 (t) ..
FI1 = m 1 u1
Fs2 Fd 2
position
Fs2 = k2 (u2−u1) . .
Fd 2 = c2 (u 2 −u1)
u 2 = u2(t)
Mass m2
p2 (t) ..
p2 (t) FI2 = m 2 u 2
x
Figure 22.1. Two-DOF, forced, damped spring-mass example system: (a) configuration, (b) DFBD.
§22.4. Matrix Equations of Motion
Consider again the two-DOF mass-spring-dashpot example system of Lecture 19. This is repro-
duced in Figure 22.1 for convenience. The physical-coordinate EOM in detailed matrix notation
are

m1 0 ü 1 c1 + c2 −c2 u̇ 1 k1 + k2 −k2 u1 p1
+ + = . (22.1)
0 m2 ü 2 −c2 c2 u̇ 2 −k2 k2 u2 p2
In compact notation,
M ü + C u̇ + K u = p. (22.2)
where M, C and K are the mass, damping and stiffness matrix, respectively, p, u, u̇ and ü are the
force, displacement, velocity and acceleration vectors, respectively. The latter four are function of
time: u = u(t), etc, but the time argument will be often omitted for brevity.
In the sequel it will be assumed that M, C and K are symmetric. Furthermore M is positive definite
(PD) whereas K is nonnegative definite (NND).
§22.4.1. Equations of Motion Using Undamped Modes

This technique attempts to reuse the modal analysis methods introduced in Lectures 19–21. Suppose
that damping is removed whence C = 0. Get the natural frequencies and mode shapes of the
unforced, undamped system governed by M ü + K u = 0, by solving the eigenproblem K Ui =
22–5
ωi2 M Ui . Normalize the vibration mode shapes Ui ⇒ φi so that they are orthonormal with respect
to the mass matrix:
φiT M φi = δi j , (22.3)
in which δi j denotes the Kronecker delta: δi j = 1 if i = j, else zero. Let Φ be the modal
matrix constructed with the orthonormalized mode shapes φi , and denote by η the array of modal
amplitudes ηi , also called generalized coordinates. As before, assume mode superposition, so that
physical DOF are linked to modal amplitudes via
u = Φ η. (22.4)
Following the same scheme as in the previous two Lectures, the transformed EOM in generalized
coordinates are
ΦT M Φ + ΦT C Φ + ΦT KΦ = ΦT p(t). (22.5)
Define the generalized mass, damping, stiffness and forces as
Mg = ΦT M Φ, Cg = ΦT C Φ, Kg = ΦT K Φ, f = ΦT p(t). (22.6)
Of these, the generalized mass matrix Mg and the generalized stiffness matrix Kg were introduced
in Lecture 20. If Φ is built by stacking mass-orthonormalized vibration modes as columns, it was
shown there that
Mg = I, Kg = diag[ωi2 ]. (22.7)
That is, Mg reduces to the identity matrix while Kg becomes a diagonal matrix with squared
frequencies stacked along its diagonal. The generalized forces f(t), also called modal forces, were
introduced in Lecture 21. The new term in (22.5) is the generalized damping matrix Cg = ΦT C Φ,
also called the modal damping matrix in the literature.
Substituting (22.7) into (22.6) we arrive at the modal EOM for the damped system:
η̈(t) + Cg η̇(t) + diag(ωi2 ) η(t) = f(t). (22.8)
Here we ran into a major difficulty: matrix Cg generally will not be diagonal. If that happens, the
modal EOM (22.8) will not decouple. We seem to have taken a promising path, but hit a dead end.
§22.4.2. Three Ways Out

There are three ways out1 of the dead end:
• Diagonalization. Keep working with (22.8), but make Cg diagonal through some artifice.
• Complex Eigensystem. Set up and solve a different (augmented) eigenproblem that diagonal-
izes two matrices that comprise M, C and K as submatrices. The name comes from the fact
that it generally leads to frequencies and mode shapes that are complex.
1 In the textbook of Craig-Kurdila the first two approaches listed above are called mode superposition through real modes
of the undamped system, and mode superposition through complex modes of the damped system, respectively. See their
Section 10.1.10 for details. Those mouthfuls are abbreviated above to more easily remembered labels.
22–6
§22.5 DIAGONALIZATION BY MODAL DAMPING
• Direct Time Integration or DTI. Integrate directly (22.2) numerically in time without passing
to modal coordinates.
Each approach has strengths and weaknesses. (Obviously, else we would mention just one.)
Diagonalization allows straightforward reuse of undamped frequencies and mode shapes, which
are fairly easy to obtain with standard eigensolution software. The uncoupled modal equations
often have physical interpretation and can be validated with experiments. Only real arithmetic is
necessary.
The down side is that we dont solve the original EOM (22.2) directly, so some form of approximation
is inevitable. This is counteracted by the fact that structural damping is often difficult to quantify
since it can come from many sources, as discussed earlier. Thus the approximation in solving (22.2)
may be tolerable in view of modeling uncertainties. This is particularly true if damping is light.
There are certain problems, however, in which diagonalization cannot adequately represent damping
effects within engineering accuracy. Three of them are:
(1) Structures with localized damper devices (e.g., shock absorbers, piezoelectric dampers)
(2) Structure-media interaction (e.g., building foundations, tunnels, aeroelasticity, marine struc-
tures, surface ships)
(3) Active control systems
In those situations one of the two remaining approaches must be taken.
The complex eigensystem method is mathematically irreproachable and can solve (22.2) without
additional approximations. No assumptions as to light versus heavy damping are needed. But
it involves a substantial amount of preparatory work since the EOM must be transformed to the
so-called state space form. For a large number of DOF, solving complex eigensystems is unwieldy.
Physical interpretation of complex frequencies and modes is less immediate and may require sub-
stantial expertise in math as well as engineering experience. But in the three scenarios listed above
it should be preferred to diagonalization.
Direct time integration (DTI) has the advantage of being completely general. In fact it can handle not
only the linear EOM (22.2), but also nonlinear systems, as well as non-viscous damping mechanisms.
No transformations to mode coordinates is required, and no complex arithmetic emerges. The main
disadvantage is that it requires substantial expertise in computational handling of ODE, which is a
hairy topic onto itself. Since DTI can only handle fully numerically-specified models, the approach
is not particularly useful during preliminary design stages when many free parameters float around.
Because the last two approaches lie outside the scope of an introductory course (they are usually
taught at the graduate level), our choice is easy: diagonalization it is.
§22.5. Diagonalization by Modal Damping

In this approach the generalized damping matrix Cg = ΦT C Φ is assumed to be diagonal from the
start by using modal damping factors
Cg = diag[ 2ξi ωi ] Mi i = 1, 2, . . . n. (22.9)
in which n is the number of DOF and ξi denotes the damping factor for the i th natural mode. As
a result, the modal EOM (22.8) decouple, reducing to n canonical second-order equations in the
22–7
modal amplitudes
η̈i + 2ξi ωi η̇i + ωi2 η = f i (t), i = 1, 2, . . . n. (22.10)
These equations can be solved using the methods described in Lecture 21. The solutions can be
superposed via the mode decomposition assumption u(t) = Φ η(t) to get the physical response.
The method is straightforward. Two technical difficulties remain. First, how are modal damping
factors picked? Second, what is the error incurred by the decoupling assumption?
§22.5.1. Damping Factor Guessing

One time honored approach is educated guessing. This is of course the only possibility if (a) little is
known about the damping level and sources and (b) there are no experimental results; for example
the structure only exists on paper. The structural engineer then makes recourse to experience with
similar systems. With an air of authority she says: “No problema. Let’s assign 1% to modes 1
through 5, 2% to modes 6 through 20, and 4% to all higher ones.” Done.
There is some method in the madness. First, damping factors of well constructed structures are
typically small compared to unity: 1 to 5% is typical. Second, damping generally increases with
frequency, the reason being that more hysteresis cycles take place within a fixed time interval.
Third, the effect of light damping is not significant on the response unless the force (or base motion)
happens to excite a resonance; so the difference between, say, 1% and 2% may be well within
modeling uncertainty.
Note that this procedure directly constructs a diagonal Cg rather than C. If it is desirable to obtain
C a posteriori, one may use the following relations,the second of which assumes Mg = I:

n
−T −1
C=Φ Cg Φ = M Φ Cg Φ M = T
2ξi φi (Mφi ) (Mφi )T . (22.11)
i=1
Taking account of mode orthogonality gives the useful relation
φiT Cφi = 2ξi ωi . (22.12)
§22.5.2. Energy Equivalent Damping Factor

This technique is applicable if C is available. It can be shown that the power dissipated by viscous
damping in an unforced system governed by the physical-coordinates EOM (22.6) with p = 0 is
D= 1
2
u C u̇. (22.13)
The kinetic and potential energies of the system are
T = 1
2
u̇ M u̇, V = 1
2
u K u, (22.14)
Now suppose that the structure is moving in one of the undamped modes, say φi , The damping
factor can be chosen so that the energy dissipated in that mode matches that of the full discrete
system. This leads to the Rayleigh Quotient rule. (Topic to be developed further, since it is not in
C-K.)
22–8
§22.8 DAMPED MATRIX EOM
§22.6. Diagonalization by Rayleigh Damping
This diagonalization procedure is widely used for civil structures, especially for seismic response
calculations. The viscous damping matrix C is directly defined as a linear combination of the mass
and stiffness matrix:
C = a0 M + a1 K. (22.15)
in which a0 and a1 are selected constants (with appropriate physical dimensions). The method
is also called proportional damping in the literature. Applying the modal matrix congruential
transformation to this C results in
Cg = ΦT C Φ = diag[Cr ] = diag[a0 + a1 ωi2 ] = diag[2ξi ωi ], (22.16)
whence the effective modal damping factor is

1 a0
ξi = + a1 ωi . (22.17)
2 ωi
Choosing the damping factor for two modes of different frequencies and solving (22.17) for a0 and
a1 yields C from (22.15). In practice the stiffness proportional term is more physically relevant
(because, as previously noted, damping usually increases with frequency). Thus structural engineers
tend to adjust the choice of those frequencies so that ξi is roughly minimized for the lowest frequency
mode.
§22.7. Remainder of Lecture
The reminder of this lecture just presents the spring-mass-dashpot example sytem of previous
lectures. For various damping levels, this is solved exactly (as a coupled system, by direct numerical
integration) and also modally, with a particular damping-matrix-diagonalization method. The
results are compared for various damping levels.
§22.8. Damped matrix EOM
Consider again the two-DOF mass-spring-dashpot example system of Lecture 19. This is repro-
duced in Figure 22.1 for convenience. The physical-coordinate EOM in detailed matrix notation
are

m1 0 ü 1 c1 + c2 −c2 u̇ 1 k1 + k2 −k2 u1 p1
+ + = . (22.18)
0 m2 ü 2 −c2 c2 u̇ 2 −k2 k2 u2 p2
Suppose m 1 = 2, m 2 = 1, k1 = 6, k2 = 3, c2 = c1 = c, p1 = 0, and p2 = F2 cos t, in

which c (a dashpot coefficient), (excitation frequency of harmonic force), and F2 (harmonic
force amplitude) are for the moment kept as free parameters. Then

2 0 ü 1 2c −c u̇ 1 9 −3 u1 0
+ + = . (22.19)
0 1 ü 2 −c c u̇ 2 −3 3 u2 F2 cos t
22–9
;;
(a) (b)
k1 c1
Static equilibrium
position .
Fs1 = k1 u1 Fd1 = c1 u 1
u 1 = u1(t)
Mass m1
p1 (t) p1 (t) ..
FI1 = m 1 u1
Fs2 Fd 2
position
Fs2 = k2 (u2−u1) . .
Fd 2 = c2 (u 2 −u1)
u 2 = u2(t)
Mass m2
p2 (t) ..
p2 (t) FI2 = m 2 u 2
x
Figure 22.2. Two-DOF, forced, damped spring-mass example system: (a) configuration, (b) DFBD.
The data matrices and vectors are

2 0 2c −c 9 −3 0
M= , C= , K= , p= . (22.20)
0 1 −c c −3 3 F2 cos t
The undamped circular natural frequencies and associated modal matrix (orthonormalized with
respect to the mass matrix) were obtained in Lecture 20 and are reproduced here for convenience:
 1 
√ √1
3  6 3  0.4082 0.5773
ω1 = , ω2 = 6, Φ = [ φ1 φ2 ] =
2 2
= . (22.21)
2 √2 − √1 0.8165 −0.5773
6 3
With the foregoing choice of Φ, the generalized matrices and forces are
 c c 
3 − √
1 0 3 2
Mg = ΦT M Φ = = I, Cg = ΦT C Φ =  ,
0 1 − √ c 5c
2

3 2 2 (22.22)
√
3/2 0
Kg = Φ K Φ =
T
= diag[3/2, 6], f(t) = Φ p(t) = F2 cos t
T 6
0 6 1
√
3
The modal EOM in terms of the amplitudes η1 (t) and η2 (t) collected in vector η, are
η̈(t) + Cg η̇(t) + diag[3/2, 6] η(t) = f(t). (22.23)
Since Cg is not diagonal if c > 0, the modal EOM (22.23) are coupled through that matrix.
22–10
§22.10 COMPARISON OF EXACT VERSUS DIAGONALIZED RESPONSES
§22.9. Diagonalization By Energy Balance, aka RQ Disgonalization

We will investigate the following “Rayleigh quotient” (RQ) diagonalization technique†

C1R Q 0 2cφ1T Cφ1 φ2T Cφ2 5c
CgR Q = diag[C1R Q , C2R Q ] = RQ , C1R Q =
T
, C RQ
2 = = T
= .
0C2 φ1 φ1 5 φ2 φ2 2
(22.24)
This relation is obtained by equating the energy dissipated by damping over one cycle. Replacing
Cg with CgR Q gives the “RQ modal EOM”

2c 5c
η̈(t) + diag , η̇(t) + diag[3/2, 6] η(t) = f(t). (22.25)
5 2
These now decouple to

2c 3 2 5c 1
η̈1 (t) + η̇1 (t) + η1 (t) = √ F2 cos t, η̈2 (t) + η̇2 (t) + 6η2 (t) = − √ F2 cos t.
5 2 6 2 3
(22.26)
The effective modal damping factors are
2c 5c
ξ1R Q = = 0.1633 c, ξ2R Q = = 0.5103 c. (22.27)
5(2ω1 ) 2(2ω2 )
Ther IC required for solving the ODE (22.26) are, from Lecture 20, η(0) = ΦT M u(0) and
η̇(0) = ΦT M u̇(0), and replacing numbers:
 1   2(u 1 (0) + u 2 (0)) 
√ √2 √
η10 2 0 u (0)
η0 = = 6 6  1
=  2(u (0) −6 u (0))  ,
η20 1
√ −√ 1 0 1 u 2 (0) 1 √ 2
3 3 3 (22.28)
 1   2(u̇ 1 (0) + u̇ 2 (0)) 
√ √2 √
η̇10 6  2 0 u̇ 1 (0)
η̇0 = = 6 =  2(u̇ (0) −6 u̇ (0))  ,
η̇20 1
√ −√ 1 0 1 u̇ 2 (0) 1 √ 2
3 3 3
Once the solutions η1R Q (t) and η2R Q are obtained thay can be combined via the modal matrix as
u R Q (t) = Φ η R Q to produce the response in physical coordinates.
§22.10. Comparison of Exact Versus Diagonalized Responses

This will be demonstrated in class using Mathematica from a laptop if there is time, else it will be
skipped. The damping level will be varied as well as the IC and forcing function.
† Do not confuse with the Rayleigh damping, also known as proportional damping, approach. Same name, different
matrices.
22–11
23
Stability Of
Structures:
Basic Concepts
23–1
Lecture 23: STABILITY OF STRUCTURES: BASIC CONCEPTS
TABLE OF CONTENTS
Page
§23.2. Terminology 23–3
§23.3. Testing Stability 23–3
§23.3.1. Stability of Static Equilibrium . . . . . . . . . . . . 23–3
§23.3.2. Stability of Dynamic Equilibrium . . . . . . . . . . 23–7
§23.4. Static Stability Loss 23–7
§23.4.1. Buckling Or Snapping? . . . . . . . . . . . . . . 23–7
§23.4.2. Response Diagrams . . . . . . . . . . . . . . . 23–8
§23.4.3. Primary Equilibrium Path and the Design Critical Load . . . 23–9
§23.4.4. Stability Models . . . . . . . . . . . . . . . . 23–9
§23.4.5. Stability Equations Derivation . . . . . . . . . . . . 23–10
§23.5. Exact Versus Linearized Stability Analysis 23–11
§23.5.1. Example 1: The HCR Column: Geometrically Exact Analysis . 23–11
§23.5.2. Example 1: The HCR Column: LPB Analysis . . . . . . 23–12
§23.5.3. Example 2: The PCR Column: Geometrically Exact Analysis . 23–13
§23.5.4. Example 2: The PCR Column - LPB Analysis . . . . . . 23–15
23–2
§23.3 TESTING STABILITY
This Lecture presents basic concepts on structural stability, describes procedures for testing stabil-
ity, classifies models and analysis methods, and concludes by contrasting exact versus linearized
determination of critical loads.
The term stability has both informal and formal meanings. As regards the former, the American
Heritage Dictionary lists the following three. 1. Resistance to sudden change, dislodgment, or
overthrow. 2a. Constancy of character or purpose: tenacity; steadfastness. 2b. Reliability;
dependability. Related verb: to stabilize. Related adjective: stable. Antonyms: stability loss,
instability, to destabilize, unstable.
The formal meaning is found in engineering and sciences, concerning stability of systems.1 Broadly
speaking, structural stability can be defined as the power to recover equilibrium. It is an essential
requirement for all structures. Jennings2 provides the following historical sketch:
“Masonry structures generally become more stable with increasing dead weight. However
when iron and steel became available in quantity, elastic buckling due to loss of stability
of slender members appeared as a particular hazard.”
§23.2. Terminology
For referential convenience, Tables 23.1 and 23.2 provide a dictionary of terms that will be used
in the four Lectures on stability. Particularly important are the often used configuration and state,
which are related but should not be confused. “Configuration,” which is generically denoted by
C, is qualitative: it describes the disposition or arrangement of system components. On the other
hand, “state” is quantitative; it is precisely determined by a set of state variables called degrees of
freedom, or DOF. This set is finite in discrete models and infinite in continuous models.
In this and following Lectures we shall only use physical kinematic variables, such as displacements
and rotations, as DOF. We shall generically denote the DOF variables by u. In discrete models, u
is simply the vector that collects all DOF, and is called the state vector. In continuous models, u is
a function of the spatial (position) coordinates.
§23.3. Testing Stability

The stability of a mechanical system, and of structures in particular, can be tested (experimentally
or analytically) by observing how it reacts when external disturbances are applied. Here we have
to distinguish between statics and dynamics.
1 As defined in Table 23.1, a system is a functionally related group of components forming or regarded as a collective
entity. This definition uses “component” as a generic term that embodies “element” or “part,” which connote simplicity,
as well as “subsystem,” which connotes complexity. In this course we shall be concerned about mechanical systems
governed by Newtonian mechanics, with focus on structures.
2 A. Jennings, Structures: From Theory to Practice, Taylor and Francis, London, 2004, Chapter 7.
23–3
Table 23.1. Terminology Related to Mechanical Systems
Term Definition
System A functionally related set of components regarded as a physical entity.

Configuration The relative disposition or arrangement of system components.
State The condition of the system as regards its form, structure or constitution.
Degrees of Freedom A set of state variables that uniquely characterizes the state. Abbr.: DOF.
Kinematic DOF A DOF that is directly linked to the system geometry; e.g., a displacement.
Input The set of all actions that can influence the state of a system, or component.
Output The set of all quantities that characterize the state of a system, or component.
Model A mathematical idealization of a physical system.
Discrete Model A model with a finite # of DOF. Often expressed as vector equations.
Continuous Model A model with an infinite # of DOF. Often expressed as a ODE or PDE.
Event A change in the state variables produced by an agent.
Behavior A pattern of events.
Reference State A state of the system adopted as base or origin to measure relative changes.
Often the same as the undeformed state (cf. Table 23.2).
Motion The change in system geometry, as measured from a reference state.
Kinematics The study of system motion independently of force agents.
Kinetics The study of forces as action agents, and their effect on the system state.
Kinematic Constraint Any condition that restricts the system motion, usually expressed in terms of
kinematic DOF. Also simply called constraint.
Environment A set of entities that do not belong to the system, but can influence its behavior.
Open System A system that is influenced by entities outside the system (its environment).
Closed System A system that is not affected by entities outside the system.
Interaction The mutual effect of a system component, or group of such components,
on other components.
Forces The action agents through which effects are transmitted between system
components, or between environment entities and system components.
Internal Forces Forces that act between system components.
External Forces Forces that act between environment entities and system components.
Constraint Force A force manifested by removing a constraint while keeping it enforced.
Reaction Force A constraint force that is an external force.
Applied Load An external force specified as data. Also simply called load.
Some of these terms are useful in a more general context; for example active control.
§23.3.1. Stability of Static Equilibrium
For simplicity we will assume that the structure under study is elastic, since memory and historical
effects such as plasticity or creep introduce additional complications such as path dependence,
which are beyond our scope. The applied forces are characterized by a load factor λ, also called a
load parameter or load multiplier. This value scales a set of reference loads to provide the actual
applied loads.
Setting λ = 0 means that the structure is unloaded and takes up an equilibrium configuration C0
called the undeformed state. Furthermore, we assume that this state is stable in the sense defined
later. As λ is monotonically varied away from 0 the structure deforms and assumes equilibrium
23–4
§23.3 TESTING STABILITY
Table 23.2. Terminology Related to Static Stability Analysis
Term Definition
Reference Loads A set of applied loads taken as reference for application of a load factor.
Load Factor A scalar, denoted by λ, which scales reference loads to get the actual
applied loads. Also called load parameter and load multiplier.
System Response Values of the DOF, or subset thereof, expressed as function of the load factor,
or of the load level if only one load is applied. Also simply called response.
Equilibrium State A state in which internal and external forces are in equilibrium.
The associated configuration is called an equilibrium configuration.
Undeformed State The equilibrium state under zero applied loads, or, equivalently, λ = 0.
The associated configuration is called an undeformed configuration.
Equilibrium Response A system response in which all states are equilibrium states.
State Space A RCC frame with a DOF subset as axes.
Response Space A RCC frame with the load factor as one axis, and a DOF subset as the others.
Response Plot A display of the system response in response space.
Equilibrium Path An equilibrium response viewed in response space.
Perturbation An externally imposed disturbance of an equilibrium state while actual
loads are kept fixed. It may involve application of forces or motions.
Allowed Perturbation A perturbation that satisfies kinematic constraints. Also called admissible
perturbation, and (in the sense of variational calculus) virtual variation.
Stability The ability of a system to recover an equilibrium state upon being
disturbed by any of the allowed perturbations.
Instability The inability of a system to recover an equilibrium state upon being
disturbed by at least one allowed perturbation.
Stable Qualifier for an equilibrium state, or configuration, at which stability holds.
Unstable Qualifier for an equilibrium state, or configuration, at which instability occurs.
Neutrally Stable Qualifier for an equilibrium state, or configuration, at which transition
between stability and instability occurs.
Critical A qualifier that flags the occurrence of neutral stability. Applicable to state,
configuration, load, and load factor. For example: critical load.
Critical Point In a equilibrium response plot, a location where a critical state occurs.
Bifurcation Point A critical point at which two or more equilibrium paths cross.
Limit Point A critical point at which the load factor reaches a maximum or minimum.
Buckling Name used by structural engineers for the occurrence of a bifurcation point.
Snapping Name used by structural engineers for the occurrence of a limit point.
Also called snap-through, snap buckling, and snap-through buckling.
configurations C(λ). These are assumed to be (I) continuously dependent on λ and (II) stable for
sufficienttly small |λ|.
How is stability tested? Freeze λ at a value, say λd , in which d connotes “deformed.” The associated
equilibrium configuration is Cd = C(λd ). Apply a perturbation to Cd , and remove it. What sort
of perturbation? Any action that may disturb the state, for example a tiny load or a small imposed
motion. It must meet two conditions: any kinematic constraints (for example structural supports)
are satisfied, and the applied loads are kept fixed. Such perturbations are qualified as allowed or
admissible, as noted in Table 23.2.
23–5
Motion oscillates about

equilibrium configuration S: Stable*
Apply an
allowed or decays toward it
perturbation
Transition
Equilibrium Subsequent between stable N: Neutrally stable
configuration motion outcome and unstable
Motion is either unbounded, or

oscillates about, or decays toward, U: Unstable
another equilibrium configuration
* Strictly speaking, S requires stability
for all possible admissible perturbations
Figure 23.1. Stability test outcomes.
Applying and removing an allowed perturbation will trigger subsequent motion of the system.
Three possible outcomes are sketched in Figure 23.1.
S: Stable For all admissible perturbations, the structure either returns to the tested configu-
ration Cd or executes bounded oscillations about it. If so, the equilibrium is called
stable.
U: Unstable If for at least one admissible perturbation the structure moves to (decays to, or
oscillates about) another configuration, or “takes off” in an unbounded motion,
the equilibrium is unstable.
N: Neutral The transition from stable to unstable occurs at a value λcr , which is called the
critical load factor. The configuration Ccr = C(λcr ) at the critical load factor is
said to be in neutral equilibrium. The quantitative determination of this transition
is a key objective of the stability analysis.
The foregoing classification has gaps and leaves some details unanswered.
First, speaking about “moving” or “returning” introduces time into the picture. Indeed the concept
of stability is necessarily dynamic in nature3 There is a before: the act of applying the perturbation
to the frozen configuration, and an after: what happens upon removing it. Many practical methods
to assess critical loads, however, factor out the time dimension as long as certain conditions4 are
verified. Those are know as static criteria.
Second, the concept of “perturbation” as “small imposed change” is imprecise. How small is a
“tiny load” or a “slight deflection”? The idea is made more mathematically precise later when we
3 For example, Bazant and Cedolin in Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, Dover,
2003, comment on p. 144: “Failure of structures is a dynamical process, and so it is obviouly more realistic to approach
buckling and instability from a dynamical point of view.”
4 E.g., conservative loading: applied loads derive from a potential.
23–6
§23.4 STATIC STABILITY LOSS
introduce linearized stability, also called “stability in the small.” This is a natural consequence of
assuming infinitesimal configuration changes.
§23.3.2. Stability of Dynamic Equilibrium
Stability of motion is a more general topic that includes the static case as a particular one. (As
previously noted, the concept of stability is essentially dynamic in nature.)
Suppose that a mechanical system is moving in a predictable manner. For example, a bridge
oscillates under wind, an airplane is flying a predefined trajectory under automatic pilot, a satellite
orbits the Earth, the Earth orbits the Sun. What is the sensitivity of such a motion to changes of
parameters such as initial conditions? If the system includes stochastic or chaotic elements, like
turbulence, the analysis will require probabilistic methods.
To make such problems mathematically tractable it is common to restrict the kind of motions in such
a way that a bounded reference motion can be readily defined. For example, a bounded periodic
motion of structure oscillating under periodic excitation. Departures under parametric changes are
studied. Transition to unbounded or unpredictable motion is taken as a sign of instability.
An important application of this concept are vibrations of structures that interact with external
or internal fluid flows: bridges, buildings, airplanes, fluid pipes. The steady speed of the flow
may be taken as parameter. At a certain flow speed, increasing oscillations may be triggered:
this is called flutter. Or a non-oscillatory unbounded motion happens: this is called divergence. A
famous example of flutter in a civil structure was the collapse of the newly opened Tacoma-Narrows
suspension bridge near Seattle in 1940 under a moderate wind speed of about 40 mph.
Modeling and analysis of dynamic instability is covered in other courses, primarily at the graduate
level because it requires fancier mathematical tools and heavier use of complex analysis. In this
course we will consider only static stability. Moreover, the class of problems, methods and examples
will be severely constrained so as to fit within four lectures.
§23.4. Static Stability Loss

As just noted, we restrict attention to stability of static equilibrium. Two behavioral assumptions
are introduced:
Linear Elasticity. The structural material is, and remains, linearly elastic. Displacements and
rotations, however, are not necessarily small.
Conservative Loading. The applied loads are conservative, that is, derivable from a potential.
For example, gravity and hydrostatic loads are conservative. On the other hand, aerodynamic and
propulsion loads (wind gusts on a bridge, rocket thrust, etc) are often nonconservative.
The main reason for the second restriction is that loss of stability under nonconservative loads is
inherently dynamic in nature, and thus lies beyond our scope.
§23.4.1. Buckling Or Snapping?
Under the foregoing restrictions, two types of instabilities may occur:
Bifurcation. Structural engineers use the more familiar name buckling for this one. The structure
reaches a bifurcation point, at which two or more equilibrium paths intersect. What happens after
the bifurcation point is traversed is called post-buckling behavior.
23–7
Load or load (a) Load or load (b)

parameter parameter Terms in red are
Limit point those in common
Equilibrium (snapping) use by structural
path engineers
L Equilibrium paths
B
Initial linear
response
Bifurcation
Representative point (buckling) Representative
R deflection R deflection
Reference state Reference state
Figure 23.2. Graphical representation of static equilibrium paths and their critical points: (a) a response path with no
critical points; (b) multiple response paths showing occurrence of two types: bifurcation point (B) and limit point (L).
Snapping. Structural engineers use the term snap-through or snap buckling for this one. The
structure reaches a limit point at which the load, or the loading parameter, reaches a maximum.
What happens after the limit point is traversed is called post-snapping behavior.
Bifurcation points and limit points are instances of critical points. The importance of critical points
in static stability analysis stems from the following property:
Transition from stability to instability can only occur at critical points
Reaching a critical point may lead to immediate destruction (collapse) of the structure. This depends
on its post-buckling or post-snapping behavior, and nature of the material. For some scenarios the
knowledge of such behavior is important since immediate collapse may lead to loss of life.
On the other hand, there are some configurations where the structure keeps resisting significant
— or even increasing — loads after traversing a critical point. Such “load-sustaining” designs are
obviously preferable from a safety standpoint.
§23.4.2. Response Diagrams

To illustrate the occurrence of static instability as well as critical points we will often display load-
deflection response diagrams.5 This is a plot of equilibrium configurations taken by a structure
as a load, or loading parameter, is gradually and continuously varied. The load, or the loading
parameter λ, is plotted along the vertical axis while a judiciously chosen representative deflection,
which could be an angle, is plotted along the horizontal axes. A common convention is to take zero
deflection at zero load. This defines the reference state, labeled as point R in such plots.
A continuous set of equilibrium configurations forms an equilibrium path. Such paths are illustrated
in Figure 23.2. The plot in Figure 23.2(a) shows a response path with no critical points. On the
5 Students should be familiar with this visualization technique if they have done tension or torsion mechanical tests.
23–8
§23.4 STATIC STABILITY LOSS
Load or Load or
load factor (a) (b)
load factor
L L
B B
Critical load, Critical load,
or load factor, Primary or load factor, Primary
of interest equilibrium Representative of interest equilibrium Representative
path deflection path deflection
R R
Reference state Reference state
Figure 23.3. For design purposes, a critical load is that associated with the critical point first encountered when
traversing the primary equilibium path from the reference state. In (a) the critical load occurs at bifurcation point B
beacuse the limit point occurs later. In (b) the critical load occurs at the limit point L, because bifurcation occurs later.
other hand, that in Figure 23.2(b) depicts the occurrence of two critical points: one bifurcation and
one limit point. Those points are labeled as B and L, respectively, in the Figure.
§23.4.3. Primary Equilibrium Path and the Design Critical Load

A complex structure in general will exhibit multiple critical points, with a mixture of bifurcation
and limit points. An important question for design engineers is:
Which critical point should be chosen as determining the critical load to be

used for estimating the safety factor against instability?
Most textbooks say: pick the one associated with the lowest6 critical load or load factor. That is
fine if post-buckling behavior is not considered. For a more comprehensive answer, it is useful to
introduce the following definition:
The primary equilibrium path is the one that passes through the reference state (which is often the
same as the undeformed or unloaded state).
Then we define the design critical load as follows: the one located on the primary equilibrium path
that is nearest to the reference state. This choice makes engineering sense since most structures are
designed to operate on the primary equilibrium path while in service. The definition is illustrated
in Figure 23.3.
Note that in the case illustrated by Figure 23.3(b). the limit point L defines the design critical load
because it is encountered first while traversing the primary equilibrium path starting from R. It does
not matter that the bifurcation point B occurs at a lower load factor unless post-critical behavior is
important in design, which is rarely the case.
6 If the load factor can take either sign, as happens in some types of instability scenarios, lowest means in the sense of
absolute value, i.e., closest to zero.
23–9
§23.4.4. Stability Models

Stability models of actual structures fall into two categories:
Continuous. Such models have an infinite number of degrees of freedom (DOF). They lead to
ordinary or partial differential equations (ODEs or PDEs) in space, from which stability equations
may be derived by perturbation techniques. Obtaining nontrivial solutions of the perturbed equations
generally leads to trascendental eigenproblems, even if the underlying model is linear.
Discrete. These models have a finite number of DOF in space. Where do these come from?
Often they emerge as discrete approximations to the underlying continuum models. Two common
discretization techniques are:
(1) Lumped parameter models, in which the flexibility of the structure is localized at a finite
number of places. One common model of this type for columns: joint-hinged rigid struts
supported by extensional or torsional springs at the joints.
(2) Finite element models that include the so-called geometric stiffness effects.
Stability equations for discrete models may be constructed using various devices. For lumped
parameter models one may resort to either perturbed equilibrium equations built via FBDs, or to
energy methods. For FEM models only energy methods are practical. All techniques eventually
lead to matrix stability equations that take the form of an algebraic eigenproblem.
Both continuous and lumped-parameter discrete models for column problems are presented later
in this Lecture, and in the following two Lectures.
§23.4.5. Stability Equations Derivation

The equations that determine critical points are called characteristic equations in the applied math-
ematics literature.7 In structural engineering the names stability equations and buckling equations
are common. Two methods are favored for deriving those equations.
Equilibrium Method. The equilibrium equations of the structure are established in a perturbed
equilibrium configuration. This is usually done with Free Body Diagrams (FBD). The resulting
equations are examined for the occurrence of nontrivial solutions that produce disturbed equilibrium.
configurations. Those are obtained by disturbing original equilibrium positions through admissible
buckling modes. If those equilibrium equations are linearized for small perturbations, one obtains
an algebraic eigenproblem. The eigenvalues give values of critical loads while eigenvectors yield
the buckling mode shapes.
Energy Method. The total potential energy of the system is established in terms of the degrees
of freedom. The Jacobian matrix of the potential energy function, taken with respect to those
degrees of freedom, is established and tested for positive definiteness as the load parameter (or set
of parameters) is varied. Loss of such property occurs at critical points. These may be in turn
categorized into bifurcation and limit points according to a subsequent eigenvector analysis.
The energy method is more general for structures subject to conservative loading. It has two
major practical advantages: (1) merges naturally with the FEM formulation and so it can be
efficiently implemented in general-purpose codes, and (2) requires no a priori assumptions as
7 Often this term is restricted to the determinantal form of the stability eigenproblem. This is the equation whose roots
give the eigenvalues, which can be interpreted physically as critical loads, or as critical load parameters.
23–10
§23.5 EXACT VERSUS LINEARIZED STABILITY ANALYSIS
(a) (b) (c)

P = λ Pref λ
vA = L sin θ
A P = λ Pref
A Equilibrium path
A' of tilted column
θ
L Bifurcation
point
rigid MB = k θ L λcr = k /(LPref )
k Primary equilibrium
;;;;
path: vertical
B (untilted) column
B θ
P Stable
Unstable
Figure 23.4. Geometrically exact analysis of the hinged cantilevered rigid (HRC) column: (a) untilted
column, (b) tilted column with nonzero θ , (c) equilibrium paths intersecting at a bifurcation point.
regards admissible perturbations. But since energy methods are not covered at the undergraduate
level, only equilibrium methods will be presented here. Such techniques are necessarily restricted
to simple 1D problems amenable to FBDs, but those are sufficient to gain basic understanding of
stability analysis.
§23.5. Exact Versus Linearized Stability Analysis
As noted above, we will use only the equilibrium method to set up stability equations. Its key feature
is that FBDs must take the perturbed configuration into account. For certain simple problems it
is possible to establish the equations using the exact geometry of the deflected structure. Such
analyses will be called geometrically exact. A couple of examples follow to illustrate exact versus
linearized results.
§23.5.1. Example 1: The HCR Column: Geometrically Exact Analysis

Consider the configuration depicted in Figure 23.4(a). A rigid strut of length L stabilized by a
torsional spring of stiffness k > 0 is axially loaded by a vertical dead load P = λPr e f , in which
Pr e f is a reference load and λ a dimensionless load parameter. The load remains vertical as the
column tilts. (Observe that k has the physical dimension of force × length, i.e. of a moment.) This
configuration will be called a hinged cantilevered rigid column, or HCR column for brevity. The
definition P = λk/L renders λ dimensionless, which is convenient for result presentation. As state
parameter (and only DOF) we pick the tilt angle θ as most appropriate for the ensuing analysis.
For sufficiently small P the column remains vertical as in Figure 23.4(a), with θ = 0. The only
possible buckled shape is the tilted column shown in Figure 23.4(b). This Figure depicts the FBD
required to analyze equilibrium of the tilted column. Notice that θ is not assumed small. Taking
moments with respect to the hinge B as sketched in the figure, we obtain the following equilibrium
23–11
equation in terms of λ and θ:
k θ = P v A = λ Pr e f L sin θ ⇒ k θ − λ Pr e f L sin θ = 0. (23.1)
The equation on the right has two equilibrium solutions:
k θ
θ = 0 for any λ, λ= . (23.2)
Pr e f L sin θ
These pertain to the untilted (θ = 0) and tilted (θ = 0) equilibrium paths, respectively. Since
lim(θ/ sin θ) → 1 as θ → 0, those paths intersect when
k k
λcr = , or Pcr = . (23.3)
Pr e f L L
The two paths are plotted in Figure 23.4(c). The intersection (23.3) characterizes a bifurcation
point B. The diagram shows four branches emanating from B. Three are stable (drawn as full lines)
and one is unstable (drawn as dashed line). Note that the applied load may rise beyond the critical
Pcr = λcr Pr e f = k/L by moving to a tilted configuration. It is not difficult to show that the
maximum load occurs if θ → 180◦ , for which P → ∞; this is a consequence of the assumption
that the column is rigid, and that it may fully rotate by that amount about the hinge without being
impeded, say, by hitting the ground.
§23.5.2. Example 1: The HCR Column: LPB Analysis

The geometrically exact analysis that leads to (23.2) has the advantage of providing a complete
solution. In particular, it shows what happens after the bifurcation point B is traversed. For this
particular configuration the structure maintains load-bearing capabilities while tilted, which is the
hallmark of a safe design.
But for more complicated cases this approach becomes impractical because it involves solving
systems of nonlinear algebraic or differential equations. Even for the Euler column presented in
Lecture 25, a geometrically exact analysis leads to elliptic functions.
Often the engineer is interested only in the critical load. This is especially true in preliminary design
scenarios, when the main objective is to assess safety factors against buckling. If so, it is more
practical to work with a linearized version of the problem. Technically the full technical name is
linearized prebuckling (LPB) analysis. This approach relies on the following assumptions:
• Deformations prior to buckling are neglected. Consequently the analysis can be carried out in
the reference configuration (undeformed) geometry.
• Perturbations of the reference configuration are restricted to infinitesimal displacements and
rotations.
• The structure remains linearly elastic up to buckling.
• Both structure and loading do not exhibit any imperfections.
• The critical state is a bifurcation point.
23–12
We apply these rules to the HCR column of Figure 23.4(a). The equilibrium equation (23.1) is
linearized by assuming an infinitesimal tilt angle θ << 1, whence sin θ ≈ θ and that expression
becomes

k − λ Pr e f L θ = 0. (23.4)
This is the stability equation. Since the product of two numbers is zero, at least one must be zero:
θ = 0 or k − λ Pr e f /L = 0, or both. The solution θ = 0 reproduces the untilted configuration. For
buckling to occur, we must have θ = 0. If so, the expression in parenthesis must vanish, which
requires
k k
λ = λcr = or Pcr = λcr Pr e f = . (23.5)
Pr e f L L
This reproduces the critical load given in (23.3).
Note that LPB analysis only provides the critical load. It does not give any information on post-
buckling behavior. If this is necessary, a more comprehensive analysis, such as the geometrically
exact one carried out in the previous subsection, is required.
It should also be noted that if the loss of stability is by snap buckling, it cannot be obtained by LPB
analysis. The main reason is that deformations prior to buckling are essential in the determination
of limit points, whereas the first LPB assumption listed above explicitly precludes those.
§23.5.3. Example 2: The PCR Column: Geometrically Exact Analysis
Consider next the configuration pictured in Figure 23.5(a). This one differs from the HRC column
in the type of stabilizing spring. A rigid strut of length L is hinged at B and supports a vertical
load P = λPr e f at end A. The load remains vertical as the column tilts. The column is propped
by an extensional spring of stiffness k attached to A. This configuration will be called a propped
cantilevered rigid column; or PCR column for short. As before, the only DOF is the tilt angle θ.
For the geometrically exact analysis is it important to know what happens to the spring as the
column tilts. One possible assumption is that it remains horizontal, as pictured in Figure 23.5(b).
If so, the FBD in the tilted configuration will be as shown in Figure 23.6(a). Taking moments about
B yields the equilibrium condition
λPr e f sin θ = k L sin θ cos θ. (23.6)
(Important: do not cancel out sin θ from both sides of (23.6) yet — that will cause the primary
equilibrium path θ = 0 to be lost as a solution.) This equation has two equilibrium solutions
kL
θ = 0 for any λ , λ= cos θ. (23.7)
Pr e f
These solutions yield the vertical- and tilted-column equilibrium paths, respectively, which are
plotted in Figure 23.6(b) on the λ versus θ plane. The two paths intersect at θ = 0 and λ = λcr =
k L/Pr e f , which is a bifurcation point. Consequently the critical load is given by
kL
λcr = or Pcr = λcr Pr e f = k L . (23.8)
Pr e f
23–13
;; ;
(a) (b)
P = λ Pref vA = L sin θ
;; ;
C P = λ Pref
A A
C' uA = L (1−cos θ)
;; ;
k
θ A'
L
rigid L
;; B B
P
Figure 23.5. Geometrically exact analysis of a propped rigid cantilever (PRC) column with extensional
spring remaining horizontal: (a) untilted column, (b) tilted column.
(a) (b) λ
vA = L sin θ
P = λ Pref
A
A' Bifurcation
FA = k vA point Equilibrium path
θ Primary equilibrium of tilted column
path: vertical
(untilted) column
L cos θ
λ cr = k L/Pref
B −90 +90
θ
Stable
P Unstable
Figure 23.6. Geometrically exact analysis of the hinged cantilevered rigid (HRC) column: (a) untilted
column, (b) tilted column with nonzero θ , (c) equilibrium paths intersecting at a bifurcation point.
Of the four branches that emanate from the bifurcation point B, only one (the θ = 0 path for λ < λcr )
is stable. Once B is reached, the tilted column supports only a decreasing load P, which vanishes
at θ = ±90◦ . Consequently this configuration is poor from the standpoint of post-buckling safety.
Another reasonable assumption is that the spring attachment point to the wall, called C in Fig-
ure 23.7(a), stays fixed. The distance AC is parametrized with respect to the column length as ηL,
in which η is dimensionless. If the column tilts, the spring also tilts as pictured in Figure 23.7(b).
The geometrically exact FBD for this case is shown in Figure 23.7(c). As can be observed, it is
considerably more involved than for the spring-stays-horizontal case. We quote only the final result
23–14
(a) (b)
; ;
ηL P = λ Pref
vA = L sin θ
P = λ Pref
; ;
A A
A' uA = L (1− cos θ)
C C
k
θ
L
;;
rigid L
;;
B B
P
(c)
ϕ = angle A'CA
positive CW vA = L sin θ P = λ Pref
C A
A' uA = L (1− cos θ)
FA = k ds
ds : tilting θ
spring
elongation
B
P
Figure 23.7. Geometrically exact analysis of a propped rigid cantilever (PRC) column with wall-
attached extensional spring: (a) untilted column; (b) tilted column; (c) FBD for tilted equilibrium.
for the two equilibrium path equations:

k L η cos θ + sin θ
θ =0 for any λ, λ= . (23.9)
Pr e f η + sin θ
As in previous cases, the first solution corresponds to the untilted column whereas the second one
pertains to the tilted one. Figure 23.8 shows response plots on the λ versus θ plane for the tilted
column, that is, the second solution in (23.9), for the four cases η = 14 , η = 12 η = 1 and η = 5,
drawn with Pr e f = 1, k = 1, and L = 1.
Although the bifurcation points stay in the same location: λ = 1 and θ = 0, the post-buckling
response is no longer symmetric for with respect to θ. That deviation is most conspicuous when
η < 1, since if so the spring tilting has a highly noticeable effect if θ < 0. The sharp drop of λ
towards −∞ occurs when the spring and the tilted column are nearly aligned. Comparing the plot
for η = 5 with that in Figure 23.6(b), it is clear that as η >> 1 the response approaaches that of the
spring-stays-horizontal tilted-column in (23.7), as may be expected. This may be mathematically
proven by taking the limit of the second of (23.9) as η → ∞.
23–15
4 λ 2 λ
L
3
Note occurrence of η = 1/4 1.5 η = 1/2
limit points
2
B
Note occurrence of 1
limit points
B 0.5
1
L
−60 −40 −20 R 20 40 60 80
−60 −40 −20 R 20 40
angle θ (deg) −0.5 angle θ (deg)
−1
−1
B 1 B
1
λ λ
0.75 0.8
0.5
η=1 0.6 η=5
0.25
0.4
−60 −40 −20 R 20 40 60 80
−0.25
angle θ (deg) 0.2
−0.5
angle θ (deg)
−0.75
−60 −40 −20 R 20 40 60 80
−1
Figure 23.8. Geometrically exact analysis of a PRC column with wall-attached extensional spring: λ vs./ θ
response diagrams for tilted column and four values of η [defined in Figure 23.7(a)], with Pr e f = k = L = 1.
Note: the untilted-column equilibrium path θ = 0 is not plotted beyond B to reduce clutter.
§23.5.4. Example 2: The PCR Column - LPB Analysis

To linearize the PRC problem, assume again that θ is so small that sin θ ≈ θ and cos θ ≈ 1. Then
the equilibrium equations of the two foregoing cases (horizontal spring and wall-attached spring)
collapse to
(λ Pr e f − k L) θ = 0. (23.10)
The two solutions of (23.10) represent the equilibrium paths θ = 0 and λ = k L/Pr e f , which pertain
to the untilted (vertical) and tilted column, respectively. Consequently
kL
λcr = , or Pcr = λcr Pr e f = k L . (23.11)
Pr e f
This result is independent of assumptions on how the spring wall-attachment point behaves after
the column buckles. This is to be expected since linearization filters out that information. Once
again, the LPB analysis provides no information on post-buckling behavior.
23–16
24
Stability Of
Structures:
Discrete Models
24–1
Lecture 24: STABILITY OF STRUCTURES: DISCRETE MODELS
TABLE OF CONTENTS
Page
§24.2. Example 1: A Cantilevered Two-Strut Column 24–3
§24.2.1. Equilibrium Analysis . . . . . . . . . . . . . . . 24–3
§24.2.2. Critical Loads . . . . . . . . . . . . . . . . . 24–4
§24.2.3. Buckling Mode Shapes . . . . . . . . . . . . . . 24–5
§24.3. Example 2: A Pinned-Pinned Three-Strut Column 24–6
§24.3.1. Equilibrium Analysis . . . . . . . . . . . . . . 24–6
§24.3.2. Critical Loads . . . . . . . . . . . . . . . . . 24–7
§24.3.3. Buckling Mode Shapes . . . . . . . . . . . . . . 24–8
§24.4. Example 3: Strut Propped by Inclined Springs 24–8
§24.4.1. Equilibrium Analysis . . . . . . . . . . . . . . . 24–8
§24.4.2. Critical Load . . . . . . . . . . . . . . . . . 24–9
§24.4.3. Optimal Rise Angle . . . . . . . . . . . . . . . 24–9
24–2
§24.2 EXAMPLE 1: A CANTILEVERED TWO-STRUT COLUMN
This Lecture describes the formulation determination of critical loads of a discrete structrural system
by the LPB approach as a matrix algebraic eigenproblem. For a system with n degrees of freedom
(DOF), the eigensystem fits the following form
Ax = λBx (24.1)
Here A and B are square matrices of order n, the n eigenvalues λi are the critical load factors, and the
eigenvectors xi are the associated buckling shapes. Generally only the smallest eigenvalue, which will
be denoted by λ1 , will be of interest for engineering design.
Matrices A and B may be obtained by either equilibrium analysis, or by energy methods, as discussed
in the preceding Lecture. If the equilibrium approach is used, the matrices A and B are not necessarily
symmetric, but may be symmetrized if necessary. If the energy metod is used, the symmetry of both
matrices is guaranteed from the start.
In this course we will cover only the equilibrium method, which is based on doing Free Body Diagrams
(FBD) on a slightly perturbed configuration. This technique is illustrated through several examples.
In all cases the system is a column consisting of one or more rigid links (called “struts” by structural
engineers) linked at joints, and stabilized by extensional or torsional springs.
§24.2. Example 1: A Cantilevered Two-Strut Column

The problem investigated here in detail is the buckling of the two-hinged-cantilever column depicted
in Figure 24.1(a).1 Rigid links of equal length L/2 are connected at the joints A, B and at the bottom
C by frictionless hinges. The column is propped by two extensional springs attached at A and B with
the stiffnesses shown in the figure.
Data: Length L and spring stiffness k. Required: Set up the stability equations as a matrix eigen-
problem, find critical loads in terms of the data and display buckling shapes as separate diagrams.
Draw an arbitrary, but admissible2 buckling shape as in Figure 24.1(b). This shape is completely
defined by the two lateral deflections v A and v B of points A and B, respectively, positive to the right.
Note: The deflections are highly exaggerated in the figure for visibility; they are actually infinitesimally
small.
§24.2.1. Equilibrium Analysis

Figure 24.2(c1,c2,c3) shows free body diagrams (FBD) of links AB, BC and of the whole column
ABC, respectively. To facilitate visualization applied forces are pictured in black, spring forces in
blue, whereas internal reactions are shown in red.
The two lateral deflections v A and v B are taken as independent degrees of freedom (DOF). Consequently,
two equilibrium equations are required obtained from FBDs are required. Three combinations of two
FBDs are possible: AB and BC, AB and whole column, or BC and whole column. In all FBDs
translational force equilibrium holds identically by construction. The remaining moment equilibrium
conditions become part of the matrix stability equation. Since the three FBDs are linearly dependent
1 Adapted from E. P. Popov, Engineering Mechanics of Solids, Prentice Hall, 2nd ed., 1999.
2 Admissible in the sense stated in Table 23.2 of the previous Lecture. That is, the assumed buckling shape must satisfy all
kinematic compatibility and boundary conditions.
24–3
(a) P (b) vA
P
k
A A'
A
v deflections
L/2 + to the right
rigid
vB
3k/2 B
L B B'
L/2
rigid
C C 2 DOF: vA , v B
Figure 24.1. Cantilevered two-strut column. (a): structure; (b): deflected shape defined by two DOF v A and v B ;
P
vA P
A vA
A' FA A
(c1) A' FA
FBD of Applied forces (c3)
link AB in black, spring
vB forces in blue,
reactions in red
FA +FB B B' FB vB
B FB
P P B'
(c2) vB FA +FB FBD of whole
B B' column
FBD of FA +FB
link BC C
Restoring spring P
FA +FB forces + to the left Note split of
C
AF = k vA , FB = 3kv
B /2 spring force at B
P
Figure 24.2. Cantilevered two-strut column. (c1,c2,c3): FBDs of link AB, link BC and whole column,
respectively. FBD color convention: applied forces in black, spring foces in blue, reactions in red.
[check the figure: merging (c1) and (c2) gives (c3)], the results of the analysis should be exactly the
same, although the stability matrix equations will be different. For the ensuing derivation we pick (c1)
and (c3).
For link AB, taking moments with respect to the displaced hinge point B’, positive CW, gives

M B = P(v A − v B ) − FA ( 12 L) = P(v A − v B ) − 12 Lkv A = 0. (24.2)
For the whole column ABC, taking moments with respect to the hinge C, positive CW, gives

MC = Pv A − FA L − FB ( 12 L) = Pv A − k Lv A − 32 k( 12 L)v B = 0. (24.3)
24–4
§24.2 EXAMPLE 1: A CANTILEVERED TWO-STRUT COLUMN
Pcr = Pcr1 = 0.25 kL Pcr2 = 1.50 kL
1 1
A' A'
B' −1 2/3
B'
Figure 24.3. Buckling mode shapes for the cantilevered two-strut column.
§24.2.2. Critical Loads

Combining (24.2) and (24.3) in a matrix equation gives the stability equation

P − 12 k L −P vA 0
= , or Av = 0. (24.4)
P − k L − 34 k L vB 0
Matrix A is called the stability matrix or characteristic matrix. The system (24.4) has two kinds of
solutions. The trivial solution v = 0, or v A = v B = 0, corresponds to the undeflected (vertical)
column, which holds for any P. To find the critical loads we must have v = 0, that is, at least one of the
{v A , v B } must be nonzero. Consequently A must be singular or, equivalently, have zero determinant.
Since the entries of A depend on P, this is an algebraic eigenproblem in P.
The eigenvalues are the solution of the characteristic equation, which is obtained by setting the deter-
minant of A to zero:
det(A) = 0, or P 2 − 74 k L P + 38 k 2 L 2 = 0. (24.5)
Solving this quadratic equation for P gives the two critical loads:
Pcr 1 = 14 k L = 0.25 k L , Pcr 2 = 32 k L = 1.50 k L . (24.6)
The smallest one is the critical load:
Pcr = Pcr 1 = 14 k L = 0.25 k L . (24.7)
§24.2.3. Buckling Mode Shapes

Replacing Pcr 1 into (24.4) gives the following equation to determine the eigenvector v1 :

k L −3 −3 v A1 0 v A1 −1
= ⇒ v1 = = c1 (24.8)
4 −1 −1 v B1 0 v B1 1
in which c1 is an arbitrary nonzero scaling factor.

Likewise, replacing Pcr 2 into (24.4) gives the following equation to determine the eigenvector v2 :

k L 2 −3 v A2 0 v A2 3/2
= ⇒ v2 = = c2 (24.9)
4 4 −6 v B2 0 v B2 1
24–5
(a) P (b) P>Pcr (c) P

A A' A
vB
(1/3) (2FB +FC ) =
rigid
L/3 (1/3) k (2vB +vC )
k
B FB = k vB
B' B'
v deflections vC
rigid
L/3 + to the right
k
C FC = k vC
rigid
C' C'
L/3 (1/3) (FB +2FC ) =

(1/3) k (vB +2vC )
D D D
Figure 24.4. Stability analysis of three-strut column propped by extensional springs: (a) reference
configuration; (b) admissible buckled shape showing realistic geometry change (struts do not change length);
(c) linearized configuration: B and C displace by small horizontal movements whereas A stays fixed.
in which c2 is an arbitrary nonzero scaling factor.

A common normalization condition for eigenvectors is to make their largest component equal to +1.
This is done for (24.8) and (24.9) by taking c1 = 1 or c1 = −1 (either one works), and c2 = 2/3,
respectively. The two eigenvectors normalized as per that criterion are

1 1
v1 = for Pcr 1 , v2 = for Pcr 2 . (24.10)
−1 2/3
The buckling shapes defined by (24.10) are plotted in Figure 24.3.
§24.3. Example 2: A Pinned-Pinned Three-Strut Column
This is a slight varint of one treated in Jennings (op. cit. in footnote 1). Figure 24.4(a) shows an axially
loaded column built with three rigid struts of equal length. The struts are pinned at the joints. It is
stabilized by two extensional spring of stiffness k working as elastic supports at the internal joints B
and C. Determine the critical loads of this configuration by linearized equilibrium analysis.

Figure 24.4(b) depicts how the structure might realistically displace upon buckling. Because the struts
are rigid and cannot change length, joints A, B and C will necessarily move downward while the springs
tilt. Such a configuration would be correct for geometrically exact analysis.
The linearized version, which assumed infinitesimal displacements, is shown in Figure 24.4(c). Here
B and C move horizontally by v B and v D , respectively, whereas A does not move at all. The horizontal
reactions R A and R D shown in Figure 24.4(c) are obtained in terms of FB and FD by statics on taking
moments with respect to D and A, respectively, in the linearized deformed configuration.
24–6
§24.3 EXAMPLE 2: A PINNED-PINNED THREE-STRUT COLUMN
vB
P P
A FB = k vB
B'
vB
RA = (1/3) (2FB +FC ) vC
= (1/3) k (2vB +vC )
FB = k vB FC = k vC
B' C'
vC
RD = (1/3) (FB +2FC )
= (1/3) k (vB +2vC )
FC = k vC D
C'
P P
Figure 24.5. Free Body Diagrams for the 3-strut column of Figure 24.4.
Pcr1 = kL/9 Pcr2 = kL/3

A A
Buckling mode #1 Buckling mode #2
(antisymmetric) (symmetric)
B 1 B 1
−1 C C 1
D D
Figure 24.6. Critical loads and buckling mode shapes for the 3-strut column of Figure 24.4.
The two FBDs used to form the stability equations are shown in Figure 24.5. We take moments with
respect to C’ in the left FBD and with respect to B’ in the right FBD:
L 2L L 2L
MC = −P vC − FB + R A = 0, M B = −P v B − FC + R D = 0. (24.11)
3 3 3 3
§24.3.2. Critical Loads

Substituting R A = (2FB + FC )/3, R D = (FB + 2FC )/3, FB = k v B and FC = k vC , and casting
(24.11) in matrix form we obtain the stability equations

kL 2 1 vB vB
=P , (24.12)
9 1 2 vC vc
This is a 2 × 2 matrix eigenproblem with P as the eigenvariable. The characteristic equation is

kL 2 1 1 0 k2 L 2 4k P L
C(P) = det −P = − + P 2. (24.13)
9 1 2 0 1 27 9
24–7
δAC
(a)
P (b) P ϕ
A vA vA (+)
A' (c)
FAC FAD vA (+)
ϕ k k ϕ ϕ
δAD
rigid
rigid
L L
90o −ϕ 90o −ϕ
C ϕ ϕ D ϕ ϕ D
B
C
B
Figure 24.7. Rigid strut propped by two inclined springs.
This is a quadratic polynomial in P. The two roots of C(P) = 0 give the critical loads:
kL kL
Pcr 1 = , Pcr 2 = . (24.14)
9 3
An alternative eigenvalue problem is obtained on premultiplying by the inverse of the LHS matrix in
(24.12). This gives

3P 2 −1 vB vB
= , (24.15)
kL −1 2 vc vC
an eigenproblem which gives the same critical loads (24.14). This is the form derived by Jennings
(except for scale factors)
§24.3.3. Buckling Mode Shapes
The eigenvectors corresponding to the critical loads (24.14), normalized to 1 as largest entry, are

−1 1
v1 = , v2 = . (24.16)
1 1
These are plotted in Figure 24.6. Note that the mode shape corresponding to the lowest critical load,
Pcr 1 = k L/9, is antisymmetric. This is in sharp contrast to the Euler column buckling treated in
Lecture 25, which buckles symmetrically.
§24.4. Example 3: Strut Propped by Inclined Springs
The rigid column of Figure 24.7(a) is propped by two symmetrically located extensional springs AC
and AD of stiffness k. The spring rise angles with respect to ground are ϕ, as indicated. Both springs
can take tension and compression. The structure moves only in the plane of the figure. Take the lateral
displacement v A (positive to the right: →) as only degree of freedom (DOF). Only the critical load is
of interest; thus assume from the start that displacements and rotations are very small.
24–8
§24.4 EXAMPLE 3: STRUT PROPPED BY INCLINED SPRINGS

Displace A by v A , positive to the right. For infinitesimal displacements, spring AC elongates by
v A cos ϕ whereas spring AD shortens by the same amount. The resulting spring forces are FAC =
k v A cos ϕ along AC (tension if v A > 0) and FAD = FAD along AD (compression if v A > 0), acting as
pictured in Figure 24.7(b,c). Their resultant is a force of magnitude FA = FAC cos ϕ + FAD cos ϕ =
2 k v A cos2 ϕ, acting horizontally at A in opposition to v A .
Taking moments with respect to B (positive CW) yields the characteristic equation

M B = 0 ⇒ Pv A − FA L = Pv A − 2k L cos2 ϕ v A = (P − 2k L cos2 ϕ) v A = 0. (24.17)
There is no need to put this in matrix form, since the systems has only one DOF.
§24.4.2. Critical Load

For buckling v A = 0, so the factor of v A in the characteristic equation (24.17) must vanish. This gives
Pcr = 2k L cos2 ϕ. (24.18)
As a function of ϕ, Pcr is maximized for ϕ = 0◦ , as may be expected. But making ϕ → 0◦ has two
disadvantages:
1. If springs must be attached to ground level (aligned with the strut base), their length goes to ∞.
2. The spring stiffness k generally drops as the spring length increases.
So in practice there is an nonzero optimal rise angle, in terms of getting the highest critical load under
constraints on spring fabrication. One case, which gives a best rise of about 35◦ , is studied below.
§24.4.3. Optimal Rise Angle

Suppose that the propping springs are actually taut cables of elastic modulus E and cross section area
A. The equivalent spring stiffness of each is k ≈ E A/L s , in which L s = L/ sin ϕ is the length of each
cable. If the rise angle ϕ decreases, L s increases and k goes down, counteracting the beneficial effect
of increasing cos2 ϕ. Replacing the equivalent k into (24.18) we get
Pcr = 2E A cos2 ϕ sin ϕ, (24.19)
from which the column height L has disappeared. Maxima and minima of this Pcr with respect to ϕ
occur at d Pcr /dϕ = E A cos ϕ(cos2 ϕ − 2 sin2 ϕ) = E A cos ϕ 3 cos (2ϕ) √ − 1 = 0. The smallest
positive root of this equation , which satisfies cos(2ϕ) = 1/3 or cos ϕ = 2/3, gives the optimal rise
angle

2 4
ϕbest = arccos ≈ 35.26◦ ⇒ Pcr,max = √ E A ≈ 0.7698 E A. (24.20)
3 3 3
24–9
25
Stability Of
Structures:
Continuous Models
25–1
Lecture 25: STABILITY OF STRUCTURES: CONTINUOUS MODELS
TABLE OF CONTENTS
Page
§25.2. Example 1: The Euler Column 25–3
§25.2.1. Critical Load Analysis . . . . . . . . . . . . . . 25–4
§25.2.2. Linearization Limitations . . . . . . . . . . . . . 25–5
§25.2.3. Reformulation as Eigenproblem . . . . . . . . . . . 25–6
§25.3. Example 2: The Fixed-Pinned Column 25–6
§25.3.2. Reformulation as Eigenproblem . . . . . . . . . . . 25–8
§25.4. Example 3: Elastically Restrained Column 25–9
§25.4.1. Problem Description . . . . . . . . . . . . . . 25–9
§25.4.3. Equivalent Spring Constant . . . . . . . . . . . . 25–11
25–2
§25.2 EXAMPLE 1: THE EULER COLUMN
This Lecture focuses on buckling analysis of continuous models of elastic systems. The analysis
process leads to ordinary or partial differential equations. Those can be solved in closed form only
for simple 1D structures, such as prismatic columns treated by linearized prebuckling (LPB).
Despite limitations in terms of obtaining analytical solutions, such problems are important for
displaying key intrinsic features of continuous models. For example, the associated eigenproblem
posseses an infinite (but denumerable) set of eigenvalues, meaning that there is an infinite number
of critical loads. (This should not be surprising, since continuous models have an infinite number
of degrees of freedom.) Of these, the lowest buckling load is of primary interest to designers.
We illustrate the treatment of continuous models via several examples that involve prismatic elastic
columns with simple boundary conditions. As a result, all of the examples in this Lecture have
closed form solutions.
§25.2. Example 1: The Euler Column
(a) P (b) (c)

P P
A A
A y y y
x x x
X'
X X' X
X
v(x) v(x)
L P M(x) = −P v
elastic
constant EI assumed (+ as drawn)

buckled shape
B B
P
Figure 25.1. The Euler column: (a) problem definition; (b) FBD of whole column assuming a
buckled shape v(x) with its amplitude exaggerated for visibility; (c) FBD at distance x from top.
The Euler column is shown in Figure 25.1(a). It is a homogeneous, prismatric, elastic column
pinned (hinged, simply supported) at both ends A and B, and subjected to axial load P. The elastic
modulus is E. Reference axes are chosen as follows: x is longitudinal, with origin at end A; z is
normal to the plane of the figure and selected so that I = Izz is the minimum second order moment
of inertia of the column cross section about z; finally, y lies in the plane of the figure. For example,
if the cross section is a solid rectangle of width b and thickness t < b, the minimum I = Izz is
b t 3 /12, and y, z will be aligned with the short and long cross-section dimensions, respectively.
25–3
It is emphasized that the minimum moment of inertia, called here I = Izz as noted above, is the one
that determines the buckling load.1 With the axes chosen as shown, the column will buckle in the
x y plane. Figure 25.1(b) pictures a buckled shape that satisfies the end boundary conditions and is
C 1 continuous. Such a curve is called a kinematically admissible buckling mode shape.
The mode shape illustrated by Figure 25.1(b) is mathematically defined by the deflection curve
v(x), whose determination is part of the problem.2 This deflection is infinitesimal because of the
linearized prebuckling (LPB) assumption. In that Figure (and others that follow) it is drawn with
exaggerated amplification for visibility convenience.
§25.2.1. Critical Load Analysis

Figure 25.1(b) also shows reactions at B, which here reduce to just P. There are no end moment
reactions because the column is hinged at both A and B. Furthermore, taking moments with respect
to A and B one shows that both lateral reactions (along y) vanish.
Next, do a FBD of a segment AX, with X located at a distance x from A, as shown in Figure 25.1(c).
The displaced X is labeled X’, and the distance from X to X’ is v(x). At this cross section we
will have a bending moment Mz (x), which is positive as drawn in the figure (reason: a +Mz (x)
compresses the “beam top surface”, which by convention lies on the +y side.) Taking moments
with respect to X’ , equilibrium requires Mz (x) + P v(x) = 0. But according to beam deflection
theory, Mz (x) = E Izz v (x) = E I v (x). Replacing gives
E I v (x) + P v(x) = 0. (25.1)
This is a homogeneous, second order, linear ODE in the unknown deflection v(x). For convenience
in reducing (25.1) to standard (canonical) form introduce

def P
λ = + , (25.2)
EI
Dividing (25.1) through by E I and substituting (25.2) we get the canonical form
v (x) + λ2 v(x) = 0, (25.3)
Its general solution is

v(x) = A cos λx + B sin λx, (25.4)
in which A and B are determined from the kinematic boundary conditions v(0) = 0 and v(L) = 0.
The first one requires that A = 0, whence (25.4) reduces to
v(x) = B sin λx. (25.5)
1 If the section happen to have equal bending inertia in all directions, as in the case of a circular or square section, buckling
can occur in any direction. The actual buckling plane is determined by imperfections, a topic not treated in this course.
2 The notation choice is not accidental, for v(x) can be interpreted as a lateral beam deflection. That symbol was used in
previous lectures that dealt with beam deflections.
25–4
§25.2 EXAMPLE 1: THE EULER COLUMN
π 2E I
Pcr2 = 4 π 2E I
2
Pcr3 = 9 π 2E I
2
Pcr = Pcr1 =
L2 L L
L/3
L/2
L L/3
L/2
L/3
Mode shape 1 Mode shape 2 Mode shape 3

(n=1) (n=2) (n=3)
Figure 25.2. First three buckling mode shapes for the Euler column.
This is called a characteristic stability equation or simply characteristic equation; the names
buckling equation and critical load equation are also used. It has two possible solution types:

B = 0 ⇒ v(x) = 0: the column remains straight
(25.6)
B = 0 ⇒ v(x) = 0: the column buckles with shape sin λx
These two solution types are designated as trivial and nontrivial, respectively, in the literature.
The critical loads are the values of P at which nontrivial solutions are possible. These are determined
by applying the second end condition: v(L) = 0. Since B = 0 we must have sin λL = 0.
Its solutions are λL = nπ, in which n = 1, 2, . . . is a positive integer. Square both sides for
convenience: λ2 L 2 = n 2 π 2 , replace λ2 by P/(E I ), solve for P and tag those loads as critical:
n2 π 2 E I
Pcr,n = , n = 1, 2, . . . (25.7)
L2
As can be observed, there is an infinite number of critical loads. These are the nontrivial solutions
of the characteristic equation (25.5). The lowest one is associated with n = 1:
π2E I
Pcr,1 = Pcr = . (25.8)
L2
This is called the Euler critical load.

The buckling mode shapes associated with the set of critical loads (25.7) are
n πx
vcr,n (x) = B sin . (25.9)
L
If n = 1 the mode shape is a half sine wave, similar to that drawn in Figure 25.1(c). If n = 2 we get
a complete sine wave, if n = 3, a one-and-a-half sine wave, etc. Those buckling shapes are drawn
for n = 1, 2, 3 in Figure 25.2.
25–5
§25.2.2. Linearization Limitations

The LPB assumption does yields the critical loads, but also brings about the following inconsisten-
cies, some of which have an air of paradox.
• The foregoing solution says nothing about the amplitude of the buckling modes (25.8) in terms
of the data. This is a consequence of the LPB assumptions, which led to a simple linear ODE
(25.3), but filters out that information. Determination of this postbuckling behavior is quite
involved, since it requires solving a nonlinear ODE in terms of elliptic functions. Such analysis
shows that the buckled column continues to take up on increasing load, albeit very slowly, as
long as it remains elastic.

• Differentiating (25.5) three times gives a nonzero transverse shear load Vy (x) = E I B vcr 1 =
−E I (π /L )B cos(π x/L). Evaluation at x = 0 or x = L gives a nonzero transverse shear
3 3
force whereas Figure 25.1(b) shows zero horizontal reactions there. A similar inconsistency
arises if one differentiates four times. This gives a nonzero lateral load p(x) but no such load
exists. This “paradox” indicates that the sine-wave result must be corrected once the buckling
mode amplitude becomes finite.
• The analysis does not take into account imperfections such an initially crooked column, or
axial load eccentricity. That topic is beyond the level of this course, but a basic result is used
in the justification of the “Southwell plot” experimental procedure described in Lecture 26.
§25.2.3. Reformulation as Eigenproblem

Lecture 24 emphasized that a discrete stability problem handled through the LPB assumptions may
be presented as a matrix eigenproblem. This is also the case for the linearized continuous problem,
but it requires reformulation. For the Euler column, the two end boundary conditions v(0) = 0 and
v(L) = 0 may be written conjointly as

1 0 A 0
= . (25.10)
cos λL sin λL B 0
√
This is a matrix eigenproblem in λ = P/(E I ), but note that (unlike the discrete case) this
unknown does not appear linearly. For a nontrivial solution, meaning that A and B are not both
zero, the determinant of the matrix on the left side must vanish. This leads to the condition
sin λL = 0, (25.11)
which agrees with that previously found from the characteristic equation (25.5). Mathematicians
call this a trascendental eigenproblem since the characteristic equation (25.11) is a trascendental
(non-polynomial) equation. Observe that this reformulation does not gain much here. But it
becomes useful in more complicated problems, such as the next two examples.
§25.3. Example 2: The Fixed-Pinned Column

As second example we consider the configuration shown in Figure 25.3(a). The buckled column
is drawn in Figure 25.3(b) in a deflected position. The most notable difference with respect tp the
Euler column of Figure 25.1(a) is the presence of additional reaction forces: the lateral reactions
R A and R B , and the fixed end moment M B . These are positive as shown in Figure 25.3(b).
25–6
§25.3 EXAMPLE 2: THE FIXED-PINNED COLUMN
(a) P (b) (c)

P P
RA A RA A
A y y y
x x x
X'
X X' X RA
X
v(x) v(x)
L P M(x)
elastic assumed (+ as drawn)
constant EI
buckled shape
MB
B
B RB = R A
P
Figure 25.3. Column fixed at one end and simply supported (hinged, pinned) at the other.

Equilibrium of y forces in Figure 25.3(b) gives R A = R B . Equilibrium of moments taken with
respect to either A or B yields M B = R A /L. Next consider the FBD of the portion AX shown in
Figure 25.3(c). Equilibrium with respect to X’ gives the non-homogeneous ODE
MB x
E I v (x) + Pv(x) = R A x = . (25.12)
L
The LHS is the same as in (25.1), but the RHS is no longer zero. Dividing through by E I and
def
setting λ2 = P/(E I ) gives
MB x λ2 M B x
v (x) + λ v(x) =
2
= . (25.13)
EI L PL
The general solution is the sum of the homogeneous solution v H (x) = A cos λx + B sin λx and
the particular solution v P (x) = λ2 M B x/(λ2 P L) = M B x/(P L):
MB x
v(x) = A cos λx + B sin λx + . (25.14)
PL
The three kinematic BCs are v A = v(0) = 0, v B (0) = v(L) = 0 and v B = v (L) = 0. These
provide three equations:
MB MB
A = 0, A cos λL + B sin λL + = 0, −A λ sin λL + B λ cos λL − = 0. (25.15)
P PL
Solving these equations simultaneously one obtains the characteristic equation
tan λL = λL . (25.16)
25–7
ClearAll[x,P,A,B,MB,EI,λ,λλ ]; c={A,B,MB};
vG=A*Cos[λ*x]+B*Sin[λ*x]+MB*x/(λλ*L);
ODE=D[vG,{x,2}]+λλ*vG-MB*x/(L);
chk=Simplify[ODE/.{λλ->(P/EI),λ->Sqrt[P/EI]}];
Print["check vG=",chk]; vA=vG/.x->0;
vB=Simplify[vG/.x->L]; vpB=Simplify[D[vG,x]/.x->L];
Print["vA=",vA," vB=",vB," v'B=",vpB];
S={Coefficient[vA,c],Coefficient[vB,c],Coefficient[vpB,c]};
S=Simplify[S/.λλ->P/EI];
Print["Stab matrix=",S//MatrixForm];
Sdet=Simplify[Det[S]]; Print["Sdet=",Sdet];
Figure 25.4. Mathematica script to produce the characteristic equation for Example 2.
The smallest root of this transcendental equation, to 4 places, is λL = 4.493. The corresponding
critical load is
20.19 E I
Pcr = . (25.17)
L2
Since 20.19 ≈ 2.05 π 2 this critical load is approximately twice that of the Euler column. Therefore,
fixing one end has substantially increased the critical load.
§25.3.2. Reformulation as Eigenproblem

The three kinematic BC: v A = v(0) = 0, v B = v(L) = 0, and v B = v (L) = 0, written out in
(25.15), can be combined in matrix form as
 
1 0 0
A 0
 cos λL λL 1 
 sin P  B = 0 , (25.18)
−λ sin λL λ cos λL P1L MB 0
The matrix that appears in the LHS is the stability matrix, say S. Its determinant is
1
det S = sin λL − λL cos λL . (25.19)
PL
Since P and L are not zero, this determinant vanishes only if sin λL − λL cos λL = 0, or
tan λL = λL. This is the same characteristic equation (25.16) found previously. The proce-
dure can be conveniently carried out symbolically using a computer algebra system (CAS). For
example, Figure 25.4 shows a Mathematica script that produces (25.19).
25–8
§25.4 EXAMPLE 3: ELASTICALLY RESTRAINED COLUMN
(a) (b) (c) (d) P

P P P
RA RA
A y y,v A
x x M(x)= −P v+RA x
X
RA beam-column
L P
v(x)
k
MB = − k θB C elastic beam
B
B
RA
x P
−θB
Figure 25.5. Column with elastic restraint: a torsional spring, at end B.
§25.4. Example 3: Elastically Restrained Column
The title configuration is of interest because it covers the three boundary condition cases tested in
the Column Buckling Lab.3 as special cases.
§25.4.1. Problem Description
Consider the modification of the classical Euler column pictured in Figure 25.5(a). The column
is simply supported (pinned) at both ends A and B, and axially loaded by P. The rotation at B is
further restrained by a torsional spring with stiffness k.
Column AB has length L and constant flexural rigidity E I , in which I ≡ Izz is the minimum
moment of inertia of the cross section that controls buckling. (For a rectangular cross section of
width b and thickness t < b, I = Izz = b t 3 /12.) The column is simply supported (pinned) at
both ends A and B, and axially loaded by P. The rotation at B is further restrained by a torsional
spring with stiffness k. Select x as shown and consider the buckling shape v(x) sketched in Figure
A.1(b). The end rotation at B is θ B = v (L), positive CCW. The spring at B applies a restoring end
moment M B = −k θ B . For convenience in obtaining dimensionless equations we define
EI
k=β (25.20)
L
where β is a numerical coefficient. If β = 0 the problem reduces to that of the classical Euler
column, whereas if β → ∞ we obtain the case of a column simply supported at A and fixed
(clamped) at B.
3 Prior to Fall 2011, that used to be a formal Experimental Lab. Since then it has become a “Lab-Homework” or “Labwork,”
meaning that test results are used as part of an assigned Homework exercise rather than a Lab Report.
25–9

The FBD of the complete column is shown in Figure 25.5(b). Taking moments with respect to B
one finds that
R A = M B /L = −k θ B /L = −β E I θ B /L 2 . (25.21)
Now cut the column at section X at distance x as sketched in Figure 25.5(c). Moment equilibrium
with respect to X yields the second order differential equation
x E I θB x
E I v + P v = R A x = −kθ B = −β , (25.22)
L L2
which results from equating the bending moment M(x) = E I v to −Pv + R A x. Dividing through
by E I and calling λ2 = P/E I produces the canonical form
θB x
v + λ2 v = −β . (25.23)
L2
The general solution of (25.23) is the sum of the homogeneous and particular solutions:
θB x
v(x) = C1 sin λx + C2 cos λx − β . (25.24)
λ2 L 2
Since v(0) = 0, C2 = 0. The rotation is θ(x) = v (x) = C1 λ cos λx − βθ B /(λ2 L 2 ). Evaluating
this at x = L yields θ B = θ(L) = C1 λ cos λL − βθ B /(λ2 L 2 ), from which we can solve for the end
rotation:
λ3 L 2
θ B = C1 cos λL . (25.25)
β + λ2 L 2
Inserting this into the solution (25.24) with C2 = 0 gives v(x) = C1 sin λx − C1 βλx cos λL/(β +
λ2 L 2 ). Applying now the second boundary condition: v(L) = 0, furnishes the stability equation

βλL
v(L) = C1 sin λL − cos λL = 0. (25.26)
β + λ2 L 2
For buckling to occur, C1 = 0, and the expression in parentheses in (25.26) must vanish. Calling
α = λL, which is also a dimensionless variable, we obtain the trascendental equation
αβ
tan α = (25.27)
α2+β
For a given β ≥ 0 we seek the solution αcr > 0 of (25.27) closest to zero. If β = 0 there is no
closed form solution and it is better to proceed numerically, using for example a Newton solver. It
is easily shown that for β = [0, ∞], π ≤ αcr < 4.5, so α ≈ 4 is a good start point for a Newton
iteration. Such a solver is implemented in Mathematica in the built-in function FindRoot. For
example, the statements
beta=100; Print["alphacr=",alpha/.FindRoot[Tan[alpha]==
alpha*beta/(alpha^2+beta), {alpha,4}]];
25–10
§25.4 EXAMPLE 3: ELASTICALLY RESTRAINED COLUMN
return αcr = 4.4494 as the desired numerical solution of (25.27) for β = 100. Here is a table for
selected values of β:
β 0 1 3 10 100 1000 10000 ∞

αcr π 3.4056 3.7264 4.1323 4.4494 4.4889 4.4930 4.4934
L e /L 1.0000 0.9224 0.8431 0.7602 0.7061 0.7000 0.6992 0.6991
The critical load Pcr and the effective length L e (tabulated above) are
αcr
2
EI π2E I π
Pcr = = , Le = L (25.28)
L2 L 2e αcr
This solution also applies to the buckling of a column AB rigidly connected at B to an elastic
beam BC, as illustrated in Figure 1(d). This is actually one of the configurations to be tested in
Lab 4. All that is needed is to work out the appropriate value of k = M BBC /θ BBC obtained from a
moment-deflection analysis of beam BC subjected to an applied end moment M BBC = M B .
§25.4.3. Equivalent Spring Constant

Suppose the restraining elastic beam BC in Figure A.1 has length L r = L BC , modulus Er and
moment of inertia Ir about z. Assume a simply support condition at C. Under an end moment
M B applied as pictured in Figure 25.5(d), the restraining beam deflects by4 vr = −M B xr (L r2 −
xr2 )/(6Er Ir L r ), in which xr is the distance from C. The end rotation is θ B = (dvr /d x)|xr →L r =
M B L r /(3Er Ir ), whence the equivalent torsional spring stiffness is k = M B /θ B is 3Er Ir /L r .
Comparing with (1) shows that β = 3(L/L r )(Er Ir )/(E I ). If the restraining beam is fabricated of
the same material and has the same cross section dimensions of the beam-column, as is the case in
Lab 4, Er = E, Ir = I , and β = 3L/L r .
4 E. M. Popov, Engineering Mechanics of Solids, Prentice Hall, NJ, 1999, case 5 of Table 10, p. 853.
25–11
26
Stability Of
Structures:
Additional Topics
26–1
Lecture 26: STABILITY OF STRUCTURES: ADDITIONAL TOPICS
TABLE OF CONTENTS
Page
§26.2. Unified Buckling Formula 26–3
§26.2.1. Effective Length . . . . . . . . . . . . . . . . 26–3
§26.2.2. Critical Load in Terms of Slenderness Ratio . . . . . . . 26–3
§26.3. Failure Mode: Buckling Versus Yield 26–4
§26.3.1. Failure Envelopes . . . . . . . . . . . . . . . . 26–4
§26.3.2. Long Versus Short Columns . . . . . . . . . . . . 26–4
§26.4. The Southwell Plot 26–6
§26.4.1. Effect of Imperfections . . . . . . . . . . . . . . 26–6
§26.4.2. Southwell Plot Example . . . . . . . . . . . . . 26–7
§26.5. The ITL Buckling Demo 26–9
§26.5.1. Module Description . . . . . . . . . . . . . . . 26–10
§26.5.2. Experimental Procedure . . . . . . . . . . . . . 26–10
§26.5.3. Specimen Dimensions . . . . . . . . . . . . . . 26–12
§26.5.4. Restraint Beam Attachment Details . . . . . . . . . 26–12
§26.5.5. Miscellaneous Reminders . . . . . . . . . . . . . 26–13
§26.5.6. Comparing with Analytical Buckling Loads . . . . . . . 26–14
26–2
§26.2 UNIFIED BUCKLING FORMULA
This Lecture covers additional topics that are useful for Homework 10 and the Final Exam, as well
as the “Labwork” demo on Friday November 30. The data collected from the lab is to be used for
Homework Exercise 10.6.
§26.2. Unified Buckling Formula
The bucking load formulas for prismatic elastic columns depends on the end conditions but contains
many common ingredients. For design purposes it is convenient to have a unified formula that makes
it look like the buckling load of an “effective” Euler column.
§26.2.1. Effective Length

The effect of different end conditions on the critical load of a column can be unified by expressing
it in the form
π2E I
Pcr = 2 (26.1)
Lef f
Here L e f f denotes the effective length of the column. This length turns out to have a simple physical
interpretation: distance between the inflexion points of the buckling curve associated with Pcr .
P
P P P
L
L eff =0.7 L
L eff = L L eff =2L L L L eff =L/2
fictitious
continuation
about fixed end
pinned-pinned free-fixed pinned-fixed fixed-fixed

(Euler column) (cantilever)
Figure 26.1. Effective lengths of columns with different end conditions.
Values of L e f f for three common support conditions are given in Figure 26.1. (The value of 0.7 L
listed for the pinned-fixed case is correct to 3 places; a more accurate value is 0.6992 L but 0.7 is
easier to remember.) The unified definition (26.1) allows us to extend results derived for the Euler
column to other cases, simply by replacing the appropriate effective length in the formula.
26–3
§26.2.2. Critical Load in Terms of Slenderness Ratio

In formula (26.1), I (the √
minimum second moment of inertia of the cross section) may be replaced
by A r 2 , in which r = + I /A is the radius of gyration of the column cross section. Examples of
r for three column cross sections:
• is a b × h rectangle with h ≤ b, I = b h 3 /12, A = b h, then r 2 = I /A = h 2 /12
If the section√
and r = 12 h/ 3 ≈ 0.2887 h.
• If a solid circle of radius R, I = π R 4 /4, A = π R 2 , then r 2 = I /A = R 2 /4 and r = 12 R.
• mean radius R and thickness t << R, I ≈ π R 3 t, A = 2π Rt,
If a thin-wall circular tube of √
r 2 = I /A = R / 2 and r = R/ 2 = 0.707 R.
Substituting I = A r 2 in (26.1) gives
π2 E I π2 E A r2 π2 E A
Pcr = = = , (26.2)
L 2e f f L 2e f f s2
in which s = L e f f /r is called the slenderness ratio of the column.
§26.3. Failure Mode: Buckling Versus Yield

A key question in column design is: will be column fail first by yield or elastic buckling? In other
words, which is the failure mode?
The average axial stress at the critical load Pcr expressed as (26.2) is
Pcr π2 E
σcr = = 2 , (26.3)
A s
This is called the critical stress. A key design question is: how does the critical stress compare with
the yield stress σY of the column material? Evidently if σcr > σY the critical load formula is not
valid, since the derivation is based on the assumption that the column is linearly elastic at buckling.
To visualize the limitation of the buckling formula, we will use a graphical interpretation of (26.3).
§26.3.1. Failure Envelopes

Introduce two dimensionless ratios for the column material:
def σcr def E
σ̄cr = , Ē = , (26.4)
σY σY
Dividing both sides of (26.3) by σY and introducing the ratios (26.4) gives a dimensionless version
π 2 Ē
σ̄cr = , (26.5)
s2
This is graphed in Figure 26.2 for the following materials: structural steel (E = 210 GPa, σY = 210
MPa, Ē = E/σY = 1000), aluminum alloy (E = 70 GPa, σY = 280 MPa, Ē = E/σY = 250), and
fir wood (E = 12.6 GPa, σY = 35 MPa in compression, Ē = E/σ y = 360). The curves delimit
what are called failure envelopes.
26–4
§26.3 FAILURE MODE: BUCKLING VERSUS YIELD
1.2 Failure by yield

(short columns)
1.0 Failure by buckling
(long columns)
σcr = σ cr /σ Y
0.8
0.6
Fir Wood Steel
0.4 Aluminum
0.2
0.0
0 50 100 150 200
Slenderness ratio s = Leff /r
Figure 26.2. Column failure modes: buckling versus yield plotted in terms of the
slenderness ratio. Solid lines (“Euler hyperbolas”) represent failure by buckling; dashed
line failure by yield; black dots mark the critical slenderness ratios.
§26.3.2. Long Versus Short Columns

The following three failure cases may be distinguished according to the slenderness of the column.
• Long column. If the slenderness exceeds a critical value scr defined below, the column will
fail by elastic buckling, in which case σcr < σY or σ̄cr < 1. Failure occurs at one of the solid
curves, which are known as Euler hyperbolas in the literature.
• Short column. If the slenderness is less than than the critical value scr defined below, it either
fails by yield at σ = P/A = σY , or σ̄cr = 1 and the failure occurs at the dashed line.
• Goldilocks column. Not too short, not too long. If s = scr , where scr 2
= π 2 Ē = π 2 E/σY
the column will fail simultaneously by buckling and yield. That transition is marked by the
black dot in Figure 26.2. It represents the most efficient use of the material, so it is an optimal
design in that particular sense.
Remark 26.1. Most Mechanics of Materials books show failure diagrams by plotting (26.3) directly, that is,
without dividing through bt σY . As a result the horizontal scale is dimensionless (the slenderness ratio) but
the vertical scale (in stress) is not. The resulting graphs depend on column material as well as physical units
chosen.
Example 26.1. A pinned-pinned steel column with E = 210 GPa and σY = 210 MPa has a pin-to-pin length
of L = 5 m = 5000 mm and a b × h solid rectangular cross section with b = 0.12 m = 120 mm, and h = 0.08
m = 80 mm. Will the column fail first by yield or elastic buckling?
Solution. The critical Euler buckling load is Pcr = π 2 E I /L 2 since L e f f = L for the pinned-pinned case.
The minimum second moment of inertia is I = bh 3 /12 because h < b. Replace and divide by A = bh to
get σcr = Pcr /A = π 2 E h 2 /(12L 2 ) = 44.8 N/mm2 = 44.8 MPa. Compare to yield: σcr < σY = 210 MPa.
Thus the column will fail first by buckling .
√ s = L e f f /r , where L e f f = L and r = I /A = h /12. A

2 2
Alternatively one can check the√slenderness ratio:
quick computation gives s = L 12/ h = 5000 12/80 ≈ 216, which is way into the “long column” range as
can be quickly checked from Figure 26.2.
26–5
Example 26.2. A fixed-fixed steel column (same E and σY as in the previous example) of length L = 6 m has
a solid circular cross section of unknown radius R. Find (1) R in mm so the column fails simultaneosuly by
yield and by elastic buckling, and (2) which maximum load Pmax the so-designed column can support if the
safety factor against both buckling and yield is s F = 4.
Solution. For (1), equate σcr = σY and solve for R. Details: L e f f = L/2, I = (π/4)R √
4
, A = π R2,
r 2 = √I /A = R 2 /4. σcr = π 2 E (R 2 /4)/(L 2 /4) = π 2 E R 2 /L 2 = σY , whence R = σY /E L/π =
L/(π Ē) = L/99.34 = 6000/99.34, hence R = 60.4 mm . The area of this optimal design is
A = π R2 = 11461 mm . The failure load is Pcr = σY A = 2.407 × 1012 N. Divide this by 4 to get
2
Pmax = 6.02 × 1011 N .
§26.4. The Southwell Plot
The Southwell plot1 provides a clever graphical method for nondestructive critical-load testing
of columns (as well as other structural components that may fail by buckling). This method is
particularly useful for field tests of structures or structural components that are likely to be damaged
by being taken up to and beyond the critical load, such as reinforced concrete columns or advanced
composite materials.
§26.4.1. Effect of Imperfections

For the method to work well, the column should be subjected to eccentric axial loads, or to tiny
lateral loads that “trigger” measurable lateral deflections prior to buckling. Let P be the applied
axial load P < Pcr . Let vm (P) denote a measured lateral deflection at load level P. Southwell
observed that the following approximation holds as P approaches Pcr :
a
vm ≈ (26.6)
Pcr − 1
P
where a is a constant. Equation (26.6) may be transformed to
vm
vm = Pcr −a (26.7)
P
If vm /P and vm are plotted along the x and y axis, respectively, (26.7) becomes the equation of a
straight line: y = mx − a, in which the slope m = Pcr .
This observation furnishes a simple but effective experimental method: from measurements at
several axial loads P < Pcr one obtains vm and vm /P as data points. These are plotted along the
1 R. V. Southwell, On the Analysis of Experimental Observations in Problems of Elastic Stability, Proc. Roy. Soc. London,
Series A, 135, pp. 601–616, April 1932. Bio note: Sir Richard Southwell was honored by a nobiliary title in 1948
because of his contributions to the British WWII effort. He was a developer of “relaxation” methods, which were used
to solve systems of hundreds of equations prior to computers. (Before digital computers appeared ca. 1951, the term
“computers” meant humans working on numerical calculations. To solve large systems of equations by relaxation, dozens
of “computers” — mostly women during WWII — were gathered in a large room; each did part of the computations
on desk calculators, receiving and passing results to neighbors.) One important application was the prediction of the
Normandy weather for D-day: June 6, 1944. It worked.
26–6
§26.4 THE SOUTHWELL PLOT
vm
experimental data points
~ Pcr
slope =
vm
P
Figure 26.3. The Southwell plot.
vertical and horizontal axes, respectively, as illustrated in Figure 26.3. This is the Southwell plot.
A straight line is fitted to these points. Its slope estimates Pcr . The value of the constant a is of
little importance.
This technique has the important advantage of being non-destructive because P need not exceed
the critical load Pcr . For this reason it is often used in aerospace structures fabricated of expensive
materials such as composites.
If (26.6) or (26.7) were to hold exactly, the points would fall neatly on a straight line, as Southwell
discovered. This will not be generally the case in experiments because of various factors, which
are discussed in §26.5.
§26.4.2. Southwell Plot Example

This example illustrates an actual Southwell plot produced with data collected in the Fall 2010
semester in the lab demo described in §26.5. The example pertains to the pinned-pinned end
support case of that lab. Four load eccentricities of 1.5 through 6.0 mm were used. Tray-load
versus lateral deflection results are collected in Table 26.1.
The offsets listed in Table 26.1 were not part of the collected data. They were selected afterwards
by trial and error so the lower left portion of the 4 Southwell plots look reasonable. If an offset is
too low or too high, that part of the plot may look like an “outlier” region. (Those portions of the
plots are associated with lower applied loads; as a result vm /P is more sensitive to deflection data.)
The Mathematica script used to produce the Southwell plots is listed in Figure 26.4. The plots are
shown in Figure 26.5 (axis labels were added by Adobe Illustrator). Since the 4 slopes look fairly
consistent, the best-fit red line shown in Figure 26.6 was chosen by “eyeballing” using the pen tool
of Adobe Illustrator. That slope is approximately (see the red dashed lines):
50 mm 50 mm
= = 51.02 N. (26.8)
(1.09 − 0.11) mm/N 0.98 mm/N
Thus the predicted buckling load is
PcrSouthwell = 51.02 N. (26.9)
26–7
Table 26.1 - Experimental Data for Pinned-Pinned Column (Fall 2010)

(Converted from Excel spreasheets to TEX table format)
e = 1.5 mm (1 notch) e = 3 mm (2 notches) e = 4.5 mm (3 notches) e = 6 mm (4 notches)

Offset=1.7 mm† Offset = 4 mm Offset = 4.5 mm Offset = 7.5 mm
TLoad(N)‡ Def(mm) TLoad(N) Def(mm) TLoad(N) Def(mm) TLoad(N) Def(mm)
4 2.1 4 5 3 5 3 8
8 2.8 8 7 6 6 6 9
12 3.5 12 8 9 7 9 10
16 4.2 16 9 12 8 12 11
20 5.6 20 11 15 9.5 15 13
24 7.5 24 15 18 11.5 18 15
28 11.2 28 21.5 21 13.5 21 18
30 14.5 30 25.5 24 16.5 24 21
32 19.8 32 34 27 22.5 27 28
34 30.0 34 48 30 31.5 30 39
36 53.1 35 61 33 48.5 33 59
34 59
† Offsets are chosen by trial and error so lower left portion of the S-plots look reasonable.
‡ TLoad means “tray load”. Actual load on tested columns is (4/3)× tray load.
The analytical critical load is given by Euler’s formula for the pinned-pinned case: Pcr = π 2 E I /L 2 ,
where I is the minimum second moment of inertia of the cross section. The Mathematica script of
Figure 26.7, which uses the geometric and material data for the lab stated in §26.5, computes
PcrEuler = 50.54 N. (26.10)
This corresponds to a critical “tray load” of (3/4) × 50.54 = 37.9 N. As can be seen the agreement
for this BC case is quite good, since the two values differ by only about 1%.
26–8
§26.5 THE ITL BUCKLING DEMO
<<Graphics`MultipleListPlot`;
(* Southwell plots for Pinned-Pinned column - Fall 2010 lab *)
offs1=1.7; offs2=4; offs3=4.5; offs4=7.5;

PPdata1={{4,2.1},{8,2.8},{12,3.5},{16,4.2},{20,5.6},{24,7.5},{28,11.2},
{30,14.5},{32,19.8},{34,30},{36,53.1}};
PPdata2={{4,5},{8,7},{12,8},{16,9},{20,11},{24,15},{28,21.5},
{30,25.5},{32,34},{34,48},{35,51}};
PPdata3={{3,5},{6,6},{9,7},{12,8},{15,9.5},{18,11.5},{21,13.5},
{24,16.5},{27,22.5},{30,31.5},{33,48.5},{34,59}};
PPdata4={{3,8},{6,9},{9,10},{12,11},{15,13},{18,15},{21,18},
{24,21},{27,28},{30,39},{33,59}};
PPSouth1=Table[N[{(PPdata1[[i,2]]-offs1)/((4/3)*PPdata1[[i,1]]),
PPdata1[[i,2]]-offs1}],{i,1,Length[PPdata1]}];
MultipleListPlot[PPSouth1,PPSouth2,PPSouth3,PPSouth4,
PlotJoined->True,Frame->True];
Figure 26.4. Mathematica script to produce Southwell plots for the data of Table 26.1.
50
40
Deflection in mm
30
20
10
0
0 0.2 0.4 0.6 0.8 1 1.2
Deflection over column axial load in mm/N
Figure 26.5. Southwell plots for pinned-pinned column tested with 4 eccentricities.
50
40
Deflection in mm
30
20 "Eyeballed" best-fit
10
0
0 0.2 0.4 0.6 0.8 1 1.2
Deflection over column axial load in mm/N
Figure 26.6. Red line is “eyeballed best fit” to plots of Figure 26.5.
26–9
Critical load of pinned-pinned test column (Euler column)
Em=190000*Nw/mm^2; L=600*mm; t=1.67*mm; w=25*mm;

Izz=w*t^3/12; Pcr=N[Pi^2]*Em*Izz/L^2; Ptray=Pcr*3/4;
Print["Pcr=",Pcr," Ptray=",Ptray];
Pcr=50.542 Nw Ptray=37.90 Nw
Figure 26.7. Critical load for pinned-pinned column computed using the Euler buckling formula.
§26.5. The ITL Buckling Demo

This section describes an experimental lab on column buckling that will be carried out on Friday
December 2 during lab hours. Prior to the Fall 2010 offering, this used to be Lab 3, a formal
experimental lab scheduled similarly to the first two ones. That it was converted to a ”labwork” or
“homelab”, meaning that experiments are done during the two regular lab sections. The collected
data is collectively incorporated into Homework Exercise 10.6, which deals with the Southwell plot
described in §26.4.
§26.5.1. Module Description

The apparatus being used is an off-the-shelf module produced by a British company by the name of
Hi-Tech Ltd.2 As befits its European provenance, relevant dimensions and weights are provided in
SI units. The Beam-Column Buckling Module is sketched in Figure 26.8. It consists of a mobile
steel frame that supports a load arm and end conditions for the beam-column. A high-strength steel
beam-column is provided with a nominal Young’s Modulus of 200–205 GPa. Pinned and clamped
conditions can be simulated at the top of the supplied beam-column with different adapters. The
boundary condition at the bottom of the beam can be varied from pinned to clamped through the
use of a restraint beam that provides adjustable levels of torsional stiffness at this point.
The end fixtures clamped to the beam-column also provide off-axis notches for applying eccentric
loads. These notches are spaced 1.5 mm apart. Typically, an eccentric load is applied by offsetting
the knife-edge contacts an equal number of notches off-center at the top and bottom of the beam-
column.
The ratio of the load arm lengths from the pivot point to the beam-column and from the pivot point
to the load tray is 3:4 as shown in Figure 26.8. Remember to convert the loads applied at the load
tray to those resulting at the beam-column appropriately. During loading, some friction between
the load arm and frame may develop. The effect of this friction on the deflection response can be
alleviated by gently tapping on the frame after loads are applied.
§26.5.2. Experimental Procedure

The experiment procedure consists of making load-displacement measurements for several combi-
nations of boundary conditions and eccentric loads. Initially, data will be collected for a nominally
pinned-pinned beam-column using two or three different eccentric load configurations. (The exact
number will be decided upon during the labs.) This data will be used to explore the benefits of
using the Southwell plot.
2 Not a very accurate name, but the apparatus is rugged and serves its purpose.
26–10
3L L
load arm knife edge

counterweight stop
beam-column specimen
of high-strength steel load tray
ruler
beam clamps
frame restraint beam
Figure 26.8. ITLL Beam-Column Buckling Module (aka Hi-Plan Module).
Following this, the BC at the beam-column base will be adjusted to a nominally clamped state and
then to an elastically restrained state.
Do not forget to make careful measurements of the beam-column dimensions.
BC Case 1: Pinned-Pinned Response

For this case the restraining beam should be detached from the beam-column. Assemble the beam-
column with zero eccentricity and with the ruler centered along the length of the beam-column as
demonstrated in the demos. (Use a tiny 2 N weight on the load tray to keep the beam-column in
place during assembly.)
Check that the lower knife-edge is located directly beneath the upper knife-edge by measuring the
distances from these boundaries to the right hand side of the frame. (How would you model errors
in this alignment?) Align the ruler with a reference point on the beam to ”zero” (don’t forget your
2 N load) the measurement. Measure the lateral deflections induced by applied loads varying from
”zero” to the critical load. (You may need to help the beam-column buckle to the right as described
in the demos.) Since the deflections of the beam-column will increase rapidly as the critical load
is approached, reduce the loading increments near this load. It is difficult to resolve deflections to
better than ±0.25 mm, so obtaining the larger deflection data near failure is important. You should
have at least 10 load increments. Once you are satisfied with the repeatability of your data for this
configuration, adjust the load eccentricity to 3 mm (2 notches) and repeat the load-deflection data
acquisition procedure. You may need to shift the ruler or adjust your displacement measurements
accordingly.
Repeat the procedure yet again for 6 mm of load eccentricity. Use a spreadsheet to plot the load-
deflection data and to generate Southwell plots of deflection vs. compliance. Fit lines to the
three data series in the Southwell plots to derive values for the critical loads for all three cases.
26–11
How does the apparent critical load in the load-deflection plots and Southwell analysis compare to
the theoretical prediction? (When measuring the geometry of the beam-column, consider which
dimension(s) are most critical.)
BC Case 2: Fixed-Pinned Response

Return the knife-edges to the on-axis notches (no eccentricity). Attach the restraint beam to the
bottom of the beam-column and position the restraint beam clamps as close as possible to the bottom
of the beam-column. Be sure to tighten the clamps to both the beam and the frame. Attachment
details are shown in Figure 26.9 on next page.
This will provide a nominally fixed boundary condition. In this test, we would again like to measure
the lateral deflection where it is greatest. Consider the buckled deflection profile you would expect
for these boundary conditions and adjust the vertical position of the ruler accordingly. You can
proceed to buckle the beam-column to check your decision. Explain how you selected this ruler
position in your lab report. Once the system is configured, repeat the pinned-pinned procedure for
collecting and analyzing the load-deflection data. Be sure to do the three eccentricities: 0, 3 and 6
mm.
How do the measured values for critical load, derived from Southwell’s plot, compare with the
theory for a fixed-pinned beam-column? How “fixed” does the lower boundary condition appear
to be?
BC Case 3: Elastically Restrained-Pinned Response

Loosen both restraint beam clamps and slide the left-most one completely off of the restraint beam.
Position the remaining clamp 500 mm from the beam column and position the top grips of the
restraining clamp to provide about 0.5 mm of clearance to allow lateral motion of the restraint
beam. Approximate where the maximum beam-column lateral deflection will occur and position
the ruler accordingly. Again, explain in the lab report how this was done. Repeat the data collection
and analysis procedure one last time for this case. Be sure to do the three eccentricities: 0, 3 and 6
mm. Attachment details are shown Figure 26.9 below.
How would you model the BC provided by the restraint beam? (The appropriate theory is provided
in Lecture 18.) When this is done, how do the measured values for critical load, derived from
Southwell’s plot, compare with the theory?
§26.5.3. Specimen Dimensions

From measurements taken on 11/11/03, the beam-column specimen (and restraint beam) in the
buckling module seems to be nominally b = 1” = 25.4 mm wide and t = 1/16” = 1.58 mm thick. A
caliper measurement gives a thickness of about 68/1000 in = 1.7 mm, a bit larger than 1/16”. [This
is curious because all apparatus dimensions are supposed to be metric.] Note that the buckling
load Pcr is very sensitive to the thickness t because Izz = bt 3 /12.) Please recheck these specimen
dimensions.
The span between knife-edges is 60 cm = 600 mm. The elastic modulus E of the high strength
steel used in the specimen and restraint beam has a nominal range of 200 to 205 GPa.
26–12
Tighten clamp screws

to preclude rotations
Place aluminum spacer

below restraint beam
and under clamp grips
Figure 26.9. Attachment detail for fixed-pinned case.
Keep this clamp sufficiently loose

so it acts roughly as a hinge, allowing
rotation but little vertical motion
Slide this clamp
out of the way
distance 50 cm = 500 mm
Figure 26.10. Attachment detail for elastically restrained-pinned case.
§26.5.4. Restraint Beam Attachment Details
Fixed-Pinned Case. Slide the aluminum spacer (in blue plastic box) below restraining beam and
under grip of clamps, as sketched in Figure 26.9. Tighten clamps enough to prevent rotation of the
beam-column specimen, but do not overtighten. During this operation, please check that upper end
of beam-column stays aligned with knife-edge. If misalignment occurs, loosen one side slightly as
necessary to keep clamps even.
Elastically Restrained-Pinned Case. Slide one clamp out of the way. (Also do not continue to
use the aluminum spacer, save it in the blue plastic box.) Place the other clamp 50 cm (500 mm)
away from the beam column as sketched in Figure 26.10. Attach the bottom clamp screw to set the
distance. Keep the upper screws sufficiently loose to allow rotational motion as well as axial slip,
but not so excessively loose as to allow significant vertical motions. This configuration justifies
idealization of the clamp as a hinge BC in the theoretical analysis.
§26.5.5. Miscellaneous Reminders
• If something is wrong with the equipment, or you find that items are missing, do not proceed.
Report problem to instructor or TAs.
• At least one student in each subsection should bring a notebook to record measurements.
Transcribe data later to a spreadsheet. It is recommended to bring also a pocket metric ruler.
• For dimensions of the apparatus, beam column specimen and restraint beam, you may check
Instruction Manual HST/12, “Column Buckling Failure.” This is kept with the module. All
dimensions given there are in mm. But it is equally as easy to measure those yourself.
26–13
• Once finished with all Friday demos, disssamble the apparatus and place all items and tools
in the blue plastic box. Do not forget to insert the two restraint-beam attachment screws back
in the beam specimen to avoid loss.
§26.5.6. Comparing with Analytical Buckling Loads

The elastically restrained case is worked out in §25.4 of Lecture 25. This solution includes the
pinned-pinned and fixed-pinned columns as special cases by making the torsional spring constant
zero and infinite, respectively.
26–14

Plan Stress PDF

Uploaded by

Copyright:

Available Formats

Plan Stress PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Plan Stress PDF

Uploaded by

Copyright:

Available Formats

6

§6.2. Thin Plate in Plate Stress

In-plane internal forces y + sign conventions

dx dy In-plane stresses In-plane body forces

In-plane strains In-plane displacements

§6.3. 2D Stress Transformations

§6.3.1. Why Are Stress Transformations Important?

§6.3.2. Method of Equations

σnn = σx x cos2 θ + σ yy sin2 θ + 2 τx y sin θ cos θ,

Couple of checks are useful to verify these equations. If θ = 0◦ ,

σnn = σx x , σtt = σ yy , τnt = τx y , (6.2)

σx x + σ yy = σnn + σtt . (6.4)

§6.3.3. Double Angle Version

§6.3.4. Principal Stresses, Planes, Directions, Angles

This is satisfied for θ = θ p if

To evaluate (6.8) it is convenient to go through the following staged sequence:

σ1 = σav + R, σ2 = σav − R. (6.10)

§6.3.5. Maximum Shear Stresses

§6.3.6. Principal Stress Element

2θ2 = 36.88 +180 = 216.88

§6.3.7. Mohr’s Circle

§6.4. What Happens in 3D?

§6.4.1. Including the Plane Stress Thickness Dimension

σ1 = 110, σ2 = 10, σ3 = 0. (6.15)

The principal stresses in 3D are the eigenvalues of the 3D stress matrix

Overall All possible

Yellow-filled circle is the inplane one

§6.4.2. 3D Mohr Circles

§6.4.3. Overall Maximum Shear

§6.4.4. Plane Stress Revisited

§6.4.5. The Sphere Paradox

§17.3. Free Vibrations of Undamped SDOF Oscillator

u(0) = u 0 , u̇(0) = v0 . (17.3)

§17.3.1. Equation of Motion Of Undamped Oscillator

§17.3.2. Undamped Response in Terms of Complex Exponentials

λ2 + ωn2 = 0, ⇒ λ1,2 = ±i ωn , (17.7)

§17.3.3. Undamped Response in Terms of Trigonometric Functions

e±iθ = cos θ ± i sin θ, (17.9)

the solution (17.8) becomes

u(t) = (C1 + C2 ) cos ωn t + i(C1 − C2 ) sin ωn t. (17.10)

u(t) = A1 cos ωn t + A2 sin ωn t. (17.11)

u(0) = u 0 = A1 , u̇(0) = v0 = A2 ωn , (17.12)

Replacing into (17.11) yields

u(t) = U cos(ωn t − α). (17.14)

§17.3.4. Effect of Initial Conditions

u(t) = u 0 cos ωn t. (17.16)

Frequency f n is expressed in cycles per second or Hertz, abbreviated to Hz (1 Hz = cycle/s). The

§17.4. Free Vibrations of Viscous-Damped SDOF Oscillator

§17.4.1. Equation of Motion of Damped Oscillator

coefficient. We thus convert (17.18) to the canonical form of the EOM:

§17.4.2. Damped Response in Terms of Trigonometric Functions

u(t) = A eλt . (17.21)

which is quadratic in λ. Its roots are given by

Figure 17.4. Response of underdamped (and critically damped) spring-dashpot-mass oscillator

u(t) = e−ξ ωn t (A1 cos ωd t + A2 sin ωd t) . (17.27)

This equation may in turn be rewritten in the phased form

u(t) = U e−ξ ωn t cos(ωd t − α) (17.29)

This kind of response is illustrated in Figure 17.4.

§17.4.4. Critically Damped Case: ξ = 1