J Eng Math

DOI 10.1007/s10665-011-9478-0

A theoretical approach to the interaction between buckling

and resonance instabilities
Alberto Carpinteri · Marco Paggi

Received: 10 April 2010 / Accepted: 9 May 2011

© Springer Science+Business Media B.V. 2011

Abstract This article deals with the interaction between buckling and resonance instabilities of mechanical sys-
tems. Taking into account the effect of geometric nonlinearity in the equations of motion through the geometric
stiffness matrix, the problem is reduced to a generalized eigenproblem where both the loading multiplier and the
natural frequency of the system are unknown. According to this approach, all of the forms of instabilities inter-
mediate between those of pure buckling and pure forced resonance can be investigated. Numerous examples are
analyzed, including discrete mechanical systems with one to n degrees of freedom, continuous mechanical systems,
such as oscillating deflected beams subjected to a compressive axial load, as well as oscillating beams subjected
to lateral–torsional buckling. A general finite element procedure is also outlined, with the possibility to apply the
proposed approach to any general bi- or tri-dimensional framed structure. The proposed results provide a new
insight in the interpretation of coupled phenomena such as flutter instability of long-span or high-rise structures.

Keywords Buckling · Continuous mechanical systems · Discrete mechanical systems · Finite elements ·
Flutter · Resonance

1 Introduction

Buckling, resonance, and flutter are the main forms of instabilities of the elastic equilibrium of structural systems.
In buckling instability, by removing the hypothesis of small displacements, by means of which deformed structural
configuration can be distinguished from the undeformed one, it is possible to show that the solution of an elastic
problem can represent a condition of stable, neutral, or unstable equilibrium, depending on the magnitude of the
applied load. As is well known, buckling is usually observed in slender structural elements subjected to a compres-
sive stress field, such as columns of buildings, machine shafts, struts of trusses, thin arches, and shells. In some

Electronic pre-print of this manuscript is available at arXiv:0802.0756v1.

A. Carpinteri · M. Paggi (B)

Department of Structural and Geotechnical Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, Turin 10129, Italy
e-mail: [email protected]
A. Carpinteri
e-mail: [email protected]

123
 A. Carpinteri, M. Paggi

cases, elastic instability can also take place for special loading and geometrical conditions, as, for example, in the
lateral-torsional buckling of slender beams [1, Chap. 17].
The phenomenon of resonance is also particularly important in structural engineering. It represents a form of
dynamic instability, which occurs when an external periodic frequency matches one of the natural frequencies of
vibration of the mechanical system. In this case, therefore, structural design deals with the determination of such
dangerous natural frequencies according to modal analysis [2–4].
Finally, the phenomenon of flutter is a form of aeroelastic instability observed in long-span or high-rise structures
subjected to wind loads, such as towers, tall buildings [5] and suspended [6–8] or cable-stayed bridges [9]. In this
case, the instability is attributed to motion-induced or self-excited forces, which are loads induced or influenced
by the deformation of the structure itself [10,11]. Forces originated in this way are coupled to deformations. This
feedback mechanism can give rise to an amplification effect on the initial deformation, leading to premature fail-
ure of the structure. The well-known dramatic Tacoma Narrows bridge disaster of 1940 is a famous example of
this catastrophic interaction, and it is still very much in the public eye today. In this field, which is not at present
completely understood, aeroelastic instability is often considered as the result of the interaction between buckling
(static) and resonance (dynamic) instabilities. However, only a few theoretical formulations have been proposed
for modeling aerodynamic forces and, in most investigations, empirical models are set up in which the parameters
related to the fluid–structure interaction are established by experiments [12].
In the present study, we deal with the phenomenon of interaction between buckling and resonance instabilities.
A state-of-the-art survey of the existing literature shows that this problem has been mainly addressed in the field of
multi-parameter stability theory [13–20], where the conditions for stability of a mechanical system are studied with
reference to a perturbation of the problem parameters. In this framework, instability domains of some continuous
oscillatory systems subjected to buckling loads were provided [17,18]. However, to make the problem analytically
treatable, the analysis was mainly limited to certain continuous mechanical systems such as deflected beams and
beams experiencing lateral-torsional buckling. More importantly, the concept of the geometric stiffness matrix,
typical of structural engineering approaches, was neglected in these contributions.
In our treatment, we will adopt the classical notions of vibration and stability of conservative systems. Stability
problems for conservative systems appear in the study of elastic structures under the action of potential forces like
stationary loads. The determination of the critical loads corresponding to bifurcation of the elastic equilibrium is
required in the design of elastic structures. Determination of frequencies and modes of vibration is also a typical
requirement in the design of structural elements and in many cases modification of frequencies and modes by
changing design parameters is necessary to avoid resonance and noise. In the presence of several parameters, the
graph of natural frequencies of vibration in the parameter space represents a boundary of the stability domain. The
analysis of such graphs is the main focus of the present article, where we consider the influence of stationary loads
on the natural frequencies of vibration of elastic conservative discrete or continuous mechanical systems. Analysis
of these stability boundaries, and therefore the study of the interaction between the aforementioned elementary
forms of instability, allows changing design parameters to either modify the natural frequencies of the system, or
to increase the values of the critical buckling loads.
Special focus will be given to the role played by the geometric stiffness matrix, which contributes to the reduction
of the elastic stiffness matrix due to the effect of the geometric nonlinearity. According to the proposed formulation,
we will demonstrate that the interaction between buckling and resonance leads to a generalized eigenvalue problem
where both the buckling loads and the natural frequencies of the system are unknown and represent eigenvalues.
This methodology was outlined in the valuable textbook [2, Chap. 11], although it was restricted to oscillating
deflected beams under compressive axial loads. In this article, it will be shown that this methodology can be applied
to a generic mechanical system—discrete or continuous. Hence, we will consider discrete mechanical systems with
one to n degrees of freedom and continuous mechanical systems such as oscillating deflected beams and beams
showing lateral-torsional coupled deformations. Moreover, a general procedure is also established in the framework
of the finite element method. This approach will permit to inspect all the forms of structural instability intermediate
between pure buckling and pure resonance. These limit cases are instead observed either when the dynamics of the
system is neglected, or when the external buckling forces are equal to zero.

