applied

sciences
Article
Study of Nonlinear Aerodynamic Self-Excited Force in Flutter
Bifurcation and Limit Cycle Oscillation of Long-Span
Suspension Bridge
Jieshan Liu 1,2,3 , Fan Wang 1,2, * and Yang Yang 1,2

1 School of Mechanics and Construction Engineering, Jinan University, Guangzhou 510632, China;
[email protected] (J.L.); [email protected] (Y.Y.)
2 Key Laboratory of Disaster Forecast and Control in Engineering, Ministry of Education,
Guangzhou 510632, China
3 CCCC Fourth Harbor Engineering Institute Co., Ltd., Guangzhou 510230, China
* Correspondence: [email protected]

Abstract: This article establishes a nonlinear flutter system for a long-span suspension bridge, aiming
to analyze its supercritical flutter response under the influence of nonlinear aerodynamic self-excited
force. By fitting the experimental discrete values of flutter derivatives using the least squares method,
a polynomial function of flutter derivatives with respect to reduced wind speed is obtained. Flutter
critical value is determined by the linear matrix eigenvalues of a state-space equation. The occurrence
of a supercritical Hopf bifurcation in the nonlinear system is determined by the Jacobian matrix
eigenvalues of the state-space equation and the system’s vibrational response at the critical state. The
vibrational response of the supercritical state is obtained through Runge–Kutta integration, revealing
the presence of stable limit cycle oscillation (LCO) and unstable limit cycle oscillation in the system,
and through analyzing the relationship between the LCO amplitude and wind speed. Considering
cubic nonlinear damping and stiffness, the effects of different factors on the nonlinear flutter system
are analyzed.

Citation: Liu, J.; Wang, F.; Yang, Y. Keywords: flutter bifurcation; nonlinear aerodynamic self-excited force; nonlinear damping; nonlinear
Study of Nonlinear Aerodynamic stiffness; supercritical Hopf bifurcation; limit cycle oscillation (LCO)
Self-Excited Force in Flutter
Bifurcation and Limit Cycle
Oscillation of Long-Span Suspension
Bridge. Appl. Sci. 2023, 13, 10272.
1. Introduction
https://doi.org/10.3390/
app131810272 Since the issue of bridge flutter was put forward, numerous studies have been con-
ducted on aerodynamic vibrations in bridges. Flutter, as a form of aerodynamic instability,
Academic Editor: Jean-Jacques
is the most destructive type of vibration induced by wind. When wind speed is low,
Sinou
structural vibrations caused by the interaction between the structure and the flow field
Received: 1 July 2023 are dissipated by system damping, keeping the bridge in a stationary state. However,
Revised: 7 September 2023 when wind speed exceeds a critical value, the motion energy generated by the interac-
Accepted: 10 September 2023 tion between the structure and the flow field cannot be completely dissipated by system
Published: 13 September 2023 damping, resulting in flutter. Flutter is a self-excited vibration phenomenon of elastic
structures induced by their interaction with fluid in the flow field. This study investigates
the nonlinear flutter response when wind speed exceeds the critical value. The occurrence
of self-limiting flutter, known as limit cycle oscillations in bridges, has been indicated by
Copyright: © 2023 by the authors.
the Tacoma Narrows Bridge flutter event and recent research and experiments. This type of
Licensee MDPI, Basel, Switzerland.
limit cycle oscillation, often referred to as “hard flutter,” differs from linear divergence-type
This article is an open access article
flutter theories and does not always lead to divergence [1–5]. Notably, nonlinearity in the
distributed under the terms and
aerodynamic self-excited forces of bluff body bridges is significant [6–9], and the nonlinear
conditions of the Creative Commons
phenomenon stands out under angle of attack conditions [10,11]. Moreover, there exists
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
a significant potential for advancing our assessment of the vibrational performance of
4.0/).
bridges through machine learning algorithms based on artificial neural networks. These

Appl. Sci. 2023, 13, 10272. https://doi.org/10.3390/app131810272 https://www.mdpi.com/journal/applsci

Appl. Sci. 2023, 13, 10272 2 of 19

algorithms employ geometric information and dynamic parameters as inputs to construct

vibrational models of bridges. The integration of computational fluid dynamics (CFD) with
machine learning models is used to predict wind loads, and it has demonstrated satisfactory
results in forecasting and evaluating the aerodynamic elastic responses of bridges [12–14].
This underscores the widespread attention garnered by the study of nonlinear flutter in
long-span bridges. Although bridge design still does not allow for situations beyond the
flutter critical wind speed, studying flutter behavior beyond the critical state contributes to
a deeper understanding of the structural responses of bridges under flutter conditions.
The flutter instability process caused by aeroelastic effects is generally analyzed using
the aerodynamic self-excitation theory proposed by Scanlan et al. [15]. Traditional studies
have established linear flutter equations based on this theory, but the limitations of linear
systems are evident. They can only solve the flutter critical state of the aerodynamic self-
excited system. For example, Náprstek et al. proposed stability conditions and unified
linear variables for bluff bridge deck sections under unsteady flow [16,17]. However, as the
span of suspension bridges increases, the vibration response of bridges after flutter becomes
increasingly important. Therefore, it is crucial to establish nonlinear motion equations for
self-excited flutter in suspension bridges.
In the conventional Scanlan aerodynamic self-excited force model, the interaction
between the bridge cross-section and aerodynamic forces is represented by a linear multi-
plication of flutter derivatives and motion components. The higher-order terms of flutter
derivatives are disregarded. However, the aerodynamic self-excited force is inherently non-
linear. Considering the critical post-flutter of bridges, it becomes essential to formulate the
appropriate equations to describe this nonlinear motion. Numerous scholars have delved
into the study of nonlinear flutter. Building upon Scanlan’s research, they have demon-
strated a strong correlation between the flutter derivatives related to torsional motion and
nonlinearity. They have determined a functional relationship between nonlinear flutter
derivatives, reduced wind speed, and flutter amplitudes [2,11]. Furthermore, through wind
tunnel experiments with box girder models, they have proposed a nonlinear aeroelastic self-
excited force model, coupling bridge lift and torsion [18]. Additionally, investigations into
complex nonlinear fluid–structure interaction phenomena in turbulent flow around bluff
bodies have been conducted through wind tunnel tests and numerical simulations [19–21].
It has been found that the presence of a strong shear layer flow separation contributes
significantly to the nonlinearity of aerodynamic self-excited forces, manifested in variations
in the amplitude of the first harmonic and the appearance of higher harmonics. As for
the nonlinear motion system of bridge flutter, researchers have focused on the critical
bifurcation state and post-bifurcation behavior. Research has mainly centered on analyzing
various mechanical and aerodynamic characteristics (such as mechanical damping ratio,
natural frequency, initial angle of attack, and flutter derivatives) and their influence on the
motion behavior of a given cross-section after the critical flutter state [22–24]. Moreover,
the critical state between static and flutter behavior of nonlinear systems is represented
by Hopf bifurcations [25]. Measures to suppress the occurrence of flutter motion require
attention, primarily in terms of optimizing and controlling the placement of tuned mass
dampers [26,27]. However, previous studies have often concentrated on specific local
aspects of nonlinear motion systems, such as establishing wind tunnel models correspond-
ing to bridge structures to analyze flutter motion states, exploring the acquisition and
identification of high-order flutter derivatives in nonlinear aerodynamic self-excited forces,
and investigating effective energy dissipation and vibration suppression measures post
flutter. There remains a relatively limited amount of research focused on establishing
nonlinear motion equations for bridge flutter systems and studying critical flutter states
and post-flutter motion behaviors.
Based on the status of current research, it is evident that nonlinear aerodynamic self-
excitation exists, and the discrete values of the nonlinear components of flutter derivatives
can be obtained through wind tunnel section model tests. Therefore, this study considers
the effect of nonlinear aerodynamic self-excited force on bridges, establishes nonlinear
Appl. Sci. 2023, 13, 10272 3 of 19

flutter equations for long-span suspension bridges, and fits the function expression of
flutter derivatives with respect to reduced wind speed based on the discrete values of
flutter derivatives obtained from the literature. Subsequently, the nonlinear bifurcation
and limit cycle oscillation phenomena are investigated based on the motion equations,
analyzing and determining the bifurcation types and studying the properties of LCO, as
well as predicting the development of the relationship between LCO and wind speed.
Finally, on this basis, the effects of nonlinear damping and structural nonlinear stiffness on
flutter LCO are further considered.

2. Framework of Nonlinear Flutter Analysis

In this section, motion equations for a long-span suspension bridge under the influence
of nonlinear aerodynamic self-excited force are established. Flutter equations are obtained
by fitting the discrete values of flutter derivatives in the nonlinear aerodynamic force model
using the least squares method. The framework also considers the nonlinear stiffness and
nonlinear damping of the structure, and four computational analysis cases are set up.

2.1. Flutter Motion Equations for Suspension Bridges

The current research on flutter utilizes the motion expression proposed by Scanlan [28].
This model considers bridges as Euler–Bernoulli beams, thereby neglecting shear defor-
mations. Furthermore, it assumes that the plane (cross-section), which is perpendicular
to the beam’s central axis before deformation, remains planar after deformation, and that
regardless of pre- and post-deformation processes, the plane of the cross-section remains
perpendicular to the beam’s axis. Because the centroid of the bridge’s cross-section co-
incides with its shear center, the vibration equation for the beam, derived based on the
Euler–Bernoulli beam theory and the method of separation of variables, can be expressed in
terms of its natural modes, mass, and stiffness. The first-mode response of the bridge con-
tributes significantly to its vibration, resulting in a two-degree-of-freedom cross-sectional
motion model consisting of a first-order vertical bending and torsional frequencies of the
bridge, along with their corresponding masses and stiffness.
As shown in Figure 1, the cross-section of the bridge deck of the Nansha Bridge over
the Nizhou Channel is allowed to undergo motion in the vertical bending and torsional
degrees of freedom. It experiences aerodynamic self-excitation forces in a uniform flow
field, including aerodynamic lift and aerodynamic torque, acting at the centroid of the
cross-section [15,28]. Therefore, flutter motion equations can be expressed as follows. On
the left side of the equation lies the motion equation representing the characteristics of the
bridge, while on the right side lies the aeroelastic self-exciting force exerted on the bridge:
. . !
.. . h Bα h
mh h + 2ξh mh ωh h + mh ω2h h = ρU 2 B KH1∗ + KH2∗ + K2 H3∗ α + K2 H4∗ (1)
U U B
. . !
.. . h Bα h
I α + 2ξ α Iωα α + Iωα2 α 2 2
= ρU B KA1∗ + KA2∗ + K2 A3∗ α + K2 A4∗ (2)
U U B

