STABILITY INVESTIGATIONS OF NON-CONSERVATIVE DYNAMIC SYSTEMS

By

Mahmoud Nabil Abdullatif

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Mechanical Engineering – Doctor of Philosophy

2020

STABILITY INVESTIGATIONS OF NON-CONSERVATIVE DYNAMIC

ABSTRACT

SYSTEMS

By

Mahmoud Nabil Abdullatif

Elastic stability continues to be a subject of considerable interest due to the increasing de-

mands of lightweight structures. It is well known that structures may lose stability through

divergence or flutter and the nature of the instability depends on the type of loading. A

structure may lose stability through divergence when the loading is conservative but may

lose stability through divergence or flutter when the loading is non-conservative.

In this

work, we investigate stability transitions in structures due to: damping in the presence of

non-conservative loading, terminal dynamic moment with and without intermediate support,

and external flow with dynamic moment. The role of damping on non-conservative systems,

which was first investigated by Ziegler, is revisited here. It is shown that increasing the level

of damping can cause an unstable non-conservative system to become stable, then unstable,

then stable again at the same value of the non-conservative load. This sequence of stabil-

ity transitions, which has not been reported in the literature, is found to exist for several

non-conservative systems, including articulated linkages with follower end forces and fluid-

conveying pipes. The stability transitions were investigated by applying the Routh-Hurwitz

criterion twice: once for the characteristic polynomial and the second time for the poly-

nomial that guarantees the existence of a second-order auxiliary polynomial in the Routh

array. Nonconservative loading in elastic structures has previously been limited to follower

forces and forces generated by a fluid jet exiting a fluid-conveying pipe. A new type of

non-conservative loading is introduced here in the form of a terminal dynamic moment. The

terminal moment is dynamic in the sense that it is proportional to the slope or curvature

of a point along the structure. Irrespective of whether the moment is slope-dependent or

curvature-dependent, stability is lost in a beam through divergence when the constant of

proportionality is positive, and through flutter when the constant of proportionality is neg-

ative. Some of the theoretical investigations are supported by experiments with a cantilever

beam. In non-conservative systems, the introduction of an intermediate support is known

to result in stability transitions between flutter and divergence. This result was confirmed

for the new type of non-conservative loading, namely, the terminal dynamic moment. For

a cantilever beam with terminal dynamic moment, a rich set of stability transitions were

observed for the first time. This includes multiple stability transitions between divergence

and flutter and between different modes of flutter. In some cases, the transitions are ac-

companied by an abrupt change in the critical load, which is commonly referred to as the

“jump phenomenon” in the literature. External flow is known to result in non-conservative

loading and, therefore, stability transitions are investigated in a hinged beam in external

flow. A dynamic moment, proportional to the curvature of the beam at some point along

its length, is applied at the hinge. The combined effect of non-conservative loading due to

both dynamic moment and external flow is investigated with the objective of developing a

mechanism for energy harvesting from flow fields.

Copyright By
MAHMOUD NABIL ABDULLATIF
2020

I dedicate this work to my parents, who have raised me to be the person I am today, and
to my beloved wife, Sarah, for nursing me with love and affection and always being by my

side through this tremendous journey.

iv

ACKNOWLEDGEMENTS

Thanks to Almighty ALLAH for giving me the strength and ability to learn, understand,
and complete this work.

I would like to thank my advisor, Ranjan Mukherjee, for his guidance and support.

I am thankful for Aren Hellum for helping me see new avenues of progress when I got stuck.

A tremendous thank you to my wife, Sarah, for always being by my side: helping me get
through the lows and celebrating the highs with me. I could not have done this without you.

I am also thankful for my wonderful colleagues, Amer Allafi, Sanders Aspelund, Connor Boss,
Sheryl Chau, Krissy Kamensky, and Nilay Kant, for their constant source of camaraderie
and fruitful discussions.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 1

1.2.1
1.2.2

Introduction

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Stability Transitions Induced by Damping . . . . . . . . . . . . . . .
Stability Transitions Induced by Intermediate Support
. . . . . . . .
1.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Research Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2 Divergence and Flutter Modes of Instability . . . . . . . . . .
2.1 Discrete System Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Continuous System Example . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 3 Stability Transitions Induced by Damping . . . . . . . . . . . .
3.1 Two Degree-Of-Freedom System with a Follower Force . . . . . . . . . . . .
3.1.1
Effect of γ1 on the Stability of the System . . . . . . . . . . . . . . .
3.1.2 Effect of γ2 on the Stability of the System . . . . . . . . . . . . . . .
3.2 Two Degree-Of-Freedom Fluid-Conveying Pipe . . . . . . . . . . . . . . . . .
3.2.1
Effect of Relative Damping on the Stability of the System . . . . . .
3.2.2 Effect of Absolute Damping on the Stability of the System: . . . . . .
. . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Three-Link System with a Follower Force . . . . . . . . . . . . . . . .
3.3.2 Three-Link Fluid-Conveying Pipe . . . . . . . . . . . . . . . . . . . .
3.4 Cantilevered Fluid-Conveying Pipe . . . . . . . . . . . . . . . . . . . . . . .

3.3 Three Degree-Of-Freedom Systems

Chapter 4 Non-conservative Effects of Dynamic Terminal Moment

. . .
4.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Galerkin method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

x

1
1
2
2
4
5
6
7

8
8
12

15
15
18
21
28
29
31
32
32
34
37

41
41
43
43
45

4.3

Instability Investigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Dynamic Moment Proportional to Positive Slope
. . . . . . . . . . .
4.3.2 Dynamic Moment Proportional to Negative Slope . . . . . . . . . . .
4.3.3 Dynamic Moment Proportional to Positive and Negative Curvature .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Hardware Description . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 Model Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.3 Flutter Oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.4 Experimental Investigations

5.1.1 Mathematical Model

Chapter 5 Stability Transitions Induced by an Intermediate Support . .
5.1 Cantilever Beam with Terminal Dynamic Moment . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Terminal Dynamic Moment Proportional to Slope . . . . . . . . . . . . . . .
5.2.1 Model Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 Effect of Intermediate Support on Critical Stability . . . . . . . . . .
. . . . . . . . . . . .
5.3.1 Effect of Intermediate Support on Critical Stability . . . . . . . . . .
5.4 Cantilever Beam with a Follower End Load . . . . . . . . . . . . . . . . . . .
5.4.1 Mathematical Model
. . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.3 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.4 Effect of Intermediate Support on Critical Stability . . . . . . . . . .

5.3 Terminal Dynamic Moment Proportional to Curvature

6.1.1 Mathematical Model
6.1.2

Chapter 6 Non-conservative Behavior of Hinged Beam with a Dynamic
Moment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1 Hinged Beam with Dynamic Moment . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
Instability Investigation . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Effect of External Flow on Flutter Instability . . . . . . . . . . . . . . . . .
6.2.1 Mathematical Model
. . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.2 Flutter Instability and Nature of Oscillations . . . . . . . . . . . . . .
6.3 Fluid-Immersed Hinged Beam Affixed to a Rigid Body . . . . . . . . . . . .
6.3.1 Mathematical model
. . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.2 Method of Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.3 Flutter Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46
46
48
51
52
52
54
54

57
57
57
59
60
61
65
67
74
74
75
77
78

81
81
81
83
85
85
87
90
90
93
97

Chapter 7 Conclusion

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.1 Divergence and Flutter Modes of Instability . . . . . . . . . . . . . . . . . . 104
7.2 Stability Transitions Induced by Damping . . . . . . . . . . . . . . . . . . . 104
7.3 Non-conservative Effects of Dynamic Terminal Moment . . . . . . . . . . . . 106
7.4 Stability Transitions Induced by an Intermediate Support . . . . . . . . . . . 107
7.5 Non-conservative Behavior of Hinged Beam with a Dynamic Moment
. . . . 109

vii

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

viii

LIST OF TABLES

Table 4.1:

Natural frequencies of the beam in Fig.4.8: analytical, numerical, and
experimental values. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 4.2:

Analytical, numerical and experimentally observed values of frequency
and proportionality constant at the flutter instability point for the
. . . . . . . . . . . . . . . . . . . . .
experimental setup in Fig.4.8.

Table 5.1:

Critical values of C s for C s > 0 that results in divergence. . . . . . .

Table 5.2:

Critical values of C s for C s < 0 that results in flutter.

. . . . . . . .

Table 6.1:

Assumed values for the non-dimensional parameters used in the anal-
ysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

55

61

61

95

ix

LIST OF FIGURES

Figure 2.1: Two-link fluid conveying pipe with torsional joint stiffnesses.

. . . .

8

Figure 2.2: The locus of the system roots as the non-dimensional fluid velocity
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

varied.

Figure 2.3:

(a) Fluter instability for uf = 5.007. (b) Divergence instability for
ud = 7.434. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.4: Nonlinear simulation of θ for different values of u.

. . . . . . . . . .

Figure 2.5: Cantilevered beam with end loads . . . . . . . . . . . . . . . . . . .

Figure 2.6:

Figure 2.7:

(a) Cantilevered beam with axial end loading (b) variation of the
first two modal frequencies of a cantilevered beam as the axial load-
ing increases (c) locus of the first two modal frequencies illustrating
divergence instability.
. . . . . . . . . . . . . . . . . . . . . . . . . .

(a) Cantilevered beam with follower end loading (b) variation of the
first two modal frequencies of a cantilevered beam as the follower load-
ing increases (c) locus of the first two modal frequencies illustrating
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
flutter instability.

10

10

11

12

14

14

Figure 3.1: A two-link system with a follower force.

. . . . . . . . . . . . . . . .

16

Figure 3.2: Absolute instability curve in γ1-p space for the two-link system in
Fig.3.1 when γ2 = 1.0. For p ∈ [0, p∗), there is only one real root of
γ1. For p ≥ p∗, there are three real roots of γ1 but only one of them
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
is positive.

Figure 3.3: Two-link system in Fig.3.1 with γ1 = 1.0: (a) Absolute instability
curve in γ2-p space showing the S-D-S effect of damping (b) locus
of the roots that determine system stability as γ2 increases for the
particular case of p = 2.08.
. . . . . . . . . . . . . . . . . . . . . . .

Figure 3.4: Contour plot for ∆ = 0 in γ2-p plane for the two-link system with
follower force: regions with positive and negative values of ∆ are
. . . . . . . . . . . . . . . . .
denoted by (+) and (−), respectively.

20

22

24

x

Figure 3.5: Flutter instability curves in γ1-p space for the two-link system with
follower force for: (a) γ2 = 0.02 and p = 1, (b) γ2 = 0.02 and p = 2.05,
(c) γ2 = 0.02 and p = 5. . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.6: The locus of the roots that determine system stability as γ1 increases
for the particular cases: (a) γ2 = 0.02 and p = 1, (b) γ2 = 0.02 and
p = 2.05, (c) γ2 = 0.02 and p = 5.
. . . . . . . . . . . . . . . . . . .

Figure 3.7: Flutter instability curves in p-γ1 space for the two-link system with
. . . . . . . . . . . . . .

follower force and range of γ2 ∈ [0.01, 0.10]

Figure 3.8: Contour plot for ∆ = 0 in γ1-p plane for the two-link system with
follower force: regions with positive and negative values of ∆ are
. . . . . . . . . . . . . . . . .
denoted by (+) and (−), respectively.

Figure 3.9: Flutter instability curves in p-γ2 space for the two-link system with
. . . . . . . . . . . . . .

follower force and range of γ1 ∈ [0.25, 2.00].

Figure 3.10: Flutter instability curves in γ1-γ2 space for the two-link system with
. . . . . . . . . . . . . . . . . . . . . . .

follower force p ∈ [0.5, 2.0].
(a) Relative damping and (b) absolute damping is introduced in
the fluid-conveying pipe considered by Benjamin [1] with equal link
lengths ℓ and equal torsional joint stiffness k. . . . . . . . . . . . . .

Figure 3.11:

Figure 3.12: Contour plot for ∆ = 0 in β-u plane for the two-link fluid-conveying
pipe with relative damping: regions with positive and negative values
. . . . . . . . . . . .
of ∆ are denoted by (+) and (−), respectively.

Figure 3.13: Absolute instability curves in γ-u space for the two-link pipe with
relative damping for different values of β; the S-D-S effect of damping
is seen for β = 0.04 and 0.08. . . . . . . . . . . . . . . . . . . . . . .

Figure 3.14: Absolute instability curves in γ-u space for the two-link pipe with
absolute damping showing D, S-D, and S-D-S effects for (a) β = 1.50,
(b) β = 0.70, and (c) β = 0.10, respectively. For (b) and (c), the
critical flow velocity ucr for γ = 0 lies outside the range of u; these
values are 2.245 for β = 0.70 and 5.007 for β = 0.10. . . . . . . . . .

Figure 3.15: A three-link system with a follower force.

. . . . . . . . . . . . . . .

24

25

26

26

27

27

28

30

30

32

33

Figure 3.16: Absolute instability curves in γ3-p space for the three-link system with
follower force showing D, D-S, S-D-S, and S effects of damping for
(a) γc = 0.01, (b) γc = 2.0, (c) γc = 4.0, and (d) γc = 6.0, respectively. 34

xi

Figure 3.17:

(a) Flutter instability curves and absolute instability curve in γ-u
space for the three-link fluid-conveying pipe in Example 3.6 for β =
0.10 (b) locus of all roots of the system as u increases for the particular
case of γ = 0.2.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.18:

(a) Flutter instability curves in γ-u space for the three-link fluid-
conveying pipe in Example 3.7 for β = 0.10 (b) A magnified view
of the region of intersection of the flutter instability curves is shown
. . . . . . . . . . . . . . .
along with the absolute instability curve.

Figure 3.19: Absolute instability curves in c-u space for the cantilever fluid-conveying
pipe with β = 0.1 for: (a) α = 0.86 (b) α ∈ [0.2, 1.0]. The S-D-S
. . . . . . . . . . . . . . . . .
effect of damping is seen for α = 0.86.

35

36

39

Figure 3.20: Absolute instability curves in cs-u space for the cantilever fluid-conveying

pipe with β = 0.35 showing the D-S effect of damping. . . . . . . . .

40

Figure 4.1: A flexible cantilever beam with a tip mass and a terminal dynamic
moment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

Figure 4.2: Terminal moment is proportional to the positive slope of the beam:
(a) Variation in the first natural frequency with C s (C s ≥ 0) for four
different values of α. Analytical and numerical results are plotted
using solid lines and dashed lines, respectively.
(b) Analytical re-
sults showing variation of the second, third, fourth, and fifth natural
frequencies with C s for α = 0.25. . . . . . . . . . . . . . . . . . . . .

Figure 4.3: Terminal moment is proportional to the negative slope of the beam:
(a) Variation in the first four natural frequencies with C s (C s ≤
0) for α = 0.1 shows the first flutter instability mode. Analytical
and numerical results are plotted using solid lines and dashed lines,
respectively. (b) Analytical results showing the first flutter instability
mode for five different values of α in range 1: [0.02, 0.63], including
the case shown in (a); the critical stability point is shown using the
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
• sign.

Figure 4.4: Terminal moment is proportional to the negative slope of the beam:
Analytical results show the (a) second flutter instability mode for
range 2: [0.64, 0.95], (b) Switching from first to second flutter insta-
bility mode for α = 0.6394 and C
s = −52.5.
In both figures, the
. . . . . . . . . . .
critical stability point is shown using the • sign.

∗

47

48

49

xii

Figure 4.5: Terminal moment is proportional to the negative slope of the beam:
Critical stability curve for α ∈ [0.02, 0.95]. The critical stability curve
(solid line) was obtained numerically using a ten-mode approxima-
tion; the analytical results are shown for discrete values of α using
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
the • sign.

Figure 4.6: Terminal moment is proportional to the positive curvature of the
beam: Variation in the first natural frequency with C c (C c ≥ 0)
for three different values of α. Analytical and numerical results are
plotted using solid lines and dashed lines, respectively. . . . . . . . .

Figure 4.7: Terminal moment is proportional to the negative curvature of the
beam: Critical stability curve for α ∈ [0.01, 0.88]. The critical sta-
bility curve (solid line) was obtained numerically using a ten-mode
approximation; the analytical results are shown for discrete values of
α using the • sign. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 4.8: Cantilever beam with tip mass: the tip mass is comprised of a motor
for application of a dynamic moment and an air bearing to avoid
excitation of torsional modes. . . . . . . . . . . . . . . . . . . . . . .

Figure 4.9: Analytical and numerical results showing the first flutter instability
mode for the Euler-Bernoulli beam corresponding to the experimen-
tal setup in Fig.4.8. The analytical results are shown using a solid
line; the numerical results are shown using a dashed line but are
indistinguishable from the analytical results.
. . . . . . . . . . . . .

Figure 4.10: A comparison of normalized, non-dimensional mode shapes of the
cantilever in Fig.4.8: solid and dashed lines show the mode shapes
. . . . . . . . . .
obtained analytically and from experimental data.

Figure 5.1: A flexible cantilever beam with an intermediate support subjected to
. . . . . . . . . . . . . . . . . . . . . .

a terminal dynamic moment.

Figure 5.2: Critical stability curve for the beam with intermediate support: ¯Cs >
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.25.

Figure 5.3: Critical stability curve for the beam with intermediate support: ¯Cs >
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.50.

Figure 5.4: Critical stability curve for the beam with intermediate support: ¯Cs >
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.75.

Figure 5.5: Critical stability curve for the beam with intermediate support: ¯Cs <
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.25.

50

51

52

52

55

56

58

62

62

63

64

xiii

Figure 5.6: Critical stability curve for the beam with intermediate support: ¯Cs <
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.50.

Figure 5.7: Critical stability curve for the beam with intermediate support: ¯Cs <
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 and κ = 0.75.

Figure 5.8: Terminal moment is proportional to the positive curvature: (a) C cr is
plotted on a fine mesh grid of κ and α values; divergence instability
is shown by dark points on the grid (b) Top view of (a). . . . . . . .

Figure 5.9: Case 1: (a) Critical stability curve - divergence instability at α = 0.1
and flutter instability for α ∈ [0.101, 0.9] (b) Second flutter instability
mode at α = 0.513 showing veering of the locus of the first natural
frequency.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.10:

(a) Case 2: Critical stability curve - divergence instability for κ ∈
(0.47, 0.54] ∪ (0.72, 0.78] and flutter instability for other values of κ
(b) Case 3: Critical stability curve - divergence instability for α ∈
(0.631, 0.642]∪ (0.673, 0.900] and flutter instability for other values of
α. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.11: Case 4: Locus of the first six natural frequencies as C c is increased
from zero for α = 0.25 and (a) κ = 0.249, (b) κ = 0.251. . . . . . . .

Figure 5.12: Case 5: (a) Critical stability curve - divergence instability for α ∈
(0.62, 0.65], first-mode flutter instability for α ∈ (0.65, 0.715] and
second-mode flutter instability for α ∈ (0.715, 0.78] (b) Second-mode
flutter instability at α = 0.716 showing veering of the locus of the
first natural frequency.
. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.13: Terminal moment is proportional to the negative curvature: (a) C cr
is plotted on a fine mesh grid of κ and α values; divergence instability
is shown by dark points on the grid (b) Top view of (a). . . . . . . .

Figure 5.14: A flexible cantilever beam with an intermediate support subjected to
a follower force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5.15:

(a) Variation in the critical follower force σ as the location of the
intermediate support κ is varied; both analytical (solid line) and nu-
merical (dots) results are shown (b) Variation in the critical value of
σ for 7 different values of κ; stability is lost through flutter for κ =
{0.40, 0.45, 0.49, 0.50} and through divergence for κ = {0.51, 0.52, 0.53}
(c) Loci of the first two natural frequencies for κ = {0.5, 0.5001} show-
ing jump and veering phenomena; mode of instability changes from
flutter to divergence and the critical value of σ jumps from 71.8 to
80.6 when κ changes from 0.5 to 0.5001. . . . . . . . . . . . . . . . .

64

65

68

69

70

71

72

73

74

79

xiv

Figure 6.1: A flexible hinged-beam with dynamic moment

. . . . . . . . . . . .

82

Figure 6.2: Critical stability curve for the hinged beam: C versus α for the case
where k = 0.5 and C > 0. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6.3: Critical stability curve for the hinged beam: C versus α for the case
where k = 0.5 and C < 0. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6.4: The hinged-beam in Fig.6.1 immersed in axial flow . . . . . . . . . .

Figure 6.5: Root locus of the first two fundamental frequencies for the hinged
beam without external flow: C < 0 and k = α = 0.5. . . . . . . . . .

Figure 6.6: Argand diagram for k = α = 0.5, U e = 0.7, and β = 0.9; the critical
. . . . . . . . . . . . . . . . . . . . . . . . . .

C is denoted by C cr.

Figure 6.7: Root locus of the second mode complex conjugate frequency for the
hinged beam subjected to external folw as C varies negatively for
k = α = 0.5, U e = 0.7, and β = 0.9.
. . . . . . . . . . . . . . . . . .

Figure 6.8: Fluid submersible tail mechanism affixed to a rigid body.

. . . . . .

Figure 6.9: Free-body diagrams of the tail mechanism and the rigid body. . . . .

Figure 6.10: Critical stability curves in C−ue space: (a) a range of α ∈ (0.10−0.20]
(b) for α = 0.20 the efficiency, η, and the critical frequency, Ωcr, are
denoted as (η, Ωcr) for five different points along the critical stability
curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6.11: Critical stability curves in C−ue space: (a) a range of α ∈ [0.22−0.38]
with an increment of 0.02 (b) for α = 0.22 the efficiency, η, and the
critical frequency, Ωcr, are denoted as (η, Ωcr) for three different points
along the critical stability curve.
. . . . . . . . . . . . . . . . . . . .

84

85

86

88

88

89

90

91

98

99

Figure 6.12: Critical stability curves in C−ue space: (a) critical stability curves are
depicted for eight different values of α ∈ [0.40, 0.70] (b) for α = 0.58
the efficiency, η, and the critical frequency, Ωcr, are denoted as (η, Ωcr)
for three different points along the critical stability curve.

. . . . . . 100

Figure 6.13: Critical stability curves in C−ue space: (a) critical stability curves are
depicted for four different values of α, namely, α = 0.80, 0.84, 0.90, 0.92
(b) for α = 0.90 the efficiency, η, and the critical frequency, Ωcr, are
denoted as (η, Ωcr) for three different points along the critical stabil-
ity curve; The darkness region of the critical stability curve depicts
the region of negative thrust. . . . . . . . . . . . . . . . . . . . . . . 101

xv

Figure 6.14: Region 1: critical stability curves in C−ue space: (a) critical stability
curves are depicted for seven different values of α ∈ [0.18, 0.36] (b) for
α = 0.30 the efficiency, η, and the critical frequency, Ωcr, are denoted
as (η, Ωcr) for three different points along the critical stability curve. 102

Figure 6.15: Region 2: critical stability curves in C − ue space: (a) critical sta-
bility curves are depicted for three different values of α, namely,
α = 0.50, 0.52, 0.58 (b) for α = 0.62 the efficiency, η, and the critical
frequency, Ωcr, are denoted as (η, Ωcr) for three different points along
the critical stability curve.

. . . . . . . . . . . . . . . . . . . . . . . 102

Figure 6.16: Region 3: critical stability curves in C−ue space: (a) critical stability
curves are depicted for five different values of α ∈ [0.66, 0.98] (b) for
α = 0.98 the efficiency, η, and the critical frequency, Ωcr, are denoted
as (η, Ωcr) for three different points along the critical stability curve. 103

xvi

Chapter 1

Introduction

1.1 Background

The dynamics and stability characteristics of elastic structures have been studied extensively

since the early work by Prandtl [2] and Reut [3]. Depending on the nature of the boundary

loading, elastic systems exhibit divergence and flutter modes of instability. In the absence

of damping, a loss of stability through divergence occurs when the lowest frequency of a

system passes through the origin and the property of the equilibrium shifts from stable to

unstable.

In contrast, flutter instability occurs when two successive frequencies approach

each other, coexist, and become complex conjugates. In many instances, the system reaches

a limit cycle post-flutter and converts the gained energy into a steady periodic oscillation.

According to Herrmann and Bungay [4], a flexible system, in the absence of damping, loses

its stability solely through divergence when the boundary loading is conservative in nature.

On the other hand, the same system may lose stability through divergence or flutter when

subjected to a non-conservative boundary load. The critical load for divergence can always

be determined through the static Euler method; on the contrary, the critical load for flutter

can be determined through the kinetic method alone.

A non-conservative boundary load can be applied to both discrete and continuous sys-

tems; the common examples of the discrete system are Ziegler’s double pendulum [5] and

Benjamin’s articulated fluid conveying pipe [6]. The most common example of the contin-

uous system is a cantilevered column subjected to a follower loading. Several variations of

the column problem are discussed in the literature: a Beck’s column in which a cantilever

is subjected to tangential follower boundary loading; a Pfluger’s column, which is a Beck’s

column with a concentrated end mass; a Leipholz’s column, which is a beam subjected to

uniformly distributed follower loading along the column span; and finally, a Hauger’s column,

1

where a linearly increasing distributed follower loading is subjected to a cantilevered beam.

A survey on the dynamic stability of these columns can be found in [7]. Another typical

example of a continuous non-conservative system is a fluid conveying pipe, which is exten-

sively reviewed by Paidoussis and Li [8]. Several techniques were addressed in the literature

to tackle the stability problems. Early work [3,9,10] focused on the static instability aspects

of a column under follower loading. More recently, Kounadis [11] utilized nonlinear dynamic

analysis to reexamine the mechanism of stability on the problem considered by Herrmann

and Bungay [4]. Qiu and Nemate-Nasser [12] reported experimental results on instability

modes of an articulated cantilever subjected an impinging air-jet at its free end. A compre-

hensive discussion of the literature on non-conservative systems, including both theoretical

and experimental investigations, can be found in the review paper by Elishakoff [13].

1.2 Literature Review

In this work we investigate stability transitions in structures due to damping in the presence

of non-conservative loading, due to terminal dynamic moment with and without intermediate

support, and due to external flow with a dynamic moment. Since the terminal dynamic

moment is introduced here for the first time, we review the literature on stability transitions

induced by damping and intermediate support.

1.2.1 Stability Transitions Induced by Damping

A counter-intuitive property of some non-conservative systems is that the application of

damping can render the system unstable. This property is frequently referred to as Ziegler’s

paradox [14] after a classical finding that the addition of viscous damping to the joints of a

double pendulum subjected to a follower end force destabilizes the system. Since Ziegler’s