123
Interaction between buckling and resonance

l l

l l

k ϕ m ϕ N
N

ϕ m ϕ k

l cos ϕ l cosϕ l cosϕ l cosϕ

Fig. 1 Schematic of the first one-degree of freedom system Fig. 2 Schematic of the second one-degree of freedom system

Finally, we will discuss on the possibility to apply the proposed approach to the analysis of flutter instability of
cable-suspended bridges. This will represent a novelty of our approach with respect to the models available in the
literature. In fact, the use of the geometric stiffness matrix may provide the proper link between the multi-parameter
stability theory, typical of rational mechanics, and the bridge engineering approach. Specific structural applications
will be the subject of forthcoming articles.

2 Discrete mechanical systems

2.1 Discrete mechanical systems with one degree of freedom

Let us consider the mechanical system shown in Fig. 1, consisting of two rigid rods connected by an elastic hinge of
rotational stiffness k and constrained at one end by a hinge and at the other by a roller support. A mass m is placed
at the intermediate elastic hinge and the system is loaded by a horizontal axial force N . Considering the absolute
rotation ϕ of the two arms as the generalized coordinate, the total potential energy, W , and the kinetic energy, T ,
of the whole system are given by
1
W (ϕ) = k(2ϕ)2 − 2Nl(1 − cos ϕ),
2
 2  2 (1)
1 d 1 d 1
T (ϕ̇) = m (l sin ϕ) + m (l − l cos ϕ) = ml 2 ϕ̇ 2 .
2 dt 2 dt 2
The equation of motion can be determined by writing Lagrange’s equation:
 
∂ ∂T ∂T ∂W
− =− . (2)
∂t ∂ ϕ̇ ∂ϕ ∂ϕ
In the present case, this yields
ml 2 ϕ̈ = −4kϕ + 2Nl sin ϕ, (3)
which can be suitably linearized about ϕ = 0:
ml 2 ϕ̈ = −4kϕ + 2Nlϕ. (4)
iωt
Looking for the solution to Eq. 4 in the general form ϕ = ϕ0 e , where ω denotes the natural angular frequency
of the vibrating system, we obtain the following equation:
 
4k − 2Nl − ω2 ml 2 ϕ0 = 0. (5)
A nontrivial solution to Eq. 5 exists if and only if the term in brackets is equal to zero. This critical condition
corresponding to the bifurcation of the equilibrium establishes a one-to-one relationship between the applied axial
force, N , and the angular frequency, ω:

123
 A. Carpinteri, M. Paggi

2k ml 2
N= − ω . (6)
l 2
Moreover, Eq. 6 admits two important limit conditions for, respectively, N = 0 and m = 0. In the former  case,
Eq. 6 gives the natural angular frequency of the system according to pure modal analysis, that is ω1 = 4k/ml 2 .
In the latter, the pure critical Eulerian load is obtained, that is N1 = 2k/l. Dividing Eq. 6 by N1 , we obtain the
following relationship between N and ω in a nondimensional form:
 2  
ω N
+ = 1. (7)
ω1 N1
As a second example, let us consider the mechanical system shown in Fig. 2, consisting in two rigid rods on
three supports, of which the intermediate one is assumed to be elastically compliant with stiffness k. As in the
previous case, a mass m is placed at the intermediate hinge and the system is loaded by a horizontal axial force N .
Considering the absolute rotation ϕ of the two arms as the generalized coordinate, the total potential energy, W ,
and the kinetic energy, T , of the whole system are
1
W (ϕ) = k(l sin ϕ)2 − 2Nl(1 − cos ϕ),
2 (8)
1
T (ϕ̇) = ml 2 ϕ̇ 2 .
2
Following the procedure discussed above, we determine the equation of motion by employing Lagrange’s equa-
tion (2):
ml 2 ϕ̈ = −l sin ϕ(kl cos ϕ − 2N ), (9)
which can be suitably linearized about ϕ = 0
ml 2 ϕ̈ = −lϕ(kl − 2N ). (10)
Looking for the solution to Eq. 10 in the general form ϕ = ϕ0 eiωt , where ω denotes the natural angular frequency
of the vibrating system, we obtain the following condition:
 
kl 2 − 2Nl − ω2 ml 2 ϕ0 = 0. (11)

As in the previous example, by setting equal to zero the term in brackets, we obtain a one-to-one relationship
between the applied axial force, N , and the angular frequency, ω:
kl ml 2
N= − ω . (12)
2 2
This equation admits two important limit conditions for, respectively, N = 0 and m = 0. In the former case, Eq. 12
√
gives the natural angular frequency of the system according to pure modal analysis, that is ω1 = k/m. In the
latter, the pure critical Eulerian load for buckling instability is obtained, that is N1 = kl/2. Dividing Eq. 12 by N1 ,
we obtain the same relationship between the nondimensional terms N /N1 and (ω/ω1 )2 as in the previous example
(see Eq. 7).
A graphical representation of the condition (7) in Fig. 3 shows that the resonance frequency is a decreasing func-
tion of the compressive axial load. This demonstrates, for the mechanical systems with a single degree of freedom,
that the condition of bifurcation of the equilibrium can be reached for a compressive axial force, N , lower than the
Eulerian buckling load, N1 , provided that the system is subjected to an external excitation with frequency ω given
by Eq. 7. Conversely, failure due to resonance can take place for ω < ω1 , provided that the system is loaded by an
axial force N given by Eq. 7.
Finally, the issue of stability or instability of the mechanical system in the correspondence of the bifurcation
point can be discussed as in the static case, i.e., by evaluating the higher order derivatives of the total potential
energy W .

123
Interaction between buckling and resonance

l l l

k k
N
x1 m m
x2

Fig. 3 Nondimensional frequency versus nondimensional axial Fig. 4 Schematic of the two-degrees of freedom system
force for the single degree of freedom systems

2.2 Discrete mechanical systems with n degrees of freedom