where m and I are mass and moment of inertia per unit length of the bridge deck, respec-
. ..
tively. h, h, and h represent the displacement, velocity, and acceleration of the vertical
. ..
bending degree of freedom, while α, α, and α represent the angular displacement, angular
velocity, and angular acceleration of the torsional degree of freedom. ξ h and ξ α are the
damping ratios of the suspension bridge in the vertical bending and torsional modes,
respectively. ωh and ωα are the first-order circular frequencies of vertical bending and
torsion, respectively. On the right-hand side of Equation (1) is the aerodynamic lift, and on
the right-hand side of Equation (2) is the aerodynamic torque. ρ is air density, U is fluid
flow velocity, B is the width of the cross-section of the stiffened beam, and K is reduced
frequency, defined as K = Bω/U. Hi∗ and Ai∗ are eight flutter derivatives, which are
 Appl. Sci. 2023, 13, x FOR PEER REVIEW 4 of 20

Appl. Sci. 2023, 13, 10272 𝑈 is fluid flow velocity, 𝐵 is the width of the cross-section of the stiffened beam, 4and of 19𝐾
∗ ∗
is reduced frequency, defined as 𝐾 = 𝐵𝜔/𝑈 . 𝐻 and 𝐴 are eight flutter derivatives,
which are dimensionless quantities related to the reduced frequency or reduced wind
dimensionless
speed and canquantities
be expressedrelated to the reduced frequency or reduced wind speed and can
as functions.
be expressed as functions.
𝐻 ∗ , 𝐴∗ = 𝑓 𝑉 = 𝑓 𝐾 (3)
Hi∗ , Ai∗ = f (Vr ) = f (K ) (3)
where 𝑉 is the reduced wind speed, defined as 𝑉 = 𝑈/𝐵𝑓 . This study focuses on the
where Vr is the reduced wind speed, defined as Vr = U/B f . This study focuses on the
NanshaBridge
Nansha Bridgeover
overthetheNizhou
NizhouChannel,
Channel,which whichisisa asuspension
suspensionbridge
bridgewith
witha amain
mainspan
span
of 1688 m and a bridge deck width of 41.7 m. The relevant design parameters
of 1688 m and a bridge deck width of 41.7 m. The relevant design parameters of the bridge of the bridge
arelisted
are listedininTable
Table1 1[29].
[29].

Figure1.1.Motion
Figure Motionmodel
modelofofbridge
bridgedeck
decksection.
section.

Table 1. Parameter values of the suspension bridge.

Table 1. Parameter values of the suspension bridge.
Parameter Value
Parameter Value
𝑚 65,649 (kg/m)
m 65,649 (kg/m)
𝐼 7,808,150 (kg·m)
I 7,808,150 (kg·m)
B 𝐵 41.7(m)
41.7 (m)
ξ h𝜉 0.005
0.005
ξ α𝜉 0.005
0.005
0.45497 (rad/s)
ωh
𝜔 0.45497 (rad/s)
ωα 1.33285 (rad/s)
𝜔 1.33285 (rad/s)

2.2. Cubic Damping and Cubic Torsional Stiffness

2.2. Cubic Damping and Cubic Torsional Stiffness
Taking into account the structural nonlinearity of the bridge’s stiffness and the pos-
Taking into account the structural nonlinearity of the bridge’s stiffness and the pos-
sibility of obtaining strong nonlinear damping using tuned mass dampers, the vibration
sibility of obtaining strong nonlinear damping using tuned mass dampers, the vibration
equations of the abovementioned flutter motion Equations (1) and (2) are expressed in the
equations of the abovementioned flutter motion Equations (1) and (2) are expressed in the
following form:
following form: .. .
m h h + Ch h + k h h = F ( L ) (4)
𝑚 ℎ+𝐶 ℎ +𝑘 ℎ=𝐹 𝐿 (4)
.. .