publication in 1952, investigations of the destabilizing effect of damping have continued to

the present day, [15–20]. The effect of both internal (structural) and external (viscous)

damping has been investigated in systems subjected to follower end forces; this includes

2

fluid-conveying pipes where the conveyed fluid contributes to an additional damping-like

velocity-dependent force in addition to the follower end-force.

Research interest in the general phenomenon of destabilization by damping has been ro-

bust. Bottema [21] investigated the conditions for damping to induce instability in Ziegler’s

work [14] and Bolotin [22] provided some general results for the stability of two-DOF non-

conservative systems. Herrmann and Jong [23] established that a relationship exists be-

tween the damped and undamped critical loads of a system. Nemat-Nasser and coworkers

provided extensions to multi-dof [24] and continuous [25] systems and destabilization by

damping in a continuous system were observed experimentally by Gregory and Paidoussis in

fluid-conveying pipes [26]. Twenty-one years after Ziegler’s finding, the destabilization-by-

damping phenomenon was known to be sufficiently common that Done [27] felt compelled

to demonstrate mathematically that there exist damping configurations which are never

destabilizing. This result was supported by observations published earlier [22], [28], and

corroborated in other non-conservative systems such as fluid-conveying pipes with different

boundary conditions [29], curved pipes [30], tapered columns [31], panel flutter [32], and

fluid-conveying shells [33]. Although the destabilizing effect implies divergent behavior of

trajectories in the neighborhood of the equilibrium; Kounadis [34] used nonlinear analysis

to show that the trajectories, despite diverging from the equilibrium, may remain bounded.

Several researchers have provided a physical explanation of the destabilizing effect of

damping. For example, Semler et al. [15] used energy methods to investigate discrete systems

with viscous damping and continuous systems with Kelvin-Voigt (structural) damping. For

discrete systems, however, it was shown that the phase difference plays a vital role in the

determination of the mode that extracts energy from the follower force and the magnitude

of the critical load. For small damping, Gallina [16] used eigenvalue analysis to show that

a positive derivative of an eigenvalue with respect to the damping coefficient results in

instability. Sugiyama and Langthjem [17] used energy methods to investigate Beck’s column

and concluded that the follower force could only do work in each cycle if the column vibration

3

has a traveling wave component; additionally, the wavenumbers at the point of application of

the force determines the rate of energy growth. Kirillov and Seyranian [35] used perturbation

methods to derive necessary and sufficient conditions for the matrix of velocity-dependent

forces to guarantee asymptotic stability.

1.2.2 Stability Transitions Induced by Intermediate Support

The transition between flutter and divergence in non-conservative systems was found to

occur due to the introduction of an intermediate support.

In some cases, the transition

is accompanied by an abrupt variation in the critical non-conservative loading, commonly

referred to as the “jump phenomenon”. This phenomenon was reported earlier by Zorii and

Chernukha [36] and followed by a series of papers by Kounadis [37–39] in which regions

in the parameter space related to divergence and flutter were identified. In related work,

Kounadis and Economou [40] spotted that the transition from flutter to divergence and vice

versa is associated with the jump phenomenon. Elishakoff and Hollkamp [41] equally report

this unforeseen variation in the critical non-conservative loading.

A later work by Elishakoff and Lottati [42], based on an exact solution for a simply

supported and clamped columns resting on intermediate support, reported that the jump

phenomenon in the earlier work could not be observed. Sundararajan [11], on the other hand,

noticed that the addition of elastic support might destabilize the undamped continuous

non-conservative system. There has been a surge of interest in examining the effect of

intermediate support, taking into account some of the previous stability investigations on

various aspects of the flexible systems. Examples are two degree-of-freedom flexible systems

subjected to follower loading [43], flexible fluid conveying pipes [29, 44], and rotating shafts

with axial loading [45]. Although most of the studies dealt with undamped systems, where

stability investigations were carried out using both analytical [42, 46–48] and numerical [41,

44, 49] methods, studies involving damping have also been delineated [50, 51]. The influence

of intermediate support on the dynamic stability remains a fascinating topic for researchers,

4

even to this day, for further exploration. Some of the recent studies include: the transverse

vibrations for an axially moving beam supported by rollers [52], the influence of elastic

end support on flutter and buckling instability of the Beck’s column in the presence of an

arbitrary number of weak sections [53], and buckling behavior of beams and columns with

additional intermediate flexible supports [54].

1.3 Motivation

The pioneering works on structural stability of non-conservative systems was the motivation

behind this work. Among these works, noteworthy is the counterintuitive role of damping,

whereby the addition of damping causes destabilization of the structure. A non-conservative

system can exhibit multiple stability transitions [15], such that increasing the damping for

a given non-conservative force can stabilize an unstable system, which then becomes un-

stable again as damping is further increased. There has been little or no efforts on in-

vestigation of multiple stability transitions and this may be attributed to several factors.

Many efforts in the field consider the damping levels to be too small to produce the phe-

nomenon [17, 18, 23, 24, 35, 55]. The appearance of additional stability transitions is also

somewhat subtle, depending on the method of analysis. In at least one case [29], the phe-

nomenon was overlooked by the authors, or was simply outside the scope of their interest.

This relative lack of investigation stimulated the first part of this work where a procedure is

developed in the context of examples, to study multiple stability transitions. Three damping-

induced stability transitions has been observed, a phenomenon that has not been reported

heretofore. The majority of the results discussed in the literature also pertain to the com-

monly known non-conservative loads, namely, follower forces and forces induced by the fluid

jet in fluid-conveying pipes. A new type of non-conservative loading is introduced in beam-

like structures in the form of a terminal dynamic moment. The moment is proportional to

the slope or the curvature of a point along the beam length. Thereupon, the stability char-

acteristics of a cantilever beam subjected to the dynamic moment are also investigated. An

5

extension of this work, inspired by the work done by Lee [56], include the role of intermediate

support on structural stability. This investigation reveals a rich set of stability transitions

between divergence and flutter and different modes of flutter that have not been reported

earlier.

1.4 Research Contributions

This work provides additional insight into the role of system parameters on the stability

characteristics of structures. The study begins with the investigation of the role of damping

and additional stability transitions are shown to be both possible and common. A non-

conservative system that is unstable for zero damping is stabilized, destabilized, and resta-

bilized as the damping is increased. The observation of this phenomenon is made possible

by utilizing the Routh-Hurwitz procedure [57] twice: once for the characteristic polynomial

and a second time for the auxiliary polynomial. Although other works have employed the

Routh procedure to assess the stability of the system [15], it has not ben used to determine

the number of marginal stability points.

Stability characteristics of a cantilevered beam with a tip mass subjected to a terminal

non-conservative moment was also studied. The magnitude of the bending moment varies

with the slope or the curvature of the beam at some point along its length. To the best of

our knowledge, the stability analysis of such systems has not appeared in the literature, and

this work is the first to show both theoretical and experimental investigations with this new

type of non-conservative load.

Flutter instability, which is the main topic of investigation in this work, is explored due

to its ability to produce a sustained oscillation. Such oscillations are ideal for engineering

applications like underwater propulsion [58,59] and energy harvesting [60]. With this in mind,

the instability effect of the non-conservative dynamic moment is investigated in the presence

of external flow. Galerkin approximation and analytical solutions are used to determine the

range of parameters for which flutter oscillations can be obtained.

6

1.5 Thesis Outline

This work is organized as follows: Chapter 2 provides a short discussion on the two types of

structural instabilities using examples of discrete and continuous systems. Chapter 3 provides

new results on the role of damping on system stability; the stability transitions of different

non-conservative systems are examined through examples. Chapter 4, explores a new way

to destabilize a cantilever beam by using a non-conservative dynamic moment. Chapter

5 investigates the role of the intermediate support on stability of a structure subjected to

the dynamic moment. In Chapter 6 the instability of a hinged beam subjected to a non-

conservative moment is examined in the presence of external flow Finally, concluding remarks

are provided in Chapter 7.

7

Chapter 2

Divergence and Flutter Modes of Instability

Elastic structures can lose stability through divergence or flutter. In the absence of damping,

a loss of stability through divergence occurs when the lowest frequency of the system passes

through the origin and the property of the equilibrium shifts from stable to unstable. In

contrast, flutter instability occurs when two successive frequencies approach each other,

coexist, and become complex conjugates. While the flutter mode of instability can only

occur for non-conservative loading, the divergence mode of instability can occur for both

conservative and non-conservative loadings. In general, both discrete and continuous systems

can lose stability through flutter or divergence and the load at which stability is lost is known

as the critical load. In this chapter we consider two examples: an articulated fluid conveying

pipe and a cantilever beam subjected to terminal loading, to illustrate the both types of

instabilities in a discrete system and a continuous system, respectively.

2.1 Discrete System Example

Y

U

ℓ2

k2

φ

ℓ1

θ

k1

X

Figure 2.1: Two-link fluid conveying pipe with torsional joint stiffnesses.

8

Consider the two-link fluid-conveying tube shown in Fig.2.1. The non-dimensional equa-

tions of motion in absence of gravity force, following Benjamin’s derivation [6], can be ex-

pressed as

where

−u2αβ sin(θ − φ) + (1 + κ)θ − φ + uα2β ˙θ + 2uαβ cos(θ − φ) ˙φ

1
(3α sin(θ − φ) ˙φ2 + 2α2(3 + α)¨θ + 3α cos(θ − φ) ¨φ) = 0
2
α cos(θ − φ)¨θ −
α sin(θ − φ) ˙φ2 + ¨φ + βu ˙φ + φ − θ = 0

3
2

+

3
2

(2.1)

β , 3M

(M + m)

, α , ℓ1
ℓ2

,

κ , k1
k2

Using small angle assumption, the nonlinear equation of motion can be linearized to

(α3 + 3α2)¨θ +

3
2

α ¨φ + βu(α2 ˙θ + 2α ˙φ) − αβu2(θ − φ) + (1 + κ)θ − φ = 0
α¨θ + ¨φ + βu ˙φ + φ − θ = 0

3
2

In matrix form, Eq.(2.2) can be written as



(α3 + 3α2)

3
2α

3
2α

1





¨θ

¨φ

 +

βuα2 2βuα

0

βu



˙θ

˙φ



 +

(1 + κ) − αβu2 αβu2 − 1

−1

1





θ

φ

 =

0

0

(2.2)

 (2.3)

Assuming k1 = k2 = k ⇒ κ = 1, ℓ1 = ℓ2 = ℓ ⇒ α = 1, β = 0.1, and starting from u = 0
results in purely imaginary roots - see Fig.2.2. As the non-dimensional fluid velocity u is

varied from zero, the purely imaginary roots turn-into complex conjugate roots due to the

damping effect that arises from internal fluid flow. By increasing u further, the second mode

turns back towards the imaginary axis while the first mode keeps moving away from the

axis. At uf = 5.007, the system loses its stability as the second mode crosses the imaginary

axis. Since this instability is associated with damping, it is commonly known by single-

9

mode flutter [61]. At the same time, increasing u further render a second loss of stability at

ud = 7.434 where the system exhibits divergence type of instabilities as shown in Fig.2.2.

Im(ω)

2

1.5

1

0.5

uf = 5.007

ud = 7.434

ud = 7.434

-3

-2

-1

1

2

Re(ω)

u = 0.000

-0.5

-1

-1.5

-2

uf = 5.007

Figure 2.2: The locus of the system roots as the non-dimensional fluid velocity varied.

The linear system in Eq.(2.3) is simulated for an arbitrary initial conditions. The plot of θ

is shown in Fig.2.3(a) for uf = 5.007 and in Fig.2.3(b) for ud = 7.434.

15.0

θ

0.0

-15.0

0.0

8000

θ

0.0

40.0

time
(a)

100.0

0.0

20.0

time
(b)

40.0

Figure 2.3: (a) Fluter instability for uf = 5.007. (b) Divergence instability for ud = 7.434.

According to Kounadis [11], the nonlinear system will not destabilize globally as the

critical parameter exceeds its critical value. This can be seen by simulating the nonlinear

10

equation of motion in Eq.(2.1). The results, shown in Fig.2.4, indicate that the system

exhibits a limit cycle at critical stability point as well as post flutter. The system might

exhibit instability at much higher values of u, which is not shown here.

u < uf

u = uf

time

time

u > uf

time

Figure 2.4: Nonlinear simulation of θ for different values of u.

The results presented in Figs.2.2-2.3 reflects two different dynamic stability perspectives,

namely, local (linear) and global (nonlinear). It should be noted that the subsequent chapters

are mainly going to focus on the local dynamic stability behaviors of slender structures.

Ultimately, the main objective is to foretell the flutter instability points of these structures

by estimating the critical values of the system parameters.

11

2.2 Continuous System Example

v

v(u, τ )

σ

η

Figure 2.5: Cantilevered beam with end loads

u

The flutter and divergence instabilities in continuous system can be investigated using a

straight cantilever beam subjected to both conservative and non-conservative loadings. These

forces are an axial end load η and follower force σ applied to beam free-end, respectively, as

shown in Fig.2.5. The motion of the beam is assumed to be planar and the transverse non-

dimensional displacement of the beam is denoted by v(u, τ ). Assuming small deformations,

the non-dimensional equation of motion can be expressed as [62]:

v ′′′′(u, τ ) + (σ + η)v ′′(u, τ ) + ¨v(u, τ ) = 0

(2.4)

where the geometric and natural boundary conditions correspond to cantilevered beam are

v(0, t) = 0,

v ′(0, τ ) = 0,

v ′′(1, τ ) = 0,

v ′′′(1, τ ) + ηv ′(1, τ ) = 0

(2.5)

The dynamic instabilities of the system can be examined through numerical analysis, in par-

ticular, using the Galerkin method [63], where the equation of motion Eq.(2.4) is multiplied

by a weight function and integrated over the length of the beam. Assuming the solution to

be of the form

v(u, τ ) =

NXi=1

ai(τ ) φi(u)

(2.6)

where N is the number of terms in the approximation, φi(u), i = 1,· · · , N, are functions

12

that satisfy the boundary conditions in Eq.(2.5), and ai(τ ), i = 1,· · · , N, are coefficients
that will be determined. Substituting Eq.(2.6) into Eq.(2.4) and integrating over the length,

yields

Ka + (σ + η)Ga + M¨a = 0

where the (i, j)-th elements of the K, G and M matrices are defined as

j (1) +Z 1
j(1) −Z 1

0

0

φi(u)φj(u) du

Kij , −ηφi(1)φ′′′
Gij , φi(1)φ′

Mij , Z 1

0

φ′′
i (u)φ′′

j (u) du

φ′
i(u)φ′

j(u) du

The eigenvalue problem can be performed for a given σ or η, where the first two natural

frequencies can be depicted a function of the unknown loading. Figure 2.6(b) shows the

first two natural frequencies as a function of the axial end load η for zero follower force,

i.e. σ = 0. The result reveals a decreasing nature in the natural frequencies as η increases,

and estimates the end loading at which the beam buckles to be ηcr = 2.48. Figure 2.7(b)

shows the same first two natural frequencies as a function of the follower force σ for η = 0.

Note that as σ increases, the first modal frequency increases while the second modal nat-

ural frequency decreases. These two frequencies merge, resulting in flutter instability, at

a critical value of the follower force given by σcr = 20.05. To demonstrate this, the root

loci of the system are depicted in Fig.2.6(c) and Fig.2.7(c), respectively. A loss of stability

through divergence occurs when the lowest natural frequency of an undamped system passes

through the origin as the equilibrium of the system shifts from stable to unstable, as shown in

Fig.2.6(c). In contrast, in the absence of damping, flutter instability occurs when two succes-

sive natural frequencies approach each other, coexist, and then become complex conjugate

– see Fig.2.7(c). According to Paidoussis [61], this instability is known by coupled-mode

fluttering via Hamiltonian-Hopf bifurcation.

13

v

v

v(u, τ )

1
(a)

u

η

η

ω

30

20

10

ω1

ω2

ηcr = 2.48

-15

-10

-5

5

10

15

η

Im(ω)

x

iω2(η = 0)

(b)

ηcr = 2.48

root locus

x

iω1(η = 0)

x

iω1(η = 0)

x

iω2(η = 0)

(c)

Re(ω)

Figure 2.6: (a) Cantilevered beam with axial end loading (b) variation of the first two modal
frequencies of a cantilevered beam as the axial loading increases (c) locus of the first two
modal frequencies illustrating divergence instability.

v

v

v(u, τ )

1
(a)

u

σ

σ

ω

30

20

10

ω2

ω1

σcr = 20.05

-15

-10

-5

5

10

15

σ

Im(ω)

x

iω2(σ = 0)

(b)

σcr = 20.05

x

iω1(σ = 0)

x

iω1(σ = 0)

root locus

x

iω2(σ = 0)

(c)

Re(ω)

Figure 2.7: (a) Cantilevered beam with follower end loading (b) variation of the first two
modal frequencies of a cantilevered beam as the follower loading increases (c) locus of the
first two modal frequencies illustrating flutter instability.

14

Chapter 3

Stability Transitions Induced by Damping

Chapter 2 introduced the two types of instabilities in elastic structures, namely, divergence

and flutter. Damping plays an important role in flutter instability as we have seen - double

mode flutter - in the absence of damping is replaced by single mode flutter in the pres-

ence of damping. A counter-intuitive property of damping is that it can render some non-

conservative systems unstable. This property is frequently referred to as Ziegler’s paradox

after the classical finding that the addition of damping destabilizes a stable double pendulum

subjected to a follower end force. Since this finding, investigations of the phenomenon have

continued to the present day. For example, It has been shown that some non-conservative

systems may exhibit up to two stability transitions when the level of damping is increased

while the non-conservative load is healed fixed. This chapter focuses on additional stabil-

ity transitions associated with damping; it is shown that a system which is unstable for

zero damping can be stabilized, then destabilized, and finally restabilized as the level of

damping is increased. This action is illustrated with the help of examples which show var-

ious taxonomies of damping-induced stabilization and destabilization behavior in various

non-conservative systems.

3.1 Two Degree-Of-Freedom System with a Follower Force

Consider the two-link system shown in Fig.3.1, which is subjected to the follower force P .

Constrained to move in the horizontal plane, the two links are comprised of point masses m1

and m2, connected by massless rods of equal length ℓ. The rotational joints have torsional

stiffness k1 and k2 and their damping coefficients are c1 and c2. The stability characteristics

15

of the system can be investigated from the linearized equations of motion1:

(m1 + m2)ℓ2 ¨θ1 + m2ℓ2 ¨θ2 + (c1 + c2) ˙θ1 − c2

˙θ2 + (k1 + k2 − P ℓ)θ1 − (k2 − P ℓ)θ2 = 0
m2ℓ2 ¨θ1 + m2ℓ2 ¨θ2 − c2 ˙θ1 + c2 ˙θ2 − k2θ1 + k2θ2 = 0

(3.1)

Y

P

m2

ℓ

k2, c2

m1

θ2

ℓ

θ1

k1, c1

X

Figure 3.1: A two-link system with a follower force.

For a comparison with the results in [15], we assume m1 = 2m2 = 2m, k1 = k2 = k, and

scale the time, damping and force variables as follows:

to obtain the non-dimensional equations [15]:

τ =pk/mℓ2 t,




1 1

3 1

¨θ1

¨θ2

γ1 =p1/kmℓ2 c1,

 +

γ1 + γ2 −γ2
γ2
−γ2

γ2 =p1/kmℓ2 c2,
 +

2 − p p − 1
−1

1



p = (ℓ/k) P

(3.2)



θ1

θ2

 =

0

0



(3.3)

˙θ1

˙θ2



In Eq.(3.3), the (˙) and (¨) denote the first and second derivatives with respect to the non-

1The equations of motion presented here are a special case of the equations presented in [64].

16

dimensional time variable τ . Assuming the solution to be of the form



θ1

θ2

we get the characteristic polynomial

 =

Θ1

Θ2

 esτ

Det

3s2 + (γ1 + γ2)s + 2 − p s2 − γ2s + p − 1
s2 + γ2s + 1

s2 − γ2s − 1

 = 0

⇒ 2s4 + (γ1 + 6γ2)s3 + (7 + γ1γ2 − 2p)s2 + (γ1 + γ2)s + 1 = 0

The Routh array [57] can be constructed as follows:

s4

s3

s2

s1

s0

2

(7 + γ1γ2 − 2p) 1

(γ1 + 6γ2)

(γ1 + γ2)

1

ǫ21

ǫ11

1

where ǫ21 and ǫ11 are given by the expressions

ǫ21 =

ǫ11 =

(γ1 + 6γ2)(7 + γ1γ2 − 2p) − 2(γ1 + γ2)
(γ1 + 6γ2)(7 + γ1γ2 − 2p)(γ1 + γ2) − 2(γ1 + γ2)2 − (γ1 + 6γ2)2

γ1 + 6γ2

(γ1 + 6γ2)(7 + γ1γ2 − 2p) − 2(γ1 + γ2)

(3.4)

(3.5)

Flutter instability is characterized by a pair of complex conjugate roots with zero real part.

The characteristic polynomial in Eq.(3.4) will consequently have a second-order auxilliary

polynomial [57] as a factor. This corresponds to the third case of the Routh-Hurwitz criterion

17

[57], and necessitates

ǫ11 = 0

⇒ (γ1 + 6γ2)(7 + γ1γ2 − 2p)(γ1 + γ2) − 2(γ1 + γ2)2 − (γ1 + 6γ2)2 = 0

(3.6)

To ensure that the roots are purely imaginary, we additionally need ǫ21 > 0. The critical

frequency can then be obtained as follows

(ǫ21s2 + 1) = 0 ⇒ ωcr =r 1

ǫ21

=s

(γ1 + 6γ2)

(γ1 + 6γ2)(7 + γ1γ2 − 2p) − 2(γ1 + γ2)

(3.7)

Remark 3.1: An auxilliary polynomial, by definition, is symmetric with respect to both

the real and imaginary axes [57]. This implies that a second-order auxilliary polynomial

must be comprised of either a pair of roots on the imaginary axis located symmetrically with

respect to the real axis, or a pair of roots on the real axis located symmetrically with respect

to the imaginary axis. The first case is associated with the onset of flutter instability and

will be the subject of investigation of this paper. The second case will not be pursued further

as it corresponds to a system that is already unstable in divergence, and not a system at the

onset of instability.

We now investigate two separate cases that are discussed below.

3.1.1

Effect of γ1 on the Stability of the System

To this end, we rewrite the left-hand side of Eq.(3.6) as a polynomial in γ1:

α1 γ3

1 + α2 γ2

1 + α3 γ1 + α4 = 0

α1 , γ2, α2 , 7γ2

2 + 4 − 2p, α3 , γ2(6γ2

2 + 33 − 14p), α4 , 4γ2

2(1 − 3p)

(3.8)

(3.9)

The coefficients of the polynomial, αi, i = 1, 2, 3, 4, are functions of p and γ2; here, both of

18

them are assumed to be positive. The discriminant of the cubic equation in Eq.(3.8), namely

∆ = 18α1α2α3α4 − 4α3

2α4 + α2

2α2

3 − 4α1α3

3 − 27α2

1α2
4

(3.10)

determine the number and types of roots of γ1 [65]. In particular, Eq.(3.8) has

• three distinct real roots if ∆ > 0.

• a multiple root and all its roots are real if ∆ = 0.

• has one real root and two complex conjugate roots if ∆ < 0.

For a real root to be an admissible solution of γ1, it has to be positive. On the other hand, the

number of roots with positive real parts will be equal to the number of sign changes [57] in

the first column of the Routh array constructed using the coefficients of the cubic polynomial

in Eq.(3.8), namely

γ3
1

γ2
1

γ1
1

γ0
1

γ2

(7γ2

2 + 4 − 2p)

γ2(6γ2

2 + 33 − 14p)
4γ2
2(1 − 3p)

µ11

4γ2

2(1 − 3p)

µ11 , γ2(6γ2

2 + 33 − 14p) −

4γ3
2(1 − 3p)
(7γ2
2 + 4 − 2p)

where

(3.11)

The number of positive real roots can therefore be deduced from the sign of ∆ in Eq.(3.10)

and the number of sign changes in the first column of the Routh array in Eq.(3.11). This is

illustrated with the help of a numerical example.

Example 3.1: Consider the case where γ2 = 1.0. From Eqs.(3.9) and (3.10), we have

α1 = 1, α2 = 11 − 2p, α3 = 39 − 14p, α4 = 4(1 − 3p)

⇒ ∆ = 25(16p4 − 64p3 − 104p2 + 1008p − 1763)

19

When p is assumed positive, it can be shown ∆ = 0 for p = p∗ ≈ 3.245. For p = p∗, we
can use Eq.(3.8) to compute γ1 = γ ∗
1 ≈ 2.69. For p ∈ [0, p∗), ∆ < 0 and therefore there is
only one real root of γ1. Using the first column of the Routh array in Eq.(3.11), it can be
shown that there will be only one sign change for p ∈ (0.¯3, p∗). This implies that the critical
load pcr = 0.¯3 for γ1 = 0 and it increases to pcr = p∗ as γ1 increases to γ ∗
1. For p ≥ p∗,
∆ ≥ 0 and therefore there are three real roots of γ1. However, only one of these three roots
will be positive; this can be concluded from the single sign change in the first column of the

Routh array in Eq.(3.11). These results, which can be verified from the plot of pcr in Fig.3.2,

implies that γ1 has a purely stabilizing (S) effect when γ2 = 1.0. The positive sign of ǫ21

in Eq.(3.5) can be used to ascertain that the nature of the instability is flutter, the critical

frequency can be determined using Eq.(3.7).

1
γ

6

4

2

0

0

stable

γ ∗
1

≈ 2.69

unstable

0.¯3

p∗ ≈ 3.245

2

4

6

p

Figure 3.2: Absolute instability curve in γ1-p space for the two-link system in Fig.3.1 when
γ2 = 1.0. For p ∈ [0, p∗), there is only one real root of γ1. For p ≥ p∗, there are three real
roots of γ1 but only one of them is positive.

We now introduce the following definitions:

Flutter instability curve: The flutter instability curve is the locus of points in the pa-

rameter space for which a pair of complex conjugate roots maintain their position on the

imaginary axis.

20

Absolute instability curve: The absolute instability curve is the locus of points in the

parameter space for which a pair of complex conjugate roots lie on the imaginary axis while

all other roots lie in the open left-half plane.

Remark 3.2: The curve in Fig.3.2 is the flutter instability curve corresponding to the first

mode of the system in the γ1-p space; it is also the absolute instability curve.

3.1.2 Effect of γ2 on the Stability of the System