Let us consider the mechanical system with two degrees of freedom shown in Fig. 4, consisting of three rigid rods
connected by two elastic hinges of rotational rigidity k, and constrained at one end by a hinge and at the other by a
roller support. A mass m is placed at the intermediate elastic hinges and the system is loaded by a horizontal axial
force N . Assuming the vertical displacements x1 and x2 of the elastic hinges as the generalized coordinates, the
total potential energy, W , and the kinetic energy, T , of the whole system are given by
   

1 x1 x2 − x1 2 x2 x2 − x1 2
W (x1 , x2 ) = k arcsin − arcsin + arcsin + arcsin
2 l l l l
      
x1 x2 x2 − x1
− Nl 3 − cos arcsin − cos arcsin − cos arcsin , (13)
l l l
 
1 1 1 1 2x1 ẋ1 x2 ẋ2 x2 ẋ1 x1 ẋ2 2
T (ẋ1 , ẋ2 ) = m ẋ12 + m ẋ12 x12 + m ẋ22 + m + − − .
2 2 2 2 l l l l
Performing a Taylor series expansion of Eq. 13 about the origin, and assuming x1 /l < 1/10 and x2 /l < 1/10,
we obtain
k   N  
W (x1 , x2 ) ∼ = 2 5x12 + 5x22 − 8x1 x2 − x12 + x22 − x1 x2 ,
2l l
1 1
∼ m ẋ + m ẋ .
T (ẋ1 , ẋ2 ) = 2 2
(14)
2 1 2 2
The equations of motion are identified by considering Lagrange’s equations:
 
∂ ∂T ∂T ∂W
− =− , i = 1, 2. (15)
∂t ∂ x˙i ∂ xi ∂ xi
In matrix form, they are
⎛ ⎞ ⎛ ⎞
   5k 4k   2 1    
m 0 ẍ1 ⎜ 2 − 2 ⎟ x ⎜ − ⎟ x 0
+ ⎝ l4k 5kl ⎠ 1 − N ⎝ l 1 2l ⎠ 1 = (16)
0 m ẍ2 x2 − x 0
− 2 2
l l2 l l
The same result could be arrived at by considering the nonlinear expressions of W and T , writing Lagrange’s
equations and then linearizing them. In this article, we have preferred to work with the linearized expressions of
W and T to simplify the algebra involved in the differentiation. Looking for the solution to Eq. 16 in the general

123
 A. Carpinteri, M. Paggi

form q = q0 eiωt , where ω denotes the natural angular frequency of the system, we obtain the following equation,
written in symbolic form:
 
−ω2 [M] + [K ] − N [K g ] q0 = 0, (17)
where [M], [K ] and [K G ] denote, respectively, the mass matrix, the elastic stiffness matrix and the geometric
stiffness matrix of the mechanical system. Their expressions can be simply obtained by comparing Eq. 17 with
Eq. 16.
A nontrivial solution to Eq. 17 exists if and only if the determinant of the resultant coefficient matrix of the vector
q0 vanishes. This yields the following generalized eigenvalue problem:
 
det [K ] − N [K g ] − ω2 [M] = 0, (18)

where N and ω2 represent the eigenvalues. For this example, Eq. 18 provides the following relationships between the
eigenvalues ω and N :
k N1
ω12 = 2
− , (19a)
ml ml
k 3N2
ω22 = 3 2 − . (19b)
ml ml
As limit cases, if m = 0, then we obtain the Eulerian buckling loads (bl):
k
N1bl = , (20a)
l
3k
N2bl = , (20b)
l
whereas, if N = 0, then we obtain the natural frequencies (nf) of the system:

k
ω1 =
nf
, (21a)
ml 2

9k
ω2nf = . (21b)
ml 2
System (17) yields the eigenvectors corresponding, respectively, to the eigenfrequencies (19a) and (19b) as
functions of N :
⎛ ⎞ ⎛ ⎞
  4k/l − N1   4k/l − N2
x1 x x1
= ⎝ 6k/l − 3N1 2 ⎠ , = ⎝ 14k/l − 5N2 ⎠ (22)
x2 x2
1 1
Dividing Eqs. 19a and 19b by (ω1nf )2 , we derive the following nondimensional relationships between the eigen-
values:
 2  
ω1 N1
= 1− , (23a)
ω1nf N1bl
 2  2  
ω2 ω2nf N2bl N2
= − bl . (23b)
ω1nf ω1nf N1 N1bl
In analogy with the results for the single degree of freedom systems, a graphical representation of Eqs. 23a and
23b is provided in Fig. 5. We notice that both the eigenfrequencies are decreasing functions of the compressive
axial load. Starting from N = 0, bifurcation of the equilibrium would correspond to pure resonance instability.
Entering the diagram with a value of the nondimensional compressive axial force in the range 0 < N /N1bl < 1,
the coordinates of the points of the two curves provide the critical eigenfrequencies of the mechanical system
leading to bifurcation. Axial forces higher than N1bl in the range 1 < N /N1bl < N2bl /N1bl can only be experienced

123
Interaction between buckling and resonance

(ω nf2/ωnf1)2=9
2
(ω /ω1)
nf

l l l
1

x2
bl bl
x1 m m N
0 1 3=N 2/N 1
bl
N /N 1 k k
Fig. 5 Nondimensional frequencies versus nondimensional Fig. 6 Schematic of the second two-degrees of freedom system
axial forces for the two-degrees of freedom system in Fig. 4

if an additional constraint is introduced into the system. Moreover, we observe that the applied compressive load
influences all the eigenfrequencies. In particular, for the present example, the influence of the axial load is greater
on the highest frequency than on the lower one.
As a second example of a system with two degrees of freedom, let us examine that of Fig. 6, which consists of
three rigid rods on four supports, of which the central ones are assumed to be elastically compliant with stiffness k.
A mass m is placed at the intermediate hinges and the system is loaded by a horizontal axial force N . Assuming the
vertical displacements x1 and x2 of the elastic hinges as the generalized coordinates, the total potential energy, W ,
and the kinetic energy, T , of the whole system are given by (x1 /l < 1/10 and x2 /l < 1/10):
1     x1 
W (x1 , x2 ) = k x12 + x22 − Nl 3 − cos arcsin
2  l 
 x2  x2 − x1