Iα + 𝐼𝛼C+ α 𝐶α +𝛼 k+ α) =
α (𝑘 𝛼 F=( M𝐹 )𝑀 (5)
(5)
. .

where
where Ch𝐶 h ℎ and Cα𝐶α 𝛼represent
and represent thethenonlinear damping
nonlinear damping values of the
values offlutter system
the flutter in thein
system
the vertical
vertical bending
bending and torsional
and torsional degrees degrees
of freedom, k α (α)𝑘is structural
of freedom, 𝛼 is structural nonlinear
nonlinear stiff-
stiffness,
ness, and 𝐹 𝐿 and 𝐹 𝑀 are self-excited forces
and F ( L) and F ( M) are self-excited forces of flutter. of flutter.
Thestrong
The strong nonlinear
nonlinear damping
damping in the
in the vertical
vertical bending
bending and and torsional
torsional degrees
degrees of free-
of freedom
.
isdom is expressed
expressed as a function
as a function of verticalof vertical
velocityvelocity ℎ and velocity
h and angular angular α;velocity 𝛼; then,
then, the verticalthe
vertical
and and damping
torsional torsional damping
value ch andvalue 𝑐 and
cα are 𝑐 as
written arefollows:
written as follows:
𝐶.  ℎ = 𝑐. ℎ + 𝑐. 3 ℎ (6)
Ch h = ch1 h + ch2 h (6)
𝐶 𝛼 =𝑐 𝛼+𝑐 𝛼 (7)
.
. .3
Cα αcoefficients
In the equation, linear damping α2 α given by 𝑐 = 2𝜉 𝑚 𝜔 and 𝑐 (7)=
= cα1 α + care
2𝜉 𝐼𝜔 .
In the equation, linear damping coefficients are given by ch1 = 2ξ h mh ωh
and cα1 = 2ξ α Iωα .
Appl. Sci. 2023, 13, 10272 5 of 19

Considering the nonlinear stiffness in the torsional degree of freedom of the suspension
bridge, the bending and torsional stiffness are expressed as follows:

k h h = mh ωh2 h (8)

k α (α) = k α1 α + k α2 α3 (9)
In the equation, stiffness k h in the vertical degree of freedom is set as a linear function
of vertical displacement h. k α , k α1 represent the first-order coefficient of torsional stiffness,
while k α1 = Iωα2 , and k α2 represent the cubic term coefficient of torsional stiffness.
The cubic damping coefficient and cubic torsional coefficient are determined in Table 2.

Table 2. Values of cubic damping and stiffness terms for each case.

Cases Name ch2 cα2 kα2

Reference case 0 0 0
Case A 0 10cα1 0
Case B 0 10cα1 1000mh ωα2
Case C 10ch1 10cα1 1000mh ωα2

2.3. Expression of Nonlinear Aerodynamic Self-Excited Force

Based on the aerodynamic self-excited force model proposed by Scanlan (Equations (1) and (2)),
Gao et al. proposed a nonlinear aerodynamic self-excited force model [18]. The expressions
for the nonlinear self-excited forces F ( L) and F ( M ) in this model are as follows:
" . . ! . #
h Bα Bα h
F ( L) = ρU 2 B KH1∗ KH2∗ + KH2,02
∗
+ K2 H3∗ + K2 H3,02
∗
|α| α + K2 H4∗ + K2 H4,01
∗ ∗
|α| + K2 H4,02 α2 (10)


+
U U U B

" . . .2! . #
2
∗ h ∗ B α ∗ B α Bα h
F ( M) = ρU B 2 2
KA1 + KA2∗ + KA2,02 + KA2,03 + K2 A3∗ α + K2 A4∗ (11)
U U U2 U B

where FL represents aerodynamic lift force, FM represents aerodynamic torque, ρ is air

density, U is fluid velocity, B is the cross-sectional width of the girder, K is the reduced
frequency, K = Bω/U. Hi∗ and Ai∗ are flutter derivatives, and H2,02
∗ , H ∗ , H ∗ , H ∗ , A∗ ,
3,02 4,01 4,02 2,02
∗
A2,03 are small perturbation flutter derivatives. These flutter derivatives are dimensionless
quantities that depend on the reduced frequency or reduced wind speed. Similarly, these
flutter derivatives can be expressed as functions of the reduced wind speed or reduced
frequency. For computational convenience, in this study, these flutter derivatives are
expressed as functions of the reduced wind speed.

Hi∗ , Ai∗ = f (Vr ) = f (K ) = H, Ai1

∗ ∗
· Vr + H, Ai2 · Vr2 (12)

H2∗ ∨ A2∗ = f (Vr ) = f (K ) = Hi1

∗ ∗
∨ Ai1 ∗
· Vr + Hi2 ∗
∨ Ai2 ∗
· Vr2 + Hi3 ∗
∨ Ai3 ∗
· Vr2 + Hi4 ∗
∨ Ai4 · Vr2 (13)

Nonlinear flutter self-excited forces (10) and (11) not only include linear first-order
components but also nonlinear second-order components. In the nonlinear self-excited
force model proposed by Gao et al., the discrete results of flutter derivatives were obtained
through wind tunnel tests on a bridge deck section model [18]. In the flutter motion
equation proposed by Scanlan, the flutter derivatives in the flutter self-excited forces are
only related to the cross-sectional shape of the bridge deck panel. Therefore, it can be
assumed that the numerically obtained results of the box girder section model by Gao
et al. can be used to analyze the nonlinear self-excited forces of the Nansha Bridge over
the Nizhou Channel in this study. To facilitate a continuous solution of the flutter motion
Appl. Sci. 2023, 13, 10272 6 of 19

Appl. Sci. 2023, 13, x FOR PEER REVIEW 6 of 20

equation, the flutter derivatives are fitted as a function of the reduced wind speed using
the least squares method, as shown in the form of (12). The discrete values of the first-order
derivatives and the corresponding
flutter derivatives function graph
and the corresponding are shown
function in Figure
graph are shown 2, in
where the2,
Figure black
where
dots
the black dots represent the discrete flutter derivative values obtained from wind and
represent the discrete flutter derivative values obtained from wind tunnel tests, tunnel
the black
tests, andsolid line represents
the black the fitted flutter
solid line represents derivative
the fitted flutterfunction curve
derivative obtained
function through
curve obtained
the least squares
through the leastmethod.
squares method.

0.4 0.04

0.3 0.02
A1*

A2*
0.00
0.2
-0.02
0.1
-0.04
0.0 0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
Vr Vr

(a) Flutter derivative 𝐴∗ (b) Flutter derivative 𝐴∗

0.6 0.00

-0.02
0.4
A3*

-0.04
A4*

0.2 -0.06

0.0 -0.08
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
Vr Vr

(c) Flutter derivative 𝐴∗ (d) Flutter derivative 𝐴∗

0.0 0.3
-0.5
0.2
-1.0
H2*

-1.5
H1*

0.1
-2.0
-2.5 0.0
-3.0
-0.1
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
Vr Vr

(e) Flutter derivative 𝐻 ∗ (f) Flutter derivative 𝐻 ∗

1 0.0

0 -0.2

-0.4
-1
H4*
H3*

-0.6
-2
-0.8
-3
-1.0
0 1 2 3 4 V5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
r Vr

(g) Flutter derivative 𝐻 ∗ (h) Flutter derivative 𝐻 ∗

Figure
Figure 2.
2. Discrete
Discrete points and fitted
points and fitted curves
curves of
ofthe
thefirst-order
first-orderflutter
flutterderivatives
derivativesofof aero
aero self-excited
self-excited force.
force.
 Appl. Sci. 2023, 13, x FOR PEER REVIEW 7 of 20
Appl. Sci. 2023, 13, 10272 7 of 19

The discrete values of the second-order flutter derivatives are shown in Figure 3,
∗ ∗
whereThe thediscrete
red colorvalues of the the
represents second-order flutter
second-order derivatives
flutter are 𝐻
derivatives , , 𝐻in
shown , 𝐻 ∗, 3,
, Figure , where
∗ ∗ ∗
𝐻the , 𝐴 ,color
, red , 𝐴 ,represents
. The red the
dotssecond-order
represent the flutter
discretederivatives ∗
values obtained
H2,02 ,from∗
H3,02wind ∗
, H4,01 ∗ , A∗ ,
tunnel
, H4,02 2,02
∗ and
tests,
A2,03 . Thetheredred solid
dots line represents
represent the continuous
the discrete functionfrom
values obtained curve obtained
wind tunnel through
tests, and the
fitting.
red solid line represents the continuous function curve obtained through fitting.
1.5
0
1.0 A*2,02 A*2
-5
0.5 A*2
A*2&A*2,02

A*2&A*2,03
0.0 -10
-0.5
-15
-1.0 A*2,03
-1.5 -20
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
Vr Vr
(a) Flutter derivative 𝐴∗ (b) Flutter derivative 𝐴∗
15 6
H*2,02 H*3,02
4

H*3&H*3,02
10
H*2&H*2,02

2 H*3
5 0
H*2 -2
0
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
Vr Vr
(c) Flutter derivative 𝐻 ∗ (d) Flutter derivative 𝐻 ∗
0.5 10
8 H*4,01
H*4,01 6
0.0
H*4&H*4,01

H*4&H*4,02

-0.5
2 H*4
0
H*4 -2
-1.0
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9
Vr Vr
(e) Flutter derivative 𝐻 ∗ (f) Flutter derivative 𝐻 ∗
Figure
Figure3.3.Discrete
Discretepoints andand
points fitted curves
fitted (colored
curves in red)
(colored inof theof
red) second-order flutter derivatives
the second-order of
flutter derivatives of
aero self-excited force.
aero self-excited force.
The coefficients of the function expressions for the continuous curves of flutter deriv-
The coefficients of the function expressions for the continuous curves of flutter deriva-
atives shown in Figures 2 and 3 are listed in Tables 3 and 4.
tives shown in Figures 2 and 3 are listed in Tables 3 and 4.
Table 3. The coefficients in the curve of the fitting quadratic function of the flutter derivative.
Table 3. The coefficients in the curve of the fitting quadratic function of the flutter derivative.
𝑯∗𝒊𝟏 𝑯∗𝒊𝟐 ∗
𝑨𝒊𝟏 𝑨∗𝒊𝟐
∗ * * ∗ *
𝐻 Hi1
0.04374 Hi2
−0.05154 𝐴 0.06288Ai1 −0.00187 A*i2
∗ ∗
𝐻∗
H 0.26574
0.04374 −0.08902
−0.05154 𝐴 A∗ −0.01230
0.06288 0.01217−0.00187
1 1
𝐻3∗∗
H −0.27729
0.26574 0.02632
− 0.08902 𝐴∗A3∗ −0.01665
−0.01230 7.49839 ×0.01217
10−4
𝐻H∗4∗. −−3.76111
0.27729 0.70966
0.02632 𝐴∗ .A4∗ −0.01665 −0.13286
0.95929 7.49839 × 10−4
∗∗ ∗
H𝐻2,02.
−5.81372
3.76111 0.70966
−0.74839 𝐴∗A. 2,02 0.95929
−9.99665 1.13183−0.13286
∗∗ − ∗ −
H𝐻3,02.
5.81372
0.06882 0.74839
−0.00366 A 2,03 9.99665 1.13183
∗
H4,01 0.06882 −0.00366
∗
H4,02 6.52151 −0.82327
Appl. Sci. 2023, 13, 10272 8 of 19

Table 4. The coefficients in the curve of the fitting quartic function of the flutter derivative.

H*i1 ∨A*i1 H*i2 ∨A*i2 H*i3 ∨A*i3 H*i4 ∨A*i4

entry 1 0.29832 −0.17132 0.03441 −0.00224
entry 2 0.00214 −0.00314 3.55170 × 10−4 0

3. Flutter Critical State and Hopf Bifurcation

According to the findings of Andronov et al., if a linear expression of a system’s motion
considers the critical point as a saddle point, a node, or a focus of the critical state, then for
the original nonlinear system, the critical point is also a saddle point, a node, or a focus [30].
The determination of the critical point properties of a nonlinear system can be made by
examining the Jacobian matrix of the motion equations and the motion states in the vicinity
of the critical values. In this section, using the values of the reference case in Table 2, the
calculation methods for the critical state, bifurcation determination, and bifurcation limit
cycle of the nonlinear self-excited flutter system of the suspension bridge are presented.

3.1. Flutter Critical Wind Speed Solution

Flutter Equations (4) and (5) can be written in the form of state-space equations, i.e.,
.
h . . iT
x = f (x). Let x = h, α, h, α ; thus, the state-space equations of this nonlinear flutter
system can be written in the following form:
.
x = F(U, ω, x ) = F Linear · x + F Nonlinear (x) (14)

where F Linear is the linear coefficient matrix of the flutter equations, and F Nonlinear is the
nonlinear term of the flutter equations. Based on the analysis and calculations of the non-
linear system equations in the previous sections, the critical flutter state can be computed
using the linear part of Equation (14), and then the complete state-space equations can be
used to analyze the nonlinear flutter response. The coefficients of the linear and nonlinear
parts of Equation (14) are given as follows:
 
0 0 1 0
 0 0 0 1 
2  2 ∗
 
F Linear = ρU B K H4 2
K H3∗ KH1∗ ∗
KH2 B (15)
kh ch 
 Bmh − mh mh Umh − mh Umh


2
K A4 ∗ 2
BK A3 ∗ BKA1 ∗ 2 ∗
B KA2
I I − kIα UI UI − I
cα

 
0
 0 
. .
 
 ρU 2 B 
∗ B |α|α
2
2 ∗ 2 ∗ 2 ∗

F Nonlinear (x) =  2  (16)
 mh KH2,02 U 2 + K H3,02 |α|α + K H4,01 |α| + K H4,02 α 
.
   . 
 ρU 2 B2 ∗ B|α| ∗
.2
B2 α Bα

I KA2,02 U + KA2,03 U 2 U

The critical point of a nonlinear motion system can be obtained by solving its linear
expression. According to the Theodorsen method, as the wind speed Vr of the system
parameters changes, when the eigenvalues of the matrix (15) have positive real parts, the
corresponding wind speed VCr is the critical value of the system [31]. Therefore, the critical
flutter state can be obtained through the following steps:
1. Assume a small value of frequency ω0 and substitute it into the matrix (15). Grad-
ually increase the reduced wind speed Vr until the first pair of complex conjugate
eigenvalues of the matrix have zero real parts. Record the imaginary part Imag(λ) of
this eigenvalue as frequency ωi = Imag(λ);
2. Take ωi as the new frequency value and substitute it into the matrix (15). Again,
gradually increase the reduced wind speed Vr until the first complex eigenvalue of
Appl. Sci. 2023, 13, 10272 9 of 19

the matrix has zero real parts. Record the imaginary part Imagi (λ) of this eigenvalue
as the new flutter frequency ωi+1 = Imagi (λ).;
3. Compare ωi and ωi+1 . Repeat step 2 until |ωi+1 − ωi | approaches zero. At this point,
the frequency value is the critical flutter frequency ωC of the system, the reduced
wind speed value is the critical reduced wind speed VCr of the flutter, and the flutter
critical wind speed UC can be obtained using the relationship Vr = 2Uπ/Bω .
Through the above steps, the parameter Vr that causes the eigenvalues of the matrix
FLinear to change from negative to positive is obtained, and the flutter system reaches the
critical state. At this critical state, the matrix eigenvalues are given by:

−0.0593 + 0.4754i
 
−0.0593 − 0.4754i 
λ|Vr =7.2673 =
 0.0000 + 1.1888i 


0.0000 − 1.1888i

The flutter system’s critical flutter frequency is ωC = 1.1888 rad/s, the critical reduced
wind speed is VCr = 7.2673, and the critical flutter wind speed is UC = 57.3374 m/s.

3.2. Proof of Hopf Bifurcation

The Jacobian matrix of the state-space Equation (14) is obtained according to the rule
in Equation (17).
 
0 0 1 0
0
 .. 0.. 0.. 1.. 

JF( x) =  ∂h ∂h ∂h. ∂h. 

 (17)
 ∂h.. ∂α.. ∂h.. ∂α.. 
∂α ∂α ∂α. ∂α.
∂h ∂α ∂h ∂α

Based on the results of the critical values of the flutter, taking the reduced wind speeds
Vr = 7, Vr = Vcr = 7.2673, and Vr = 7.5, they are substituted into the Jacobian matrix to
calculate the corresponding eigenvalues, as shown below:

−0.0543 + 0.4782i
 
−0.0543 − 0.4782i 
λ|Vr =7 =
−0.0022 + 1.2001i 


−0.0022 − 1.2001i

−0.0593 + 0.4754i
 
−0.0593 − 0.4754i 
λ|Vr =7.2673 = 
 0.0000 + 1.1888i 


0.0000 − 1.1888i

−0.0638 + 0.4725i
 
−0.0638 − 0.4725i 
λ|Vr =7.5 =  0.0023 + 1.1785i 


0.0023 − 1.1785i
The obtained eigenvalues from the three sets of results above indicate that the Jacobian
matrix consistently has pairs of complex conjugate eigenvalues. Moreover, as the system
parameter Vr surpasses the flutter critical value from small to large, the Jacobian matrix
exhibits a pair of complex conjugate eigenvalues with a change in the real part from negative
to positive. This characteristic satisfies the definition of Hopf bifurcation. Therefore, the
aforementioned process demonstrates that the nonlinear flutter system experiences Hopf
bifurcation at the flutter critical state.
Furthermore, Hopf bifurcation encompasses three types of critical bifurcations. Next,
the nonlinear flutter response in the intervals on both sides of the Hopf bifurcation point
will be analyzed through calculations to further prove the Hopf bifurcation type. The
Appl. Sci.
Appl. Sci. 2023,
2023, 13,
13, xx FOR
FOR PEER
PEER REVIEW
REVIEW 10 of
10 of 20
20

Appl. Sci. 2023, 13, 10272 Furthermore, Hopf

Furthermore, Hopf bifurcation
bifurcation encompasses
encompasses three three types
types of of critical
critical bifurcations.
bifurcations.10Next, Next,
of 19
the nonlinear flutter response in the intervals on both
the nonlinear flutter response in the intervals on both sides of the Hopf bifurcation point sides of the Hopf bifurcation point
will be
will be analyzed
analyzed throughthrough calculations
calculations to to further
further prove prove the the Hopf
Hopf bifurcation
bifurcation type. type. The The
Runge–Kuttamethods
Runge–Kutta
Runge–Kutta methodswill
methods willbe
will beused
be usedtoto
used tocompute
computethe
compute themotion
the motionresponse
motion response
response of
ofof this
this
this nonlinear
nonlinear
nonlinear flut-
flut-
flutter
ter equation
ter equation
equation [32].
[32].[32].
TheThe The flutter
flutter
flutter equations,
equations,
equations, as established
as established
as established in thisin this
in this
paper, paper,
paper, have
havehave been
beenbeen formulated
formulated
formulated into
into
ainto aa first-order
first-order
first-order ordinary ordinary
ordinary
differential differential
differential
equationequationequation
represented represented
represented
as Equation as (14).
as Equation
Equation (14). Conse-
(14).
Consequently, Conse-for
quently,
quently,
the purpose for of
for thesolving
the purpose
purpose this ofequation,
of solving this
solving this equation,
theequation,
adoption of thethe
the adoption
adoption of the
of
fourth-order theRunge–Kutta
fourth-order method
fourth-order Runge–
Runge–
isKutta
Kutta method is
method
appropriate. is appropriate.
Byappropriate.
carefully selecting By carefully
By carefully selecting an
selecting
an appropriate an appropriate
time appropriate time step,
time
step, the temporal step, the temporal
the
response temporal
of the
response
response
flutter canofof the
bethe flutterthereby
flutter
obtained, can be
can be obtained,
obtained, thereby
facilitating thereby facilitating
facilitating
a more in-depth aa more
analysismore in-depth analysis
of in-depth
nonlinear analysis
vibrations. of
of
nonlinear vibrations.
•nonlinear
Wind vibrations.
speed: U = 75 m/s (beyond the critical state)
•• Wind speed:
Wind speed: 𝑈 𝑈= = 75 75 m/s
m/s (beyond
 .
(beyond the
T
the critical
critical state)
state)
Initial state
Initial state of
state of motion:
of motion: h,
motion: ℎ,h,
ℎ, ℎℎ,,α,
𝛼,α𝛼
𝛼,
. 𝛼 = = (0,
= 0,0,0.001,0
0, 0.001, 0) T (small
0,0,0.001,0 (small perturbation).
(small perturbation).
perturbation).
1
Taking
Takingwind
Taking wind
wind speed
speed
speed 𝑈
U =𝑈75=m/s, = 75
75 m/sm/s , setting
, setting
setting the system the system
the vibration vibration
system frequency frequency
vibration ωfrequency
= 1.888 rad/s, 𝜔𝜔= =
1.888 rad/s
1.888 rad/s ,, and
and usingusing the the initial
initial conditions
conditions for numerical .integration
for numerical . T
integration ℎ, ℎℎ,, 𝛼,
𝛼, 𝛼𝛼 T= =
= (0,ℎ,
 
and using the initial conditions for numerical integration h, h, α, α 0, 0.001, 0) ,
0,0,0.001,0 ,, the
0,0,0.001,0 the vibration
vibration time time history response response is is obtained
obtained as as shown
shown 1 below. Addition-
the vibration time history responsehistory is obtained as shown below. Additionally, below.theAddition-
vibration
ally,
ally, the
the vibration
vibration limit
limit cycle
cycle is
is plotted.
plotted.
limit cycle is plotted.
From Figures
From Figures 4 and and 5, 5, it
it can
can be be observed
observed that that whenwhen the the wind
wind speed
speed 𝑈 𝑈= = 75
75 𝑚/𝑠,
𝑚/𝑠, the the
From Figures 44 and 5, it can be observed that when the wind speed U= 75 m/s, the
vibration response
vibration response
response of of the
of the nonlinear
the nonlinear
nonlinear flutter flutter system
flutter system
system begins begins
begins to to diverge.
to diverge.
diverge. ThisThis indicates
This indicates
indicates that that
that at at
at
vibration
this point,
this point, thethe focus
focus of of the
the flutter
flutter system,
system, i.e., i.e., the
the zero
zero point,
point, becomes
becomes unstable.
unstable. It It becomes
becomes
this point, the focus of the flutter system, i.e., the zero point, becomes unstable. It becomes
an unstable
an unstable focus,
unstable focus,
focus, and and flutter
and flutter occurs
flutter occurs
occurs under under
under the the influence
the influence
influence of of a small
of aa small perturbation
small perturbation
perturbation α=0.001α=0.001
α=0.001
an
rad.
rad. TheThe motion
The motion trajectory
motion trajectory gradually
trajectory gradually
gradually moves moves
moves away away
away from from
from this this unstable
this unstable focus.
unstable focus. After
focus. After a certain
After aa certain
certain
rad.
period
period of of time,
of time,
time, the the flutter
the flutter
flutter tendstends towards
tends towards a
towards aa stable stable
stable state state known
state known
known as as limit cycle
limit cycle oscillation.
cycle oscillation.
period as limit oscillation.
Figure 66 represents
Figure represents the the phase
phase portraits
portraits of of Figures
Figures 44 and and 55 inin the
the velocity
velocity and
and displacement
displacement
Figure 6 represents the phase portraits of Figures 4 and 5 in the velocity and displacement
planes,illustrating
planes,
planes, illustratingthe
illustrating theoccurrence
the occurrenceofof
occurrence of stable
stable
stable limit
limit
limit cycle
cycle
cycle oscillation
oscillation
oscillation in the
in
in the the system
system
system at the
at
at the the cur-
cur-
current
rent wind
rent wind
wind speed. speed.
speed.

44
(h,,hh,α,α
,α,α)=(0,0,0.001,0)
)=(0,0,0.001,0) U=75m/s
U=75m/s
(h
hh
22 hh
h(m/s)
h(m)&&h(m/s)

00
h(m)

-2
-2

-4
-4
00 500
500 1000
1000 1500
1500 2000
2000 2500
2500
tt
Figure 4.
Figure
Figure 4. Continuous
4. Continuous vibration
Continuous vibration of of bridge
bridge flutter
flutter in
in vertical
vertical bending
bending DOF
DOF at
at wind
wind speed
speed of 75m/s
of 75
75 m/s
m/s
and initial state of motion  ℎ, .ℎ, 𝛼, .𝛼 T = 0,0,0.001,0 .
and initial
and initial state
state of
of motion
motion h,ℎ,h,
ℎ,α,
𝛼,α𝛼 = = (0,0,0.001,0
0, 0, 0.001, 0) T. .
1
0.4
0.4
(h,,hh,α,α
,α,α)=(0,0,0.001,0)
)=(0,0,0.001,0) U=75m/s
U=75m/s
0.3 (h
0.3 αα
0.2 αα
0.2
(rad/s)
(rad)&&αα(rad/s)

0.1
0.1
0.0
0.0
αα(rad)

-0.1
-0.1
-0.2
-0.2
-0.3
-0.3
-0.4
-0.4
00 500
500 1000
1000 1500
1500 2000
2000 2500
2500
tt

Figure 5. Continuous vibration of bridge flutter in torsional DOF at wind speed of 75 m/s and initial
 . . T

state of motion h, h, α, α = (0, 0, 0.001, 0)T .
1
Appl. Sci. 2023, 13, x FOR PEER REVIEW 11 of 20

Appl. Sci. 2023, 13, 10272 Figure 5. Continuous vibration of bridge flutter in torsional DOF at wind speed of 75 m/s and11
ini-
of 19
tial state of motion ℎ, ℎ, 𝛼, 𝛼 = 0,0,0.001,0 .

0.2
2

α (rad/s)
h (m/s)

0 0.0

-1

-2
-3 -2 -1 0 1 -0.2
h (m) -0.2 0.0 0.2
α (rad)
(a) Vibration LCO in vertical DOF (b) Vibration LCO in torsional DOF
Figure 6. Phase
Figure 6. Phasetrajectory diagram
trajectory diagram of LCO
of LCOin vertical andand
in vertical torsional DOF
torsional at wind
DOF speed
at wind 7575
of of
speed m/s
m/s
 . T
and initial state of motion ℎ, ℎ , 𝛼,
and initial state of motion h, h, α, α𝛼 . = 0,0,0.001,0 . T
= (0, 0, 0.001, 0) .
1

• • Wind Wind Speed:

Speed: WindWind speed
speed 𝑈U == 7575 m/sm/s (beyond
(beyond the
the critical state),
critical state),
Initial state of motion: ℎ,ℎ, 𝛼, . 𝛼 . = T 0,0,0.2,0 . T
ToInitial state ofthat
demonstrate motion:
the LCO α, α =is(0,
h, h,vibration 0, 0.2,
stable 0) repetitive,
and . it is necessary to select
2
To demonstrate that the LCO vibration is stable
initial conditions with amplitudes greater than the stable amplitude. Therefore, and repetitive, it is necessary when to the
select
wind speed 𝑈 = 75 m/s , the initial condition for the motion integration is set as
initial conditions with amplitudes greater than the stable amplitude. Therefore, when
ℎ,the wind=speed U = 75 m/s, the initial
time condition
response isfor the motion andintegration is set as
ℎ, 𝛼,. 𝛼 .  T
0,0,0.2,0 . The system’s computed shown below.
T
h,From = (0,
h, α, αFigures 0, 0.2,
7 and 8,0it) can
. Thebe system’s
observedtime that response
at the same is computed
wind speed and𝑈 shown
= 75 m/s, below.
by
2
only changing the torsional displacement α in the integration
From Figures 7 and 8, it can be observed that at the same wind speed U = 75 m/s, byinitial conditions to a value
greater than the torsional
only changing the torsional amplitude
displacementof the limit cycle
α in the oscillation,
integration the system’s
initial conditions vibration
to a value
response
greaterexhibits
than thedamping
torsionaland eventually
amplitude settles
of the limitinto a stable
cycle state after
oscillation, a certain vibration
the system’s period
of response
time. Additionally,
exhibits dampingit can be andnoticed
eventually that settles
the amplitude
into a stable of the
statelimit
aftercycle oscillation,
a certain period of
oncetime.it reaches the stable
Additionally, it can state, is the same
be noticed that theas the stable limit
amplitude of thecycle
limitoscillation amplitudes
cycle oscillation, once it
shown
reaches in Figures
the stable 4 and
state,5.isThethestable
same as limit
thecycle
stableoscillation
limit cycleobserved
oscillation in amplitudes
the steady-state shown
in Figures
vibration 4 and 5.after
response, Theastable
periodlimit cycleisoscillation
of time, observed
indeed a stable limitincycle.
the steady-state
The phase plots vibration
of
theresponse,
two degrees afterof a period
freedom of are
time, is indeedina Figure
illustrated stable limit
9, with cycle.
the The
motion phase plots of being
trajectory the two
degrees of
clockwise. Thefreedom
phase are illustrated
trajectory of the in Figure
vertical9,bending
with the degree
motion of trajectory
freedombeing startsclockwise.
at the
The phase
origin, and the trajectory
trajectoryofcurvethe vertical
indicates bending
a gradual degree of freedom
convergence startsaatclosed
towards the origin,
circularand
the The
loop. trajectory
phasecurve indicates
trajectory of the a gradual
torsionalconvergence
degree of freedomtowardsstartsa
closed
.
at 𝛼, 𝛼
circular = 0.2,0
loop. The
;
phase trajectory of the torsional degree of freedom starts
similarly, the trajectory curve gradually approaches a closed circular loop. In comparison at α, α = ( 0.2, 0 ) ; similarly, the
to trajectory
the phase plot curve gradually
shown approaches
in Figure 6, the closeda closed circular
circular looploop.
that In
thecomparison to the phase
limit cycle oscillation
trajectory converges towards is the same, proving that this closed circular loop is atrajectory
plot shown in Figure 6, the closed circular loop that the limit cycle oscillation stable
converges
limit towardsby
cycle exhibited is the
the same,
nonlinear proving that this
vibration closed
system at acircular
wind speed loop isofa 𝑈 stable
= 75limit
m/s.cycle
exhibited by the nonlinear vibration system at a wind speed of U = 75 m/s.
From the aforementioned vibration of the limit cycle oscillation, it is evident that when
wind speed exceeds the critical state, the system exhibits stable limit cycle oscillations,
thereby demonstrating that the nonlinear vibration system undergoes a supercritical Hopf
bifurcation at the critical state.
Appl. Sci.
Appl. Sci. 2023,
2023, 13,
13, 10272
xx FOR
FOR PEER
PEER REVIEW
REVIEW 12
12 of
12of 20
of 20
19

Appl. Sci. 2023, 13, x FOR PEER REVIEW 12 of 20

66
(h,hh,α,α)=(0,0,0.2,0)
(h )=(0,0,0.2,0) U=75m/s
U=75m/s
6 hh
44
(h,h,α,α)=(0,0,0.2,0) U=75m/shh
4 h
22 h

h(m/s)
& h(m/s)
2
00
h(m) & h(m/s)
h(m) &
h(m)
0
-2
-2
-2
-4
-4
-4
-6
-6
00 500
500 1000
1000 1500
1500 2000
2000 2500
2500
-6
0 500 1000 tt 1500 2000 2500
t
Figure 7.
Figure 7. Continuous
Continuous vibration of bridge flutter flutter in vertical
vertical bending
bending DOF
DOF at wind speed of 7575m/s
m/s
 . T
Figure 7. Continuous
initial state
and initial state of vibration
of motion
motion h, of bridge
ℎ, h,
ℎ, α,
.
𝛼, α𝛼 = flutter in vertical
= (0,0,0.2,0 T
0, 0, 0.2, 0). . bending DOF at wind speed of 75 m/s
and initial state of motion ℎ, ℎ, 𝛼, 𝛼 = 20,0,0.2,0 .
0.4
0.4
0.4 (h,hh,α,α)=(0,0,0.2,0)
)=(0,0,0.2,0) U=75m/s
U=75m/s
0.3 (h
0.3 α
(h,h,α,α)=(0,0,0.2,0) U=75m/sα
0.3
α αα
0.2
0.2
0.2 α
(rad/s)
& αα(rad/s)

0.1
0.1
α(rad) & α(rad/s)

0.1
0.0
0.0
(rad) &

0.0
αα(rad)

-0.1
-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3-0.3
-0.4
-0.4-0.4
0 00 500
500 500 1000
1000 1000 1500
1500 1500 2000 2000
2000 2500 2500
2500
t tt

Figure
Figure
Figure 8. Continuous
8. Continuous
8. Continuous vibration
vibration
vibration ofbridge
of bridge
of bridge flutter
flutter
flutter inintorsional
torsional
in torsional DOFDOFDOF
at atatwind
wind windspeed
speed speed
of 75of 7575
of
m/s m/sm/s
and and
ini-and ini-
initial
 . . T

tialtial
state
state ofofmotion
state motion
of motion
h, ℎ,ℎ,
ℎ,h, 𝛼,ℎα𝛼, 𝛼, 𝛼=
α, = 0,0,0,0.2,0
= (0,0,0.2,0
0, T
0.2, 0.) . .
2

4 44
0.2 0.2
0.2
3 33

2 22

1 1
1
α (rad/s)
h (m/s)

(rad/s)
αα (rad/s)
(m/s)
hh (m/s)

0 0.0
00 0.0
0.0
-1
-1
-1
-2
-2
-2
-3
-3
-3
-0.2
-4 -0.2
-0.2
-6-4 -5 -4 -3 -2 -1 0 1 2 -0.2 0.0 0.2
-4
-6 -5
-6 -4 h-3
-5 -4 -3(m)-2
-2 -1
-1 00 11 22 -0.2
-0.2 α (rad) 0.0
0.0 0.2
0.2
hh (m)
(m) (rad)
α (rad)
(a) Vibration LCO in vertical DOF
(b) Vibration LCO in torsional DOF
(a) Vibration LCO in vertical DOF (b) Vibration LCO in torsional DOF
Figure 9. Phase trajectory diagram of LCO in vertical and torsional DOF at wind speed of 75 m/s
Figure 9. Phase trajectory diagram of LCO in vertical and torsional DOF at wind speed of 75 m/s
ℎ, ℎdiagram
, 𝛼, 𝛼 =
. ofT0,0,0.2,0
 .
= (0, 0, 0.2, 0)T . and torsional DOF at wind speed of 75 m/s
and initial9.
Figure state of motion
Phase trajectory .
LCO in vertical

and initial state of motion h, h, α, α
and initial state of motion ℎ, ℎ, 𝛼, 𝛼 2 = 0,0,0.2,0 .
 From the aforementioned vibration of the limit cycle oscillation, it is evident that
when wind speed exceeds the critical state, the system exhibits stable limit cycle oscilla-
tions, thereby demonstrating that the nonlinear vibration system undergoes a supercriti-
Appl. Sci. 2023, 13, 10272 13 of 19
cal Hopf bifurcation at the critical state.

4. Analysis of Nonlinear Bifurcation and Limit Cycle Oscillation

4. Analysis of Nonlinear Bifurcation and Limit Cycle Oscillation
Using the Runge–Kutta method discussed in the previous section, the LCO ampli-
tude Using the Runge–Kutta
at a specific wind speedmethod discussedthrough
can be obtained in the previous section,
integration. the LCObyamplitude
Therefore, selecting
at a specific wind speed can be obtained through integration. Therefore, by
wind speed values within a certain range and substituting them into the reference case selecting windof
speed values within a certain range and substituting them into the reference
the nonlinear motion system, we can integrate them to obtain a set of correspondences case of the
nonlinearwind
between motion system,
speed andwethecan integrate
stable LCOthem to obtain
amplitude. a setdata
This of correspondences
can be used to between
plot the
wind speed and the stable LCO amplitude. This data can be used
bifurcation curve of the flutter, as shown in Figure 10, which illustrates to plotthe
therelationship
bifurcation
curve of the
between flutter,
wind speed asand
shown
the in Figure
stable LCO10,amplitude
which illustrates the relationship
of the torsional degree ofbetween
freedom wind
for
speed and the stable LCO amplitude of the torsional degree of freedom for this particular
this particular operating condition. In this context, the abscissa U signifies the actual wind
operating condition. In this context, the abscissa U signifies the actual wind speed, while
speed, while the ordinate A denotes the maximum amplitude of the limit cycle for the
the ordinate A denotes the maximum amplitude of the limit cycle for the torsional degree
torsional degree of freedom.
of freedom.
0.25

0.20

0.15
A (rad)

0.10

0.05

U=57.3374m/s
0.00
55 60 65 70 75 80
U (m/s)

Figure 10. Bifurcation curves for the relationship between stable LCO amplitudes and wind speed.
speed.

flutter critical
At the flutter critical state 𝑈 cr =
state U = 57.3374 m/s, the nonlinear flutter system undergoes
57.3374 m/s,
a Hopf bifurcation, and a stable LCO bifurcation curve
a Hopf bifurcation, and a stable LCO bifurcation appears on
curve appears on the
the right
right side
side of
of the
the
bifurcation point. The LCO amplitude rises with the increasing wind speed. From the
bifurcation point. The LCO amplitude rises with the increasing wind speed. From the
graph, itit can
graph, canbebeobserved
observedthat thatafter
afterthe
thewind
wind speed
speed exceeds
exceeds 𝑈U = 77.7
= 77.7 m/s,m/s,
the the trend
trend of
of the
the limit cycle amplitude shows a steep increase, and stable limit cycle vibration
limit cycle amplitude shows a steep increase, and stable limit cycle vibration cannot be cannot be
obtained at
obtained at higher
higher wind
wind speeds. Therefore, for
speeds. Therefore, for this
this operating
operating condition of the
condition of the nonlinear
nonlinear
flutter system, vibration divergence occurs when the wind speed exceeds 𝑈 = 77.7m/s.
flutter system, vibration divergence occurs when the wind speed exceeds U = 77.7 m/s.
Further research
Further research can
can bebe conducted
conducted on on the
the bifurcation curve in
bifurcation curve in Figure
Figure 10,
10, which
which may
may
reveal the presence of unstable LCO within the range of wind speeds U
reveal the presence of unstable LCO within the range of wind speeds 𝑈 < 77.7 m/s, with < 77.7 m/s, with
amplitudes larger
amplitudes larger than
than the
the stable
stable limit
limit cycle. This can
cycle. This can be
be investigated by setting
investigated by setting larger
larger
initial conditions and performing integration calculations.
initial conditions and performing integration calculations.
• Wind Speed: Wind speed U 𝑈== 75
75 m/s
𝑚/𝑠 (beyond
(beyondthe thecritical
criticalstate);
state);
Initial state of motion: ℎ, ℎ. , 𝛼, 𝛼.  T = 0,0,0.45,0 .T
Initial state of motion: h, h, α, α = (0, 0, 0.45, 0) .
By selecting a wind speed of 𝑈3= 75 m/s and initial conditions of ℎ, ℎ, 𝛼, 𝛼 =
By selecting a wind speed of U = 75 m/s and initial conditions of
0,0,0.45,0
. . T , the time response
 of the flutter system vibration can be computed.
T
h, h,As
α, αshown= (in
0, 0, 0.45, 0)11
Figures , and
the time response
12, under theofgiven
the flutter system
initial vibration
conditions, canvibration
the be computed.
am-
3
As shown in Figures 11 and 12, under the given initial conditions,
plitude of the system initially decays and eventually stabilizes into a LCO vibration re- the vibration
amplitude
sponse afterofathe system
period initially
of time. decays this
Moreover, andstable
eventually stabilizes response
LCO vibration into a LCO vibration
is consistent
response after a period of time. Moreover, this stable LCO vibration response is consistent
with the stable limit cycle response shown in Figures 4 and 5. The phase portraits of these
two figures, representing the limit cycles, are plotted in Figure 13.
• Wind Speed: Wind speed U = 75 m/s (beyond the critical state);
 . . T

Initial state of motion: h, h, α, α = (0, 0, 0.452, 0) T .
4
Appl. Sci. 2023, 13, 10272 14 of 19
Appl. Sci. 2023, 13, x FOR PEER REVIEW 14 of 20
Appl.
Appl. Sci.
Sci. 2023,
2023, 13,13, x FOR
x FOR PEER
PEER REVIEW
REVIEW 1414ofof2020

The selected initial conditions of motion are slightly greater than the previous inte-
with theinitial
gration stable values.
limit cycle
The response shown
resulting in Figures 4response
time-domain and 5. The phase
ofphase portraits
the system of these
is of
calculated and
withthe
with thestable
stablelimit
limitcycle
cycle response
response shown
shown inFigures
inareFigures4 4and
and5.5.The
The portraits
phase portraits these
of these
two
shownfigures,
in representing
Figures 14 and the limit
15. cycles, plotted in Figure 13.
two figures, representing the limit cycles, are plotted in Figure 13.
two figures, representing the limit cycles, are plotted in Figure 13.
15
1515 U=75m/s
10 (h,h,α,α)=(0,0,0.45,0) U=75m/s
1010 (h,h, α,αα,)=(0,0,0.45,0)
(h,h, α)=(0,0,0.45,0) U=75m/s
h
5 h hh
55 hh
h (m/s)

0
h (m/s)

00
h (m) & h (m/s)

-5
&&

-5-5
h (m)

-10
h (m)

-10-10
-15
-15-15
-20
-20-20
-25
-25-25 00 500
500
1000
1000
1500
1500
2000
2000
2500
2500
0 500 1000 t 1500 2000 2500
t t
Figure 11.
Figure 11. Continuous
Continuous vibration
vibration of of
bridge flutter
bridge in vertical
flutter bending
in vertical DOF atDOF
bending windatspeed 75 m/s
ofspeed
windof of 75 m/s
Figure11.
Figure 11.Continuous
and initial
Continuous vibration
vibration
state of motion  ℎof
ℎ,
of bridge
. bridge
, 𝛼,
flutter
𝛼 . =Tflutter ininvertical
0,0,0.45,0
vertical
.
bending
bending DOFatatwind
DOF windspeed of 7575
speed m/s
m/s
T
andinitial
and
and initial
initial state
ofofmotion
state
state ofmotion
motionℎ,ℎ,
ℎ,ℎh,
𝛼,, 𝛼,
h, α=0,0,0.45,0
𝛼 𝛼α, = 0,0,0.45,0 . . 0) .
= (0, 0, 0.45,
3
0.5
0.50.5 U=75m/s
0.4 (h,h,α,α)=(0,0,0.45,0) U=75m/s
0.40.4 α,α
(h,h,
(h,h, α)=(0,0,0.45,0)
α,)=(0,0,0.45,0) α
U=75m/s
0.3 α αα
0.30.3 αα
0.2
α (rad/s)

0.20.2
α (rad/s)

0.1
α (rad) & α (rad/s)

0.10.1
0.0
&&

0.00.0
α (rad)

-0.1
α (rad)

-0.1
-0.1
-0.2
-0.2
-0.2
-0.3
-0.3
-0.3
-0.4
-0.4
-0.4
-0.5
-0.5 0
-0.5 500 1000 1500 2000 2500
00 500
500 1000
1000 t 1500
1500 2000
2000 2500
2500
t t
Figure 12. Continuous vibration of bridge flutter in torsional DOF at wind speed of 75 m/s and
Figure12. 12. Continuous
Continuous vibration
vibration ofbridge
bridge flutter
of bridge in torsional
flutter DOF at wind
atspeed 75 m/s
windofspeed of and
. torsional DOF at wind speed of 75 m/s and
Figure
Figure 12.
Continuous in torsional DOF 75 m/s and
initial state of motion vibration
ℎ, ℎ, .𝛼, 𝛼 of= flutter
0,0,0.45,0 in
initial state of motion ℎ, ℎ , 𝛼, 𝛼 . =T 0,0,0.45,0 .

initial state
initial of of
state motion
motionℎ, ℎ,h,
𝛼,h, = 0,0,0.45,0
𝛼 α, α . 0) .
= (0, 0, 0.45, T
3
15
1515
0.4
10 0.40.4
1010

5 0.2
55 0.20.2
(rad/s)
h (m/s)

α α(rad/s)
h (m/s)

0
α (rad/s)
h (m/s)

0.0
00 0.00.0

-5
-5-5 -0.2
-0.2
-0.2
-10
-10
-10
-0.4
-0.4
-0.4
-15
-15
-15 -20 -15 -10 -5 0 -0.4 -0.2 0.0 0.2 0.4
-20 -15
-20 -15 -10
h -10
(m) -5-5 00 -0.4 -0.2
-0.2 0.0
h (m) -0.4 α0.0
(rad)0.20.2 0.40.4
h (m) α (rad)
α (rad)
(a) Vibration LCO in vertical DOF (b) Vibration LCO in torsional DOF
(a)Vibration
(a) VibrationLCO
LCOininvertical
verticalDOF
DOF (b)
(b)Vibration
VibrationLCO
LCOinintorsional
torsionalDOF
DOF

Figure 13. Phase trajectory diagram of LCO in vertical and torsional DOF at wind speed of 75 m/s
 . . T

and initial state of motion h, h, α, α = (0, 0, 0.45, 0)T .
3
 • Wind Speed: Wind speed 𝑈 = 75 m/s (beyond the critical state);
Initial state of motion: ℎ, ℎ, 𝛼, 𝛼 = 0,0,0.452,0 .
The selected initial conditions of motion are slightly greater than the previous inte-
Appl. Sci. 2023, 13, 10272 15 of 19
gration initial values. The resulting time-domain response of the system is calculated and
shown in Figures 14 and 15.

30
U=75m/s
(h,h,α,α)=(0,0,0.452,0) h
20
h
10

h (m) & h (m/s)

0

-10

-20

-30

0 100 200 300 400 500

t

Figure 14. Continuous

Figure 14. Continuous vibration
vibration of
of bridge
bridge flutter in vertical
flutter in vertical bending
bending DOF
DOF at
at wind
wind speed
speed of 75m/s
of 75 m/s
 . . T

and initial
and initial state
state of
of motion
motion h,ℎ,h,
ℎ,α,
𝛼,α𝛼 = = (0,0,0.452,0 T
0, 0, 0.452, 0) . .
4

1.0
(h,h,α,α)=(0,0,0.452,0) U=75m/s
0.8
α
0.6 α
0.4
α (rad) & α (rad/s)

0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
0 100 200 300 400 500
t

Figure 15.
Figure Continuous vibration of bridge flutter
15. Continuous flutter in torsional DOF at wind
wind speed
speed of 75m/s
of 75 m/s and
and
 . .𝛼 T = 0,0,0.452,0 T.
ℎ, ℎ , 𝛼,

initial state of motion
initial state of motion h, h, α, α = (0, 0, 0.452, 0) .
4

The vibration
The vibration response
response ofof the
the system
system shows
shows aa slow
slow increase
increase in in amplitude
amplitude andand max
max
speed in
speed in both
both DOFs,
DOFs, indicating
indicating aa gradual
gradual divergence
divergence of of the
theoscillations.
oscillations. The
The vibration
vibration
phase diagram
phase diagram isis shown
shown in in Figure
Figure 16.
16. AA comparison
comparison with with the
the motion
motion phase
phase diagram
diagram inin
 . T
Figure 13,
Figure 13, which
whichhad hadthe
theinitial conditions h,ℎ,h,ℎα,
initialconditions , 𝛼,α 𝛼 == (0,
. 0,0,0.45,0 under the
0, 0.45, 0) T under the same
same
3
wind speed,
wind speed, reveals
reveals the
the presence
presence of of an
an unstable
unstable limit
limit cycle
cycle oscillation
oscillation within
within the
the range
range of
of
torsional amplitudes between 0.45 and 0.452. The vibration phase trajectories
torsional amplitudes between 0.45 and 0.452. The vibration phase trajectories on either side on either
side
of ofLCO
this this LCO will gradually
will gradually movemoveawayaway
from from
it. it.
From the analysis above, at a wind speed of U = 75 m/s, the origin point of this
nonlinear flutter system is an unstable focus. Additionally, there are two LCOs in the
phase diagram. The smaller-amplitude LCO is stable, and vibration starting from the
unstable focus will eventually converge to this stable LCO. The larger amplitude limit
cycle is unstable one. Vibrations within the range between the unstable and stable LCOs
will gradually approach the stable one, while vibrations with amplitudes greater than the
unstable LCO will diverge.
In summary, the Runge–Kutta integration method can be used to solve flutter systems
with nonlinear aerodynamic self-excitation and obtain bifurcating LCOs beyond the critical
state. The stable LCO amplitude at a specific wind speed can be obtained through the
steady-state time response, while the unstable LCO requires the following calculation steps:
select a wind speed U and an initial condition (0, 0, α0 , 0) T , calculate the time response of
the nonlinear system vibrations, and continuously change the initial amplitude α0 . When
 Appl. Sci. 2023, 13, 10272 16 of 19

the response above this value converges to the stable limit cycle and the response above
this value diverges, it can be concluded that there exists an unstable limit cycle at the wind
speed U, with an amplitude of α0 . By substituting the parameter values of the reference
case into the flutter Equations (4) and (5), the stable and unstable limit cycles corresponding
Appl. Sci. 2023, 13, x FOR PEER REVIEW 16 of 20
to the wind speed can be obtained. Figure 17 shows the relationship between the limit cycle
amplitudes and wind speed.
20 0.6

15
0.4
10

0.2
5

α (rad/s)
h (m/s)

0 0.0

-5
-0.2
-10

-15 -0.4

-20
-30 -20 -10 0 -0.6
h (m) -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
α (rad)
(a) Vibration LCO in vertical DOF (b) Vibration LCO in torsional DOF
Figure
Figure16.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 16.Phase
Phasetrajectory
trajectorydiagram
diagramofofLCO
LCOininvertical
verticaland
andtorsional DOF
torsional atat
DOF wind speed
wind 75
ofof17
speed m/s
75of 20
m/s
 ℎ,.𝛼, 𝛼 .  T= 0,0,0.452,0 .
and initial state of motion ℎ,
and initial state of motion h, h, α, α = (0, 0, 0.452, 0)T .
4

From the analysis above, at a wind speed of 𝑈 = 75 m/s, the origin point of this non-
linear
0.8
flutter system is an unstable focus. Additionally, there are two LCOs in the phase
diagram. The smaller-amplitude LCO is stable, and vibration starting from the unstable
focus will eventually converge to this stable LCO. The larger amplitude limit cycle is un-
0.6 one. Vibrations within the range between the unstable and stable LCOs will gradu-
stable
ally approach the stable one, while vibrations with amplitudes greater than the unstable
A (rad)

LCO0.4will diverge.
In summary, the Runge–Kutta integration method can be used to solve flutter sys-
Bifurcation of stable LCO
tems with nonlinear aerodynamic
Bifurcation of unstable self-excitation
LCO and obtain bifurcating LCOs beyond the
0.2 state. The stable LCO amplitude at a specific wind speed can be obtained through
critical
the steady-state time response, while the unstable LCO requires the following calculation
steps:
0.0
select a wind speed 𝑈 and an initial condition 0,0, 𝛼 , 0 , calculate the time re-
sponse 55 of the nonlinear
60 65
system 70
vibrations, and75 continuously
80 change the initial amplitude
U (m/s)
𝛼 . When the response above this value converges to the stable limit cycle and the re-
Figure
Figure17.
sponse 17.Bifurcation
above curves
this value
Bifurcation for
forthe
therelationship
diverges,
curves between
betweenLCO
it can be concluded
relationship amplitudes
that
LCO and
there exists
amplitudes wind
andan speed.
unstable
wind speed. limit
cycle at the wind speed 𝑈, with an amplitude of 𝛼 . By substituting the parameter values
Taking
of theTaking
reference intocase
into account
account the parameter
into the settings for
settings(4)
flutter Equations the
forand other
the(5),
otherthree
thethreeoperating
stableoperating conditions
limitin
conditions
and unstable
Table
in Table 2, the
2, relationship
the between
relationship between the LCO
the amplitude
LCO amplitude and wind
and speed
wind after
speed
cycles corresponding to the wind speed can be obtained. Figure 17 shows the relationship the
aftersupercritical
the super-
Hopf bifurcation
critical Hopf is calculated
bifurcation using the
is calculated aforementioned method. The results are plotted
between the limit cycle amplitudes andusing
windthe aforementioned
speed. method. The results are
in Figures 18 and 19.
plotted in Figures 18 and 19.
As shown in the results of the LCO bifurcation curves, the flutter systems of all four
0.20
operating conditions undergo supercritical Hopf bifurcation at the same critical value,
leading to stable limit cycle oscillations on the side above the critical value. The stable LCO
amplitude rises with the increasing wind speed. Comparing the curves of the reference
0.15and case A, the two curves nearly overlap, but the stable LCO amplitude of case
case
A (rad)

0.10
 Appl.Sci.
Appl. Sci.2023,
2023,13,
13,xxFOR
FORPEER
PEERREVIEW
REVIEW 17 ofof 20
17 20

Appl. Sci. 2023, 13, 10272 0.8

0.8 17 of 19

0.6
0.6
A is relatively smaller than that of the reference case, and the unstable LCO amplitude

(rad)
AA(rad)
is larger than that of the reference case. This indicates that applying stronger nonlinear
0.4
damping
0.4 on the torsional degree of freedom has a limited dissipating effect on limit
cycle oscillations Bifurcation
underofofsupercritical
Bifurcation stableLCO
stable LCO flutter. Furthermore, comparing the curves of case
A and case B, Bifurcation
it can be of unstable
Bifurcation of unstable LCO
observed LCOto achieve the same LCO amplitude, case B needs
that
0.2
0.2
higher wind speeds. This suggests that considering the nonlinear stiffness of the bridge
structure reduces the stable LCO amplitude and increases the unstable LCO amplitude
under
0.0 supercritical flutter, indicating a higher stability of the bridge compared to a system
0.0
55
55 60
60 6565 70
70 75
75 80
80
without considering structuralUU(m/s)
nonlinear
(m/s) stiffness. Finally, in case C, by adding strong
nonlinear damping on the vertical bending degree of freedom, the corresponding wind
Figure 17.Bifurcation
Figure Bifurcationcurves
curvesfor forthe
therelationship
relationshipbetween
betweenLCO
LCOamplitudes
amplitudesandandwind
windspeed.
speed.
speed17. for the same LCO amplitude is the highest, demonstrating strong energy dissipation
and suppression capabilities for supercritical flutter. In conclusion, by considering various
Takinginto
Taking intoaccount
accountthe theparameter
parametersettings
settingsfor
forthe
theother
other threeoperating
operating conditions
nonlinear factors, the variation relationship between the limitthree
cycle amplitudeconditions
and wind
in
in Table
Table 2,
2, the
the relationship
relationship between
between the
the LCO
LCO amplitude
amplitude and
and wind
wind speed
speed afterthe
after thesuper-
super-
speed for the four operating conditions has been analyzed. Among them, considering
critical
critical Hopf
Hopf bifurcation
bifurcation is calculated
is calculated using
using the aforementioned
the aforementioned method. The results are
structural nonlinear stiffness and cubic damping on the verticalmethod.
bendingThe results
degree hasare
the
plotted
plotted in Figures 18 and
in Figuresinfluence
18 and 19. 19.
most significant on the LCO amplitude under supercritical flutter.
0.20
0.20

0.15
0.15
(rad)
AA(rad)

0.10
0.10

Referencecase
Reference case
0.05
0.05 CaseAA
Case
CaseBB
Case
CaseCC
Case
0.00
0.00
55
55 60
60 65
65 70
70 75
75 80
80 85
85 90
90 95
95 100
100
UU(m/s)
(m/s)

Figure18.
Figure
Figure 18. The
18.The Hopf
TheHopf
Hopf bifurcation
bifurcation
bifurcation resultsofofof
results
results the
the calculation
the examples
calculation
calculation ininTable
examples
examples Table2,2,the
in Table the relationship
2, relationship
the be-
relationship
be-
tween
tween
betweenthe
thethe amplitude
amplitude of the
of the
amplitude LCO
LCO
of the and
andand
LCO the
thethewind
wind speed.
speed.
wind speed.

0.8
0.8 Referencecase
Reference case
CaseAA
Case
CaseBB
Case
0.6
0.6 CaseCC
Case
(rad)
AA(rad)

0.4
0.4

0.2
0.2

0.0
0.0
5555 60
60 65
65 70
70 75
75 80
80 85
85 90
90 95
95 100
100
UU(m/s)
(m/s)

Figure19.
Figure
Figure 19. The
19.The Hopf
TheHopf
Hopf bifurcation
bifurcation results
results
bifurcation ofofof
results the
the calculation
calculation
the examples
examples
calculation ininTable
examples Table 2,2,the
in Table the relationship
2, relationship
the be-
be-
relationship
tween
tween the
the amplitude
amplitude of
of the
the LCO
LCO and
and the
the wind
wind speed,
speed, including
including stable
stable and
and unstable
unstable
between the amplitude of the LCO and the wind speed, including stable and unstable LCO. LCO.
LCO.

5. Conclusions
This paper considers the effects of nonlinear damping and structural stiffness on the
limit cycle response of a long-span suspension bridge’s flutter systems under supercritical
Appl. Sci. 2023, 13, 10272 18 of 19

conditions, focusing on analyzing the influence of strong nonlinear damping and cubic
stiffness on limit cycle oscillation at the flutter critical state. The nonlinear flutter self-
excited force model proposed by Gao is considered, and the flutter derivatives obtained
from wind tunnel tests are fitted using the least squares method to establish the nonlinear
flutter self-excited force equation.
The flutter equation is transformed into state-space equations, and the flutter critical
state is solved by analyzing the changes in matrix eigenvalues and obtaining the critical
wind speed value. By using the Jacobian matrix of the state-space equations, the changes in
matrix eigenvalues at the critical point are analyzed, and the supercritical Hopf bifurcation
of the flutter system is demonstrated through vibration responses. Furthermore, the results
of the four cases in Table 2 show that the flutter critical values for the systems are the same,
indicating that supercritical Hopf bifurcation occurs at the same critical value. Finally,
through a comparative analysis of the four case studies, the relationship between the
limit cycle amplitude and wind speed under supercritical conditions is analyzed, demon-
strating the effective flutter suppression and energy dissipation effects of cubic damping
and stiffness.
The above studies indicate that in terms of engineering practicality, the flutter occur-
ring in long-span suspension bridges is predictable and regular, displaying LCOs with
amplitudes increasing as wind speed rises. To mitigate the detrimental effects of flutter on
the bridge, increasing vertical bending and torsional damping of the bridge is an effective
approach. By employing tuned mass dampers to apply strong nonlinear damping in both
vertical bending and torsion motion dimensions, the amplitude of flutter oscillations in
the bridge can be significantly reduced, thereby minimizing further damage caused by
the displacement and deformation of the bridge. With the advancement of deep learning
and machine learning, our research will leverage artificial neural network algorithms to
conduct investigations, facilitating the prediction and evaluation of the vibrational response
of bridges under aerodynamic forces.

Author Contributions: Conceptualization, J.L. and F.W.; methodology, J.L. and F.W.; software, J.L.;
validation, J.L. and Y.Y.; writing—original draft preparation, J.L.; writing—review and editing, J.L.,
F.W. and Y.Y.; visualization, J.L.; supervision, F.W.; project administration, F.W.; funding acquisition,
F.W. All authors have read and agreed to the published version of the manuscript.
Funding: This research was supported by the Key Laboratory of Disaster Forecast and Control in
Engineering (Jinan University), MOE of China.
Informed Consent Statement: Not applicable.
Data Availability Statement: The data presented in this study are available on request from the
corresponding author.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Amandolese, X.; Michelin, S.; Choquel, M. Low Speed Flutter and Limit Cycle Oscillations of a Two-Degree-of-Freedom Flat Plate
in a Wind Tunnel. J. Fluids Struct. 2013, 43, 244–255. [CrossRef]
2. Gao, G.; Zhu, L.; Han, W.; Li, J. Nonlinear Post-Flutter Behavior and Self-Excited Force Model of a Twin-Side-Girder Bridge Deck.
J. Wind Eng. Ind. Aerodyn. 2018, 177, 227–241. [CrossRef]
3. Giaccu, G.F.; Caracoglia, L. A Gyroscopic Stabilizer to Improve Flutter Performance of Long-Span Cable-Supported Bridges. Eng.
Struct. 2021, 240, 112373. [CrossRef]
4. Bera, K.K.; Banerjee, A. A Consistent Dynamic Stiffness Matrix for Flutter Analysis of Bridge Decks. Comput. Struct. 2023,
286, 107107. [CrossRef]
5. Wu, L.; Woody Ju, J.; Zhang, J.; Zhang, M.; Li, Y. Vibration Phase Difference Analysis of Long-Span Suspension Bridge during
Flutter. Eng. Struct. 2023, 276, 115351. [CrossRef]
6. Náprstek, J.; Pospíšil, S.; Hračov, S. Analytical and Experimental Modelling of Non-Linear Aeroelastic Effects on Prismatic Bodies.
J. Wind Eng. Ind. Aerodyn. 2007, 95, 1315–1328. [CrossRef]
7. Král, R.; Pospíšil, S.; Náprstek, J. Wind Tunnel Experiments on Unstable Self-Excited Vibration of Sectional Girders. J. Fluids
Struct. 2014, 44, 235–250. [CrossRef]
Appl. Sci. 2023, 13, 10272 19 of 19

8. Pigolotti, L.; Mannini, C.; Bartoli, G. Experimental Study on the Flutter-Induced Motion of Two-Degree-of-Freedom Plates. J.
Fluids Struct. 2017, 75, 77–98. [CrossRef]
9. Tang, Y.; Hua, X.G.; Chen, Z.Q.; Zhou, Y. Experimental Investigation of Flutter Characteristics of Shallow Π Section at Post-Critical
Regime. J. Fluids Struct. 2019, 88, 275–291. [CrossRef]
10. Wu, B.; Chen, X.; Wang, Q.; Liao, H.; Dong, J. Characterization of Vibration Amplitude of Nonlinear Bridge Flutter from Section
Model Test to Full Bridge Estimation. J. Wind Eng. Ind. Aerodyn. 2020, 197, 104048. [CrossRef]
11. Li, Z.; Wu, B.; Liao, H.; Li, M.; Wang, Q.; Shen, H. Influence of the Initial Amplitude on the Flutter Performance of a 2D Section
and 3D Full Bridge with a Streamlined Box Girder. J. Wind Eng. Ind. Aerodyn. 2022, 222, 104916. [CrossRef]
12. Liao, H.; Mei, H.; Hu, G.; Wu, B.; Wang, Q. Machine Learning Strategy for Predicting Flutter Performance of Streamlined Box
Girders. J. Wind Eng. Ind. Aerodyn. 2021, 209, 104493. [CrossRef]
13. Mostafa, K.; Zisis, I.; Moustafa, M.A. Machine Learning Techniques in Structural Wind Engineering: A State-of-the-Art Review.
Appl. Sci. 2022, 12, 5232. [CrossRef]
14. Tinmitondé, S.; He, X.; Yan, L.; Hounye, A.H. Data-Driven Prediction of Critical Flutter Velocity of Long-Span Suspension Bridges
Using a Probabilistic Machine Learning Approach. Comput. Struct. 2023, 280, 107002. [CrossRef]
15. Scanlan, R.H.; Tomko, J.J. Airfoil and Bridge Deck Flutter Derivatives. J. Eng. Mech. Div. 1971, 97, 1717–1737. [CrossRef]
16. Náprstek, J.; Pospíšil, S. Response Types and General Stability Conditions of Linear Aero-Elastic System with Two Degrees-of-
Freedom. J. Wind Eng. Ind. Aerodyn. 2012, 111, 1–13. [CrossRef]
17. Náprstek, J.; Pospíšil, S.; Yau, J.-D. Stability of Two-Degrees-of-Freedom Aero-Elastic Models with Frequency and Time Variable
Parametric Self-Induced Forces. J. Fluids Struct. 2015, 57, 91–107. [CrossRef]
18. Gao, G.; Zhu, L.; Li, J.; Han, W.; Yao, B. A Novel Two-Degree-of-Freedom Model of Nonlinear Self-Excited Force for Coupled
Flutter Instability of Bridge Decks. J. Sound Vib. 2020, 480, 115406. [CrossRef]
19. Gao, G.; Zhu, L. Nonlinearity of Mechanical Damping and Stiffness of a Spring-Suspended Sectional Model System for Wind
Tunnel Tests. J. Sound Vib. 2015, 355, 369–391. [CrossRef]
20. Zhang, M.; Xu, F.; Ying, X. Experimental Investigations on the Nonlinear Torsional Flutter of a Bridge Deck. J. Bridge Eng. 2017,
22, 04017048. [CrossRef]
21. Tesfaye, S.; Kavrakov, I.; Morgenthal, G. Numerical Investigation of the Nonlinear Interaction between the Sinusoidal Motion-
Induced and Gust-Induced Forces Acting on Bridge Decks. J. Fluids Struct. 2022, 113, 103680. [CrossRef]
22. Zhang, M.; Xu, F.; Zhang, Z.; Ying, X. Energy Budget Analysis and Engineering Modeling of Post-Flutter Limit Cycle Oscillation
of a Bridge Deck. J. Wind Eng. Ind. Aerodyn. 2019, 188, 410–420. [CrossRef]
23. Zhang, M.; Xu, F.; Han, Y. Assessment of Wind-Induced Nonlinear Post-Critical Performance of Bridge Decks. J. Wind Eng. Ind.
Aerodyn. 2020, 203, 104251. [CrossRef]
24. Li, K.; Han, Y.; Cai, C.S.; Hu, P.; Li, C. Experimental Investigation on Post-Flutter Characteristics of a Typical Steel-Truss
Suspension Bridge Deck. J. Wind Eng. Ind. Aerodyn. 2021, 216, 104724. [CrossRef]
25. Wu, B.; Shen, H.; Liao, H.; Wang, Q.; Zhang, Y.; Li, Z. Insight into the Intrinsic Time-Varying Aerodynamic Properties of a
Truss Girder Undergoing a Flutter with Subcritical Hopf Bifurcation. Commun. Nonlinear Sci. Numer. Simul. 2022, 112, 106472.
[CrossRef]
26. Casalotti, A.; Arena, A.; Lacarbonara, W. Mitigation of Post-Flutter Oscillations in Suspension Bridges by Hysteretic Tuned Mass
Dampers. Eng. Struct. 2014, 69, 62–71. [CrossRef]
27. Zhang, M.; Xu, F. Tuned Mass Damper for Self-Excited Vibration Control: Optimization Involving Nonlinear Aeroelastic Effect. J.
Wind Eng. Ind. Aerodyn. 2022, 220, 104836. [CrossRef]
28. Scanlan, R.H. The Action of Flexible Bridges under Wind, I: Flutter Theory. J. Sound Vib. 1978, 60, 187–199. [CrossRef]
29. Hong, G. Identifying Long-Span Bridge Flutter Derivatives via the Free Vibration Method Based on the Fluent Software. Master’s
Thesis, Chang’an University, Xi’an, China, 2012.
30. Andronov, A.A.; Leontovich, E.A.; Gordon, I.I.; Maier, A.G. Qualitative Theory of Second-Order Dynamic Systems; Halsted Press:
New York, NY, USA, 1973; ISBN 0-470-03195-6.
31. Theodorsen, T. General Theory of Aerodynamic Instability and the Mechanism of Flutter; NACA Technical Report; Langley: Newport
News, VA, USA, 1935.
32. Runge, C. Ueber die numerische Auflösung von Differentialgleichungen. Math. Ann. 1895, 46, 167–178. [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.