To investigate the effect of γ2 on the stability of the system, we rewrite the left-hand side of

Eq.(3.6) as a polynomial in γ2

β1 γ3

2 + β2 γ2

2 + β3 γ2 + β4 = 0

β1 , 6γ1,

β2 , 7γ2

1 + 4(1 − 3p),

β3 , γ1(γ2

1 + 33 − 14p),

β4 , γ2

1(4 − 2p)

and consider the specific example discussed below.

Example 3.2: Consider the case where γ1 = 1.0. Using Eqs.(3.13) we get

(3.12)

(3.13)

β1 = 6,

β2 = 11 − 12p,

β3 = 34 − 14p,

β4 = 2(2 − p)

The discriminant of the cubic polynomial in Eq.(3.12) and the Routh array are as follows

∆ = 100(144p4 − 936p3 + 409p2 + 5172p − 6787)

γ3
2

γ2
2

γ1
2

γ0
2

6

(11 − 12p)

(34 − 14p)
2(2 − p)

ν11

2(2 − p)

ν11 ,

2(84p2 − 275p + 175)

(11 − 12p)

When p is assumed positive, it can be shown ∆ = 0 for p = {p∗

1, p∗

2, p∗

3} ≈ {1.964, 2.144, 4.75}.

21

1 and p∗

1, p∗

For p = p∗

2, we can use Eq.(3.12) to compute the γ2 values to be ≈ 1.052 and ≈ 0.149.
2] where ∆ ≥ 0; this implies that all three
We focus our attention on the domain p ∈ [p∗
roots of γ2 are real. Within this domain, it can be shown that there will be two sign changes in

1, 2) and three sign changes for p ∈ (2, p∗
the first column of the Routh array above for p ∈ [p∗
2].
These results can be verified from the plot of pcr vs γ2 in Fig.3.3 (a). Based on the notion

introduced in [15], γ2 has a stabilizing effect in the range [0, 0.149]; a destabilizing effect in the

range [0.149, 1.052]; and finally a stabilizing effect for γ2 > 1.052. This anomalous property,

which is attributed to three positive real roots of γ2, has not been reported in the literature to

the best of our knowledge. The locus of the roots that determine the stability of the system

are plotted in Fig.3.3 (b) for the particular case of p = 2.08. It can be seen that the system is

unstable for γ2 ∈ (0.00, 0.036)∪ (0.382, 1.907) and stable for γ2 ∈ (0.036, 0.382)∪ (1.907,∞).
Similar to Example 3.1, the positive sign of ǫ21 in Eq.(3.5) can be used to ascertain that the

nature of the instability is flutter, the critical frequency can be determined using Eq.(3.7).

The curve in Fig.3.3 is the flutter instability curve corresponding to the second mode of the

system in the γ2-p space; it is also the absolute instability curve.

Remark 3.3: Example 3.1 showed that γ1 has a purely stabilizing effect and Example 3.2

2

stable

γ2 = 1.907

γ2 ≈ 1.052

p∗

1 = 1.964

Im(s)
1.0

γ2 = 0.0

γ2 = 0.036

γ2 = 1.907

γ2 = 0.382

2
γ

1

unstable

γ2 ≈ 0.149

-0.04

γ2 = 1.907

γ2 = 0.382

p∗
2

≈ 2.144

stable

γ2 = 0.036

0

1.96

2.00

2.08

2.16

p

(a)

-1.0

(b)

Re(s)

0.04

γ2 = 0.382

γ2 = 0.036

γ2 = 0.0

Figure 3.3: Two-link system in Fig.3.1 with γ1 = 1.0: (a) Absolute instability curve in γ2-p
space showing the S-D-S effect of damping (b) locus of the roots that determine system
stability as γ2 increases for the particular case of p = 2.08.

22

showed that γ2 has a heretofore unobserved “stabilizing-destabilizing-stabilizing” (S-D-S)

effect. It should be mentioned that γ1 can also exhibit the stabilizing-destabilizing-stabilizing

(S-D-S) effect for specific values of γ2, and γ2 can have a purely stabilizing effect for specific

values of γ1.

Remark 3.4: In Example 4.2, the value of γ1 was assumed to be equal to 1.0. For γ1 ∈
{0.02, 0.05, 0.10}, γ2 exhibits the same stabilizing-destabilizing-stabilizing (S-D-S) effect. For
all three cases, there exist a range of p values for which γ2 has three positive real roots.

However, two of the three roots are very close to one another and the third root is distantly

placed. For example, for γ1 = 0.1 and p = 2.05, the roots of γ2 are 0.0026, 0.0183 and 34.1956.

The presence of the third root was overlooked in [15] for all three cases; consequently, the

authors reached the conclusion that γ2 has a stabilizing-destabilizing (S-D) effect as opposed

to the stabilizing-destabilizing-stabilizing (S-D-S) effect.

Remark 3.5: The stabilizing and destabilizing effect of damping was defined in a broader

sense in [15]. Based on this definition, the damping is said to be stabilizing when the slope of

the flutter instability curve is positive, it is said to be destabilizing when the slope is negative.

It will be shown later with the help of Example 3.7 that this definition is not suitable when

flutter instability curves of different modes intersect each other and the absolute stability

curve is comprised of segments of multiple flutter instability curves.

Overall, the role of damping in the joints on the system stability can be investigated using

the discriminant of the cubic equations in Eq.(3.8) and Eq.(3.12), respectively. For instance,

the contour plot of ∆, Eq.(3.10), in γ2-p plane is depicted in Fig.3.4 to show the positive

and the negative regions of ∆. A negative value of ∆ implies that there can be at most one

real root, and therefore at most one admissible (positive and real) solution; this implies that

γ1 will certainly not exhibit the stabilizing-destabilizing-stabilizing effect. A positive value

of ∆ implies that there exists three real roots of γ1. However, all of these roots may not

be admissible solutions; therefore, γ1 may or may not exhibit the stabilizing-destabilizing-

stabilizing effect.

23

6

p

6

p

(+)

(+)

0

0

(−)

γ2

2

0

0

(0.02, 5.00)

(+)

(−)

(0.02, 2.05)

(0.02, 1.00)

(+)

γ2

0.1

Figure 3.4: Contour plot for ∆ = 0 in γ2-p plane for the two-link system with follower force:
regions with positive and negative values of ∆ are denoted by (+) and (−), respectively.

It can be seen from Fig.3.4 that ∆ is positive for γ2 = 0.02 and p = 1.00, 2.05, and 5.00,

respectively. Therefore, it may be possible to obtain three admissible solutions for γ1.

0.04

stable

1
γ

9.0

400

stable

stable

1
γ

unstable

1
γ

(1, 0.008)

unstable

)
1
1
.
0
,
5
0
.
2
(

0.00

0.2

p

(a)

1.8

0.0

2.04

(2.05, 3.74)

(2.05, 1.01)

p

(b)

280

2.09

4.8

unstable

(5, 300)

p

(c)

6.0

Figure 3.5: Flutter instability curves in γ1-p space for the two-link system with follower force
for: (a) γ2 = 0.02 and p = 1, (b) γ2 = 0.02 and p = 2.05, (c) γ2 = 0.02 and p = 5.

Figure 3.5 (a) and (c) show that γ1 has purely stabilizing effects since they have a single

stability transition. Figure 3.5 (b) reveals a “stabilizing-destabilizing-stabilizing” (S-D-S)

effect due to three stability transitions. These results can also be verified from the number

of sign changes in the first column of the Routh array in Eq.(3.11). For (a) and (c), a single

sign change is obtained. However, three sign changes are observed in (b). The locus of the

eigenfrequencies of the system are plotted in Fig.3.6: γ1 is increased from zero for the cases

24

shown in Fig.3.5.

Im(s)

1.5

γ1 = 0

Im(s)

γ1 = 0

γ1 = 0.008

1.0

γ1 = 0

0.5

-0.01

0.01

Re(s)

-0.04

0.02

Re(s)

-0.004

-1.5

-1.0

(a)

-0.5

γ1 = 0.11

γ1 = 3.74
γ1 = 1.01

(b)

Im(s)

2.0

0.5

γ1 = 0

0.003

Re(s)

-0.5

γ1 = 300

-2.0

(c)

Figure 3.6: The locus of the roots that determine system stability as γ1 increases for the
particular cases: (a) γ2 = 0.02 and p = 1, (b) γ2 = 0.02 and p = 2.05, (c) γ2 = 0.02 and
p = 5.

It is worth noting that for the case where the system exhibited the (S-D-S) effect of

damping (γ2 = 0.02, p = 2.05), the system is originally unstable. As γ1 is increased,

the system is stabilized when the fundamental complex conjugate pair of roots cross the

imaginary axis. When γ1 is increased further, the system is destabilized by the second

pair of complex conjugate roots which cross the imaginary axis to move from the left-half

plane to the right-half plane. For even higher values of γ1 the system regains stability at

γ1 = 3.74. This interesting behavior, where the S-D-S effect is produced by two pairs of

complex conjugate roots, is not discussed in the literature.

To complete the investigation, the flutter instability curve in p-γ1 space is plotted for a

range of γ2 values in the range [0.01, 0.10] - see Fig.3.7. The figure discloses that at higher

values γ2 (γ2 > 0.06), γ1 loses its ability to produce the three stability transitions (stabilizing-

destabilizing-stabilizing). On contrary, for lower values of γ2 the system does exhibit those

three stability transitions.

Similarly, the role of damping in the second link on the system stability is investigated

through the contour plot of the discriminant - see Fig.3.8. Consider the horizontal dotted

25

5

1
γ

0
2.00

1
0
.
0
=
2
γ

2
0
.
0
=
2
γ

3
0.0
=
γ2

0.04
=

γ2

=

γ2

0.05

0.06
0.07
γ 2

=

=

γ 2

γ 2

0.0 8
γ 2

=

=

0.0 9
=

γ 2

0.1 0

p

2.25

Figure 3.7: Flutter instability curves in p-γ1 space for the two-link system with follower force
and range of γ2 ∈ [0.01, 0.10]

line in Fig.3.8 where p = 2. For γ1 ∈ [0.25, 2.00], the discriminant ∆ is positive for four
values of γ1 and negative for the other four values. A plot of the flutter instability curves in

p-γ2 space for γ1 ∈ [0.25, 2.00], show that for γ1 > 1 the discriminant ∆ < 0. Consequently,
Eq.(3.12) has one admissible solution, which also can be verified from a single sign change

in the first column of Routh array. At the same time, for γ1 < 1, the discriminant ∆ > 0,

which may result in three admissible solution for (3.12). The first column of the array shows

three sign changes and accordingly it is possible to obtain three admissible solutions for γ2.

Figure 3.9 shows both cases where Eq.(3.12) has single and three stability transitions for a

5

(−)

(+)

p

5
2
.
0

0
5
.
0

5
7
.
0

0
0
.
1

5
2
.
1

0
5
.
1

5
7
.
1

0
0
.
2

2

(+)

0

0

(−)

γ1

5

Figure 3.8: Contour plot for ∆ = 0 in γ1-p plane for the two-link system with follower force:
regions with positive and negative values of ∆ are denoted by (+) and (−), respectively.

26

range of γ1 ∈ [0.25, 2.00].
5

5
2
.
0
=
1
γ

0
0.5
=
γ1

0.75
=

γ1

1.00

=

γ1

1.25

=

γ 1

=

γ 1

1.5 0

1

γ

1

γ

5

1 .7

0

2 . 0

=

=

p

5

2
γ

0

0

Figure 3.9: Flutter instability curves in p-γ2 space for the two-link system with follower force
and range of γ1 ∈ [0.25, 2.00].

20

2
γ

p = 0.50

p = 0.75
p = 1.00
p = 1.25
p = 1.50
p = 1.75
p = 2.00

0

0

0.2

γ1

stable

γ2 = 6.50

γ2 = 0.55

1

Figure 3.10: Flutter instability curves in γ1-γ2 space for the two-link system with follower
force p ∈ [0.5, 2.0].

As a final analysis, the effect of γ1 and γ2 on the stability of the system is examined for a

range of p values in [0.5, 2.0]; this is shown with the help of Fig.3.10. The system is initially

unstable for γ1 = γ2 = 0 and p ∈ [0.5, 2.0]. Increasing γ1 in the range [0, 1.1], the system
gains stability based upon the value of p. In the meantime, increasing γ2 caused the system

losing its stability. However, Increasing γ2 beyond a certain value (6.5 for example) for the

range of p the two-link system regains stability. To this end, this work extensively studied

the two-link articulated rod subjected to follower end force and highlighted the condition on

27

which stabilizing-destabilizing-stabilizing transition in the system stability occurs due to the

application of damping.

3.2 Two Degree-Of-Freedom Fluid-Conveying Pipe

Y

U

Y

U

ℓ

k, c2

ℓ

k, c2

ℓ

θ1

θ2

k, c1

(a)

ℓ

θ1

θ2

k, c1

(b)

X

X

Figure 3.11: (a) Relative damping and (b) absolute damping is introduced in the fluid-
conveying pipe considered by Benjamin [1] with equal link lengths ℓ and equal torsional
joint stiffness k.

The two-link articulated fluid conveying pipe studied by Benjamin [1] is investigated here

under the assumptions that the link lengths are equal, the torsional stiffness of the joints are

the same, and the motion of the links are restricted to the horizontal plane. Two different

models of damping, namely, relative damping and absolute damping, are introduced in the

system; this is illustrated with the help of Figs.3.11 (a) and (b). We use the same scaling of

variables in [1] together with the following scaling for the damping variable

γ1 = c1s

3

(M + m)kℓ3 ,

γ2 = c2s

3

(M + m)kℓ3

where M and m denote the mass per unit length of the fluid and pipe, respectively. We

investigate the effect of damping on the stability of the system for the two damping models

separately.

28

3.2.1

Effect of Relative Damping on the Stability of the System

Example 3.3: Consider the case where γ1 = 0 and γ2 = γ. The non-dimensional equations

of motion are given as follows:



4

3/2

3/2

1



¨θ1

¨θ2



 +

βu + γ 2βu − γ
βu + γ
−γ



˙θ1

˙θ2



 +

2 − βu2 βu2 − 1

−1

1





θ1

θ2

 =

0

0

 (3.14)

where θ1 and θ2 denote the absolute joint angles of the two links as shown in Fig.3.11, u

denotes the non-dimensional fluid velocity, and β is the mass fraction

β , 3M/(M + m)

Using the procedure outlined in Section 3.1, the characteristic polynomial can be obtained
as follows:

1.75s4 + (8γ + 2βu)s3 + (9 + 4γβu + β2u2

− 2.5βu2)s2 + (γ + 5βu − β2u3)s + 1 = 0

(3.15)

Constructing the Routh array from Eq.(3.15), the conditions for flutter instability can be

expressed by the following cubic equation in γ:

α1γ3 + α2γ2 + α3γ + α4 = 0

α1 , 128βu

α2 , βu2(−128β2u2 + 704β − 80) + 25
α3 , βu(−64β3u4 + 80β2u4 + 328β2u2 − 694βu2 + 1314)
α4 , β2u2(−8β3u4 + 13β2u4 + 40β2u2 − 102βu2 + 169)

(3.16)

The discriminant of the cubic polynomial ∆ can be evaluated using Eq.(3.10). The contour
plot in Fig.3.12 shows the regions where ∆ is positive and negative. A negative value of

29

10.0

u

10.0

(+)

u

(+)

(−)

(+)

(0.08, 9.00)

(0.04, 8.00)

(−)

(0.08, 5.75)

(−)

(−)

β

0.0

0.00

(+)

1.00

0.0

0.00

β

0.10

Figure 3.12: Contour plot for ∆ = 0 in β-u plane for the two-link fluid-conveying pipe with
relative damping: regions with positive and negative values of ∆ are denoted by (+) and
(−), respectively.

∆ implies that there can be at most one real root, and therefore at most one admissible

(positive and real) solution; this implies that γ will certainly not exhibit the stabilizing-

destabilizing-stabilizing effect. A positive value of ∆ implies that there exists three real

roots of γ. However, all of these roots may not be admissible solutions; therefore, γ may or

may not exhibit the stabilizing-destabilizing-stabilizing effect. It can be seen from Fig.3.12

that ∆ is positive for β = 0.08 and u = 5.75; therefore, it not surprising that γ has three

admissible solutions - see Fig.3.13.

3.5

2.5

1.0

γ

0.0

2.0

0
.
1
=
β

5
.
0
=
β

2
.
0
=
β

5
7
.
5
=
u

8
0
.
0
=
β

1.6

5
7
.
5
=
u

unstable

stable

0
0
.
8
=
u

4
0
.
0
=
β

4.0

6.0

u

8.0

0.0

5.6

5.7

5.8

Figure 3.13: Absolute instability curves in γ-u space for the two-link pipe with relative
damping for different values of β; the S-D-S effect of damping is seen for β = 0.04 and 0.08.

30

The flutter instability curves for different values of β are shown in Fig.3.13. Among the

different β values considered in Fig.3.13, γ has a stabilizing-destabilizing-stabilizing effect

only for β = 0.04 and 0.08. This can be verified from the three intersections of the critical

stability curve for β = 0.04 with the vertical line u = 8.00 and the critical stability curve

for β = 0.08 with the vertical line u = 5.75. From Fig.3.12 it can be verified that the points

(β, u) = (0.04, 8.00) and (0.08, 5.75) lie in the region where ∆ > 0.

Remark 3.6: From the results in Fig.3.13 it can be inferred that the damping parameter

γ has a stabilizing-destabilizing-stabilizing (S-D-S) effect for small values of β.

3.2.2 Effect of Absolute Damping on the Stability of the System:

Example 3.4: Consider the case where γ1 = γ2 = γ. The non-dimensional equations of

motion are given as follows:



4

3/2

3/2

1



¨θ1

¨θ2



 +

βu + γ

2βu

0

βu + γ



˙θ1

˙θ2



 +

2 − βu2 βu2 − 1

−1

1





θ1

θ2

 =

0

0

 (3.17)

The characteristic polynomial is given by the expression:

1.75s4 + (5γ + 2βu)s3 + (9 + γ2 + 2γβu + β2u2 − 2.5βu2)s2
+ (3γ + 5βu − γβu2 − β2u3)s + 1 = 0

(3.18)

Constructing the Routh array from Eq.(3.18), the condition for flutter instability can be

described by a quartic equation of γ; this implies that γ may have up to four admissible

solutions. However, a plot of the critical stability curves for different values of β and u

indicate that γ has one, two, or three admissible (positive and real) solutions. Three specific

cases are shown in Fig.3.14; it can be seen that γ has a purely destabilizing (S) effect when

the quartic equation has one admissible solution - see Fig.3.14(a), a stabilizing-destabilizing

(S-D) effect when the quartic equation has two admissible solutions - see Fig.3.14(b), and

31

1.0

1.6

1.2

unstable

unstable

γ

0.5

γ

0.8

γ

0.6

unstable

stable

0.0

1.67

u

(a)

stable

stable

1.73

0.0

2.28

0.0

2.33

5.18

u

(b)

5.23

u

(c)

Figure 3.14: Absolute instability curves in γ-u space for the two-link pipe with absolute
damping showing D, S-D, and S-D-S effects for (a) β = 1.50, (b) β = 0.70, and (c) β = 0.10,
respectively. For (b) and (c), the critical flow velocity ucr for γ = 0 lies outside the range of
u; these values are 2.245 for β = 0.70 and 5.007 for β = 0.10.

a stabilizing-destabilizing-stabilizing (S-D-S) effect when the quartic equation has three ad-

missible solutions - see Fig.3.14(c).

3.3 Three Degree-Of-Freedom Systems

3.3.1 Three-Link System with a Follower Force

Consider the three-link system with the follower force P , shown in Fig.3.15. Each link is

comprised of a point mass m and a massless rod of length ℓ and is constrained to move in the

horizontal plane. The rotational joints have the same torsional stiffness k but their damping

coefficients are different, equal to c1, c2 and c3, respectively. The stability characteristics of

the system can be investigated from the non-dimensional linearized equations of motion:



3 2 1

2 2 1

1 1 1



¨θ1

¨θ2

¨θ3



+





0

γ1 + γ2 −γ2
−γ2
0

γ2 + γ3 −γ3
γ3
−γ3

2p

2 − p −1 − p
−1
0

2 − p −1 + p
−1

1





θ1

θ2

θ3

=





0

0

0

(3.19)







˙θ1

˙θ2

˙θ3



+



32

θ1

ℓ

k, c1

m

ℓ

θ2

Y

m

ℓ

P

k, c2

m

θ3

k, c3

X

Figure 3.15: A three-link system with a follower force.

where p is the non-dimensional follower force, γi, i = 1, 2, 3, are the non-dimensional damping

coefficients, and τ is the non-dimensional time variable; the non-dimensional variables are

obtained using the scaling expressions provided in Eq.(3.2). Assuming the solution to be of

the form θi = Θiesτ , i = 1, 2, 3, we get a sixth-order characteristic polynomial in s. We now

consider the following example.

Example 3.5: Consider the case where γ1 = γ2 , γc = constant. Constructing the Routh

array and equating the coefficient of s1 equal to zero, we get a quintic polynomial in γ3 that

describes the flutter instability curve; the coefficients of the polynomial are functions of γc

and p. The quintic polynomial has five roots but numerical solutions based on different values

of γc and p yield one, two, or three admissible (positive and real) solutions. Four specific

cases are shown in Fig.3.16; it can be seen that γ3 can have a purely destabilizing (D) or

a purely stabilizing (S) effect when the quintic equation has one admissible solution - see

Fig.3.16(a) and Fig.3.16(d), a destabilizing-stabilizing (D-S) effect when the quintic equation

has two admissible solutions - see Fig.3.16(b), and a stabilizing-destabilizing-stabilizing (S-

D-S) effect when the quintic equation has three admissible solutions - see Fig.3.16(c).

33

0.4

9.0

stable

3
γ

0.2

unstable

3
γ

4.0

2.0

1.8

unstable

stable

3
γ

1.0

1.0

3
γ

unstable

stable

unstable

stable

0.0

0.1

p
(a)

0.8

0.0

0.75

1.05

0.0

1.27

p
(b)

p
(c)

0.2

1.32

1.3

1.8

p
(d)

Figure 3.16: Absolute instability curves in γ3-p space for the three-link system with follower
force showing D, D-S, S-D-S, and S effects of damping for (a) γc = 0.01, (b) γc = 2.0, (c)
γc = 4.0, and (d) γc = 6.0, respectively.

3.3.2 Three-Link Fluid-Conveying Pipe

The three-link articulated fluid conveying pipe studied by Benjamin [1] is investigated here

for the case where the link lengths are equal, the torsional stiffness of the joints are the same,

and the motion of the links are restricted to the horizontal plane. Similar to the two-link

pipe shown in Fig.3.11 (a), the configuration of the system is described using absolute joint

angles θi, i = 1, 2, 3. Relative damping is introduced in the system and the non-dimensional

damping coefficients of the three joints are assumed to be γi, i = 1, 2, 3. The non-dimensional

equations are:



7

4.5 1.5

4.5

4

1.5

1.5 1.5

1



¨θ1

¨θ2

¨θ3



+





βu + γ1 + γ2

−γ2
0

+



34

2βu

2βu − γ2
βu + γ2 + γ3 2βu − γ3
βu + γ3

−γ3

−βu2 + 2 −1
−βu2 − 1
2

0

−1

βu2
βu2 − 1

1







˙θ1

˙θ2

˙θ3

θ1

θ2

θ3




=



0

0

0

(3.20)



Example 3.6: Consider the case where γ1 = γ2 = 0.10 and γ3 = γ. The condition for

flutter instability is a quintic polynomial in γ whose coefficients are functions of β and u.

The quintic polynomial has five roots but numerical solutions based on β = 0.10 yielded a

maximum of three admissible (positive and real) solutions. There are two flutter instability

curves - see Fig.3.17 (a); the one to the left is also the absolute instability curve. It can be

seen from this figure that γ has three admissible solutions for both u = 12.0 and 16.5. For

u = 12.0, two of the three roots (γ = 4.63, 5.78) lie on the absolute instability curve. The

system is unstable for γ ∈ (0, 4.63), stable for γ ∈ (4.63, 5.78), and unstable for γ > 5.78;
consequently, the system exhibits a stabilizing-destabilizing (S-D) effect. For u = 16.5,

none of the roots lie on the absolute instability curve; therefore, the system is unstable and

damping has no effect on the stability property of the system.

The locus of the roots are plotted in Fig.3.17 (b) for γ = 0.20. Consistent with Fig.3.17

(a), the locus is shown for u ∈ [4, 20]. At the starting point, where u = 4.0, there are three
pairs of complex conjugate roots in the left-half plane. At u = 4.49, one pair of roots crosses

the imaginary axis into the right-half plane; this identifies the absolute instability curve,

which is also a flutter instability curve. When u is increased to 12.44, a second pair of roots

6

4

2

γ

stable
stable

γ = 5.78

γ = 4.63

absolute instability curve

unstable

u = 4.0

3

u = 12.44

u = 4.49

flutter instability curves

-0.8

0.8

unstable

(4.49, 0.20)

0

4

8

u = 12.44

5
.
6
1
=
u

16

0
.
2
1
=
u

u

(a)

20

-3

(b)

Figure 3.17: (a) Flutter instability curves and absolute instability curve in γ-u space for
the three-link fluid-conveying pipe in Example 3.6 for β = 0.10 (b) locus of all roots of the
system as u increases for the particular case of γ = 0.2.

35

cross the imaginary axis into the right-half plane; this identifies the second flutter instability

curve that intersects the γ = 0.20 line at u = 12.44.

Example 3.7: Consider the case where γ2 = 0.0, γ1 = 0.5 and γ3 = γ. The conditions

for flutter instability results is an octic polynomial in γ whose coefficients are functions

of β and u. The octic polynomial has eight roots but numerical solutions with β = 0.10

yielded a maximum of four admissible solutions. The flutter instability curves are shown

in Fig.3.18(a); it can be seen from this figure that γ has four roots for u = 12.0. Similar

to Example 3.6, there are two flutter instability curves which correspond to the second and

third modes of the system. Unlike Example 3.6, however, the two flutter instability curves

12

8

γ

γ = 11.03

stable

unstable

flutter instability curve
(third mode)

γ

4

flutter instability curve
(second mode)

γ = 3.48

unstable

γ = 2.67

0

4

8

16

20

12
u

(a)

4

3

2

1

0

4

γ = 3.48

γ = 2.67

stable

e

v

r

t e
y

u

c

o l u
s
b ili t

b
a

a
s t

i n

unstable

11
u

(b)

unstable

γ = 2.705

γ = 2.068

u = 12.0

18

Figure 3.18: (a) Flutter instability curves in γ-u space for the three-link fluid-conveying pipe
in Example 3.7 for β = 0.10 (b) A magnified view of the region of intersection of the flutter
instability curves is shown along with the absolute instability curve.

intersect each other; consequently, the absolute instability curve is piecewise continuous and

comprised of segments of both flutter instability curves - see Fig.3.18(b). Similar results

have been reported in the literature earlier [66] where nonlinear analysis was used to show

the intersection of two curves in the bifurcation diagram.

For u = 12.0, three of the four roots (γ = 2.67, 3.48 and 11.03) lie on the absolute

instability curve - see Figs.3.18 (a) and (b). The system is unstable for γ ∈ (0, 2.67), stable
for γ ∈ (2.67, 3.48), unstable for γ ∈ (3.48, 11.03) and stable for γ ∈ (11.03,∞); consequently,
the system exhibits a stabilizing-destabilizing-stabilizing (S-D-S) effect.

36

Remark 3.7: The absolute instability curve in Example 3.7 has a positive slope and two

negative slopes in the range γ ∈ [2.068, 2.705] - see Fig.3.18 (b). This is at odds with the
definition of “destabilizing” and “stabilizing” damping proposed in [15]. We instead propose

the following definition.

Destabilizing (Stabilizing) damping: The addition of damping is “destabilizing” (“sta-

bilizing”) if it causes the system to become unstable (stable) at some value of the non-

conservative forcing variable (u, in Example 3.7). This definition may be overly-reductive

but the emphasis on the non-conservative forcing variable rather than the local slope is

needed to cover all observed cases of stability transitions.

Remark 3.8: The destabilizing effect of damping was first described by Ziegler [14] relative

to an undamped system. The above definition is valid for changes in stability properties

between any two non-negative levels of damping.

3.4 Cantilevered Fluid-Conveying Pipe

For the sake of completeness, we consider a continuous system, namely, a cantilever fluid-

conveying pipe with a damping force applied at a fixed point along the length of the pipe. In

the absence of gravity, the non-dimensional equation of motion of the pipe and its boundary

conditions are:

∂4y(x, t)

∂x4 + u2 ∂2y(x, t)

∂x2 + 2upβ

∂2y(x, t)

∂x∂t

+ cs

∂5y(x, t)

∂4x∂t

+ c

∂y(α, t)

∂t

+

∂2y(x, t)

∂t2

= 0 (3.21)

y(0, t) = 0,

∂y
∂x

(0, t) = 0,

∂2y
∂x2 (1, t) = 0,

∂3y
∂x3 (1, t) = 0

where u is the non-dimensional fluid velocity, β is the ratio of the mass of the fluid to the

mass of the fluid and pipe, cs is the non-dimensional structural damping coefficient, c is the

non-dimensional external damping coefficient, and α ∈ [0, 1] defines the point of application
of the external damping force. Using Galerkin approximation [63], we assume the solution

37

to be of the form

y(x, t) =

NXi=1

ai(t) φi(x)

(3.22)

where N is the number of terms in the approximation, φi(x), i = 1,· · · , N, are functions
that satisfy the boundary conditions, and ai(t), i = 1,· · · , N, are coefficients that will be
determined. Substituting Eq.(3.22) into Eq.(3.21) and integrating over the length of the

pipe, we get

M ¨a + C ˙a + Ka = 0

where the (i, j)-th elements of the M, C and K matrices are defined as

φiφj dx

0

Mij , Z 1
Cij , 2upβZ 1
Kij , Z 1

φ′′
i φ′′

0

0

φiφ′

j dx + csZ 1
j dx + u2(cid:20)φi(1)φ′

0

j(1) −Z 1

0

φ′′
i φ′′

j dx + c φi(α)φj(α)

φ′
iφ′

j dx(cid:21)

and (.)′ denotes the first spatial derivative of (.).

Example 3.8: Consider the fluid-conveying pipe with external damping only, i.e., cs =

0, c 6= 0. A four-mode approximation (N = 4) of the fluid-conveying pipe is obtained
using the first four mode shapes of the cantilever beam [67]. The M, C and K matrices

have dimension four and the characteristic polynomial is of order eight. Using the procedure

outlined in Section 3.1, the conditions for flutter instability can be described by a polynomial

equation in c; the order of the polynomial equation is forty-one and the coefficients are

functions of α, β and u. For α = 0.86 and β = 0.1, the polynomial equation has three

admissible solutions over a small range of u values - see Fig.3.19(a). Thus the fluid-conveying

cantilever also exhibits the stabilizing-destabilizing-stabilizing (S-D-S) effect of damping; for

u = 5.05, for example. The qualitative behavior of the system does not change when the

number of terms in the Galerlin approximation is increased to five (N = 5); this can be seen

from the plot in Fig.3.19(a).

38

25

15

c

5

0

N = 5

N = 4

unstable

stable

5.04

5.05

5.06

5.07

u
(a)

25

15

c

5

0

0
.
1
=
α

9
.
0
=
α

8
.
0
=
α

7
.
0
=
α

2
.
0
=
α

0.3

=

α

6
.
0
=
α

0 . 4

=

α

0.5
=
α

4.5

5.5

6.5
u
(b)

7.5

8.5

Figure 3.19: Absolute instability curves in c-u space for the cantilever fluid-conveying pipe
with β = 0.1 for: (a) α = 0.86 (b) α ∈ [0.2, 1.0]. The S-D-S effect of damping is seen for
α = 0.86.

We complete this example with some observations on how the stabilizing and destabilizing

effects of damping depend on the location where the external damping force is applied. The

absolute instability curves for different values of α are shown in Fig.3.19 (b) for β = 0.1. It

can be seen that the curves move towards the right as the point of application of the force

changes from the free end of the pipe (α = 1.0) to the mid-point (α = 0.5), implying that a

decrease in α in the range [0.5, 1.0] has a stabilizing (S) effect. However, further decrease in

the value of α moves the absolute instability curves back towards the left, i.e., decrease in α

in the range [0.2, 0.5] has a destabilizing (D) effect. Overall, the pipe has the best stability

characteristics when α = 0.5 and the worst stability characteristics when α = 1.0. This is

also supported by the observation that increasing damping has a stabilizing (S) effect for

α = 0.5 (u = 7.5, for example) whereas it has a destabilizing (D) effect for α = 1.0 (u = 4.5,

for example).

Example 3.9: Consider the fluid-conveying pipe with structural damping only, i.e., c =

0, cs 6= 0. Similar to the previous example, a four-mode approximation (N = 4) is used.
The conditions for flutter instability can be described by a polynomial equation in cs; the

order of the polynomial equation is sixty-four and the coefficients are functions of β and u.

For β = 0.35, the polynomial equation has two admissible solutions for cs ∈ [0.0, 0.1] - see

39

Fig.3.20. This figure indicates that the fluid-conveying cantilever exhibits a destabilizing-

stabilizing (D-S) effect for structural damping; for u = 7.7, for example. This problem was

also investigated in [15], where it was claimed that structural damping has a destabilizing

effect for β > 0.29. While this is a correct statement for very small values of damping

(approx. cs < 0.27), slightly higher values of damping have a stabilizing effect.

0.10

stable

0.06

s
c

0.02

0.00

0.027

7.6

unstable

u = 7.7

8.0

u

8.4

8.8

Figure 3.20: Absolute instability curves in cs-u space for the cantilever fluid-conveying pipe
with β = 0.35 showing the D-S effect of damping.

40

Chapter 4

Non-conservative Effects of Dynamic Terminal Moment

In this chapter, we introduce a new type of non-conservative load in the form of a terminal

dynamic moment. For a flexible cantilevered beam the dynamic moment is applied at the

free end and its magnitude varies with the slope or the curvature of the beam at some points

along its length. The stability analysis of such a system has not appeared in the literature

and the non-conservative effects of this type of loading is shown here using both theory

and experiment. This chapter is organized as follows. In Section 1, we provide a problem

formulation where the equation of motion of a cantilevered beam with terminal dynamic

moment and a tip mass is presented. Section 2 highlights the procedure to solve the problem

both analytically and numerically and introduces the characteristic equation as a function of

the non-conservative load parameter. The stability investigations are introduced in Section

3 in the context of various loading frameworks and results for an experimental are presented

in Section 4.

4.1 Problem Formulation

Consider the cantilevered Euler-Bernoulli beam of length L shown in Fig.4.1; it is assumed

to have a uniform cross-sectional area A. The free end of the cantilever is subjected to a

dynamic bending moment M. This bending moment is the reaction of a torque produced by

an actuator mounted at the free end. The mass and mass moment of inertia of the actuator

about the free end is assumed to be m and J, respectively. The transverse displacement of

the beam is denoted by y(x, t). For small deformation, the equation of motion of the beam

and its boundary conditions are:

EI y ′′′′ + ρA ¨y = 0

(4.1)

41

y(0, t) = 0,

y ′(0, t) = 0, EIy ′′(L, t) + J ¨y ′(L, t) = M, EIy ′′′(L, t) = m¨y(L, t)

(4.2)

y

slope : y′(bx, t)
or curvature : y′′(bx, t)

y(x, t)

M

bx

x

L

m, J

x

Figure 4.1: A flexible cantilever beam with a tip mass and a terminal dynamic moment

where E, I, and ρ are the Young’s modulus of elasticity, cross-sectional area moment of

inertia, and the mass per unit volume of the beam, and y ′ and ˙y denote the spatial and time

derivatives of y(x, t). Using the following change of variables

v =

y
L

,

u =

x
L

,

τ = ts EI

ρAL4

the non-dimensional equation of motion and boundary conditions are obtained as

v ′′′′(u, τ ) + ¨v(u, τ ) = 0

(4.3)

(4.4)

v(0, τ ) = 0,

v ′(0, τ ) = 0,

v ′′(1, τ ) + ν ¨v ′(1, τ ) = M ,

v ′′′(1, τ ) = µ ¨v(1, τ )

(4.5)

where v ′ and ˙v denote the partial derivatives of v(u, τ ) with respect to u and τ , respectively,

and

M , ML
EI

,

µ , m
ρAL

,

ν , J

ρAL3

(4.6)

42

We consider two cases where the dynamic moment M is assumed to be proportional to:

(a) the slope of the beam at x =bx, and
(b) the curvature of the beam at x =bx, i.e.
M = Ms , Csy ′(bx, t)
M = Mc , Ccy ′′(bx, t)

:

:

slope-dependent

curvature-dependent

(4.7)

These relationships can be expressed using the following non-dimensional variables

M = M s , C sv ′(α, τ ),

C s , CsL/EI

M = M c , C cv ′′(α, τ ),

C c , Cc/EI

:

:

slope-dependent

curvature-dependent

,

α , bx

L

(4.8)

where α ∈ (0, 1] denotes the location of the point on the beam from where the slope or the
curvature is measured.

4.2 Solution Methods

We solve Eqs.(4.4) and (4.5) both analytically and numerically.

In addition to providing

confidence in our analytical results, the numerical method allows us to examine the structure

of the stiffness matrix and determine conditions under which the matrix loses its symmetric

property, and thereby infer the nature of the instability due to the dynamic moment. The

numerical model is also useful for investigating the effect of damping, when present, on the

stability characteristics of the system [68].

4.2.1 Analytical Solution

Using variable separation, we substitute v(u, τ ) = U(u)T (τ ) to get solutions of the form

T (τ ) = A cos ωτ + B sin ωτ,

U(u) = P1eβu + P2e−βu + P3eiβu + P4e−iβu

ω , β2

(4.9)

43

where A and B are constants that can be obtained from initial conditions, and P1 through P4

are constants that can be obtained from the boundary conditions in Eq.(4.5). The boundary

conditions result in the following transcendental characteristic equations:

[βµ(cos β − cosh β) + sin β − sinh β](cid:2)β(sin β + sinh β) + C s(cos αβ − cosh αβ)
+β4ν(cos β − cosh β)] − [cos β + cosh β − βµ(sin β − sinh β)](cid:2) C s(sin αβ + sinh αβ)
[βµ(sin β − sinh β) − (cos β + cosh β)](cid:2) C c(cos αβ + cosh αβ) − (cos β + cosh β)
+β3ν(sin β + sinh β)] − [sin β − sinh β + βµ(cos β − cosh β)](cid:2) C c(sin αβ + sinh αβ)