− cos arcsin − cos arcsin
l l
  N   (24)
∼ 1
= k x1 + x2 −
2 2
x1 + x2 − x1 x2 ,
2 2
2 l
∼ 1 m ẋ 2 + 1 m ẋ 2 .
T (ẋ1 , ẋ2 ) =
2 1 2 2
In this case, the linearized Lagrange’s equations (15) yield the following matrix form:
⎛ ⎞
      2 1    
m 0 ẍ1 k0 x1 ⎜ − ⎟ x 0
+ − N ⎝ l 1 2l ⎠ 1 = . (25)
0 m ẍ2 0k x2 − x 2 0
l l
Looking for the solution to Eq. 25 in the general form q = q eiωt , where ω denotes the natural angular frequency
0
of the system, we obtain the following equation, written in symbolic form:
 
−ω2 [M] + [K ] − N [K g ] q0 = 0, (26)

123
 A. Carpinteri, M. Paggi

where [M], [K ], and [K g ] denote, respectively, the mass matrix, the elastic stiffness matrix, and the geometric
stiffness matrix of the mechanical system. As can be readily seen, the geometric stiffness matrix for this problem
is the same as that of the previous example.
A nontrivial solution to Eq. 26 exists if and only if the determinant of the resultant coefficient matrix of the vector
q0 is equal to zero. This yields the following generalized eigenvalue problem:
 
det [K ] − N [K g ] − ω2 [M] = 0, (27)

where N and ω2 are the eigenvalues of the system. For this example, Eq. 27 provides the following relationships
between the eigenvalues:
k N1
ω12 = −3 , (28a)
m ml
k N2
ω22 = − . (28b)
m ml
As limit cases, if m = 0, then we obtain the Eulerian buckling loads:
1
N1bl = kl, (29a)
3
N2bl = kl, (29b)
whereas, if N = 0, then the natural frequencies of the system coincide:

k
ω1 = ω2 =
nf nf
. (30)
m
System (26) yields the eigenvectors corresponding, respectively, to the eigenfrequencies (28a) and (28b), as
functions of the axial force:
⎛ ⎞ ⎛ ⎞
  N1 /l   N2 /l
x1 x
= ⎝ 5N1 /l − 2k ⎠ , 1
= ⎝ 3N2 /l − 2k ⎠ (31)
x2 x2
1 1

Dividing Eqs. 28a and 28b by (ω1nf )2 , we obtain the following nondimensional relationships between the eigen-
values:
 2  
ω1 N1
= 1− , (32a)
ω1nf N1bl
 2  
ω2 N1bl N2
= 1 − bl . (32b)
ω1nf N2 N1bl
A graphical representation of Eqs. 32a and 32b is provided in Fig. 7. Also in this case, both the frequencies are
decreasing functions of the compressive axial load. However, for the present example, the influence of the axial
load is greater on the lower frequency of the system than on the higher one.

3 Continuous mechanical systems

3.1 Oscillations of deflected beams under compressive axial loads

Let us consider a slender elastic beam of constant cross-section, inextensible and not deformable in shear, though
deformable in bending, constrained at one end by a hinge and at the other by a roller support, loaded by an axial
force, N (see Fig. 8). In this case, with the purpose of analyzing the free flexural oscillations of the beam, the

123
Interaction between buckling and resonance

1
2
(ω /ωnf1)

0 1 N bl2/N bl1=3
N /N bl1
Fig. 7 Nondimensional frequencies versus nondimensional Fig. 8 Undeformed and deformed configurations of a deflected
axial forces for the two-degrees of freedom system in Fig. 6 beam under compressive axial force

differential equation of the elastic line with second-order effects can be written by replacing the distributed load
with the force of inertia:
∂ 4v ∂ 2v ∂ 2v
EI + N = −μ , (33)
∂z 4 ∂z 2 ∂t 2
where EI denotes the flexural rigidity of the beam, and μ is its linear density (mass per unit length). Equation 33
can be rewritten in the following form:
∂ 4v 2∂ v
2 μ ∂ 2v
+ β = − , (34)
∂z 4 ∂z 2 EI ∂t 2
where we have set β 2 = N /EI.
Equation 34 is an equation with separable variables, the solution being represented as the product of two different
functions, each one depending on a single variable:
v(z, t) = η(z) f (t). (35)
Introducing Eq. 35 into Eq. 34 and dividing the resulting equation by η f leads to
d2 f d4 η d2 η
+ β2 2
2 EI dz 4 dz = ω2 ,
− dt = (36)
f μ η
where ω2 represents a positive constant, the left- and the right-hand sides of Eq. 36 being at the most functions of
the time t and the coordinate z, respectively. From Eq. 36, there follow two ordinary differential equations:
d2 f
+ ω2 f = 0, (37a)
dt 2
d4 η 2d η
2
+ β − α 4 η = 0, (37b)
dz 4 dz 2
with

4 μω
2
α= . (38)
EI
While Eq. 37a is the equation of the harmonic oscillator, with the well-known general solution
f (t) = A cos ωt + B sin ωt, (39)

123
 A. Carpinteri, M. Paggi

Equation 37b has the following general solution:

η(z) = Ceλ1 z + Deλ2 z + Ee−λ1 z + Fe−λ2 z , (40)
where λ1 and λ2 are functions of α and β:
 