−β(cos β + cosh β) + β4ν(sin β + sinh β)] = 0

(4.10a)

−(sin β + sinh β) − β3ν(cos β − cosh β)] = 0

(4.10b)

where Eqs.(4.10a) and (4.10b) correspond to the two cases where the terminal dynamic

moment is slope-dependent and curvature-dependent, respectively. By solving Eq.(4.10),

we get the non-dimensional natural frequencies ω = β2. The constants P1 through P4 are

obtained from the boundary conditions as follows:

P1 = 0, P2 = 1.0, P3 = 0, P4 =

where

γs

γc

: M = M s is slope-dependent

: M = M c is curvature-dependent

γs , C s(sin αβ + sinh αβ) − β [cos β + cosh β − β3ν(sin β + sinh β)]
C s(cos αβ − cosh αβ) + β [sin β + sinh β + β3ν(cos β − cosh β)]
C c(cos αβ + cosh αβ) − (cos β + cosh β) + β3ν(sin β + sinh β)
C c(sin αβ + sinh αβ) − (sin β + sinh β) − β3ν(cos β − cosh β)

γc , −

(4.11a)

(4.11b)

44

Each solution of β obtained from Eq.(4.10) results in a unique mode shape given by the

expression

U(u) =

(cos βu − cosh βu) + γs(sin βu − sinh βu)
(cos βu − cosh βu) + γc(sin βu − sinh βu)

: M = M s is slope-dependent

(4.12)

: M = M c is curvature-dependent

It is important to note that the terminal moment, independent of whether it is slope- or

curvature-dependent,

is a boundary condition. Therefore, a boundary-value problem is

solved to obtain the characteristic equation and determine critical stability points.

4.2.2 Galerkin method

To obtain a numerical solution, the Galerkin method [63] is used; the solution to Eq.(4.4) is

assumed to be of the form

v(u, τ ) =

NXi=1

ai(τ ) φi(u)

(4.13)

where N is the number of terms in the approximation, φi(u), i = 1,· · · , N, are assumed
modes that satisfy the geometric boundary conditions, and ai(τ ), i = 1,· · · , N, are the
modal amplitudes. Substituting Eq.(4.13) into Eq.(4.4) and integrating over the length, we

get

M¨a + Ka = 0

(4.14)

where the (i, j)-th elements of the M and K matrices are defined as

0

Mij ,Z 1
Z 1
Z 1

Kij ,

0

0

φiφj du + µ φi(1)φj(1) + ν φ′

i(1)φ′

j(1)

φ′′
i φ′′

φ′′
i φ′′

j du − C s φ′
j du − C c φ′

i(1)φ′

j(α)

i(1)φ′′

j (α)

: M = M s is slope-dependent

: M = M c is curvature-dependent

(4.15)

and where φ′ denotes the derivative of φ with respect to u. The second and third terms in the

expression of Mij denote the contribution of the mass and rotor inertia of the actuator to the

45

mass matrix. The second term in the expression of Kij denotes the geometric stiffness asso-

ciated with the terminal dynamic moment. When the dynamic moment is slope-dependent,

it is clear that the stiffness matrix K will be symmetric if α = 1 and asymmetric otherwise.

When the dynamic moment is curvature-dependent, the stiffness matrix K will be asymmet-

ric. A nonconservative system is associated with an asymmetric stiffness matrix [56], and

therefore the cantilever beam may lose stability through flutter for both cases where the

terminal moment is slope- or curvature-dependent.

4.3

Instability Investigation

The analytical and numerical methods are both used to investigate the instabilities of the

cantilever beam due to the terminal dynamic moment. For the analytical method, the non-

dimensional natural frequencies ωi = β2

i , i = 1, 2,· · ·, are computed from Eq.(4.10) and
plotted as a function of C s or C c (depending on whether the moment is slope- or curvature-

dependent), for discrete values of α. For the numerical method, the natural frequencies are

obtained by solving the eigenvalue problem associated with Eq.(4.14). The numerical method

is based on a ten-mode approximation (N = 10) and the behavior of the first few natural

frequencies is investigated. The values of the mass and inertia ratios, defined in Eq.(4.6), are

chosen as µ = 7.03, ν = 0.0381. We investigate the two cases where the dynamic moment is

slope-dependent and curvature-dependent, and consider positive and negative values of the

proportionality constant for each of these cases.

4.3.1 Dynamic Moment Proportional to Positive Slope

The results for positive values of C s are shown in Fig.4.2 (a). The first natural frequency of

the beam with a mass at the tip (µ = 7.03, ν = 0.038) is equal to 0.639 in the absence of the

terminal dynamic moment (Cs = 0). When C s is positive, the system loses stability through

divergence: as the value of C s is increased, the first natural frequency reduces to zero. The

1These values were chosen to be identical to those in our experimental setup, described in Section 4.4.

46

0.639

0.7

1
ω

α = 1.00

0.0

0.0

α

=

0

.

7

5

1.5

α = 0.50

Cs
(a)

α = 0.25

120

ω5

ω

ω4

ω3
ω2

20

C s = 111.4

3.0

4.5

0

Cs
(b)

120

Figure 4.2: Terminal moment is proportional to the positive slope of the beam: (a) Variation
in the first natural frequency with C s (Cs ≥ 0) for four different values of α. Analytical and
numerical results are plotted using solid lines and dashed lines, respectively. (b) Analytical
results showing variation of the second, third, fourth, and fifth natural frequencies with C s
for α = 0.25.

critical value of C s, denoted by C

∗

s, is different for different values of α but the overall trend

is the same. The analytical and numerical methods provide similar results, which are shown

for four different values of α, α = {0.25, 0.50, 0.75, 1.00}. It should be noted that the K
matrix is asymmetric when α 6= 1 and the system loses stability through divergence. If C s
is increased beyond its first critical value, some of the higher natural frequencies become

complex for α 6= 1. Plots of some of the higher natural frequencies are shown in Fig.4.2
(b) for α = 0.25; it can be seen that ω3 and ω4 become equal for C s = 111.4, they assume

complex values (not shown) for higher values of C s.

In general, if an undamped system loses stability through flutter due to variation of a

system parameter, two consecutive natural frequencies approach each other, assume real

equal values at the critical value of the parameter, and assume complex conjugate values for

higher values of the parameter. This is described as coupled-mode flutter via a Hamiltonian-

Hopf bifurcation [69] in the literature. This motivates us to adopt the following definition

for the ease of discussion in the next section.

n-th Flutter Instability Mode: The system loses stability through the n-th flutter in-

stability mode if the n-th and (n + 1)-th natural frequencies of the system are equal at the

47

lowest critical value of a system parameter, i.e., the n-th and (n + 1)-th natural frequencies

exhibit coupled-mode flutter via a Hamiltonian-Hopf bifurcation.

4.3.2 Dynamic Moment Proportional to Negative Slope

When C s is negative, the system loses stability through flutter. The results are discussed

for two ranges of α, described below:

range 1: α ∈ [0.02, 0.63] first flutter instability mode
range 2: α ∈ [0.64, 0.95]

second flutter instability mode

C s = −153.8

10

2
ω

,
1
ω

α = 0.3

α = 0.5

α = 0.2

-160

0

0

α = 0.6

α = 0.1

-22

Cs
(b)

ω1

C

∗

s = −17.8

Cs
(a)

ω4

ω3

ω2

70

ω

0

0

Figure 4.3: Terminal moment is proportional to the negative slope of the beam: (a) Variation
in the first four natural frequencies with Cs (C s ≤ 0) for α = 0.1 shows the first flutter
instability mode. Analytical and numerical results are plotted using solid lines and dashed
lines, respectively. (b) Analytical results showing the first flutter instability mode for five
different values of α in range 1:
[0.02, 0.63], including the case shown in (a); the critical
stability point is shown using the • sign.
The first flutter instability mode is observed for range 1: α ∈ [0.02, 0.63] and the results are
shown in Fig.4.3. It was not possible to find a critical value of C s for α < 0.02. For a specific

value of α in range 1, namely α = 0.1, the first four natural frequencies are plotted in Fig.4.3

(a) as the value of C s is decreased from zero. It is seen that the system loses stability for

∗

C s = C

s = −17.8: the first two natural frequencies become equal at this value of C s, they
assume complex values for more negative values of C s and are not shown. The analytical and

numerical results provide almost identical plots for all four natural frequencies. Analytical

48

results for variation of the first two natural frequencies as a function of C s are provided

in Fig.4.3 (b) for five different values of α in range 1:

[0.02, 0.63]; each plot is terminated

when the two natural frequencies become equal (shown by the • sign) at the point of flutter
instability. It is seen that C
s decreases as α is increased from 0.1 to 0.3 but increases when α

∗

is increased further. This can be viewed as a “destabilizing-stabilizing” effect of the location

of the point from where the slope of the beam is measured; further discussion on this topic

appears later in this section.

24

3
ω

,
2
ω

10

0

α = 0.94

25

ω

0.80

=

α

α = 0.70

α = 0.64

α = 0.65

-50

0

0

Cs
(a)

C

∗

s = −52.5

ω3

ω2

ω1

Cs
(b)

-60

Figure 4.4: Terminal moment is proportional to the negative slope of the beam: Analytical
results show the (a) second flutter instability mode for range 2: [0.64, 0.95], (b) Switching
s = −52.5. In both figures,
from first to second flutter instability mode for α = 0.6394 and C
the critical stability point is shown using the • sign.

∗

The second flutter instability mode is observed for higher values of α; the results are

provided in Fig.4.4 (a). Variation of the natural frequencies ω2, ω3 as a function of C s are

provided in Fig.4.4 (a) for five different values of α in range 2:

[0.64, 0.95]; each plot is

terminated when the frequencies become equal (shown by the • sign) at the point of flutter
instability. Similar to the observation in Fig.4.3 (b), the “destabilizing-stabilizing” effect of

α can be observed in Fig.4.4 (a): the value of C

∗

s decreases as α is increased from 0.64 to

0.80 but increases when α is increased further. For α > 0.95, it becomes increasingly difficult

to determine the mode of flutter instability. For α = 0.6394, which lies in between ranges 1

and 2, the mode of flutter instability switches from the first to the second as C s is decreased.

For C

∗

s = −52.5, ω2 = ω3 and ω1 veers off resulting in the second flutter instability mode -

49

see Fig.4.4 (b). This veering phenomenon has been reported earlier in the literature [43].

To summarize the results presented in this section and to further investigate the “destabilizing-

stabilizing” effect of α, we plot the critical stability curve in Fig.4.5 for α ∈ [0.02, 0.95], which
spans ranges 1 and 2. The numerical results (solid line) match well with the analytical re-

sults. The critical stability curve is comprised of two curves for the two ranges of α, each

having a “catenary-like” shape. It is possible to find a value of C s (horizontal line) which

intersects both “catenaries” twice. For example, the horizontal line C s = −11 intersects the
catenaries in ranges 1 and 2 twice. This implies that for C s = −11, increasing the value of α
results in a “destabilizing-stabilizing-destabilizing-stabilizing” effect2, where the first desta-

bilization occurs in the first flutter instability mode and the second destabilization occurs in

the second flutter instability mode.

-80

s
C

-11

-10

0.0

α = 0.02

α = 0.575

α = 0.775

range 1

range 2

unstable

unstable

stable

stable

0.5
α

1.0

Figure 4.5: Terminal moment is proportional to the negative slope of the beam: Critical
stability curve for α ∈ [0.02, 0.95]. The critical stability curve (solid line) was obtained
numerically using a ten-mode approximation; the analytical results are shown for discrete
values of α using the • sign.

2Similar multiple stability transitions have been observed in non-conservative systems for variation in the
level of damping starting with the work by Semler et al. [15]; a survey of the literature along with some new
results can be found in [68].

50

4.3.3 Dynamic Moment Proportional to Positive and Negative Curvature

In this section, we discuss the case where the terminal moment is proportional to the curva-

ture of the beam. For positive values of C c, both analytical and numerical results are shown

in Fig.4.6 for three different values of α, α = {0.25, 0.50, 0.75}. Similar to Fig.4.2 (a), the

0.7

1
ω

0.638

0.0

0.0

α

=

0.25

α

=

α

=

0.50

0.75

Cc

1.0

Figure 4.6: Terminal moment is proportional to the positive curvature of the beam: Variation
in the first natural frequency with C c (C c ≥ 0) for three different values of α. Analytical
and numerical results are plotted using solid lines and dashed lines, respectively.

first natural frequency of the beam with a tip mass (µ = 7.03, ν = 0.038) is equal to 0.638

in the absence of the terminal moment (C c = 0). Also, similar to the results in Section 4.3.1

(moment is proportional to the positive slope), the system loses stability through divergence.

However, unlike the results in Section 4.3.1, the critical value of C c is almost identical, equal

to unity, for different values of α. For negative values of C c, the system loses stability through

flutter. Similar to the results in Section 4.3.2 (moment is proportional to the negative slope),

higher values of α are associated with higher modes of flutter instability but there are four

ranges, as described below, and shown in Fig.4.7:

range 1: α ∈ [0.01, 0.29] first flutter instability mode
range 2: α ∈ [0.30, 0.74]
range 3: α ∈ [0.75, 0.84]
range 4: α ∈ [0.85, 0.88]

third flutter instability mode

fourth flutter instability mode

second flutter instability mode

51

α = 0.01

α = 0.30

α = 0.75

range 1

range 2

5
8
0

.

8
8
0

.

3

4

e
g
n
a
r

e
g
n
a
r

-100

c

C

0
0.0

α

0.9

Figure 4.7: Terminal moment is proportional to the negative curvature of the beam: Critical
stability curve for α ∈ [0.01, 0.88]. The critical stability curve (solid line) was obtained
numerically using a ten-mode approximation; the analytical results are shown for discrete
values of α using the • sign.

4.4 Experimental Investigations

4.4.1 Hardware Description

Laboratory experiments were performed using a steel (E = 200 GPa, ρ = 8000 kg/m3)

cantilever beam with tip mass; the experimental setup is shown in Fig.4.8. The beam has
a uniform rectangular cross-section with height h = 0.05 m and thickness t = 5 × 10−4 m.

Top view

y

z

bx

strain gage

strain gage

t

motor housing

Side view

L

h

smooth surface

x

motor

x

air

bearing