−β 2 ± β 4 + 4α 4
λ1,2 = . (41)
2
As in the modal analysis, the constants A and B can be determined on the basis of the initial conditions, while the
constants C, D, E, and F can be determined by imposing the boundary conditions. As will be shown in the sequel,
for a given value of β, the parameters ω and α can be determined by solving a generalized eigenvalue problem
resulting from the imposition of the boundary conditions. From the mathematical point of view, this eigenvalue
problem is analogous to that shown for the discrete systems. On the other hand, since we are considering a contin-
uous mechanical system having infinite degrees of freedom, we shall obtain an infinite number of eigenvalues ωi
and αi , just as also an infinite number of eigenfunctions f i and ηi . The complete integral of the differential equation
(33) may therefore be given the following form, according to the principle of superposition:
∞
v(z, t) = ηi (z) f i (t), (42)
i=1
with
f i (t) = Ai cos ωi t + Bi sin ωi t, (43a)
λ1i z λ2i z −λ1i z −λ2i z
η(z) = Ci e + Di e + Ei e + Fi e . (43b)
It is important to remark that the eigenfunctions ηi are still orthonormal functions, as in the classical modal analysis
(see the mathematical demonstration reported in the Appendix). This permits us to determine the constants Ai and
Bi in Eq. 43a via the initial conditions (see also [1, p. 315]):
v(z = 0) = v0 (z), (44a)
∂v
(z = 0) = v̇0 (z). (44b)
∂t
As regards the boundary conditions, let us consider as an example a beam supported at both ends, of length l:
⎧ ⎛ ⎞⎛ ⎞ ⎛ ⎞
⎪ η(0) = 0, 1 1 1 1 C 0
⎪
⎪
⎨ η (0) = 0, ⎜ λ2 λ22 λ12 λ2 ⎟ ⎜ D ⎟ ⎜ 0 ⎟
2 ⎟ ⎜ ⎟ ⎜
⎜ ⎟
⇒ ⎜ λ1l λ −λ −λ ⎟⎜ ⎟ = ⎜ ⎟. (45)
⎪
⎪ η(l) = 0, ⎝ e 1 e 2 l e 1 l e 2 l ⎠ ⎝ E ⎠ ⎝ 0⎠
⎪
⎩ 
η (l) = 0, λ21 eλ1 l λ22 eλ2 l λ21 e−λ1 l λ22 e−λ2 l F 0
For a nontrivial solution to the system in Eq. 45, the determinant of the coefficient matrix has to vanish. The
resulting eigenequation determines, for each given value of the parameter β, the eigenvalues αi of the system.
Finally, the corresponding natural eigenfrequencies ωi can be obtained by inverting Eq. 38.
As an illustrative example, the first three nondimensional frequencies of the simply supported beam shown in
Fig. 8 are shown in Fig. 9 as functions of the applied nondimensional axial force. Parameters, ωi and Ni denote,
respectively, the ith frequency of the system determined according to modal analysis and the ith buckling load
determined according to the Euler’s formula. In close analogy with the discrete mechanical systems, the curves in
the (ω/ω1nf )2 vs. N /N1bl plane are represented by straight lines. Also in this case, the coordinates of the points along
these lines provide the critical conditions leading to the system instability in terms of frequency of the excitation
and magnitude of the applied compressive axial force.

3.2 Oscillations and lateral-torsional buckling of beams

Let us consider a beam of thin rectangular cross section, constrained at the ends so that rotation about the longitu-
dinal axis Z is prevented. Let this beam be subjected to uniform bending by means of the application at the ends

123
Interaction between buckling and resonance

(ω /ω )2= 81
3 1
(ω/ω )2
1

2
(ω2/ω1) = 16

1
0 1 4 =N2/N1 9 =N3/N1
N/N
1

Fig. 9 Nondimensional frequencies versus nondimensional Fig. 10 Schematic of a beam subjected to lateral-torsional buck-
axial force for the continuous system in Fig. 8 ling. Sketch of the beam in the YZ plane (a), lateral deflection
of the beam in the XZ plane (b), the rotation du/dz (c), and the
rotation ϕz (d)

of two moments m contained in the plane Y Z of greater flexural stiffness. The sketch of the beam in the YZ plane
is shown in Fig. 10a. The lateral deflection of the beam in the XZ plane is shown in Fig. 10b. The rotations du/dz
and ϕz are illustrated, respectively, in Fig. 10c, d.
Considering a deformed configuration of the beam, with deflection thereof in the X Z plane of smaller flexural
rigidity, and simultaneous torsion about the Z axis (see Fig. 10), bending–torsional out-of-plane vibrations of the
beam are described by the following partial differential equations:

∂ 4u ∂ 2 ϕz ∂ 2u
EI y + m 2 = −μ 2 ,
∂z 4 ∂z ∂t
(46)
∂ 2 ϕz ∂ 2u ∂ 2 ϕz
− GIt 2 + m 2 = −μρ 2 2 ,
∂z ∂z ∂t
where u(z, t) and ϕz (z, t) are, respectively, the out-of-plane deflection and the twist angle of the beam cross section;
√
EI y and GIt are the bending and torsional rigidities; μ is the mass of the beam per unit length, and ρ = I P /A is
the polar radius of inertia of the beam cross-section.
A solution to the system (46) can be found in the following variable-separable form [21]:

u(z, t) = U (t)η(z), (47)

ϕz (z, t) =
(t)ψ(z), (48)

where the functions η(z) and ψ(z) are such that the boundary conditions η(0) = η(l) = η (0) = η (l) = ψ(0) =
nπ z
ψ(l) = 0 are satisfied. According to Bolotin [21], we can assume η(z) = ψ(z) = sin , with n being a natural
l
number. In this case, we obtain the following matrix form:

123
 A. Carpinteri, M. Paggi

⎛ ⎞ ⎛ ⎞
   n4π 4   n2π 2    
μ 0 ⎜ EI 0 ⎟ U ⎜ 0
l2 ⎟
Ü y U 0
+⎝ l4 ⎠ −m⎝ ⎠ = , (49)
0 μρ 2
¨ n2π 2
n2π 2
0
0 GIt 2 0
l l2
which can be symbolically rewritten as
 
[M] q̈ + [K ] q − m K g q = 0, (50)
 
where q = (U,
)T . The mass matrix [M], the elastic stiffness matrix [K ], and the geometric stiffness matrix K g
in Eq. 50 can be defined in comparison with Eq. 49. Looking for a general solution in the form q = q0 eiωt , we
obtain:
   
[K ] − m K g − ω2 [M] q0 = 0. (51)
Looking for nontrivial solutions, we derive the following generalized eigenvalue problem:
   
det [K ] − m K g − ω2 [M] = 0, (52)

where m and ω2 are the eigenvalues of the system.