Figure 4.8: Cantilever beam with tip mass: the tip mass is comprised of a motor for appli-
cation of a dynamic moment and an air bearing to avoid excitation of torsional modes.

52

The length and cross-sectional area of the beam are

L = 0.56 m, A = h t = 2.5 × 10−5 m2

The tip mass houses a motor3 that is used for applying the dynamic moment at the free

end of the beam. The relatively large mass of the motor can cause torsional deformation

and/or excite the torsional modes of vibration of the thin beam; a pair of air bearings is

therefore used to support the tip mass and allow it to slide freely over a smooth surface during

bending motion of the beam. The tip mass, which includes the motor, motor housing, and

air bearings, and the ratio of the tip mass to the mass of the beam are given below

m = 0.787 kg, µ , m
ρAL

= 7.03

The mass moment of inertia of the motor, motor housing, and air bearings, and the corre-

sponding non-dimensional constant are given below

J = 0.00135 kgm2,

ν ,

J
ρAL3 = 0.038

In our experiments, the terminal moment was proportional to the curvature of the beam,

which was measured using a pair of strain gages located at a distance of

bx = 0.078 m ⇒ α = bx

L

= 0.14

A P3 strain gage measurement unit4 was used to condition the strain gage signals and

output a voltage proportional to the beam curvature. A dSpace DS1104 board and the

Matlab/Simulink environment was used for closed-loop excitation of the beam using motor

3The motor used in our experiment is a product of Micromo, model number 3863H0124CR. It was
controlled using an analog servo drive; the servo drive is a product of Advanced Motion Controls, model
number 12A8.

4The P3 is a product of Vishay Precision Group.

53

torque as the input and curvature measurement as the output.

4.4.2 Model Verification

The first two natural frequencies of the beam were evaluated experimentally. A sinusoidal

torque was applied to the beam using the motor and the frequency was gradually increased.

The resonant behavior of the beam was used to identify the natural frequencies; they are

listed as ω1 and ω2 in Table 4.1 to distinguish them from the non-dimensional natural fre-

quencies ω. The experimentally determined frequencies match well with those obtained

analytically and numerically (ten-mode Galerkin approximation).

Table 4.1: Natural frequencies of the beam in Fig.4.8: analytical, numerical, and experimen-
tal values.

Natural Frequencies
ω1 (rad/s) ω2 rad/s
20.685
20.688
21.5

1.470
1.470
1.6

Analytical
Numerical
Experimental

4.4.3 Flutter Oscillations

For our experimental setup with α = 0.14, the cantilever beam with tip mass was found

to lose stability through the first flutter instability mode for Cc < 0. The experiments

were conducted in the following manner. Initially, the beam was at rest in its undeformed

configuration. The value of C was reduced below zero (made more negative), incrementally.

The value of C was held constant for some time before it was reduced again. As the value

of C approached the critical value, the beam could be seen vibrating with a very small

amplitude. This amplitude would neither grow nor decay. At the critical value, the amplitude

of vibration increased suddenly and the beam started undergoing limit cycle oscillations.

Although the post-flutter behavior is not the subject of investigation of this work, it is well

known that system nonlinearities lead to limit cycle oscillations at and beyond the point of

flutter instability - see [34, 70, 71], for example.

54

The experimental results pertain to nonlinear limit cycle oscillations at or slightly be-

yond the point of instability; nevertheless, they show reasonable match with analytical results

obtained using a linear model. The analytical and numerical (ten-mode Galerkin approxi-

mation) results based on the Euler-Bernoulli beam model are shown in Fig.4.9; they indicate

is C ∗

that the critical frequency is ω ∗ = 7.34 rad/s or f ∗ = 1.167 Hz and the critical value of Cc
c = −0.2635. These values match reasonably well with the values obtained from experi-
ments - see Table 4.2. The experimental flutter frequency was obtained from the strain gage

voltage signal, which is not presented here.