As limit cases, if μ = 0, then we obtain the critical bending moments given by Prandtl’s formula [1, Chap. 17,
p. 559]:
nπ 
m nc = EI y GIt , (53)
l
whereas, if m = 0, then we obtain the natural flexural and torsional eigenfrequencies of the beam:

 nπ 2 EI
y
ωn =
flex
,
l μ

nπ GIt
ωntors = . (54)
ρl μ
Considering a rectangular beam with a depth-to-span ratio of 1/3 and with a thickness to depth ratio of 1/10,
the evolution of the first two flexural and torsional eigenfrequencies of the system are shown in Fig. 11 as functions
of the applied bending moment. In this case, the curves in the nondimensional plane (ω/ω1flex )2 vs. m/m 1c are
no longer straight lines. This fact can be ascribed to the coupling between torsional and flexural vibrations of the
beam. Moreover, when m is increased from zero (pure resonance instability) up to the critical bending moment
computed according to Parandtl’s formula (pure buckling instability), m 1c , we note that the resonance frequen-
cies related to flexural oscillations progressively decrease from ω1flex down to zero in correspondence of the critical
bending moments given by Prandtl’s formula. Conversely, the resonance frequencies related to torsional oscillations
increase. From the mathematical point of view, this is because the sum of the squares of the two eigenfrequencies
for a given value of n is constant when the applied moment m is varied.

4 Finite elements procedure

When the mechanical system cannot be reduced to the schemes previously analyzed, it is possible to apply the finite
element method [22,23]. According to this approach, the equations of motion for an elastic system with a finite
number of degrees of freedom can be expressed in matrix form, taking also into account the effect of the geometric
nonlinearity through the geometric stiffness matrix [24,25].
For the sake of generality, let the elastic domain V be divided into subdomains Ve , and let each element contain
m nodal points, each one having g degrees of freedom. In compact form, the displacement vector field defined on
the element Ve may be represented as
ηe g,1 = [ηe ]g,g×m δe g×m,1 (55)

123
Interaction between buckling and resonance

Fig. 11 Nondimensional 90
flexural and torsional 80 n=2 (torsional)
frequencies versus
nondimensional bending 70

moments for the continuous 60

2
system in Fig. 10

(ω/ωflex)
50

1
40

30
n=1 (torsional)
20

10 n=2 (flexural)
n=1 (flexural)
0
0 0.5 1 1.5 2
m/m
1c

where [ηe ] of dimensions (g, g × m) is the matrix collecting the shape functions and δe is the nodal displacements
vector.
The deformation characteristic vector is obtained by differentiation:
qe d,1 = [∂]d,g ηe g,1 , (56)
where [∂] is the kinematic operator relating strains to displacements, whereas d denotes the dimension of the strain
characteristic vector, e.g., d = 3 for a beam in plane or d = 6 for a beam in space. Hence, introducing Eq. 55 into
Eq. 56, we obtain:
qe = [∂] [ηe ] δe = [Be ] δe (57)
where the matrix [Be ] is calculated by derivation of the shape functions.
According to these definitions, we can obtain the following matrix equation for each finite element (see
[1, Chap. 11] for more details):
  
[Me ] δ¨e + [K e ] − K ge δe = 0, (58)
 
where [Me ] , [K e ], and K ge denote, respectively, the local mass matrix, the local elastic stiffness matrix, and the
local geometric stiffness matrix of the finite element.
As usual, the local mass and elastic stiffness matrices are given by [22,23]:

[Me ] = [ηe ]T [μ] [ηe ] dV, (59a)
Ve

[K e ] = [Be ]T [H ] [Be ] dV. (59b)
Ve
The local geometric stiffness matrix can be computed as follows [24]:

 
K ge = [G e ]T [Se ] [G e ] dV, (60)
Ve
where the matrix [Se ] is related to the components of the stress field for a 3D continuum (g = 3):
⎛ ⎞
σx [I3 ] τx y [I3 ] τx z [I3 ]
[Se ] = ⎝ τx y [I3 ] σ y [I3 ] τ yz [I3 ] ⎠ , (61)
τx z [I3 ] τ yz [I3 ] σz [I3 ]
where [I3 ] denotes a unit matrix with dimensions (3
 × 3). The matrix [G e ] is related to the first derivative of the
shape functions through the differential operator ∂ :
 
[G e ] = ∂ [ηe ], (62)

123
 A. Carpinteri, M. Paggi

where
⎛ ⎞
∂
[I
⎜ ∂x 3 ⎟ ]
  ⎜
⎜ ∂ [I ] ⎟
⎟
∂ =⎜ 3 ⎟. (63)
⎜ ∂y ⎟
⎝ ∂ ⎠
[I3 ]
∂z
According to this formulation, we note that the geometric stiffness matrix is a function of the stress components
through the matrix [Se ]. In the case of a compressive stress field, the geometric stiffness terms become negative and
reduce the corresponding elements of the local elastic stiffness matrix, just as shown for the discrete mechanical
systems. We also remark that this formulation is quite general, since the information related to the finite element
topology is simply included in the matrix [ηe ] which collects the shape functions and in the differential operator
[∂] (see [1, Chap. 17] for more details).
By performing the usual operations of rotation, expansion, and assemblage of the mass, elastic stiffness, and
geometric stiffness matrices of the element, Eq. 58 can be written in global form:
  
[M] δ̈ + [K ] − λ K g δ = 0. (64)

Looking for the solution to Eq. 64 in the general form δ = δ0 eiωt , where ω is the natural frequency of the system,
we can formulate the generalized eigenproblem as in the cases discussed above:
   
det [K e ] − λ K g − ω2 [Me ] = 0. (65)

Therefore, the numerical procedure for the determination of the frequency-loading multiplier diagram consists of
the following steps.
1. For a given loading configuration defined by the loading multiplier λ, determine the stress field according to a
linear elastic stress analysis.
2. Compute the local mass matrix, the local elastic stiffness matrix and the local geometric stiffness matrix for
each finite element.
3. Perform the rotation, expansion, and assembly operations to obtain the global matrices.
4. Solve the generalized eigenvalue problem of Eq. 65 and find the eigenfrequencies of the system, ωi2 , with
i = 1, . . . , g × n.
5. Iterate the above-described procedure for different values of λ.

5 Discussion and conclusions

The problems of elastic instability (buckling) and dynamic instability (resonance) have been the subject of extensive
investigation and have received a large attention from the structural mechanics community. Nonetheless, the study
of the interaction between these elementary forms of instability is still an open problem.
The phenomenon of flutter instability of the Tacoma Narrows Bridge occurred on November 7, 1940, can be
reinterpreted as the result of such a catastrophic interaction. This cable-suspended bridge was solidly built, with
girders of carbon steel anchored in huge blocks of concrete and was the first of its type to employ plate girders to
support the roadbed. While in the earlier designs any wind would simply pass through the truss, in the new design
of the 1940, the wind would be diverted above and below the structure. Shortly after construction, it was discovered
that the bridge would sway and buckle dangerously in windy conditions. This resonance was flexural, meaning the
bridge buckled along its length, with the roadbed alternately raised and depressed in certain locations. However,
the failure of the bridge occurred when a never-before-seen twisting mode occurred.
The report to the Federal Works Agency [26] excluded the phenomenon of pure forced resonance as the actual
reason of instability: “…it is very improbable that resonance with alternating vortices plays an important role in

123
Interaction between buckling and resonance

the oscillations of suspension bridges. First, it was found that there is no sharp correlation between wind velocity
and oscillation frequency such as is required in case of resonance with vortices whose frequency depends on the
wind velocity….” A new theory for the interpretation of these complex aerodynamic instabilities was developed by
Scanlan [27,28] and then elaborated by various researchers [29–32]. Basically, the so-called flutter theory considers
the following equation of motion for the mechanical system in the finite element framework [9,30]:

[M] δ̈ + [C] δ̇ + [K ] δ = Fmi + Fmd , (66)

where Fmi and Fmd are, respectively, the motion-independent wind force vector and the motion-dependent aero-
elastic force vector. A special attention is given to the structural damping, which is included in the equations of
motion through the damping matrix [C]. The motion-dependent force vector is then put in relationship with the
nodal displacements of the system, δ, and the nodal velocities, δ̇, according to the flutter derivative matrices, [K ∗ ]
and [C ∗ ], that are empirically determined in the wind tunnel by using section models of the bridge [33]. As a result,
the problem becomes highly nonlinear, and the flutter velocity, Ucr , and the flutter frequency, ωcr , can be determined
from the following eigenproblem:
 
1  ∗ 1  ∗
det [K ] − ρUcr K − ωcr [M] + ωcr [C] − ρUcr ωcr C
2 2
= 0, (67)
2 2
where ρ is the air density and 1/2ρUcr2 is the wind pressure.
It is important to remark that the eigenproblem resulting from the classical approach to flutter instability shares
most of the features of the generalized eigenproblem that we have analyzed in the present study. In fact, in both
cases, two eigenvalues have to be determined from the eigenequation. However, the flutter theory gives prominence
to the role played by the structural damping, although this is generally less than 1% (see e.g., [30]). Moreover,
the value of mechanical damping seems to represent a sort of free parameter in the model. In fact, this parameter
is usually assumed, like in [30], rather than experimentally evaluated. Hence, the classical approaches to flutter
instability used in the field of bridge engineering do not provide a single pair of critical frequency and wind speed,
but rather a number of pairs, each one for a given value of the assumed structural damping. This suggests that
also flutter instability can be mathematically treated in the framework of multi-parameter stability theory, where
the stability boundary of the graph relating the flutter frequency to the critical wind speed is of main concern.
The problem of flutter instability can therefore be reconsidered according to a completely different route. The
complex wind–structure interaction can be approximated by the sum of the action of a stationary load, given by the
wind pressure acting on the bridge deck, and the action of wind gusts, characterized by their specific frequency.
Under such conditions, the interaction between buckling and resonance would apply as for the mechanical systems
addressed in the present article. The problem is now in finding the boundary of the stability domain in the graph
relating the natural frequencies to the applied stationary wind pressure. Therefore, it is the combination of the
wind speed and of the frequency of the wind gusts which is responsible for the bridge instability and not their
individual values. Moreover, according to this approach, the geometric stiffness matrix has a preeminent role and
the mechanical damping can be neglected, as usually done in most of structural engineering applications.
On this line, the collapse of the Tacoma Narrows bridge could be considered as the result of the interaction
between buckling (related to the wind pressure proportional to the square of the wind velocity) and resonance
(caused by the frequency of the wind gusts). Thus, this would give a new explanation on why the Tacoma Narrows
bridge failure took place under moderate wind velocities (wind pressure lower than the critical buckling load) and
wind gusts frequencies different from the natural frequencies of the bridge. Future developments of the present
study will regard the assessment of the proposed approach to the analysis of bridge instabilities, as well as the
comparison with the classical flutter theory on the basis of real case histories.
Finally, the role of damping can also be included in the formulation. Damping may have in fact an important
effect, as shown in [34].

Acknowledgment The support provided by the European Union to the Leonardo da Vinci project “Innovative Learning and Training
on Fracture (ILTOF)” is gratefully acknowledged.

123
 A. Carpinteri, M. Paggi

Appendix: Orthonormality of the eigenfunctions of deflected beams subjected to an axial force

As is well known, the eigenfunctions ηi of deflected beams computed according to pure modal analysis are ortho-
normal functions. It is possible to demonstrate that this fundamental property still holds when the beam is subjected
to an axial load, N , as that shown in Fig. 8. We may in fact write Eq. 37b for two different eigensolutions:
η 2 
j + β ηj = αjηj,
4
(68a)
ηk + β 2 ηk = αk4 ηk . (68b)
Multiplying Eq. 68a by ηk and Eq. 68b by η j , and integrating over the beam length, we obtain
l l l
ηk η
j dz +β 2
ηk ηj dz = α 4j ηk η j dz, (69a)
0 0 0
l l l
η j ηk dz + β 2 η j ηk dz = αk4 η j ηk dz. (69b)
0 0 0
Integrating by parts the left-hand sides, the foregoing equations transform as follows:
 !l  !l  l  l