20

)
s
/
d
a
r
(

2
ω

,
1
ω

0
0.00

c = −0.26
C ∗

ω∗ = 7.34

ω2

ω1

-0.15

Cc (Nm2)

-0.30

Figure 4.9: Analytical and numerical results showing the first flutter instability mode for
the Euler-Bernoulli beam corresponding to the experimental setup in Fig.4.8. The analytical
results are shown using a solid line; the numerical results are shown using a dashed line but
are indistinguishable from the analytical results.

Table 4.2: Analytical, numerical and experimentally observed values of frequency and pro-
portionality constant at the flutter instability point for the experimental setup in Fig.4.8.

Analytical
Numerical
Experimental

1.168
1.167
1.33

f ∗ (Hz) C ∗

c (Nm2)
-0.262
-0.263
-0.21

A still image of the maximum deflection of the beam was used to determine the non-

5The symbols ω∗, f ∗ and C ∗

non-dimensional quantities.

c should be distinguished from their counterparts ω∗, f

∗

and C

∗

c , which are

55

dimensional mode shape6 at or slightly above the instability point. Although this mode

shape pertains to the nonlinear oscillation of the beam, it is found to match reasonably well

with the linear mode shape computed analytically using Eq.(4.12) - see Fig.4.10. In Fig.4.10,

both mode shapes were normalized by their infinity-norm.

1.0

∞

k
v
k

v

0.0

0.0

0.5
u

1.0

Figure 4.10: A comparison of normalized, non-dimensional mode shapes of the cantilever
in Fig.4.8: solid and dashed lines show the mode shapes obtained analytically and from
experimental data.

6A large number of points were chosen along the length of the deflected beam and the Matlab function

“grabit” was used to extract the data from the image file.

56

Chapter 5

Stability Transitions Induced by an Intermediate

Support

The stability investigation in this chapter is heavily inspired by the work done by Lee [72],

where multiple stability transitions were reported due to the introduction of an intermedi-

ate elastic support. The main objective, however, is to highlight the nature of instability

transitions, in a cantilevered beam subjected to a terminal dynamic moment and with an

intermediate support. The investigation reveals a rich set of stability transitions not ob-

served heretofore; these include multiple stability transitions between divergence and flutter

and between different modes of flutter. These transitions usually involving jumps in the

critical non-conservative load and sometimes occurs with the veering phenomenon [43]. In

this chapter, the investigation is comprised of three parts; the first two parts consider the

non-conservative dynamic moment to be proportional to the slope and the curvature at a

point on the structure, respectively. In the last part, the problem of the cantilever beam

with a follower force is revisited to draw attention to the jump phenomenon that has been

overlooked in the literature [42].

5.1 Cantilever Beam with Terminal Dynamic Moment

5.1.1 Mathematical Model

Consider the Euler-Bernoulli beam of length L and cross-sectional area A, shown in Fig.5.1.

The transverse displacement of a point on the beam at a distance x from the fixed end is

denoted by y(x, t). The beam has an intermediate pinned support at x = ℓ and is subjected

to a terminal dynamic moment M that is proportional to the slope or the curvature of the

57

y

slope : y′(bx, t)
or curvature : y′′(bx, t)

M

y(x, t)

ℓ

x

Figure 5.1: A flexible cantilever beam with an intermediate support subjected to a terminal
dynamic moment.

Lbx

beam at x =bx. We introduce the following change of variables
, α = bx

τ = ts EI

ρAL4 ,

κ =

u =

x
L

,

y
L

,

ℓ
L

v =

, M =

ML
EI

L

where E, I, and ρ are the Young’s modulus of elasticity, area moment of inertia, and the

mass per unit volume of the beam. The non-dimensional equation of motion is

v ′′′′(u, τ ) + ¨v(u, τ ) = 0

(5.1)

The dynamic moment is assumed to be proportional to either the slope or the curvature of

the beam at x =bx or u = α; therefore

M =

Cs

Cc

dx

dy(bx, t)
d2y(bx, t)

dx2

⇒ M s = C sv ′(α, τ )
⇒ M c = C cv ′′(α, τ )

, C s ,

CsL
EI

, C c ,

Cc
EI

To simplify the analysis, the beam is studied over two domains; the displacement of the

beam over these domains is defined as:

vL(u, τ ) = v(u, τ )

vR(u, τ ) = v(1 − u, τ )

if u ∈ [0, κ]
if u ∈ [0, 1 − κ)

(5.2)

58

The boundary conditions are

vL(0, τ ) = 0,

v ′
L(0, τ ) = 0,

R(1, τ ) = M =
v ′′



M s =
M c =

C sv ′

L(α, τ )

C sv ′

C cv ′′

R(1 − α, τ )
L(α, τ )

C cv ′′

R(1 − α, τ )

:

:

:

:

if α < κ

if α > κ

if α < κ

if α > κ

,

R (1, τ ) = 0 (5.3)
v ′′′

and continuity requires that we impose the constraints:

vL(κ, τ ) = vR(1 − κ, τ ) = 0,

L(κ, τ ) = −v ′
v ′

R(1 − κ, τ ),

L(κ, τ ) = v ′′
v ′′

R(1 − κ, τ )

(5.4)

It has already been shown [73] that for both C s 6= 0 and C c 6= 0, the system is non-
conservative and the cantilever beam may lose stability through divergence or flutter.

5.2 Terminal Dynamic Moment Proportional to Slope

To solve Eq.(5.1) together with the boundary conditions in Eq.(5.3), Galerkin approximation

[63] is used where the solution is assumed to have the following form

v(u, τ ) =

NXi=1

ai(τ ) φi(u)

(5.5)

where N is the number of terms in the approximation, φi(u), i = 1,· · · , N, are assumed
modes that satisfy the geometric boundary conditions, and ai(τ ), i = 1,· · · , N, are the modal
amplitudes. For the problem with intermediate support, we substitute Eq.(5.5) into Eq.(5.1)

and integrate over the length of the beam. After substituting the boundary conditions in

Eq.(5.3), we get

M¨a + Ka = 0

(5.6)

59

where the (i, j)-th elements of the M and K matrices are obtained as follows

0

Mij ,Z κ
Kij ,Z κ

0

φLiφLj du +Z 1−κ
Lj du +Z 1−κ

φ′′
Liφ′′

0

0

φRiφRj du

φLi and φRi, however, have the following form

φ′′
Riφ′′

Rj du −

C s φ′

Ri(0)φ′

Lj(α)

C s φ′

Ri(0)φ′

Rj(α)

:

:

if α < κ

if α > κ

(5.7)

φLi(uL) = cosh (qiuL) − cos (qiuL) + λi[sin (qiuL) − sinh (qiuL)]

λi ,

cos (κqi) − cosh (κqi)
sin (κqi) − sinh (κqi)

φRi(uR) = ηi[cosh (qi uR) + cos (qi uR)] − µi[sin (qi uR) + sinh (qi uR)]
[sin ([1 − κ]qi) + sinh ([1 − κ]qi)]
[cos ([1 − κ]qi) cosh ([1 − κ]qi) + 1]

ηi , [cos (κqi) cosh (κqi) − 1]
[sin (κqi) − sinh (κqi)] ×
cos ([1 − κ]qi) + cosh ([1 − κ]qi)
sin ([1 − κ]qi) + sinh ([1 − κ]qi)

µi ,

The qi’s in the above equations are obtained by solving

[cos (κq) sinh (κq) − cosh (κq) sin (κq)] × [cos ([1 − κ]q) cosh ([1 − κ]q) + 1]+
[cos ([1 − κ]q) sinh ([1 − κ]q)− cosh ([1 − κ]q) sin ([1 − κ]q)][cos (κq) cosh ([1 − κ]q)− 1] = 0

5.2.1 Model Verification

Previous equations are used to investigate the accuracy of the Galerkin approximation by

comparing results obtained using the analytical solution for the beam without intermediate

support. These results are used to gain confidence in our Galerkin approximation for the

beam with intermediate support, for which analytical results are not readily available. For

the beam without intermediate support, the accuracy of the approximation is verified for

N = 4, 5, 6, 7, 8, by comparing the critical value of the proportionality constant ¯Cs for α =

60

0.25, 0.5, 0.75. The results are shown in Tables 5.1 and 5.2 for positive and negative values

of the proportionality constant ¯Cs. It should be noted that ¯Cs > 0 results in divergence and

¯Cs < 0 results in flutter.

Table 5.1: Critical values of Cs for C s > 0 that results in divergence.

α

Analytical N = 4 N = 5 N = 6 N = 7 N = 8

0.25

0.5

0.75

4.05

2

1.35

3.89

2.04

1.31

4.1

1.99

1.33

3.96

1.99

1.35

4

2.01

1.35

4.03

2.01

1.35

Table 5.2: Critical values of C s for C s < 0 that results in flutter.

α

Analytical N = 4 N = 5 N = 6 N = 7 N = 8

0.25

0.5

-7.15

-13.5

-7.51

-6.82

-7.26

-7.13

-7.03

-11.89

-14.38

-14.38

-13.1

-13.1

0.75

-14.65

-21.67

-15.95

-13.83

-13.33

-13.7

The results in Tables 5.1 and 5.2 reveal that N = 6 provides a sufficiently accurate solution.

For this reason, a six-mode based approximation (N = 6) is used in the subsequent analysis.

5.2.2 Effect of Intermediate Support on Critical Stability

In this section we investigate instability of the beam with intermediate support at three

different locations along the length of the beam, namely, κ = 0.25, 0.50, 0.75. For each of

these cases, the instability is investigated for positive and negative values of C s. We rely on

a six-mode Galerkin approximation of the system described by Eqs.(5.6) and (5.7).

Dynamic Moment Proportional to Positive Slope

We first study the case where the moment is proportional to the positive slope, i.e., C s >

0. The results of the first case, κ = 0.25, are presented in Fig.5.2. Unlike the beam

61

with no intermediate support where stability is lost through divergence alone, here stability

is lost through both divergence and flutter. For α ∈ [0.1, 0.13] ∪ [0.17, 0.25] ∪ [0.83, 1.0],
stability is lost through flutter; divergence instability occurs for the intermediate range of

α ∈ (0.25.0.83). The six-mode approximation fails to capture loss of stability for the ranges
of α ∈ (0.13, 0.17) ∪ (0.8, 0.83).

3
1
.
0

7
1
.
0

100

M
I
F

h
t
5

s
C

M
I
F

h
t
4

α = 0.25

α = 0.8

Divergence

3
8
.
0

M
I
F

d
r
3

0
0.1

α

1.0

Figure 5.2: Critical stability curve for the beam with intermediate support: ¯Cs > 0 and
κ = 0.25.

7
1
0

.

2
3
0

.

50

8
4
0

.

1
5
0

.

9
6
0

.

3
7
0

.

8
8
0

.

s
C

M
I
F

t
s
1

M
I
F

h
t
4

M
I
F

d
n
2

M
I
F

d
n
2

M
I
F

h
t
4

M
I
F

h
t
5

0
0.1

α

1.0

Figure 5.3: Critical stability curve for the beam with intermediate support: ¯Cs > 0 and
κ = 0.50.

For κ = 0.50, the results are provided in Fig.5.3. Similar to the results shown in Fig.5.2

for κ = 0.25, stability is lost through both flutter and divergence. However, divergence

62

instability occurs for a very narrow intermediate range of α ∈ (0.48.0.51). For α ∈ [0.1, 0.48]∪
[0.51, 0.73] ∪ [0.88, 1.0], stability is lost through flutter and the mode of flutter instability
switches randomly as α changes. The six-mode approximation, again, fails to capture loss

of stability for the range of α ∈ (0.73, 0.88).

40

s
C

0
0.1

5
3

.

0

5
4

.

0

4
5

.

0

4
7

.

0

M
I
F

t
s
1

M
I
F

h
t
5

n
o
i
g
e
R
e
c
n
e
g
r
e
v
i
D

M
I
F

d
n
2

α

1.0

Figure 5.4: Critical stability curve for the beam with intermediate support: ¯Cs > 0 and
κ = 0.75.

For κ = 0.75, stability investigations indicated that the loss of stability through flutter

for α ∈ [0.1, 0.54] and through divergence for α ∈ [0.54, 0.74]. The six-mode approximation
failed to capture loss of stability for higher values of α.

A comparison of the results presented in Figs.5.2, 5.3 and 5.4 with those for the beam

with no intermediate support, Section 4.3.1, indicates that stability is lost through both

flutter and divergence when ¯Cs > 0. Furthermore, the range of α values that lead to flutter

and divergence instability and the mode of flutter instability strongly depend on the location

of the intermediate support.

Dynamic Moment Proportional to Negative Slope

We complete this section by investigating the case where the moment is proportional to the

negative slope, i.e., C s < 0. For α varying between 0.1 and 1.0, we again consider the three

cases where the intermediate support locations are κ = 0.25, 0.50, 0.75. Similar to the cases

63

with ¯Cs > 0, stability is lost through both divergence and flutter. For κ = 0.25, stability

is primarily lost through different fluter instability modes as shown in Fig.5.5; divergence

instability occurs for a very small region, namely, α ∈ [0.11, 0.14]. For κ = 0.50, stability is
also lost primarily through different fluter instability modes as shown in Fig.5.6; divergence

instability occurs for a very small region, namely, α ∈ [0.25, 0.29]. For κ = 0.25 and 0.50,

1
1

.

0

4
1

.

0

4
2

.

0

-90

2
4

.

0

6
7

.

0

s
C

M
I
F

h
t
5

M
I
F

d
n
2

0
0.1

M
I
F

t
s
1

α

M
I
F

h
t
5

1.0

Figure 5.5: Critical stability curve for the beam with intermediate support: ¯Cs < 0 and
κ = 0.25.

-80

s
C

0
0.1

5
2

.

0

9
2

.

0

7
3

.

0

9
4

1
5

.

.

0

0

2
6

.

0

2
8

.

0

M
I
F

d
n
2

M
I
F

h
t
3

M
I
F

h
t
4

M
I
F

h
t
4

M
I
F

h
t
3

M
I
F

h
t
5

α

1.0

Figure 5.6: Critical stability curve for the beam with intermediate support: ¯Cs < 0 and
κ = 0.50.

the structure exhibits flutter instability for very small values of α; as α is increased, the

nature of instability switches to divergence and then switches back to flutter. In the region

64

-30

s
C

e
c
n
e
g
r
e
v
i
D

-5
0.1

7
1

.

0

6
2

.

0

2
4

.

0

6
4

.

0

2
7

4
7

.

.

0

0

M
I
F

h
t
5

M
I
F

t
s
1

M
I
F

d
n
2

0.8

e
c
n
e
g
r
e
v
i
D

M
I
F

h
t
4

α

Figure 5.7: Critical stability curve for the beam with intermediate support: ¯Cs < 0 and
κ = 0.75.

of flutter instability, the flutter instability modes switch randomly. In contrast, for κ = 0.75,

the structure looses stability alternately through divergence and flutter as the value of α

increases - see Fig.5.7.

5.3 Terminal Dynamic Moment Proportional to Curvature

In this section, we use the method of variable separation to solve Eq.(5.1). We substitute

v(u, τ ) = U(u)T (τ ) into Eq.(5.1) to get

U ′′′′
U

=

¨T
T

= ω2

(5.8)

The complete solution can then be obtained as

U(u) = C1eβu + C2e−βu + C3eiβu + C4e−iβu, T (τ ) = A cos ωτ + B sin ωτ,

β , √ω (5.9)

where C1 through C4, and A and B are constants. Following the notation in Eq.(5.2), the

shape function in Eq.(5.9) is defined over the two domains that lie to the left and right of

the intermediate support:

65

UL(u) = C1eβu + C2e−βu + C3eiβu + C4e−iβu,

UR(u) = C5eβu + C6e−βu + C7eiβu + C8e−iβu,

if u ∈ [0, κ]
if u ∈ [0, 1 − κ)

(5.10)

where C5 through C8 are four additional constants that are used to describe the shape of

the beam to the right of the intermediate support. Substitution of Eq.(5.10) into Eqs.(5.3)

and (5.4) results in eight simultaneous homogeneous algebraic equations in the unknown

constants C1 through C8. By equating the determinant of the coefficient matrix to zero, we

get the following characteristic equations for α ≤ κ:

C c cos µ1 sin(µ2 − µ1) − C c sinh(αβ) cos µ1 cos(βκ)
− C c cos(αβ) cos µ1 sinh(βκ) + C c sin(αβ) cosh µ1 cosh(βκ)
+ C c cosh µ1 sin(µ2 − µ1) + C c cosh(αβ) sin(βκ) cosh µ1
+ C c sinh(αβ) cosh µ1 cosh(βκ) − C c cosh(αβ) sinh(βκ) cosh µ1
+ C c sin(αβ) cos µ1 cosh(βκ) + C c cosh(αβ) sin(βκ) cos µ1
− C c sinh(αβ) cos(βκ) cosh µ1 + C c sinh(αβ) cos µ1 cosh(βκ)
− C c cosh(αβ) cos µ1 sinh(βκ) − C c cos(αβ) sinh(βκ) cosh µ1
− 2 cos µ1 sinh µ1 + 2 cos(βκ) sinh(βκ) − 2 sin(β) cosh µ1 cosh(βκ)
+ 2 sin µ1 cosh µ1 − 2 sin(βκ) cosh(βκ) + 2 cos µ1 cos(βκ) sinh µ1 cosh(βκ)
+ 2 cos µ1 cos(βκ) sinh(βκ) cosh µ1 = 0

and α > κ:

C c sinh µ2 sin(βκ) cos µ1(cosh µ3 + cosh β)
− 2Cc cos(βκ) sinh(βκ)(sin µ2 + sinh µ2)

+ cosh µ1[ cosh(βκ)(cid:2)C c cos µ4 + 3Cc cos(αβ) − 4β sin β(cid:3)

− 4Cc cos(µ2 − µ1) + 4β sin µ1 + 2 cos(βκ) sinh(βκ)[C c sin(µ1 − µ2)
+ cos µ1(2β − C c sinh µ2)] + C c cosh µ2[ − 4 cos µ1
+ (cos µ3 + 3 cos β) cosh(βκ) + 2 sin µ1 cos(βκ) sinh(βκ)]]
+ 2 sinh µ1[C c cos(βκ) sinh(βκ){cos(µ2 − µ1)
− sinh µ2 sin µ1 + cosh µ2 cos µ1} − 2 cos µ1(β − C c sinh µ2)]

+ cosh(βκ)[2 sin(βκ)(cid:2)C c sin µ2 + C c sinh µ2 − 2β(cid:3)

+ sinh µ1[2β(cos µ3 + cos β) − 2Cc sin(βκ) [cos(µ2 − µ1) + cosh µ2 cos µ1]
− C c sinh µ2(cos µ3 + 3 cos β)]] + 4β cos(βκ) sinh(βκ) = 0

where

µ1 , β(1 − κ), µ2 , β(1 − α), µ3 , β(1 − 2κ), µ4 , β(α − 2κ)

66

(5.11)

(5.12)

Equations (5.11) and (5.12) are functions of ω (ω = β2), κ, α, and C c, and they can

be solved for critical stability for different combination of κ and α values. To this end,

we start with C c = 0 and compute the first twelve natural frequencies ω for κ ∈ [0.1, 0.9]
and α ∈ [0.1, 0.9]. Then, we trace the change in these natural frequencies by gradually
increasing the value of C c up to the point Cc = C cr where the system loses its stability

through either divergence or flutter1. Since there is no damping in the system, stability is

lost through divergence when the fundamental frequency become equal to zero and through

coupled-mode flutter [69] when the roots associated with two consecutive natural frequencies

approach each other and become equal. The proportionality constant C c can be positive or

negative, therefore we study these two cases separately.

The analytical method is used to investigate critical stability; the numerical method

requires the computation of mode shapes and this process becomes cumbersome as the

location of the intermediate support is changed incrementally with a small step size.

In

the previous section, the numerical method was used to study critical stability for the case

where the terminal moment is proportional to the slope of the beam; however, the study was

restricted to three values of κ, namely, κ = 0.25, 0.50, 0.75.

5.3.1 Effect of Intermediate Support on Critical Stability

Moment is Proportional to Positive Curvature

We first consider the case where C c > 0, i.e., the terminal moment is proportional to the

positive curvature. In the absence of intermediate support, stability is always lost through

divergence [73]; here it will be shown that introduction of an intermediate support signifi-

cantly alters the stability characteristics. The critical value of C c depends on both κ and α;

Eq.(5.11) is used to compute Ccr for α ≤ κ and Eq.(5.12) is used to compute C cr for α > κ.
The critical stability surface is obtained by computing C cr on a fine mesh grid of κ ∈ [0.1, 0.9]
1We restricted ourselves to the locus of the first ten natural frequencies; therefore, the results reported in
the next section are based on the implicit assumption that the system does not lose stability through modes
higher than ten.

67

and α ∈ [0.1, 0.9] with a step size of 0.01; the results are shown in Fig.5.8. This figure shows
that the system loses stability through both flutter and divergence; to distinguish between

the two types of instability, divergence is shown by dark points on the grid.

60

Cc

0
0.9

0.5

α

0.9

Flutter

α

0.5

α = 0.48

0.9

α = 0.25

0.5
κ

0.1

0.1

0.1

0.1

(a)

Flutter

0.9

0.5
κ

(b)

Figure 5.8: Terminal moment is proportional to the positive curvature: (a) C cr is plotted
on a fine mesh grid of κ and α values; divergence instability is shown by dark points on the
grid (b) Top view of (a).

The critical stability surface in Fig.5.8 reveals a rich set of stability transitions between

different modes of flutter2 and between flutter and divergence. To illustrate these stability

transitions, we take a closer look at the following five cases, which are marked on Fig.5.8

(b):

Case 1: κ = 0.1, α ∈ [0.1, 0.9]
Case 2: α = 0.48, κ ∈ [0.1, 0.9]
Case 3: κ = 0.9, α ∈ [0.1, 0.9]
Case 4: α = 0.25, κ ∈ {0.249, 0.251}
Case 5: κ = 0.65, α ∈ {0.62, 0.78}

2In our earlier work [73], we defined the term “flutter instability mode” as follows: The system loses
stability through the n−th flutter instability mode if two successive natural frequencies of the system,
numbered n and (n + 1), are equal at the lowest critical value of a system parameter. This is usually referred
to as coupled-mode flutter via a Hamiltonian-Hopf bifurcation [69].

68

For Case 1, the system exhibits divergence instability at α = 0.1 and transitions immediately

to flutter instability as α is increased by 0.001; the value of C cr jumps from 1.00 to 41.87 -

see Fig.5.9 (a). In contrast to the cantilever beam with follower force, discussed in section

5.4, no veering phenomenon [43] is associated with this jump. For α > 0.100, the system

loses stability through different modes of flutter as listed below:

second flutter instability mode

fourth flutter instability mode

third flutter instability mode

range 1: α ∈ [0.101, 0.513) first flutter instability mode
range 2: α ∈ [0.513, 0.717)
range 3: α ∈ [0.717, 0.798)
range 4: α ∈ [0.798, 0.843)
range 5: α ∈ [0.843, 0.872) fifth flutter instability mode
sixth flutter instability mode
range 6: α ∈ [0.872, 0.892)
range 7: α ∈ [0.892, 0.900]

seventh flutter instability mode

It should be noted that an increase in α is associated with an ascending order in the mode of

flutter instability. Furthermore, while there is no jump in the value of C cr during the flutter-

to-flutter instability transitions, each transition is associated with the veering phenomenon.

45

41.87

3
1
5

.

0

7
1
7

.

0

8
9
7

.

0

3
4
8

.

0

2
9
8

.

0

c

C

0
0.1

unstable

C cr = 21.89

divergence, C cr = 1.0

0.5
α
(a)

2.8

ω

stable

0.872

1.4

0

0.9

C cr = 21.89

ω3

ω2

ω1

20

Cc
(b)

Figure 5.9: Case 1: (a) Critical stability curve - divergence instability at α = 0.1 and flutter
instability for α ∈ [0.101, 0.9] (b) Second flutter instability mode at α = 0.513 showing
veering of the locus of the first natural frequency.

69

For instance, at α = 0.513 the mode of flutter instability switches from the first to the second

- see Fig.5.9 (b).

For Case 2, multiple stability transitions occur between flutter and divergence; this is

different from Case 1, where a single transition between divergence and flutter is observed.

The multiple stability transitions for Case 2, including flutter to flutter transitions, are listed

below and shown in Fig.5.10 (a).

first flutter instability mode

seventh flutter instability mode

third flutter instability mode

seventh flutter instability mode

range 1: κ ∈ [0.10, 0.47]
range 2: κ ∈ (0.47, 0.54] divergence instability
range 3: κ ∈ (0.54, 0.58]
range 4: κ ∈ (0.58, 0.67]
range 5: κ ∈ (0.67, 0.72]
range 6: κ ∈ (0.72, 0.78] divergence instability
range 7: κ ∈ (0.78, 0.84]
range 8: κ ∈ (0.84, 0.90]

seventh flutter instability mode

third flutter instability mode

14.47

15

c

C

1.00

0
0.1

7
4

.

0

4
5

.

0

8
5

.

0

7
6

.

0

2
7

.

0

8
7

.

0

4
8

.

0

e
c
n
e
g
r
e
v
i

d

e
c
n
e
g
r
e
v
i

d

0.5
κ
(a)

2.8

c

C

1.0

0.9

0.1

2
2

.

0

8
2

.

0

6
3

.

0

1
5

.

0

3
6

.

0

7
6

.

0

e
c
n
e
g
r
e
v
i

d

0.64

0.9

0.5
α
(b)

Figure 5.10: (a) Case 2: Critical stability curve - divergence instability for κ ∈ (0.47, 0.54]∪
(0.72, 0.78] and flutter instability for other values of κ (b) Case 3: Critical stability curve -
divergence instability for α ∈ (0.631, 0.642] ∪ (0.673, 0.900] and flutter instability for other
values of α.

70

The multiple stability transitions between divergence and flutter occur at κ = 0.47, 0.54,

0.72, and 0.78. Among them, the stability transition at κ = 0.47 is associated with a

large jump in the value of C cr from 14.47 to 1.00; there is no jump in the value of C cr for

the other three transitions. Flutter to flutter instability transitions occur in the range of

κ ∈ (0.54, 0.72] ∪ (0.78, 0.90]; unlike Case 1 where the order of flutter instability changes
sequentially, these transitions alternate between the third and seventh flutter instability

modes. However, similar to Case 1, the critical stability curve has a “catenary-like” shape in

the regions where stability is lost through flutter. Finally, there is no veering phenomenon

associated with any of the stability transitions for Case 2.

For Case 3, the critical stability curve is shown in Fig.5.10 (b).

It can be seen that

stability is lost through divergence for α ∈ (0.63, 0.64] ∪ (0.67, 0.90], and through flutter
everywhere else. Similar to Cases 1 and 2, the critical stability curve has a “catenary-like”

shape in the flutter instability regions. The modes of flutter instability are not specified here

as they change randomly between the first five modes. Despite multiple transitions between

divergence and flutter, and flutter-to-flutter, unlike Cases 1 and 2 there is no jump in the

value of C cr.

4.5

ω

1.5

0

ω6

ω4

ω3

ω2

ω1

ω5

ω5

4.5

ω

ω6

ω4

ω3

ω2

ω1

1.5

25

0

9
4
.
0
2

C c
(a)

3
7
.
0

1

C c
(b)

Figure 5.11: Case 4: Locus of the first six natural frequencies as C c is increased from zero
for α = 0.25 and (a) κ = 0.249, (b) κ = 0.251.

71

For Case 4, the locus of the first six natural frequencies are plotted as C c is increased

from zero. These loci, shown in Figs.5.11 (a) and (b), indicate a transition from the first

flutter instability mode to the fifth flutter instability mode due to a small change in the value

of κ; importantly, the flutter-to-flutter transition is accompanied with a jump in the value

of C cr from 20.49 to 0.73.

Remark 5.1: A closer look at Fig.5.8 (b) indicates that the flutter and divergence instability

regions are not accurately demarcated by the diagonal line α = κ. For example, for α = 0.25,

the two points κ = 0.249 and κ = 0.251 lie on opposite sides of the diagonal line, yet stability

is lost through flutter for both of these cases. It has been observed that flutter-to-flutter

transitions occur with jumps in the critical value of C cr wherever the flutter instability region

extends from the region α > κ to the region α < κ. The importance of this observation

will be better appreciated in the next section where we study the case where the moment is

proportional to the negative curvature.

30

c

C

15

first-mode

flutter

e
c
n
e
g
r
e
v
i

d

0
0.62

first

jump

0.70
α
(a)

3.4

second jump

ω

second-mode

flutter

2.0

0

0.78

17.34

ω3

ω2

ω1

18

Cc
(b)

Figure 5.12: Case 5: (a) Critical stability curve - divergence instability for α ∈ (0.62, 0.65],
first-mode flutter instability for α ∈ (0.65, 0.715] and second-mode flutter instability for
α ∈ (0.715, 0.78] (b) Second-mode flutter instability at α = 0.716 showing veering of the
locus of the first natural frequency.

For Case 5, the critical stability curve is shown in Fig.5.12 (a). It shows two jumps in

the value of C cr. The value of C cr jumps from 1.0 to 9.27 during the transition between

72

divergence and flutter at α = 0.65. The value of C cr jumps from 13.29 to 17.3 during the

transition between the first and second flutter instability modes at α = 0.715. The second

jump results from the first natural frequency veering off; this is illustrated in Fig.5.12 (b).

Moment is Proportional to Negative Curvature

-30

Cc

0
0.9

0.5

α

0.9

α

0.5

0.9

0.5
κ

0.1

0.1

0.1

0.1

(a)

Flutter

0.5
κ
(b)

0.9

Figure 5.13: Terminal moment is proportional to the negative curvature: (a) C cr is plotted
on a fine mesh grid of κ and α values; divergence instability is shown by dark points on the
grid (b) Top view of (a).

We now consider the case where C c < 0, i.e., the terminal moment is proportional to the

negative curvature. In the absence of intermediate support, stability is always lost through

flutter [73]; here it will be shown that introduction of an intermediate support results in both

divergence and flutter instabilities. Similar to positive curvature case, the critical stability

surface is first obtained; the results are shown in Fig.5.13, where divergence is shown by dark

points on the grid to distinguish it from flutter. A comparison of the results in Fig.5.13 with

those in Fig.5.8 indicate several similarities and dissimilarities:

• For the case C c > 0, stability is lost through divergence only when α < κ. In contrast,

for the case C c < 0, stability is lost through divergence only when α > κ.

• For the case C c < 0, the diagonal line α = κ clearly demarcates the divergence and
flutter instability regions. In contrast, for the case C c > 0, stability can be lost through

73

flutter on both sides of the diagonal line, i.e., α > κ and α < κ.

• For both C c > 0 and C c < 0, the diagonal line α = κ is associated with a jump
in the value of C cr. For C c < 0 this jump is associated with a stability transition

between divergence and flutter. In contrast, for C c > 0, the jump can be due to a

flutter-to-flutter transition or a transition between divergence and flutter.

• For both C c > 0 and Cc < 0, the critical stability surface has an inverted dome-like
shape, or the shape of a two-dimensional catenary, when stability is lost through a

specific mode of flutter instability - see the region α > κ in Fig.5.13 (a), for example.

Thus, fixing the value of one parameter (α or κ) results in catenary-like shapes when

C cr is plotted against the other parameter - see Fig.5.9 (a), Figs.5.10 (a) and (b), and

Fig.5.12 (a).

5.4 Cantilever Beam with a Follower End Load

5.4.1 Mathematical Model

y

y(x, t)

ℓ

x

L

P

x

Figure 5.14: A flexible cantilever beam with an intermediate support subjected to a follower
force

Consider the Euler-Bernoulli cantilever beam in Fig.5.14, which is subjected to a follower

force P . The beam has length L, cross-sectional area A, and is simply supported at an

intermediate point x = ℓ. The transverse displacement of the beam is denoted by y(x, t).

74

By introducing the change of variables

v =

y
L

,

u =

x
L

,

τ = ts EI

ρAL4 ,

σ =

P L2
EI

,

κ =

ℓ
L

the non-dimensional equation of motion can be expressed as

v ′′′′(u, τ ) + σv ′′(u, τ ) + ¨v(u, τ ) = 0

(5.13)

where E, I, and ρ are the Young’s modulus of elasticity, area moment of inertia, and density

of the beam, (.)′, and ˙(.) denote the partial derivatives of (.) with respect to u and τ , respec-

tively. Similar to Section 5.1, the beam is studied over two domains and the displacement

of the beam over these domains is defined as:

vL(u, τ ) = v(u, τ )

vR(u, τ ) = v(1 − u, τ )

if u ∈ [0, κ]
if u ∈ [0, 1 − κ)

(5.14)

The above equations are identical to those in Eq.(5.2). The boundary conditions are

vL(0, τ ) = v ′

L(0, τ ) = 0,

v ′′
R(0, τ ) = v ′′′

R (0, τ ) = 0

(5.15)

and continuity requires that we impose the constraints:

vL(κ, τ ) = vR(1 − κ, τ ) = 0,

v ′
L(κ, τ ) = −v ′

R(1 − κ, τ ),

v ′′
L(κ, τ ) = v ′′

R(1 − κ, τ )

(5.16)

5.4.2 Analytical Solution

To solve Eq.(5.13) together with the boundary conditions and constraints in Eqs.(5.15) and

(5.16), we assume the solution to be of the form

v(u, τ ) = U(u) T (τ ),

(5.17)

75

Substituting Eq.(5.17) in Eq.(5.13), yields

U ′′′′
U

+ σ

U ′′
U

=

¨T
T

, ω2

(5.18)

where ω2 is constant. The complete solution can be obtained from (5.18) as follows:

T (τ ) = A sin (ωτ ) + B cos (ωτ )

U(u) = C1eiλ1u + C2e−iλ1u + C3eλ2u + C4e−λ2u

λ1 ,s(

) +r(

σ
2

σ
2

)2 + ω2,

λ2 ,s−(

) +r(

σ
2

σ
2

)2 + ω2

(5.19)

where A and B are constants that can be obtained from initial conditions, and C1 through

C4 are constants that can be obtained separately for each domain from the eight boundary

conditions and constraints in Eqs.(5.15) and (5.16). The nontrivial solution of the eight

constants result in the following transcendental characteristic equation:

hλ1λ2
hλ1λ2(cid:0)λ2

2 sin ([1 − κ]λ1) cosh ([1 − κ]λ2) − λ2

1λ2 cos ([1 − κ]λ1) sinh ([1 − κ]λ2)i×
2(cid:1) sin (κλ1) sinh (κλ2) − 2λ1λ2 {cos (κλ1) cosh (κλ2) − 1}i−

h −(cid:0)λ2
1 − λ2
2(cid:1) sin ([1 − κ]λ1) sinh ([1 − κ]λ2)+2λ2

2λ2

2i×
hλ2 sin (κλ1) cosh (κλ2) − λ1 cos (κλ1) sinh (κλ2)i = 0 (5.20)

1 cos ([1 − κ]λ1) cosh ([1 − κ]λ2)+λ4

1+λ4

1 − λ2

Equation (5.20) is a function of ω, κ, and σ and can be solved for the combination of these

parameters that result in critical stability. To this end, we start with σ = 0 and obtain the

natural frequencies ω for κ ∈ [0.1, 1.0]. Then, we trace the change in the natural frequencies
by gradually increasing the value of σ up to the point σ = σcr where the system loses its

stability through either divergence or flutter.

76

5.4.3 Numerical Solution

To be consistent with the literature where both analytical and numerical methods have been

used [41, 42], we now obtain a numerical solution using a Galerkin approximation for both

domains; the solutions are assumed to have the form:

vL(u, τ ) =

NXi=1

ai(τ ) φLi(u),

vR(u, τ ) =

NXi=1

ai(τ ) φRi(u)

(5.21)

where N is the number of terms in the approximation, φLi(u), φRi(u), i = 1,· · · , N, are
assumed modes that satisfy the boundary conditions and constraints in Eqs.(5.15) and (5.16),

and ai(τ ), i = 1,· · · , N, are the modal amplitudes. Substituting Eq.(5.21) into Eq.(5.13)
and integrating over the length of the beam, we get

M¨a + Ka = 0

(5.22)

where the (i, j)-th elements of the M and K matrices are

Mij ,Z κ
Kij ,Z κ

0

0

φLi φLj du +Z 1−κ
Lj du +Z 1−κ

φ′′
Li φ′′

0

0

φRi φRj du

(5.23)

φ′′
Ri φ′′

Rj du − σ(cid:20)Z κ

0

φ′
Li φ′

Lj du +Z 1−κ

0

φ′
Ri φ′

Rj du + φRi(0)φ′

Rj(0)(cid:21)

The assumed modes φLi were chosen to be the modes of a freely vibrating cantilever-pinned

beam, namely

φLi(u) = cosh (qiu) − cos (qiu) + λi[sin (qiu) − sinh (qiu)],

λi , cos (κqi) − cosh (κqi)
sin (κqi) − sinh (κqi)

77

Similarly, the assumed modes φRi were chosen to be the modes of a freely vibrating free-

pinned beam, namely

φRi(u) = ηi[cosh (qiu) + cos (qiu)] − µi[sin (qiu) + sinh (qiu)]

[sin ([1 − κ]qi) + sinh ([1 − κ]qi)]
[cos ([1 − κ]qi) cosh ([1 − κ]qi) + 1]

ηi , [cos (κqi) cosh (κqi) − 1]
[sin (κqi) − sinh (κqi)]
cos ([1 − κ]qi) + cosh ([1 − κ]qi)
sin ([1 − κ]qi) + sinh ([1 − κ]qi)

µi ,

In the expressions for the assumed modes above, the qi’s are obtained by solving the tran-

scendental equation

[cos (κq) sinh (κq) − cosh (κq) sin (κq)][cos ([1 − κ]q) cosh ([1 − κ]q) + 1]+
[cos ([1 − κ]q) sinh ([1 − κ]q)− cosh ([1 − κ]q) sin ([1 − κ]q)][cos (κq) cosh ([1 − κ]q)− 1] = 0

which can be obtained from Eq.(5.20) by substituting σ = 0. The last three terms in the

expression for Kij in Eq.(5.23) represents the geometric stiffness associated with the follower

force σ and the last term contributes to asymmetry. Since the stiffness matrix is asymmetric,

the system can lose stability thought flutter or divergence [72].

5.4.4 Effect of Intermediate Support on Critical Stability

Analytical and the numerical methods were used to investigate the the role of intermediate

support on the critical stability of the system. The critical stability curve obtained analyti-

cally by solving Eq.(5.20) is shown in Fig.5.15 (a). To independently verify the results, we

computed the critical value of σ for select values of κ, κ = {0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9},
by solving the eigenvalue problem for the system in Eq.(5.22) using a three-mode (N = 3)

approximation. These results, shown with the help of • symbols in Fig.5.15 (a), indicate a
perfect match with the results obtained analytically. Two observations can be made from

the results in Fig.5.15 (a). First, the system loses stability through flutter for κ ≤ 0.50

78

and through divergence for κ > 0.50. Second, the transition of stability from flutter to

divergence is accompanied by a jump in the critical value of σ:

for κ = 0.5000, stability

is lost through flutter with σcr = 71.8; for κ = 0.5001, stability is lost through divergence

with σcr = 80.6. Figures 5.15 (b) and (c) provide a closer look at the neighborhood of the

transition point κ = 0.5. It is clear from these figures that the transition of stability from

flutter to divergence is accompanied by what is known as the “veering” phenomenon, first

reported in [43]; the loci of the first two natural frequencies approach each other but veer

off just before meeting, resulting in divergence instead of flutter. It should be noted that

immediately prior to veering off, the loci of the natural frequencies meet at a point where

their slopes are unequal (κ = 0.5, σ = 71.8); this is different from the “dome-shaped” curves

created by the loci for lower values of κ - see Fig.5.15 (b).

7
.
5
4

2
.
3
5

9
.
2
6

8
.
1
7

7
.
7
7

100

Flutter

κ = 0.5001

κ = 0.5000

Divergence

60

σ = 80.6

σ = 71.8

σ

60

ω

κ increasing

8
.
1
7

29

ω

ω2

0
0
0
5
.
0
=
κ

ω1

1
0
0
5
.
0
=
κ

20

0.1

κ

(a)

0

0

1.0

71.8

σ

(b)

18

64

90

6
.
4
7

78

σcr = 80.6

σ

(c)

Figure 5.15: (a) Variation in the critical follower force σ as the location of the intermediate
support κ is varied; both analytical (solid line) and numerical (dots) results are shown (b)
Variation in the critical value of σ for 7 different values of κ; stability is lost through flutter
for κ = {0.40, 0.45, 0.49, 0.50} and through divergence for κ = {0.51, 0.52, 0.53} (c) Loci of
the first two natural frequencies for κ = {0.5, 0.5001} showing jump and veering phenomena;
mode of instability changes from flutter to divergence and the critical value of σ jumps from
71.8 to 80.6 when κ changes from 0.5 to 0.5001.

79

Remark 5.2: To the best of our knowledge, the only publication on the transition of critical

stability in a cantilever beam with a follower force due to an intermediate support is due

to Elishakoff [42]. In this work, the transition of critical stability from flutter to divergence

was correctly shown to occur at κ = 0.5; however, it was claimed that no jump phenomenon

was associated with the transition. In this work, we have used analytical methods to verify

the existence of the jump phenomenon and shown that the discontinuous jump is due to the

veering off by the loci of the natural frequencies.

80

Chapter 6

Non-conservative Behavior of Hinged Beam with a

Dynamic Moment

The stability characteristics of a hinged beam subjected to a non-conservative dynamic mo-

ment are investigated in this chapter. The dynamic moment is proportional to the curvature

of the beam at some point along its length and the analysis is extended to include interaction

of the beam with an external flow. In the absence of external flow, stability is lost through

divergence or flutter, depending on the location of the point of measurement of curvature and

the sign of the applied moment. In the presence of external flow, the mechanism of flutter in-

stability changes from double-mode flutter to single-mode flutter due to the introduction of a

damping-related term in the dynamic model. An extension of the fluid-structure interaction

problem is considered by affixing the flexible hinged beam to a rigid body. This results in

thrust applied on the rigid body due to the traveling waveform of the flexible-hinged beam at

point of flutter instability. The flutter instability conditions are determined and illustrated

through examples.

6.1 Hinged Beam with Dynamic Moment

6.1.1 Mathematical Model

Consider the hinged-beam of length L with uniform cross-sectional area A, shown in Fig.6.1.

The transverse displacement of the beam at a distance x from the hinged end is denoted

by y(x, t). The beam is subjected to a dynamic moment M at the hinged joint, which

has torsional stiffness k. The dynamic moment M is proportional to the beam curvature

at a distance bx from the hinged support, i.e., M = Cy ′′(bx, t). For small deformation, the

non-dimensional equation of motion and boundary conditions are

81

x

curvature : y′′(bx, t)

bx

x

L

y

M

k

y(x, t)

Figure 6.1: A flexible hinged-beam with dynamic moment

v ′′′′(u, τ ) + ¨v(u, τ ) = 0

v(0, τ ) = 0,

v ′′(1, τ ) = 0,

v ′′(0, τ ) = Cv ′′(α, τ ) − kv ′(0, τ )
v ′′′(1, τ ) = 0

(6.1)

(6.2)

where (.)′ and ˙(.) denote the partial derivatives of (.) with respect to u and τ , respectively,

and the non-dimensional variables are defined as follows:

v =

y
L

,

u =

x
L

,

τ = ts EI

ρAL4 , C , C

EI

,

k , kL
EI

, α = bx

L

and E, I and ρ denote the Young’s modulus, area moment of inertia, and mass per unit

volume of the beam. The stability of the system is analyzed using the Galerkin method [63],

where the solution to Eq.(6.1) is assumed to be of the form

v(u, τ ) =

NXi=1

ai(τ ) φi(u)

(6.3)

where N is the number of terms in the approximation, ai(τ ) is the i-th modal amplitudes,

and φi(u) is the i-th orthogonal eigenfunction of free vibration of the hinged-beam with the

spring. The eigenfunctions have the form

82

φi(u) =

k
qi

[cos (qiu) − cosh (qiu)] + sin (qiu) + sinh (qiu) − γ[sin (qiu) − sinh (qiu)]

(6.4)

γ , k[cosh (qi) + cosh (qi)] + β[sin (qi) − sinh (qi)]

β[sin (qi) + sinh (qi)]

where the q′

is are obtained by solving

k + cosh (q)[k cos (q) + β sin (q) − β(cos (q) sinh (q)] = 0

Substituting Eq.(6.3) into Eq.(6.1) and integrating over the length of the beam along with

satisfying Eq.(6.2), we get

M¨a + Ka = 0

(6.5)

where the (i, j)-th elements of the M and K matrices are

Mij ,Z 1
Kij ,Z 1

0

0

φi(u)φj(u) du

φ′′
i (u)φ′′

j (u) du + Cφ′

i(0)φ′′

j (α) − kφ′

i(0)φ′

j(0)

(6.6)

The second term in the expression of Kij denotes the geometric stiffness associated with

the dynamic moment that is applied at the hinge. When C 6= 0, the stiffness matrix K
will be asymmetric. This implies that the system will be non-conservative [56], and can

therefore loose stability through either flutter or divergence [4]. Since Eq.(6.5) represents a

mass-spring system, flutter oscillations will not produce traveling waves.

6.1.2

Instability Investigation

The instabilities of the hinged-beam in Eq.(6.5) is investigated using a six-mode approxima-

tion (N = 6). The point where the curvature is measured, α, was varied along the length

of the beam for both cases where C is positive and negative. The characteristic equation is

83

a function of C, k and α; here, k is assumed to be equal to 0.5. By specifying α we search

for critical stability points by varying C. Flutter instability occurs as two successive natural

frequencies approach each other and coincide before they turn into a complex conjugate pair;

on the other hand, divergence instability happens when the lowest natural frequency passes

through the origin as the equilibrium shifts form stable to unstable. The loss of stability

through divergence is always associated with the first mode but flutter instability occurs in

different flutter instability modes. For ease of discussion we recall the definition of flutter

instability mode from Chapter 3:

Flutter Instability Mode (FIM): The system looses stability through the n-th flutter

instability mode if two successive natural frequencies of the system, numbered n and (n + 1),

are equal at the lowest critical value of a system parameter.

For the case where C > 0 and α ∈ [0, 1], the results are shown in Fig.6.2. The system
loses stability through divergence for α ∈ [0, 0.4]∪[0.45, 0.47]; it loses stability through flutter
for α ∈ (0.4, 0.45) and for α > 0.47. In the region where stability is lost through flutter, the
mode of instability changes from fourth to third to second, and then jumps back to fourth

as α increases - see Fig.6.2.

45

C

0
4
.
0

5
4
.
0

7
4
.
0

1
6
.
0

1
8
.
0

n
o
i
g
e
R
e
c
n
e
g
r
e
v
i
D

n
o
i
g
e
R
M
I
F

h
t
r
u
o
F

n
o
i
g
e
R
M
I
F
d
r
i

h
T

n
o
i
g
e
R
M
I
F

d
n
o
c
e
S

0.0

0.0

Divergence Region
α

n
o
i
g
e
R
M
I
F

h
t
r
u
o
F

0.9

Figure 6.2: Critical stability curve for the hinged beam: C versus α for the case where
k = 0.5 and C > 0.

84

For C < 0, the critical stability of the system is investigated for α ∈ [0.16, 0.99]; the results
are shown in Fig.6.3. The system always loses its stability through flutter, independent of

the value of α. Similar to the case of C > 0 shown in Fig.6.2, the flutter instability mode

switches sequentially downward from the fifth mode to the first mode as α increases from

0.16 to 0.74. For α > 0.74, the mode of instability switches back to higher modes. The

critical stability curve could not be obtained for α < 0.16; this is likely due to the fact that

a six-mode approximation was used and the system loses stability in modes greater than or

equal to six.

6
1
.
0

9
1
.
0

4
2
.
0

1
3
.
0

4
4
.
0

4
7
.
0

4
8
.
0

9
9
.
0

-80

C

n
o
i
g
e
R
M
I
F

h
t
f
i

F

n
o
i
g
e
R
M
I
F

h
t
r
u
o
F

n
o
i
g
e
R
M
I
F
d
r
i

h
T

n
o
i
g
e
R
M
I
F

d
n
o
c
e
S

0.0

2.0

n
o
i
g
e
R
M
I
F

t
s
r
i
F

α

n
o
i
g
e
R
M
I
F
d
r
i

h
T

n
o
i
g
e
R
M
I
F

h
t
f
i

F

0.1

Figure 6.3: Critical stability curve for the hinged beam: C versus α for the case where
k = 0.5 and C < 0.

6.2 Effect of External Flow on Flutter Instability

6.2.1 Mathematical Model

We now consider the case where the hinged-beam, introduced in Section 6.1, is immersed in

an inviscid fluid with constant velocity Ue - see Fig.6.4. Using small deformation assumption,

the non-dimensional equation of motion and boundary conditions are [58, 69]

v ′′′′(u, τ ) + U

2

e v ′′(u, τ ) + 2pβ U e ˙v ′(u, τ ) + ¨v(u, τ ) = 0

(6.7)

85

x

curvature : y′′(bx, t)

bx

x

L

y

Ue

M

k

y(x, t)

Figure 6.4: The hinged-beam in Fig.6.1 immersed in axial flow

v(0, τ ) = 0,

v ′′(1, τ ) = 0,

v ′′(0, τ ) = Cv ′′(α, τ ) − kv ′(0, τ )
v ′′′(1, τ ) = 0

(6.8)

C ,

C
EI

,

k ,

KL
EI

A comparison of Eqs.(6.1) and (6.7) that external flow introduces two new terms in the

equation of motion. It will be shown later that these terms alter the stiffness and damping

characteristics of the system. Equations (6.7) and (6.8) are non-dimensionalized using the

following change of variables:

u =

x
L

,

v =

y
L

, α = bx

L

,

τ =

t

L2(cid:18) EI

m + Mf(cid:19)1/2

We also define the following non-dimensional variables:

U e =(cid:18)Me

EI(cid:19)1/2

UeL,

β =

Me

Me + m

(6.9)

(6.10)

where m, and Me are the mass per unit length of the beam and mass per unit length of

the external fluid that is obtained using the approximation presented in [58]. Following the

procedure outlined in Section 6.1.1, we get

M¨a + D ˙a + Ka = 0

(6.11)

86

where the (i, j)-th elements of the M, D and K matrices are

φi(u)φj(u) du

0

Mij ,Z 1
Dij , 2pβ U eZ 1
Kij ,Z 1

φ′′
i (u)φ′′

0

0

φi(u)φ′

j(u) du

2

eZ 1

0

j (u) du − U

φ′′
i (u)φ′′

j (u) du

+ Cφ′

i(0)φ′′

j (α) − kφ′

i(0)φ′

j(0) + U

2
eφi(1)φ′

j(1)

(6.12)

The term Dij represents the damping due to the external flow. Similarly, the second and

the fifth term in the expression of Kij denote the stiffness effects due to the external flow.

The third term in the expression of Kij denotes the geometric stiffness due to the dynamic

moment. Finally, the fourth term corresponds to the torsional stiffness at the rotational

joint. The stiffness matrix K is asymmetric and consequently the system may loose stability

through divergence or flutter. The main focus of this section is to highlight the effects of the

external flow on the flutter instability modes.

6.2.2 Flutter Instability and Nature of Oscillations

We investigate the effects of external flow on flutter instability. We first consider the example

form Section 6.1.2 where k = α = 0.5 and C < 0. In this case, the eigenfrequencies ω are

purely imaginary in the absence of the external flow and stability is lost through flutter. The

root locus is plotted in Fig.6.5; the roots approach each other as C decreases from zero, then

coexist at C = −2.3, and finally become complex conjugate pairs. As a result, the structure
looses stability through the first flutter instability mode - see Fig.6.3. In the presence of

external flow, the eigenfrequencies ω depend on two additional parameters, namely, the flow

velocity U e and the mass ratio β. In general, for a specific set of the system parameters, k, α,

β, U e, and C, the ω’s appear as complex conjugate pairs due to the presence of the damping.

The locus of these eigenfrequencies in the complex plane is known as the Argand diagram.

87

ω2 = 49.44(C = 0)

ωcr = 32.54
(C cr = −2.3)

)
ω
(
m
I

44.72

38.73

22.36

ω1 = 14.91(C = 0)

-10

0

10

Re(ω)

Figure 6.5: Root locus of the first two fundamental frequencies for the hinged beam without
external flow: C < 0 and k = α = 0.5.

The point on the locus where Re(ω) changes sign form positive to negative constitutes the
onset of flutter. Accordingly, the structure undergoes a sustained oscillation with a critical

frequency ω = ωcr, which is a purely imaginary. The Argand diagram for k = α = 0.5,

β = 0.9, U e = 0.7, and C < 0 is shown in Fig.6.6, where the locus of the first three natural

frequencies are presented.

Re(ω)

)
0
=
C
(
1
ω

-15

0

10

10

)
0
=
C
(
2
ω

9
7
.
1
−
=
C

)
0
=
C
(
3
ω

C cr = −1.63

Im(ω)

110

Figure 6.6: Argand diagram for k = α = 0.5, U e = 0.7, and β = 0.9; the critical C is denoted
by C cr.

The results in Fig.6.6 indicate that the structure becomes unstable for C < −1.63. For the
critical value of C = C cr = −1.63, the system undergoes sustained oscillations of the form:

v(u, τ ) = φcr(u) sin (ωcrτ )

(6.13)

88

Substituting Eq.(6.13) into Eq.(6.7), we get

[φ′′′′

cr (u) + U 2φ′′

cr(u) − ω2

crφcr(u)] sin (ωcrτ ) + [2pβUωcrφ′

cr(u)] cos (ωcrτ ) = 0

which can be expressed as

γ(u) sin (ωcrτ + η(u)) = 0

(6.14)

η(u) = tan−1 [b(u)/a(u)]

a(u) = φ′′′′

γ(u) =pa2(u) + b2(u),
b(u) = 2pβ U ωcr φ′

cr (u) + U 2 φ′′

cr φcr(u)

cr(u) − ω2
cr(u)

Equation (6.14) represents a traveling wave where the amplitude and wave speed both depend

on u, i.e., the location of the point on the beam from the hinged end. These results show that

the nature of oscillations transition from standing waves in the absence of external flow to

traveling waves in the presence of external flow. Also the mechanism of flutter instability is

different. In the absence of external flow, two successive eigenfrequencies becoming identical

at the point of instability; this is known as double mode flutter - see Fig.6.5. In contrast, in

the presence of external flow, a single complex eigenfrequency crosses the imaginary axis at

the point of instability; this is known as single mode flutter - see Fig.6.7.

)
ω
(
m
I

ωcr = 42.8i (C cr = −1.63)

ω2 = −1.34 + 49.3i (C = 0)

-1.0

-0.5

ω2 = −1.34 − 49.3i (C = 0)

30

10

-10

-30

0.5

Re(ω)

ωcr = −42.8i (C cr = −1.63)

Figure 6.7: Root locus of the second mode complex conjugate frequency for the hinged beam
subjected to external folw as C varies negatively for k = α = 0.5, U e = 0.7, and β = 0.9.

89

The fluid-structure interaction problem is extended by considering the fluid-immersed

hinged beam to be affixed to a rigid body.

6.3 Fluid-Immersed Hinged Beam Affixed to a Rigid Body

6.3.1 Mathematical model

y

L

x bx

y(x, t)

Ue

ℓ

M

M, J

k

x

Curvature : y′′(bx, t)

Figure 6.8: Fluid submersible tail mechanism affixed to a rigid body.

Consider the Euler-Bernoulli beam of length L and cross-sectional area A, shown in Fig.6.8,

which is hinged to a rigid body of a mass and mass moment of inertia M and J, respectively.

The beam is subjected to a dynamic moment M at the hinged joint, which has a torsional

stiffness k, and is submerged in a fluid with a constant relative velocity Ue. The dynamic

moment is proportional to the curvature of the beam at x =bx. The transverse displacement

of a point on the beam at a distance x from the hinged joint is denoted by y(x, t). For small

deformation in absence of the gravity, viscosity, externally imposed tension and pressurization

effect, the equation of motion of the beam can be written as follows [69]:

EI

∂4y(x, t)

∂x4 + MU 2

e

∂2y(x, t)

∂x2 + 2MUe

∂2y(x, t)

∂x∂t

+ (M + ρA)

∂2y(x, t)

∂t2

= 0

(6.15)

where E, I, and ρ are the Young’s modulus of elasticity, area moment of inertia, and the

mass per unit volume of the beam. We introduce the following change of variables

90

u =

x
L

, v =

y
L

, τ = ts EI

ρAL4 , α = bx

L

, ue =r M

EI

UeL, β =

M

M + ρA

The non-dimensional equation of motion is

v ′′′′(u, τ ) + u2

ev ′′(u, τ ) + 2u2

epβ ˙v ′(u, τ ) + ¨v(u, τ ) = 0

(6.16)

where v ′ and ˙v denote the partial derivatives of v(u, τ ) with respect to u and τ , respectively.

We now make the following assumptions:

y

V

M

M

V

ℓ

φ
θ0

x

φ
θ0

y

x

Figure 6.9: Free-body diagrams of the tail mechanism and the rigid body.

• The submersible mechanism is assumed to be traveling with constant forward velocity

through the immersing medium along the negative x-axis.

• The rigid body is symmetric about the plane containing the neutral surface of the

undeformed beam and has a small angle rotation denoted by φ - see Fig.6.9.

The above assumptions yield the following boundary conditions at x = 0:

EI

∂2y(0, t)

∂x2 + Mℓ

∂2y(0, t)

∂t2 −

EI

∂3y(0, t)

∂x3 + M

∂2y(0, t)

∂t2 −

J + Mℓ2

J + Mℓ2 −

Mℓ

J + Mℓ2 −

k
ω2
n

k
ω2
n

(cid:16)Mℓ
(cid:16)Mℓ

∂2y(bx, t)
∂2y(bx, t)

∂y(0, t)

∂x (cid:17) = 0
∂x (cid:17) = 0

∂y(0, t)

∂2y(0, t)

∂t2 + C

∂x2 + k

∂2y(0, t)

∂t2 + C

∂x2 + k

91

where ℓ is the distance from the hinged joint to the center of mass of the rigid body. These

two conditions can be derived from the free-body diagram shown in Fig.6.9 as follows:

∂2y(0, t)

∂t2 − ℓ

V = −M(
∂2λ
∂t2 + Vℓ
M = J

∂2λ
∂t2 )

where λ = φ + θ0 and

∂2λ
∂t2 =

1

J + Mℓ2 −

(cid:16)Mℓ

k
ω2
n

∂2y(0, t)

∂t2 + C

∂x2 + k

∂2y(bx, t)

∂y(0, t)

∂x (cid:17)

M = C

∂x2 − k(λ − θ0)

∂2y(bx, t)

The boundary conditions at x = L are

EI

∂2y(L, t)

∂x2

= EI

∂3y(L, t)

∂x3

= 0

Due to the second assumption, the force F does not affect the boundary conditions. To this

end, the non-dimensional form of the boundary conditions can be expressed as

v ′′(0, τ ) =

v ′′′(0, τ ) =

µψ2(γ2 + 1)
κ
ǫ + µψ2 −
Ω2
µψ
ǫ + µψ2 −

κ
Ω2

n(cid:16)µψ¨v(0, τ ) + Cv ′′(α, τ ) + κv ′(0, τ )(cid:17) − µψ¨v(0, τ )
n(cid:16)µψ¨v(0, τ ) + Cv ′′(α, τ ) + κv ′(0, τ )(cid:17) − µ¨v(0, τ )

v ′′(1, τ ) = v ′′′(1, τ ) = 0

(6.17)

where Ω is the non-dimensional natural frequencies and

92

µ =

M

(ρA + Me)L

, ψ =

ℓ
L

, C =

C
EI

, κ =

KL
EI

, ǫ =

J

(ρA + Me)L3 ; Me =

1
4

ρf πS2

In our simulations the external fluid is assumed to be water, i.e., ρf = 1000

Kg
m3 and S is the

beam height which is assumed to be 0.06m.

6.3.2 Method of Analysis

To solve Eq.(6.16) along with the boundary conditions in Eq.(6.17), we followed the proce-

dure introduced in [69] and highlighted in [74] by assuming the following separable form of

v(u, τ ):

v(u, τ ) = f (u) g(τ );

g(τ ) = eiΩτ

(6.18)

It is worth noting that Eq.(6.16) is not a self-adjoint for ue 6= 0 and therefore the eigenvalues
are in general complex values. Upon substitution of Eq.(6.18) into Eqs.(6.16)-(6.17) and

dividing through by eiΩτ we get

f ′′′′(u) + u2

ef ′′(u) + 2u2

epβiΩf ′(u) − Ω2f (u) = 0

(6.19)

f ′′(0) =

f ′′′(0) =

µψ2(γ2 + 1)
κ
ǫ + µψ2 −
Ω2
µψ
ǫ + µψ2 −

κ
Ω2

n(cid:16) − µψΩ2f (0) + Cf ′′(α) + κf ′(0)(cid:17) + µψΩ2f (0)
n(cid:16) − µψΩ2f (0) + Cf ′′(α) + κf ′(0)(cid:17) + µΩ2f (0)

f ′′(1) = f ′′′(1) = 0

(6.20)

Equation (6.19) is an ordinary differential equation with constant coefficients, as a result,

the solution of f (u) is assumed to be of the form f (u) = Aezu. Accordingly, Eq.(6.19) can

93

be written as

z4 + u2

ez2 + 2u2

epβiΩz − Ω2 = 0

(6.21)

For specific values of ue and β, the characteristic polynomial, Eq.(6.21), provides four roots

of zn, where n = 1, 2, 3, 4. It is also worth noting that zn is a function of α, µ, ψ, γ, ǫ, C, κ,

and Ω in addition to ue and β. The solution of f (u) is thereupon takes the form

f (u) = A1ez1u + A2ez2u + A3ez3u + A4ez4u

(6.22)

which requires to satisfy the boundary conditions in Eq.(6.20) and yields to



|

δ1

ζ1

δ2

ζ2

δ3

ζ3

δ4

ζ4

z2
1ez1 z2

2ez2 z2

3ez3 z2

4ez4

1ez1 z3
z3

2ez2 z3

3ez3 z3

4ez4

Z

{z

where δn, ζn, n = 1, 2, 3, 4, are defined by



A1

A2

A3

A4

=





0

0

0

0



(6.23)



}

µψΩ4 (ǫ − γ2µψ2) + Ω2z2

n (µψ2 (C (γ2 + 1) eαzn − 1) − ǫ) + k (z2

n + µψΩ2 ((γ2 + 1) ψzn − 1))

δn =

ζn =

Ω2 (Cµψz2

neαzn + ǫµΩ2 − z3

n (ǫ + µψ2)) + k (z3

n + µΩ2(ψzn − 1))

k − Ω2 (ǫ + µψ2)

k − Ω2 (ǫ + µψ2)

A non-trivial solution of Eq.(6.23), det (Z) = 0, results in a transcendental characteristic

equation that has infinite roots in Ω. For different values of ue, α, β, C, µ, ψ, γ, ǫ, κ the tran-

scendental equation can be solved numerically to get the complex frequencies. Once Ω and

zn, n = 1, 2, 3, 4, have been determined, the complete solution of Eq.(6.18) can then be

obtained by substituting Eq.(6.22) into Eq.(6.18) as

94

v(u, τ ) =

4Xn=1

AneznueiΩτ =

4Xn=1

An eRe[zn]u ei(Im[zn]u+Re[Ω]τ ) e−Im[Ω]τ

(6.24)

where the coefficients An, n = 1, 2, 3, 4, can be obtained from the nullspace of the matrix Z

in Eq.(6.23). Equation (6.24) indicates that v(u, τ ) is a product of three exponential terms

where the first term is bounded since u is bounded. The second term is also bounded due

to the imaginary exponent that results in a periodic motion. The last term, however, can

grow unbounded with time if Im[Ω] < 0. In this case, the point at which the sign of the

imaginary part of the complex frequency changes from positive to negative represents the

onset point of flutter instability. The mode at which the beam-like tail, structure, undergoes

flutter instability depends on the system parameters, namely, ue, C, α, β, µ, ψ, γ, ǫ, and κ. In

this work, however, we are going to only focus on the effect of first three parameters, that is

to say (ue, C, α), on the flutter stability; the rest of the parameters are listed in Table 6.1.

Table 6.1: Assumed values for the non-dimensional parameters used in the analysis.

Parameters Assumed Values

β

µ

ψ

γ

ǫ

κ

0.97

0.27

0.44

0.36

0.01

2.00

The force exerted by the hinged beam-like tail on the surrounding fluid can be determined

using the system parameters and Eq.(6.24) evaluated at critical stability point, that is,

(Im[Ω] = 0). At the critical stability point, Eq.(6.24) can be expressed as

v(u, τ ) =

4Xn=1

eRe[zn]u(cid:16)Re[An] cos (Im[zn]u + Re[Ω]τ ) − Im[An] sin (Im[zn]u + Re[Ω]τ )(cid:17)

95

It is important to realize that the sinusoidal terms in the above equation represent a travel-

ing waveform, where the slender beam-like tail passes the wave down along its span with an

amplitude increasing over its length. This waveform can produce a positive thrust, according

to Lighthill [75], if the phase velocity is greater than the speed of the structure relative to

external fluid. In this case

Re[Ω]
Im[zn]

>

ue√β

(6.25)

Equation (6.25) constitutes the condition under which the waveform generates a positive

thrust that produces a forward speed ue. Following [74–76] derivation, the time-average

thrust τ of the harmonic motion is

τ =

1
2

M"n(

∂y
∂t

)2 − U 2
e (

∂2y

∂x2 )2ox=L −n(

∂y
∂t

)2 − U 2
e (

∂2y

∂x2 )2ox=0#

and the average power required to provide the displacement y(x, t) is

P = UeM"n ∂y
∂t(cid:18)∂y

∂t

+ Ue

∂y

∂x(cid:19)ox=L −n∂y

∂t(cid:18) ∂y

∂t

+ Ue

∂y

∂x(cid:19)ox=0#

(6.26)

(6.27)

It is worth mentioning that Eq.(6.26) suggest a substantial thrust can be obtained if the

magnitude of

∂y
∂t

is large. At the same time, a larger magnitude of

∂y
∂t

requires large energy

input according to Eq.(6.27). The efficiency introduced in [75] is given by

η =

τ Ue
P

(6.28)

Equations (6.26-6.28) can be non-dimensionalized using the change of variables introduced

in Section 6.3.1 as follows:

τ ∗ =

τ L2
EI

, P

∗

=

P M

1

2 L3

(EI)

3
2

,

η =

τ ∗ue
P

∗

(6.29)

96

where the time-average of the non-dimensional thrust τ ∗ and power P

∗

, that are evaluated

over one cycle of the harmonic motion, should be sufficient [74]

τ ∗ =

Ωcr

Ωcr

4π Z 2π

0

∗

P

n(cid:2)β ˙v(u, τ )2 − u2
2π Z 2π

Ωcr

=

Ωcr

0

ev ′(u, τ )2(cid:3)u=1 −(cid:2)β ˙v(u, τ )2 − u2

ev ′(u, τ )2(cid:3)u=0o dτ

nhueβ ˙v(u, τ )2 − u2
−hueβ ˙v(u, τ )2 − u2

epβv ′(u, τ ) ˙v(u, τ )iu=1
epβv ′(u, τ ) ˙v(u, τ )iu=0o dτ

(6.30)

(6.31)

The results in the next section determine the critical stability points, i.e., the value of Ccr

at a given ue that produces flutter instability.

6.3.3 Flutter Instability

We determine the critical stability points, i.e., Im[Ω] = 0 in the C − ue space for a range
of α ∈ [0.1 − 0.9]. On this account, we start with C = 0 and ue = 0.01 for a specific value
of α to acquire the first ten natural frequencies Ω of the beam. These frequencies are used

as a first guess to find critical stability for the next higher value of ue; the magnitude of C

is gradually increased till one of the Ω’s satisfies the condition Im[Ω] = 0. Then, we trace

in the critical stability points to obtain the critical stability curve by gradually varying the

value of ue and compute the corresponding Ccr for the particular value of α. The critical

stability curve in the C − ue space, is depicted using the automated method proposed by
Hellum [74]. Since the proportionality constant C can be positive or negative, we study these

two cases separately.

Positive Proportionality Constant

Consider the case where C > 0, i.e., the dynamic moment M at the hinge is proportional to

the curvature and the constant of proportionality is positive. The critical stability curve is

97

obtained for each α by calculating Ccr for a range of ue ∈ (0, 20) with a step size of 0.01. The
objective of this section is to study the effect of the point of measurement of the curvature,

α, on the stability. Flutter instability is discussed for four ranges of α in which they are

categorized based on the established pattern of behavior. For instance, the critical stability

curves exhibit an arch-like shape for α ∈ (0.1, 0.38]. As α increases from 0.12 to 0.20, the
curves tend to grow in size and outstretch rightward - see Fig.6.10(a).

10

e
u

0
0.4

α = 0.12
α = 0.14
α = 0.18
α = 0.16
α = 0.20

1.4
C
(a)

11.14

10.09

6
.
0

2
.
1

8
.
1

3
.
2

11.60

(0.77, 79.70)

(0.58, 35.41)

(0.63, 1034.06)

8.33

6.54

e
u

(0.53, 135.02)

2.4

0
0.4

(0.52, 3304.67)

4.00

1.8

1.4
C
(b)

1
9

.

1

0
2

.

2

2.4

Figure 6.10: Critical stability curves in C − ue space: (a) a range of α ∈ (0.10 − 0.20] (b)
for α = 0.20 the efficiency, η, and the critical frequency, Ωcr, are denoted as (η, Ωcr) for five
different points along the critical stability curve.

The critical stability curve for α = 0.20 is detailed in Fig.6.10(b) to include Froude efficiency

η and the critical frequency Ωcr. From the figure, we infer that the required value of C to

maintain flutter instability increases as the external fluid velocity is increased for the region

of ue ∈ (0.0, 11.6]. In the same range, the efficiency improves - see the first three chosen
points in Fig.6.10(b). Beyond this region, the required value of C increases to 2.3, with a

reduction in the associated ue. Past this point, C decreases as the external fluid velocity ue

diminishes and C = 1.8 when ue reaches zero. During this portion of the curve, the efficiency

drops and is coupled with a substantial rise in the critical frequency.

98

For α ∈ [0.22− 0.38] a trend similar to that of α ∈ (0.10− 0.20] is observed - see Fig.6.11(a).
For α = 0.22, the critical stability curve is presented in detail in Fig.6.11(b), mainly to

compare with the case α = 0.20 which was presented in Fig.6.10(b). These figures indicate

that a slight change in the location of α form 0.2 to 0.22 results in a significant change in the

flutter instability characteristics. It is observed that the critical stability curve contracts to a

dome-like shape with a maximum external fluid velocity ue = 3.3. The value Ccr is reduced,

the maximum Froude efficiency is lower, and finally, the frequency flutter instability drops

drastically.

7

e
u

0

0

6

.

0

2

.

1

8

.

1

3.3

(0.45, 24.45)

2.0

e
u

(0.47, 24.76)

(0.50, 31.11)

1.7

α increasing line

3

C
(a)

1.0

0.0

0.4

1.8

C
(b)

Figure 6.11: Critical stability curves in C − ue space: (a) a range of α ∈ [0.22 − 0.38] with
an increment of 0.02 (b) for α = 0.22 the efficiency, η, and the critical frequency, Ωcr, are
denoted as (η, Ωcr) for three different points along the critical stability curve.

For α ∈ [0.40, 0.70], the critical stability curves exhibit two types of behavior, as shown in
Fig.6.12(a). The instability curves have a positive slope as ue approaches 8. Thereafter, for

α ∈ [0.40, 0.46] the critical stability curves rotate clockwise and eventually reach a zero slope.
In contrast, the critical stability curves for α ≥ 0.66 rotate counterclockwise at a higher value
of ue. For an intermediate value of α = 0.58, the result is shown in Fig.6.12(b). The Froude

efficiency and the associate flutter frequency are delineated at three different points along

the critical stability curve. It is interesting to note that the efficiency increases as both ue

99

and C increase, while the flutter frequency fluctuates. At ue = 11.68 and Ccr = 7.00, the

efficiency impressively reaches 0.98 at a relatively lower flutter frequency.

16

e
u

0

0

2
7

.

0

0
7

.

1

0
0

.

7

α = 0.70

α = 0.68

α = 0.66

12

10

(0.77, 28.64)

(0.98, 25.51)

11.68

α = 0.58

e
u

6

(0.46, 17.81)

10

0

0

10

C
(b)

α = 0.46

α = 0.44

α = 0.42

α = 0.40

C
(a)

Figure 6.12: Critical stability curves in C−ue space: (a) critical stability curves are depicted
for eight different values of α ∈ [0.40, 0.70] (b) for α = 0.58 the efficiency, η, and the critical
frequency, Ωcr, are denoted as (η, Ωcr) for three different points along the critical stability
curve.

The final examined region is α ∈ [0.80, 1.00), where the critical stability curve exhibits
multiple stability transitions1 - see shown in the Fig.6.13(a). Notably, for α = 0.90 and

C = 1.18 the beam exhibits up to five stability transitions at different values of ue - see

Fig.6.13(b). These transitions are analogous to the stability transition induced by damping

discussed in Chapter 3. For C = 1.18, the bam is unstable for ue = 0.The beam gets stabilized

at ue = 6.76, then destabilized at ue = 7.55, stabilized again at ue = 13.1, destabilized again

at ue = 15.6, and finally restabilized at ue = 18.5. It should be noted that the generated

thrust is mostly positive except in some regions where the curvature is measured closed to

the free end of the beam. The negative thrust is illustrated using a darker region on the

critical stability curve - see Figs.(6.10 - 6.13)(b).

1The area toward the left side of the critical stability curve represents the stable region where the beam
does not maintain flutter. On the contrary, the area rightward the curve represents flutter instability region.

100

α = 0.90

α = 0.92

20

e
u

α = 0.84

α = 0.80

20

16.98

e
u

10.00

6.76

0
0.0

4.5

0

0

C

(a)

2
7

.

0

8
4

.

0

8
9

.

0

(0.84, 787.3)

18.5

15.6

(0.80, 663.9)

13.1

7.55

7.12

(9E-4, 567.01)

(0.56, 601.8)

1.18

1.55

2

C

(b)

Figure 6.13: Critical stability curves in C−ue space: (a) critical stability curves are depicted
for four different values of α, namely, α = 0.80, 0.84, 0.90, 0.92 (b) for α = 0.90 the efficiency,
η, and the critical frequency, Ωcr, are denoted as (η, Ωcr) for three different points along the
critical stability curve; The darkness region of the critical stability curve depicts the region
of negative thrust.

Negative Proportionality Constant

similar to the case where the constant of proportionality is positive, we investigate the case

of a negative proportionality constant for three regions of α, described below:

region 1: α ∈ (0.1, 0.36]
region 2: α ∈ [0.5, 0.58]
region 3: α ∈ [0.66, 0.98]

In region 1 the critical stability curve resemble the first two cases when C > 0; however,

instead of growing in size, the curves shrink and change shape. For α = 0.40, the Froude

efficiency and the associated flutter frequency was computed at three different locations

along the instability curve. These points are shown in Fig.6.14(b) in which the maximum

efficiency η is less than 0.6.

101

20

8
1
.
0
=
α

α = 0.20

α = 0.24

16

α = 0.28

0
5

.

0
-

1
7

.

1
-

(0.58, 244.2)

e
u

α = 0.34

e
u

8

(0.54, 269.5)

6.00

α = 0.30

α = 0.36

0

0

-7

0
0.0

C
(a)

(0.52, 365.8)

C
(b)

-4.34

-5.5

Figure 6.14: Region 1: critical stability curves in C − ue space: (a) critical stability curves
are depicted for seven different values of α ∈ [0.18, 0.36] (b) for α = 0.30 the efficiency, η,
and the critical frequency, Ωcr, are denoted as (η, Ωcr) for three different points along the
critical stability curve.

20

e
u

0

0

α = 0.52

8
2

.

0
-

20

17.2

(0.79, 126.6)

1
8

.

1
-

6
3

.

2
-

α = 0.50

α = 0.58

11.5

e
u

(0.70, 112.1)

(0.53, 84.6)

6.7

α = 0.62

(0.49, 58.5)

-6

0

0

C

(a)

-4.07

-4

C

(b)

Figure 6.15: Region 2: critical stability curves in C − ue space: (a) critical stability curves
are depicted for three different values of α, namely, α = 0.50, 0.52, 0.58 (b) for α = 0.62 the
efficiency, η, and the critical frequency, Ωcr, are denoted as (η, Ωcr) for three different points
along the critical stability curve.

In region 2 the critical stability curves are similar to both the third and fourth cases with

C > 0. It can be seen that the system exhibits a rich set of flutter instability transitions

102

- see Fig.6.15(a). For α = 0.62 the efficiency and the critical frequency is depicted at four

locations along the flutter instability curve - see Fig.6.15(b). The figure reveals that an

external velocity of 17.2 and Ccr = −1.81 results in a positive thrust and an efficiency of
0.8. This indicates that positive thrust can be obtained using both positive and negative

proportionality constant.

In region 3, the instability curves exhibit a total of three stability transitions - see Fig.6.16(a).

In Fig.6.16(b), we show the special case of α = 0.98. The thrust is negative at some locations

on the critical stability curve for α = 0.98, this is highlighted by the darker region. Although

a positive thrust τ ∗ > 0 is needed to produce an external forward velocity and overcome

external forces such as the drag, negative thrust can be useful for backward motion. At the

two ends of the darker curve, τ ∗ = 0, operating at these points can be useful for hovering

motion and quickly changing the direction of motion.

12

e
u

0

0

5
2
2
-

.

(0.70, 840.0)

12

10.0

e
u

6.4

4.1

(3E-3, 691.2)

(4E-3, 724.0)

-1.19

-6.84

-10

0

0

-10

C

(b)

α = 0.66
α = 0.70
α = 0.80
α = 0.94

α = 0.98

C

(a)

Figure 6.16: Region 3: critical stability curves in C − ue space: (a) critical stability curves
are depicted for five different values of α ∈ [0.66, 0.98] (b) for α = 0.98 the efficiency, η, and
the critical frequency, Ωcr, are denoted as (η, Ωcr) for three different points along the critical
stability curve.

103

Chapter 7

Conclusion

7.1 Divergence and Flutter Modes of Instability

It is well-known that a non-conservative system loses stability through either divergence or

flutter as the boundary load exceeds some critical value. At the same time, a conservative

system loses its stability only through divergence. With this in mind, the two types of

structural instabilities were illustrated using two examples of discrete and continuous systems

in Chapter 2. Furthermore, a classification of flutter instability based on the existence of

damping was highlighted.

7.2 Stability Transitions Induced by Damping

Chapter 3 assessed the role of damping on the dynamic stability of non-conservative systems.

It was found that the initially unstable system, at some fixed value of the non-conservative

force, gets stabilized, then destabilized, then restabilized by increasing the amount of damp-

ing. This counter-intuitive role of damping - which we refer to as the S-D-S behavior -

was observed in every classical non-conservative system considered in this chapter. The

non-conservative systems exhibit different types of stability transition behaviors in different

regions of their damping parameter space. To the best of our knowledge, this instability

behavior has not been described in the literature. This work also highlighted the likelihood

for the behavior to exist for most non-conservative systems. The work also showed the pos-

sibility for two instability curves to cross, a phenomenon which has multiple implications

for the theoretical consideration of stabilization. To accomplish this, it becomes necessary

to delineate the stable region of the parameter space using a piecewise continuous absolute

instability curve. More critically, the term “stabilizing” must be defined in terms of the

104

behavior at a fixed value of the non-conservative forcing variable rather than by the local

slope of the instability curve as proposed in the literature.

The key to finding the multiple stability transitions, and the S-D-S behavior in particular,

was the application of the Routh-Hurwitz criterion twice, back-to-back: the first time for

the characteristic polynomial of the system and the second time for the polynomial that

guarantees the existence of a second-order auxiliary polynomial in the Routh array. For

a given system, the number of stability transitions can be directly determined for a set of

system parameters. This method works by determining the number of physically realizable

damping coefficients, regardless of their magnitude, at which the system is marginally stable.

It is, therefore, possible to demonstrate the existence of the S-D-S behavior in a system that

had been investigated earlier in the literature. However, the final restabilization occurs at

a damping level higher than the range examined in that work. The overlooked stability

transition does not always occur at very high levels of damping. This analysis has been

undertaken for each of the S-D-S transitions in this work, and similar observations have

been made.

The two degree-of-freedom linkage system subjected to a follower end force was also

investigating at greater length by varying the damping level in one of the joints while keeping

the damping in the other joint fixed. Through the analysis, it was shown that there could be

only one stability transition for some values of the fixed level of damping; for other values of

the fixed level of damping, there can be three stability transitions. The study was extended

to obtain the critical stability curves in the space of the damping parameters for different

values of the non-conservative force. These curves have similar shapes that indicate that an

initially unstable system can be made stable by increasing the damping in the first joint,

or alternately, by increasing or decreasing the damping in the second joint. It is interesting

to note that the maximum number of stability transitions ever found in the investigated

systems was three, corresponding to the S-D-S behavior. This work outlined the procedure

to determine the number of stability transitions, which is only valid for a particular set of

105

system parameters. Consequently, the existence of additional transitions cannot be ruled

out.

7.3 Non-conservative Effects of Dynamic Terminal Moment

The stability characteristics of a cantilever beam, subjected to a dynamic moment at its free

end, was investigated in Chapter 4. The applied moment is assumed to be proportional to

the slope or the curvature of the beam, measured from some point along its length. The

geometric stiffness associated with the terminal loading imparts a non-symmetric structure

to the stiffness matrix when the measurement is taken from any point other than the free

end. This suggests that the terminal loading acts as a non-conservative loading for non-

terminal measurement. Regardless of whether this non-conservative moment is proportional

to the slope or the curvature of the beam, the stability is always lost through divergence

when the proportionality constant is positive. On the contrary, for a negative proportionality

constant, the system loses stability solely through flutter.

At the divergence point, the critical value of the constant of proportionality drops as the

point of measurement shifts from the fixed end to the free end when the terminal moment

corresponded to a positive slope. However, for loading proportional to the positive curvature,

the location of the point of measurement has little to no effect on the critical values. Flutter

instability, on the other hand, is induced at higher modes as the point of measurement shifts

from the fixed end to the free end of the beam. For instance, when the terminal moment

is proportional to the negative slope, the beam loses stability in the first or second flutter

instability mode depending on the location of the point of measurement. Comparably, when

the terminal moment is proportional to the negative curvature of the beam: as the point of

measurement moves towards the free end, the mode of instability changes sequentially from

the first to the fourth flutter instability mode. An interesting observation in both flutter

instability cases is that a continuous change in the location of the point of measurement

can result in multiple stability transitions and induce different modes of instability. This

106

behavior is analogous to the multiple instability transitions introduced in Chapter 3.

Experimental validation was a setup where an electric motor was used to provide a

terminal moment to a steel cantilever beam. The motor was supported by air-bearings

to avoid torsional dynamics. Due to the simplicity of measuring the curvature utilizing a

strain gage, the terminal moment was considered to be proportional to the curvature. The

critical value of the proportionality constant, frequency, and mode shape at the critical point

matched the values predicted using a linear model.

7.4 Stability Transitions Induced by an Intermediate Support

Chapter 5 investigated the critical stability of the non-conservative cantilevered beam pre-

sented in Chapter 4, in the presence of intermediate support. The analyses were carried out

using both numerical (Galerkin approximation) and analytical methods. In order to examine

the effects of the intermediate support, both the locations of the intermediate support and

the point of the slope or curvature measurement were varied along the cantilevered beam.

For the case when the terminal loading is proportional to the slope, the critical stability

is discussed in three separate locations of the intermediate support using the numerical

method. Irrespective of whether the terminal moment is proportional to the positive slope

or the negative slope and independent of the location of support, both divergence and flutter

modes of instability are exhibited. The regions of the point of measurement of the slope

that result in divergence and flutter depend strongly on the location of the intermediate

support. For instance, when the intermediate support is close to the fixed end (conversely,

free end), divergence occurs when the point of measurement of the slope lies to the right

of the support (conversely, left of the support). However, when the intermediate support is

in the middle of the beam, divergence occurs when the point of measurement of the slope

lies near to the support, on either side. Also, unlike the case without intermediate support

presented in Chapter 4, where the flutter instability modes switch sequentially, the flutter

instability modes change randomly as the point of measurement shifts from one end to the

107

other end of the beam.

When the terminal loading is proportional to the positive curvature, five cases are dis-

cussed by keeping one parameter fixed and varying the other parameters. The results reveal

a rich set of instability transitions between flutter and divergence and flutter-to-flutter when

the intermediate support location is shifted from the fixed end to the free end of the beam.

The first case involved one stability transition from divergence to flutter that accompanied

a jump in the critical proportionally constant values. The second case involves multiple

instability transitions between divergence and flutter; however, only the the first one is ac-

companied by a jump in the critical loading. Additionally, transitions between the different

flutter instability modes were random, and therefore no veering was observed. The third case

involved multiple instability transitions, just like the first two cases, but none of the transi-

tions exhibited jumps in the critical loading. The fourth case illustrated a transition between

two random instability modes with no veering but a jump in the critical load. Finally, the

last case included two jumps in the critical loading whereby the first one is associated with a

transition between divergence and flutter, and the second one involved a transition between

two consecutive flutter modes. The flutter to flutter transition results from veering.

For the case where the terminal moment is proportional to the negative curvature, tran-

sitions between divergence and flutter occur precisely along the diagonal of the parameter

space and are accompanied by the jump phenomenon. Transitions between modes of flutter

occur randomly in the demarcated flutter region. Due to the clear demarcation of the regions

of divergence and flutter, this case bears a close resemblance to that of the cantilever beam

with follower force. Although this problem has been well-studied, it was revisited here in

this chapter to highlight the jump phenomenon associated with the instability transitions.

This jump phenomenon, which results from veering, was overlooked in the literature.

108

7.5 Non-conservative Behavior of Hinged Beam with a Dynamic

Moment

Chapter 6 presented a framework for investigating the critical stability of a hinged beam

with a non-conservative moment both in the presence and absence of external flow. The

moment, in this case, is assumed to be proportional to the curvature of the beam, measured

at some point along its length. The stability of the hinged beam was first investigated using

the Galerkin approximation. The instability analysis was carried out both in the absence

and presence of external flow; in both cases the beam dynamics is non-conservative.

In the absence of external flow, the structure lost stability through either divergence

or flutter when the moment is proportional to the positive curvature. On the other hand,

stability is lost solely through flutter when the moment is considered to be proportional to

the negative curvature. The effects of external flow were also investigated using Galerkin

approximation whereby the fluid-structure interaction resulted in flow-induced stiffness and

damping terms. The damping term results in complex eigenfrequencies that alters the mech-

anism of flutter instability from that observed in the absence of external flow. Moreover,

the standing waves generated at the critical stability point in the absence of external flow

converted to traveling waves in the presence of external flow. It is noteworthy that the ex-

istence of external flow reduces the magnitude of the critical moment. However, the nature

of instability (divergence and flutter) does not change when the external flow is introduced.

The fluid-structure interaction investigation is extended by considering the fluid-immersed

hinged beam to be affixed to a rigid body. We focus our attention to exclusively investigating

flutter instability. A thrust induced by the traveling waveform is generated at the critical

instability point. The critical stability of the structure is determined using the analytical

method within a two-parameter space of the applied non-conservative loading and the rela-

tive speed of the submersible structure. Since the proportionality constant can be positive or

negative, the onset of flutter instability is discussed separately. The critical stability curves

in the space of the two considered parameters are depicted; the results showed a similar

109

trend that can be categorized based on the curvature measurement locations. It should be

noted that specific measurement locations can result in a negative thrust, especially when

the point is taken close to the free end of the beam. Also, for some other locations, both

positive and negative constant of proportionality resulted in positive thrust.

110

BIBLIOGRAPHY

111

BIBLIOGRAPHY

[1] T. B. Benjamin, Dynamics of a system of articulated pipes conveying fluid. i. theory, in:
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences, Vol. 261, The Royal Society, 1961, pp. 457–486.

[2] L. Prandtl, Kipperscheinungen, Dissertation, Munich.

[3] V. Reut, About the theory of elastic stability, Proceedings of the Odessa Institute of

Civil and Communal Engineering, No. 1.

[4] G. Herrmann, R. W. Bungay, On the stability of elastic systems subjected to noncon-

servative forces, Journal of Applied Mechanics, ASME 31 (1964) 435–440.

[5] H. Ziegler, Die stabilit¨atskriterien der elastomechanik, Archive of Applied Mechanics 20

(1952) 49–56.

[6] T. B. Benjamin, Dynamics of a system of articulated pipes conveying fluid - i.theory,
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences 261 (1307) (1962) 457–486.

[7] M. Langthjem, Y. Sugiyama, Dynamic stability of columns subjected to follower loads:

A survey, Journal of Sound and Vibration 238 (5) (2000) 809 – 851.

[8] M. P. Paidoussis, G. X. Li, Pipes conveying fluid: A model dynamical problem 7 (1993)

137–204.

[9] E. V. Nikolai, About criterion of stability of elastic systems, Proceedings of the Odessa

Institute of Civil and Communal Engineering, No. 1.

[10] A. Pfluger, Stabilit¨atsprobleme der elastomechanik, Springer-Verlent, Berlin (1950)

P.217.

[11] A. Kounadis, Some new instability aspects for nonconservative systems under follower

loads, International Journal of Mechanical Sciences 33 (4) (1991) 297 – 311.

112

[12] Z. Qiu, S. Nemat-Nasser, Instability of an articulated cantilever induced by an impinging

airjet, International Journal of Solids and Structures 21 (2) (1985) 145 – 154.

[13] I. Elishakoff, Controversy associated with the so-called “follower forces”: Critical

overview, Applied Mechanics Reviews 58 (2005) 117.

[14] H. Ziegler, The stability criteria of elastomechanics, Archive of Applied Mechanics 20 (1)

(1952) 49–56.

[15] C. Semler, H. Alighanbari, M. Pa¨ıdoussis, A physical explanation of the destabilizing

effect of damping, Journal of applied mechanics 65 (3) (1998) 642–648.

[16] P. Gallina, Effect of damping on asymmetric systems, Journal of vibration and acoustics

125 (3) (2003) 359–364.

[17] Y. Sugiyama, M. Langthjem, Physical mechanism of the destabilizing effect of damping
in continuous non-conservative dissipative systems, International Journal of Non-Linear
Mechanics 42 (1) (2007) 132–145.

[18] O. Doar´e, Dissipation effect on local and global stability of fluid-conveying pipes, Journal

of Sound and Vibration 329 (1) (2010) 72–83.

[19] A. Luongo, F. D’Annibale, On the destabilizing effect of damping on discrete and con-
tinuous circulatory systems, Journal of Sound and Vibration 333 (24) (2014) 6723–6741.

[20] M. Tommasini, O. N. Kirillov, D. Misseroni, D. Bigoni, The destabilizing effect of
external damping: Singular flutter boundary for the pfl¨uger column with vanishing
external dissipation, Journal of the Mechanics and Physics of Solids 91 (2016) 204–215.

[21] O. Bottema, On the stability of the equilibrium of linear mechanical system, Journal of

applied mathematics and physics 6 (1955) 97–103.

[22] V. Bolotin, T. Lusher, Nonconservative Problems of Theory of Elastic Stability, Elsevier

Science & Technology, 1963.

[23] G. Herrmann, C. Jong, On the destabilizing effect of damping in nonconservative elastic

systems, Journal of Applied Mechanics 32 (3) (1965) 592–597.

113

[24] S. Nemat-Nasser, G. Herrmann, Some general considerations concerning the destabiliz-
ing effect in nonconservative systems, Journal of applied mathematics and physics 17
(1966) 305–313.

[25] S. Nemat-Nasser, S. Prasad, G. Hermann, Destabilizing effect of velocity-dependent
forces in nonconservative continuous systems., AIAA Journal 4 (7) (1966) 1276–1280.

[26] R. Gregory, M. Paidoussis, Unstable oscillation of tubular cantilevers conveying fluid. i.
theory, in: Proceedings of the Royal Society of London A: Mathematical, Physical and
Engineering Sciences, Vol. 293, The Royal Society, 1966, pp. 512–527.

[27] G. Done, Damping configurations that have a stabilizing influence on nonconservative

systems, International Journal of Solids and Structures 9 (2) (1973) 203–215.

[28] M. Pa¨ıdoussis, Dynamics of tubular cantilevers conveying fluid, Journal of mechanical

engineering sciences 612 (2) (1970) 85–103.

[29] I. Lottati, A. Kornecki, The effect of an elastic foundation and of dissipative forces on
the stability of fluid-conveying pipes, Journal of Sound and Vibration 109 (2) (1986)
327–338.

[30] R. Aithal, G. Gipson, Instability of internally damped curved pipes, Journal of engi-

neering mechanics 116 (1) (1990) 77–90.

[31] B. N. Rao, G. V. Rao, Stability of tapered cantilever columns subjected to a tip-
concentrated follower force with or without damping, Computers & Structures 37 (3)
(1990) 333–342.

[32] V. Bolotin, A. Grishko, M. Y. Panov, Effect of damping on the postcritical behaviour
of autonomous non-conservative systems, International journal of non-linear mechanics
37 (7) (2002) 1163–1179.

[33] M. Paak, M. Pa¨ıdoussis, A. Misra, Influence of steady viscous forces on the non-linear
behaviour of cantilevered circular cylindrical shells conveying fluid, International Jour-
nal of Non-Linear Mechanics 58 (2014) 167–183.

[34] A. N. Kounadis, On the paradox of the destabilizing effect of damping in non-
conservative systems, International journal of non-linear mechanics 27 (4) (1992) 597–
609.

114

[35] O. Kirillov, A. Seyranian, Stabilization and destabilization of a circulatory system by
small velocity-dependent forces, Journal of sound and vibration 283 (3) (2005) 781–800.

[36] L. Zorii, Y. A. Chernukha, Influence of support conditions on the dynamic stability of

elastic column, Prikl. Mekh 7 (2) (1971) 134–136.

[37] A. Kounadis, Static stability analysis of elastically restrained structures underfollower

forces, AIAA Journal 18 (4) (1980) 473–476.

[38] A. Kounadis, Divergence and flutter instability of elastically restrained structures under

follower forces, International Journal of Engineering Science 19 (4) (1981) 553–562.

[39] A. Kounadis, The existence of regions of divergence instability for nonconservative sys-
tems under follower forces, International Journal of Solids and Structures 19 (8) (1983)
725–733.

[40] A. Kounadis, A. P. Economou, The effects of the joint stiffness and of the constraints
on the type of instability of a frame under a follower force, Acta Mechanica 36 (3-4)
(1980) 157–168.

[41] I. Elishakoff, J. Hollkamp, Computerized symbolic solution for a nonconservative system
in which instability occurs by flutter in one range of a parameter and by divergence in
another, Computer methods in applied mechanics and engineering 62 (1) (1987) 27–46.

[42] I. Elishakoff, I. Lottati, Divergence and flutter of nonconservative systems with interme-
diate support, Computer methods in applied mechanics and engineering 66 (2) (1988)
241–250.

[43] C. Pierre, Mode localization and eigenvalue loci veering phenomena in disordered struc-

tures, Journal of Sound and Vibration 126 (3) (1988) 485 – 502.

[44] Y. Sugiyama, Y. Tanaka, T. Kishi, H. Kawagoe, Effect of a spring support on the
stability of pipes conveying fluid, Journal of Sound and Vibration 100 (2) (1985) 257–
270.

[45] H. Lee, Divergence of a rotating shaft with an intermediate support and conservative
axial loads, Computer methods in applied mechanics and engineering 110 (3-4) (1993)
317–324.

115

[46] W. Glabisz, Vibration and stability of a beam with elastic supports and concentrated
masses under conservative and nonconservative forces, Computers & structures 70 (3)
(1999) 305–313.

[47] H.-P. Lin, S. Chang, Free vibration analysis of multi-span beams with intermediate

flexible constraints, Journal of sound and vibration 281 (1-2) (2005) 155–169.

[48] H.-Y. Lin, Y.-C. Tsai, Free vibration analysis of a uniform multi-span beam carrying
multiple spring–mass systems, Journal of Sound and Vibration 302 (3) (2007) 442–456.

[49] M. De Rosa, C. Franciosi, The influence of an intermediate support on the stability be-
haviour of cantilever beams subjected to follower forces, Journal of Sound and Vibration
137 (1) (1990) 107–115.

[50] M. De Rosa, Eigenvalues’ behaviour of continuous beams on elastic supports in the
presence of damping and follower forces, International Journal of Mechanical Sciences
33 (5) (1991) 325–337.

[51] H. Lee, Effects of damping on the dynamic stability of a rod with an intermediate spring

support subjected to follower forces, Computers & structures 60 (1) (1996) 31–39.

[52] S. Park, J. Chung, Dynamic analysis of an axially moving finite-length beam with
intermediate spring supports, Journal of Sound and Vibration 333 (24) (2014) 6742–
6759.

[53] S. Caddemi, I. Cali`o, F. Cannizzaro, Influence of an elastic end support on the dynamic
stability of beck’s column with multiple weak sections, International Journal of Non-
Linear Mechanics 69 (2015) 14–28.

[54] H. Unterweger, A. Taras, S. Loschan, M. Kettler, Influence of imperfections on the
stability of beams with intermediate flexible supports, Journal of Constructional Steel
Research 136 (2017) 140–148.

[55] V. Bolotin, N. Zhinzher, Effects of damping on stability of elastic systems subjected
to nonconservative forces, International Journal of Solids and Structures 5 (9) (1969)
965–989.

[56] H. Lee, Divergence and flutter of a cantilever rod with an intermediate spring support,

International Journal of Solids and Structures 32 (10) (1995) 1371 – 1382.

116

[57] C. L. Phillips, J. M. Parr, Feedback control systems, Pearson, 2011.

[58] A. Hellum, R. Mukherjee, A. J. Hull, Flutter instability of a fluid-conveying fluid-
immersed pipe affixed to a rigid body, Journal of Fluids and Structures 27 (7) (2011)
1086 – 1096.

[59] A. Hellum, R. Mukherjee, A. B´enard, A. J. Hull, Modeling and simulation of the dy-
namics of a submersible propelled by a fluttering fluid-conveying tail, Journal of Fluids
and Structures 36 (2013) 83 – 110.

[60] A. B. Rostami, M. Armandei, Renewable energy harvesting by vortex-induced motions:
Review and benchmarking of technologies, Renewable and Sustainable Energy Reviews
70 (2017) 193–214.

[61] M. P. Paidoussis, Fluid-structure interactions: slender structures and axial flow, Vol. 1,

Academic press, 1998.

[62] A. Singh, R. Mukherjee, K. Turner, S. Shaw, Mems implementation of axial and follower

end forces, Journal of Sound and Vibration 286 (3) (2005) 637 – 644.

[63] L. Meirovitch, Computational Methods in Structural Dynamics, Sijthoff and Noordhoff,

Rockville, MD, 1980.

[64] I. G. Boruk, L. G. Lobas, On the motion of a reversible double simple pendulum with

tracking force, International Applied Mechanics 35 (7) (1999) 745–750.

[65] Z. Wang, D. R. Guo, Special functions, World Scientific, Singapore, 1999.

[66] A. Di Egidio, A. Luongo, A. Paolone, Linear and non-linear interactions between static
and dynamic bifurcations of damped planar beams, Journal of Non-Linear Mechanics
42 (1) (2007) 88–98.

[67] R. R. Craig, Structural Dynamics, John Wiley and Sons, New York, NY, 1981.

[68] M. Abdullatif, R. Mukherjee, A. Hellum, Stabilizing and destabilizing effects of damping
in non-conservative systems: Some new results, Journal of Sound and Vibration 413
(2018) 442–455.

117

[69] M. Paidoussis, Fluid-Structure Interactions: Slender Structures and Axial Flow, Vol. 2,

Academic Press, 2014.

[70] A. Luongo, F. D’Annibale, M. Ferretti, Hard loss of stability of Ziegler’s column with

nonlinear damping, Meccanica 51 (11) (2016) 2647–2663.

[71] V. Zamani, E. Kharazmi, R. Mukherjee, Asymmetric post-flutter oscillations of a can-
tilever due to a dynamic follower force, Journal of Sound and Vibration 340 (2015)
253–266.

[72] H. Lee, Divergence and flutter of a cantilever rod with an intermediate spring support,

International Journal of Solids and Structures 32 (10) (1995) 1371 – 1382.

[73] M. Abdullatif, R. Mukherjee, Divergence and flutter instabilities of a cantilever beam
subjected to a terminal dynamic moment, Journal of Sound and Vibration 455 (2019)
402–412.

[74] A. Hellum, R. Mukherjee, A. J. Hull, Flutter instability of a fluid-conveying fluid-
immersed pipe affixed to a rigid body, Journal of Fluids and Structures 27 (7) (2011)
1086–1096.

[75] M. Lighthill, Note on the swimming of slender fish, Journal of fluid Mechanics 9 (2)

(1960) 305–317.

[76] T. Wu, Hydromechanics of swimming propulsion. part 1. swimming of a two-dimensional
flexible plate at variable forward speeds in an inviscid fluid, Journal of Fluid Mechanics
46 (2) (1971) 337–355.

118