l l
      
ηk η j − ηk η j + ηk η j dz + β ηk η j 0 − β
2 2
ηk η j dz = α j ηk η j dz,
4
(70a)
0 0
0 0 0
l l l
 l  !l  l
η j ηk 0 − ηj ηk + ηj ηk dz + β η j ηk 0 − β 2
2
ηj ηk dz = αk4 η j ηk dz. (70b)
0
0 0 0
When each of the two ends of the beam is constrained by a built-in support (η = η = 0), or by a hinge (η =
η = 0), the quantities in square brackets vanish. On the other hand, when the end in z = 0 is either unconstrained
(η = η = 0), or constrained by a double rod (η = η = 0), the remaining end of the beam has to be constrained
either by a built-in support (η = η = 0), or by a simple support (η = η = 0). For both configurations, only the
terms [ηi ηk ]l0 are different from zero.
In any case, subtracting member by member, these quantities are canceled, and we have
  l
α 4j − αk4 η j ηk dz = 0, (71)
0
which leads to the orthonormality condition:
l
η j ηk dz = δi j , (72)
0
where δi j is the Kronecker delta. Thus, when the eigenvalues are distinct, the integral of the product of the corre-
sponding eigenfunctions vanishes. When, instead, the indices j and k coincide, the condition of normality reminds
us that the eigenfunctions are defined neglecting a factor of proportionality.

References

1. Carpinteri A (1997) Structural mechanics: a unified approach. Chapman & Hall, London
2. Clough RW, Penzien J (1975) Dynamics of structures. McGraw-Hill, New York
3. Carpinteri A, Pugno N (2005) Towards chaos in vibrating damaged structures. Part I: Theory and period doubling. J Appl Mech
72:511–518

123
Interaction between buckling and resonance

4. Carpinteri A, Pugno N (2005) Towards chaos in vibrating damaged structures. Part II: Parametrical investigation. J Appl Mech
72:519–526
5. Kawai H (1993) Bending and torsional vibration of tall buildings in strong wind. J Wind Eng Ind Aerodyn 50:281–288
6. Astiz MA (1998) Flutter stability of very long suspension bridges. ASCE J Bridge Eng 3:132–139
7. Scott R (2001) In the wake of Tacoma: suspension bridges and the quest for aerodynamic stability. ASCE, Reston
8. Scanlan RH (1996) Aerodynamics of cable-supported bridges. J Constr Steel Resour 39:51–68
9. Starossek U (1993) Prediction of bridge flutter through the use of finite elements. Struct Eng Rev 5:301–307
10. Simiu EH, Scanlan RH (1986) Wind effects on structures. Wiley, New York
11. Billah KY, Scanlan RH (1991) Resonance, Tacoma Narrows bridge failure, and undergraduate physics textbooks. Am J Phys
59:118–124
12. Scanlan RH (1994) Toward introduction of empiricism in the wind design of long-span bridges. In: Proceedings of the international
conference on cable stayed and suspension bridges, pp 15–28
13. Liapunov AM (1892) General problem of stability of motion. Kharkov (reproduced in Annals of Mathematics Studies 17, Princeton
University Press, Princeton 1949)
14. Bolotin VV (1963) Non-conservative problems of the theory of elastic stability. Pergamon Press, New York
15. Bolotin VV (1964) The dynamic stability of elastic systems. Holden-Day, San Francisco (Russian original: Moscow, 1956)
16. Seyranian AP (1991) Stability and catastrophes of vibrating systems depending on parameters. Technical University of Denmark,
Lyngby
17. Seyranian AP, Mailybaev AA (2003) Multiparameter stability theory with mechanical applications. World Scientific, Singapore
18. Huseyin K (1986) Multiple parameter stability theory and its applications: bifurcation, catastrophes, instabilities. Clarendon, Oxford
19. Huseyin K (2002) Dynamics, stability, and bifurcations. ASME Appl Mech Rev 55:R5–R15
20. Elfelsoufi Z, Azrar L (2005) Buckling, flutter and vibration analyses of beams by integral equation formulations. Comput Struct
83:2632–2649
21. Bolotin VV (1995) Dynamic stability of structures. In: Kounadis AN, Kratzig WB (eds) Nonlinear stability of structures: theory
and computational techniques. Springer, Wien pp 3–72
22. Bathe KJ (1982) Finite element procedures in engineering analysis. Prentice-Hall, Englewood Cliff
23. Zienkiewicz OC, Taylor RL (2005) The finite element method for solid and structural mechanics, 6th edn. Elsevier, Amsterdam
24. Rutenberg A (1982) Simplified P-delta analysis for asymmetric structures. ASCE J Struct Div 108:1995–2013
25. Przemieniecki JS (1985) Theory of matrix structural analysis. Dover, New York
26. Ammann OH, Von Karman T, Woodruff GB (1941) The failure of the Tacoma Narrows bridge. Report to the Federal Works Agency,
March 28
27. Scanlan RH, Tomko JJ (1971) Airfoil and bridge deck flutter derivatives. J Eng Mech 97:1717–1737
28. Scanlan RH (1978) The action of flexible bridges under winds, I: flutter theory. J Sound Vib 60:187–199
29. Como M, Grimaldi A, Maceri F (1985) Statical behaviour of long-span cable-stayed bridges. Int J Solids Struct 21:831–850
30. Namini A, Albrecht P, Bosch H (1992) Finite element-based flutter analysis of cable-suspended bridges. ASCE J Struct Eng
118:1509–1526
31. Brancaleoni F, Diana G (1993) The aerodynamic design of the Messina Straits bridge. J Wind Eng Ind Aerodyn 48:395–409
32. Brancaleoni F, Brotton DM (2005) The role of time integration in suspension bridge dynamics. Int J Numer Methods Eng 20:
715–732
33. Scanlan RH (1981) State-of-the-art methods for calculating flutter, vortex-induced and buffeting response of bridge structures.
Federal Highway Administration, Report No. FHWA/RD-80/050, Washington, DC
34. Holmes PJ (1978) Pipes supported at both ends cannot flutter. J Appl Mech 45:619–622

123