Nonlinear Aeroelastic in-Plane Behavior of Suspension Bridges under Steady Wind Flow

: The nonlinear aeroelastic behavior of suspension bridges, undergoing dynamical in-plain instability (galloping), is analyzed. A nonlinear continuous model of bridge is formulated, made of a visco-elastic beam and a parabolic cable, connected each other by axially rigid suspenders, continuously distributed. The structure is loaded by a uniform wind ﬂow which acts normally to the bridge plane. Both external and internal damping are accounted for the structure, according to the Kelvin-Voigt rheological model. The nonlinear aeroelastic effects are evaluated via the quasi-static theory, while structural nonlinearities are not taken into account. First, the free dynamics of the undamped bridge are addressed, and the natural modes determined. Then, the nonlinear equations ruling the dynamics of the aeroelastic system, close to the bifurcation point, are tackled by the Multiple Scale Method. This is directly applied to the partial differential equations, and provides the ﬁnite-dimensional bifurcation equations. From these latter, the limit-cycle amplitude and its stability are evaluated as function of the mean wind velocity. A case study of suspension bridge is analyzed.


Introduction
Suspension bridges are long, slender and flexible structures which are very sensitive to dynamic actions induced by wind, which potentially cause a variety of instability phenomena. The importance of wind-resistant design for these structures has been highly recognized and has led to many research works and investigations on bridge aerodynamics. Some historical views of long-span bridge aerodynamics are presented in [1,2].
The first studies were addressed on the dynamic of suspended cables, pertinent to suspension bridges but without the inclusion of the stiffening girder, see Irvine and Caughey [3] and Irvine [4]. Since the 1980s, starting with Hagerdon and Schafer [5] and Luongo et al. [6,7], there has been a large number of studies on the free and forced nonlinear vibrations of suspended cables (see the review of Rega, in [8,9]). In particular, some papers investigated wind-induced aeroelastic instabilities, as [10][11][12].
After the Tacoma Narrows Bridge failure in 1940, major efforts were made to understand the dynamic of suspension bridge, including the effects of the stiffening girder. Bleich [13] modeled the stiffening girder as a flexural beam attached to the elastic cables by inextensible hangers. Additional contributions to the continuum model were made by Pugsley [14]. Abdel-Ghaffar [15] extended the approach to include coupled vertical-torsional vibrations and the effects of cross-sectional distortion. links uniformly distributed. Both external and internal damping are accounted, both for the beam and cable, according to the Kelvin-Voigt rheological model. A steady and uniform wind flow is applied exclusively to the beam. The aeroelastic effects of the wind are evaluated via the classical quasi-static theory (a discussion about strengths and weaknesses of the method is reported in [42]). First, the linear free dynamic of the undamped model is studied; then, the nonlinear bifurcation problem is tackled by the Multiple Scale Method (MSM) which naturally reduces the infinite-dimensional system to the proper finite-dimensional Center Manifold. From the bifurcation equation, the limit cycles and their stability are evaluated as function of the mean wind velocity. The asymptotic results are validated via numerical time-integration of ordinary differential equations derived by in-space finite differences. A case study of suspension bridge is analyzed.
The paper is organized as follows. In Section 2, a continuous visco-elastic model of standard single-span suspension bridge under uniformly distributed steady wind flow, is formulated. In Section 3, the free dynamic of the undamped problem is addressed. In Section 4 the Multiple Scale Method is directly applied to the partial differential equations, to get the bifurcation equation governing the motion close to a Hopf bifurcation. In Section 5, a finite-dimensional model is formulated. In Section 6, the numerical outcomes are presented. Finally, in Section 7 the main findings of the work are summarized. An Appendix A, containing details about damping parameter identification, closes the paper.

Continuous Model
A standard single-span suspension bridge is considered, consisting of a main cable, a stiffening girder, equispaced hangers (or suspenders) and two supporting towers ( Figure 1). The main cable, modeled as uniform, elastic and perfectly flexible (no bending stiffness) is connected to the tip of towers, taken to be infinitely rigid. The girder is modeled as an Euler-Bernoulli beam, simply supported at both ends. The structure is subjected to a steady and uniform wind flow of velocity U. The dynamics of the structure are studied in the vertical plane, while the out-of-plane motions are ignored.
The main cable is modeled according the Irvine's theory [3], holding for sag-to-span ratios up to about 1:10. The cable assumes a parabolic profile under the initial dead load, and experiences vibrations in which the tangent component is statically condensed. Beam and cable are coupled by vertical hangers, here assumed to be massless, inextensible and continuously distributed along the span of the bridge (curtain of hangers) which apply equal and opposite distributed reactive forces r(s, t) to the two substructures.

Equations of Motion
Kinematics of the beam are assumed to be linear; nonlinearities are accounted for in the description of the aerodynamic forces only. The current configuration, both of the beam and of the cable, is described by the common transverse displacement field v (s, t), with s ∈ (−l/2, l/2) an abscissa along the beam axis and t the time. The equations ruling the in-plane motion, obtained eliminating the reactive forces r(s, t) between the cable and the beam equations, are: where E b I b is the flexural stiffness of the beam, E c A c is the axial stiffness of the cable, and m b and m c are the beam and cable linear mass density, respectively (and therefore the total bridge linear mass density is m := m b + m c ). Here T 0 is the prestress of the cable, assumed constant on s, and k = mg/T 0 , with g the acceleration gravity, is the prestress curvature, also assumed to be constant. The integral term accounts for the increment of tension, which depends on the time only.
In Equation (1), the non-inertial contribution of the external loads are decomposed in the damping and aerodynamic parts, as: where the subscripts α = b, c refer respectively the beam and cable. The Equations (1) are a generalization of the motion equations provided for the first time by Bleich et al. in [13], and recently re-examined in [18]. They include damping and aerodynamic loads, not present in the original formulation.

Damping Forces
A damping continuous model, accounting for both external damping (due to firm air-structure interaction) and internal damping (due to material dissipation), is used. This approach is different from that usually followed in literature, in which an 'overall' damping is attributed to the modes of the discrete system on experimental (or of literature) ground.
The contribution of the internal damping is related to two different physical mechanisms (see [37]): (a) elongation of the longitudinal fibers, both in cable and in beam, and (b) internal friction of the prestressed fibers of the cable. By adopting the Kelvin-Voigt model for mechanism (a), the elastic modulus E α (α = b, c) must be replaced by E α (1 + η α ∂ t ) , with η α an internal viscous damping coefficient, having physical dimensions of a time. By adopting the same rule for mechanisms (b), the pretension T 0 must be replaced by T 0 (1 + η c ∂ t ). In symbols: Accordingly, the time-differential damping operator is proportional to the elasto-geometric stiffness operator. Such a type of damping is usually referred as proportional, or Raleigh model, in the context of the finite-dimensional systems. The external damping is modeled by forces proportional to velocities, as: in which c eb and c ec are beam and cable external damping coefficients, from which c e := c eb + c ec is the coefficient of total external damping.
Note that the damping model introduced here is able to distinguish the beam and cable contributions. The four damping parameters can be calibrated to reproduce modal damping of as many modes.

Aerodynamic Forces
The aerodynamic part of the load is caused by the wind, which blows orthogonally to the plane of bridge, with a steady velocity U. It triggers forces on the beam and on the cable, acting both in-plane and out-of-plane of the structure. Due to the higher out-of-plane stiffness of the beam and to the smaller dimensions of the cable (which, being of circular shape, is exclusively subjected to drag forces), only in-plane forces acting on the beam are taken into account.
The quasi-static theory (see [29]) is used here to express such forces. The choice is corroborated by the discussion carried out in [28], where it is stated that: (a) the use of a nonlinear quasi-steady theory is justified in the range of low reduced frequencies, in which long-span bridges are expected to fall; and, (b) the previsions of the quasi-steady model are conservative, a few percent lower than those of a fully unsteady model. In particular, it is known that (see e.g., [43,44]), the quasi-steady assumption is valid only if the main frequency of periodic components of fluid forces, associated with vortex-shedding, is well above the dominant structural frequency, i.e., n s >> n 0 , where n 0 is the natural frequency of the crosswind vibration mode of the structure and n s = S t U/b, with S t the Struhal number and b the deck width. The occurrence of the frequency condition will be verified in the numerical application. When the inequality is not satisfied, it is needed to resort to a non-steady theory. However, while this latter is well codified in the linear range, few studies concern nonlinear dynamics (e.g., [44][45][46]).
In the framework of the quasi-steady theory, the forces depend on the structural velocityv only (not on the frequency of the oscillations), as result of the aeroelastic interaction. By expressing them in series of powers of thev U ratio, by assuming the cross-section is symmetric with respect the horizontal axis, and truncating the series to the fifth order, they read: Here, a is the air mass density, A i are dimensionless aerodynamic coefficients (depending on the shape of the cylinder cross-section), D is a transverse characteristic length, and b i := 1 2 a DA i ( i = 1, 3, 5) are dimensional coefficients, introduced to simplify the notation.

Dimensionless Equation of Motion
By introducing the following non-dimensional quantities: and parameters: the equations of motion are recast in non-dimensional form: In Equations (8), prime and dot denote differentiation with respect to non-dimensional space and time, respectively, and the tilde is dropped for convenience. It is noticed that nonlinearities are only due to aerodynamics. The dimensionless parameter Λ 2 is the well-known Irvine-Caughey parameter [3], accounting for geometric and mechanical properties of the cable. The dimensionless parameter ρ 2 (also used in [13,47,48] and others), describes the beam-to-cable stiffness ratio. For long bridges, it tends to zero, so that such bridges behave as a cable with an appended mass. However, still in this case, the ρ 2 parameter causes the occurrence of boundary layers close to the constraints and the point-loads, where the flexural behavior cannot be neglected. In this paper, attention will be devoted to bridges of medium length, in which the role of the beam is not negligible, as it will be discussed soon. A parametric study about the influence of both parameters on the free dynamic of the bridge is carried out in [18].

Linear Undamped Oscillations
First, the free dynamics of the undamped problem are addressed to evaluate the natural frequencies and modes. Although this problem has already been solved in literature (see, e.g., [13,18]), it is shortly resumed here for completeness.
The relevant equations of motion are: By separating the variables as v (s, t) = φ (s) e iωt , a space eigenvalue problem follows: whose eigensolutions (ω n , φ n ), n = 1, 2, ..., are the natural frequencies and modes, respectively. Due to the symmetry of the structure with respect s = 0, the modes are of anti-symmetric and symmetric types.

Anti-Symmetric Modes
In the anti-symmetric modes, the integral term in the equations disappears (no change of tension occurs), so that the system behaves as a string, stiffned by an in-parallel beam. Accordingly: and ω n = 2π 4π 2 n 4 ρ 2 + n 2 (12) are the nth anti-symmetric eigenpair. It is independent of Λ 2 , i.e., of the elasto-geometric properties of the cable, but depends on the parameter ρ 2 , i.e., on the stiffness ratio. It is worth noticing that the frequency of the first mode (n = 1) nearly doubles with respect to that of the string (ω n = 2π) when ρ 2 = 0.075. It is concluded that the stiffening effect of the girder is significant when ρ 2 = O 10 −2 .

Symmetric Mode
The nth symmetric natural mode is found to be: where the following definitions hold for the coefficient B i : and B N := B 0 + B 1 + B 2 is a normalization factor, such that φ n (0) = 1. In Equation (13), β 1 , β 2 are frequency-dependent wave-numbers, defined as follows: The unknown squared frequency ω 2 n is found as a root of the transcendent characteristic equation: Therefore, an explicit expression for the natural frequency of the symmetric modes cannot be provided.

Aeroelastic Stability Analysis
To investigate the behavior of the structure in the nonlinear field and close to the dynamic bifurcation, a nonlinear analysis is carried out. To solve the nonlinear problem, Equations (8), the asymptotic Multiple Scale Method (MSM) is used, by directly attacking the partial differential equations of motion. The method consists in expanding the dependent variables in a formal series of a small perturbation parameter (to be reabsorbed at the end of the procedure), as well as introducing independent time scales. Substitution of these expansions in the nonlinear equations and separation of terms of different orders, lead to linear differential equations to be solved in sequence. Here, the MSM is used as reduction method, as an alternative to the Center Manifold Method ( [49,50]). It provides a bifurcation equation ruling the dynamic on a center manifold of finite dimension (two-dimensional in the case of simple Hopf bifurcation).

Bifurcation Equation
According to the MSM, a small dimensionless parameter is introduced by a suitable ordering of variables and parameters. First, the displacement is rescaled as v (s, t) → v (s, t), so that is of the order of the amplitude of motion. Then, damping and aerodynamic coefficients are rescaled as: i.e., nonlinear coefficients are taken much larger than the linear ones (as suggested by the numerical application illustrated ahead), in order linear and nonlinear aerodynamic forces appear at the same level in the perturbation scheme. Moreover, independent time-scales t 0 := t, t 1 := t, ..., are introduced, so that v (s, t; ) = v (s, t 0 (t) , t 1 (t) , . . . ; ). Hence, by exploiting the chain rule, the following formal expansions for the time differential operators hold: (17) in which ∂ j := ∂/∂t j , j = 0, 1, . . .. Finally, the independent variable is expanded in series of as: By introducing rescaling and expansions in the equations of motion, and collecting the coefficients of the different powers of , the following perturbation equations are derived. Order 0 : Order : The generating problem, Equations (19), admits the monomodal solution: v 0 (s, t) = A (t 1 ) φ (s) e iωt 0 + c.c. (21) in which A (t 1 ) is an unknown complex modulating function, c.c. stands for complex conjugate terms, and (ω, φ (s)) is a real eigenpair of the undamped and unloaded model, determined in Section 3. When the -order problem, Equation (20), is considered, damping and aerodynamics terms appear on the right hand side of the equation, to be evaluated via Equation (21). Here, the internal damping η b would trigger point-moments at the boundaries, equal to − 2 η b ∂ 0 v 0 (±1/2); however, this last contribution is zero, because of the lower-order conditions 2 v 0 (±1/2) = 0. Equations (20) therefore read: where NRT stands for Non Resonant Terms, and Since the operator is singular (i.e., it admits a nontrivial homogeneous solution), the non-homogeneous problem cannot be solved for a generic know term. Therefore, a compatibility (or solvability) condition, known as Fredholm alternative, has to be enforced on it. This condition requires the ω-frequency forcing terms are orthogonal to the kernel of the self-adjoint operator, i.e., A detailed discussion on solvability conditions is provided in [51,52], in the framework of perturbation methods. From Equation (24), an ordinary differential equation on the t 1 −scale is obtained for the amplitude. By coming back to the true time t via backward rescaling, the equation assumes the well-known normal form for the Hopf bifurcation, i.e., where the real quantities d i are defined as follows: Finally, by letting A (t) = 1 2 a (t)e iϕ(t) , and separating real and imaginary parts in Equation (25), the real form of the bifurcation equation is derived: Equation (27)-b supplies ϕ = const, since structural nonlinearities, that would be responsible for phase modulations, have been ignored. The steady solutions a = const of Equation (27)-a are periodic solutions for the original system; the T−periodic solutions a (t + T) = a (t) are quasi-periodic solutions for the original system. Due to the normalization adopted for the natural modes, the real amplitude a describes: (a) the maximum modulus of the displacement in the anti-symmetric modes, and, (b) the displacement at s = 0 in the symmetric modes.

Linear Stability Analysis and Critical Wind Velocity
The dynamic of the bridge, close to the rest position, is ruled by the linear part of Equation (25), which admits the solution a (t) = a (0) e λt , where λ := d 0 + ud 1 . The position is asymptotically stable if λ < 0 and unstable if λ > 0. The condition of incipient instability occurs for λ = 0, that is for the wind velocity getting its critical (galloping velocity) value: at which a Hopf bifurcation occurs. Since, from Equations (26), d 0 < 0 (dissipative damping), bifurcation occurs only if d 1 > 0 (self-exciting aerodynamics, entailing A 1 < 0). Once the wind direction has been fixed, this last condition occurs only for some cross-section shape (just said, aerodynamically unstable). The critical wind velocity depends on the mode at which the bridge is oscillating, each velocity being associated to each mode. The lowest critical velocity is associated with the lowest mode (at which the internal damping is minimum), which can be either symmetric (Equation (13)) or anti-symmetric (Equation (11)). For this latter, an explicit formula can be derived, namely: The beneficial influence of the flexural stiffness should be noticed, since the critical wind velocity increases with ρ 2 . This is due to the fact that the internal damping is proportional to the stiffness, according to the Kelvin-Voigt model adopted.

Limit-Cycle Analysis
Steady non-trivial solutions of the bifurcation Equation (27)-a describe periodic oscillations occurring on limit-cycles. By lettingȧ = 0 and solving the biquadratic algebraic equation, two explicit solutions are found (to within an unessential sign): where: They describe two velocity-dependent branches (denoted by the ± index) of the bifurcation diagram a ± = a ± (u). Depending on u, two, one o zero nontrivial real solutions exist. In order to discuss the number of solutions, here the interesting case in which A 3 < 0 and A 5 > 0 is considered, which entails d 3 > 0 and d 5 < 0. This case is typical, for example, of rectangular cross-sections with a 2:1 aspect ratio (see [29,31]), according which the qualitative galloping response is plotted in Figure 2. Preliminary, it is observed that d 1 u + d 0 is negative in the subcritical u < u c range and positive in the supercritical u > u c range. Then, the sign of the discriminant is discussed. Let: be the value of velocity at which the discriminant vanishes, i.e. (u * ) = 0. Since u * < u c , it follows that: • when u < u * , it is (u) < 0, so that no real solutions a ± exist; • when u * < u < u c , it is (u) > 0; however, due d 1 u + d 0 < 0, it is (u) < d 3 , so that both the a + and a − roots are real; • when u > u c , in addition to (u) > 0, it is also d 1 u + d 0 > 0; consequently, (u) > d 3 , which entails that only a + is real.
The two branches coalesce at u = u * , at which a + = a − = u * − 2 d 5 d 3 and a turning point manifests itself.

Stability Analysis
To detect stability of the periodic solutions, the 1-dimensional dynamical system described by Equation (27)-a is perturbed, by letting a → a ± + δa, where δa is a small deviation from the steady solution solution a ± . By linearizing in the perturbation, the following variational equation is obtained: which decides about stability of a ± . From this, by using Equation (30), the following conclusions are derived: • the lower non-trivial solution a − is unstable in its domain of existence u * < u < u c ; • the upper non-trivial solution a + is stable in its domain of existence u > u * .
By summarizing (Figure 2), a subcritical Hopf bifurcation occurs at u = u c . Here, the trivial solution loses stability, and a family of unstable limit cycles of amplitude a − (u) bifurcates, whose existence is limited to a portion of the subcritical range. Indeed, at u = u * , a turning point modifies unstable to stable cycles, which extend their existence to the supercritical range. As a consequence, if the wind velocity is quasi-statically increased from below through u c , the rest position loses stability and the system jumps from zero to a finite value a + (u c ). A reverse jump occurs at u * when the wind velocity is decreased from above. An hysteresis can occur for wind velocity periodically changed. Such a behaviour is related to the opposite signs of the aerodynamic coefficients, b 3 and b 5 , causing a change of the qualitative behaviour of the system, namely: (i) at small amplitudes the cubic terms in the motion Equation (8) prevails over the quintic term; (ii) at higher amplitudes, the opposite occurs. This phenomenon is known in literature as "hard loss of stability", and observed in several different problems (see, e.g., [53], where the effects of the nonlinear damping on the post-critical behavior of the Ziegler's column are discussed).

Finite-Dimensional Model
Aimed to validate asymptotic results, a finite-dimensional model of bridge for galloping analysis is formulated. The finite-difference method is applied to transform the partial integro-differential equations (8) into a set of ordinary differential equations, to be numerically integrated for given initial conditions. The space domain (changed in [0, 1] for simplicity) is divided in N > 2 equispaced subintervals of amplitude ∆ = 1/N. The following notation is adopted: in which two external nodes, j = −1, N + 1, have been introduced. The space derivatives are expressed via the central finite differences: and the integral term via the trapezoidal rule: in which the boundary conditions f 0 , f N = 0 have been accounted for. The field equation in (8), accordingly, becomes: and the boundary equations read: Equations (37) and (38) are a set of N + 3 (ordinary and algebraic) equations for the N + 3 unknown nodal displacements v j (t), j = −1, . . . , N + 1.

Numerical Results
The behavior of the long-span suspension bridges is often well-described by the cable only. Here, however, the focus is addressed to cases in which the beam significantly influences the overall mechanics, as it happens when the dimensionless parameter ρ 2 is of order 10 −2 , (see, e.g., [13,18,48]). Therefore, a case study is considered: (a) the geometrical and mechanical characteristics are l = 195 m, E b I b = 2.4 × 10 9 Nm 2 , E c A c = 4.75 × 10 9 N, m b = 917.86 kg/m and m c = 177.6 kg/m; (b) the prestress is evaluated as T 0 = mgl 2 /8d, where d = 19.5 m is the cable sag; (c) the damping ratio is taken as ξ A = 1.5% in the antisymmetric mode, and ξ S = 1.3% in the symmetric mode (corresponding to η b = 0.0044 s, c eb = 55.88 Ns/m 2 , η c = 0.00022 s, and c ec = 6.29 Ns/m 2 , whose calibration is shown in Appendix A). A rectangular boxed cross-section of dimension 2 × 1 m (and so with an aspect ratio 2:1) is assumed for the beam. The aerodynamic dimensionless coefficients in Equation (5) are drawn from literature (obtained by wind tunnel tests in [30]), namely A 1 = −3.47, A 3 = −490.25, A 5 = 44744.21. The air mass density is a = 1.25 kg/m 3 .
This example is representative of a class of bridges, whose nondimensional parameters are: 223. All numerical values assumed by the parameters are consistent with the ordering performed in the perturbation analysis.

Linear Stability Analysis
The modal properties of the undamped model were first analyzed. The first natural mode was found to be anti-symmetric (defined by Equation (11) for n = 1) and the corresponding natural frequency is ω = 2.12 rad/s. The second natural mode was found to be symmetric (defined by Equation (13) for n = 1) and the natural frequency is ω = 3.91 rad/s (see Figure 3). The critical conditions corresponding to the antisymmetric and symmetric modes are shown in Table 1. It is observed that the bridge galloped in its first antisymmetric natural mode, in corresponding of which the critical wind velocity assumed the minimum value U c = 34.13 m/s. Since the critical velocities corresponding to the antisymmetric and symmetric modes were sufficiently far from each other, nonlinear interactions between modes, leading to multimodal galloping, were excluded.
It was observed that the quasi-steady theory was valid in the case study considered. Indeed, the natural frequency of the first crosswind vibration mode was n 0 = ω c /2π = 0.33 Hz, while the vortex-shedding frequency was n s = 1.02 Hz at the critical galloping velocity U c = 34.13 m/s, according to a Struhal number S t = 0.06 (as indicated by [54] for a rectangular cross-section with a 2:1 aspect ratio) and a deck width of 2 m. Note that the quasi-steady theory worked mainly due to the narrow width of the cross section.

Galloping Response
The postcritical behavior relevant to the case study is described by the bifurcation diagram in Figure 4. There, the (dimensionless) steady amplitude a of the anti-symmetric mode, experienced by the bridge on the limit-cycle, is plotted vs the (dimensionless) wind velocity u. It is seen that, the hard loss of stability caused amplitudes of about 1% of the length (i.e., of the order of 2 m), as far the critical wind velocity was passed. Such a bad mechanical behavior is due to the unfavorable shape of the cross-section ( [23,30]).

Validation by Finite-Dimensional Model
The asymptotic solution was validated against numerical results, relevant to the discretized version (Equations (37) and (38)) of the nonlinear continuous model (Equations (8)). The space interval was divided in N = 128 subintervals, and numerical integration carried out, via commercial routines, for given initial conditions. Regardless of these latter, after a transient had been exhausted, the motion stabilized either on a limit cycle or comes back to the trivial state. In order to compare the transient motion v j (t) to that predicted by the MSM via the amplitude a (t), an initial deflection shape coincident with the first natural mode of the undamped bridge was chosen.
Four typical time-histories are shown in Figure 5 for different wind velocities. It was observed that: (a) when u < u * (Figure 5a) the motion stabilized on the trivial solution a = 0; (b) when u * < u < u c (Figure 5b) the motion stabilized either: (I) on the trivial solution a = 0, if the initial condition was sufficiently small, or (II) on a limit cycle of amplitude a = 0, if the initial condition was sufficiently large; (c) when u > u c (Figure 5c) the motion stabilized on a limit cycle of amplitude a = 0. All the time-histories were compared with the numerical integration of the bifurcation equation Equations (27). It was seen that this latter captured at the best extent the exact (numerical) solution.
To compare the numerical to the asymptotic limit-cycles, the numerical amplitudes were extracted by the recorded time-histories for different wind velocities, and bullets superimposed to the asymptotic bifurcation diagram (see Figure 5d), showing an excellent accordance. Note the power of the perturbation methods in providing accurate predictive formulae (as for the galloping velocity, Equation (29) and for the amplitudes of the limit-cycles, Equation (30)).

Conclusions
The nonlinear aeroelastic behavior of suspension bridges, undergoing dynamic in-plain instability (galloping), has been analyzed. A nonlinear continuous model of bridge has been formulated, made of a visco-elastic beam and a parabolic cable, connected each other by axially rigid suspenders, uniformly distributed. A damping continuous model, accounting for both external damping and internal damping (according to the Kelvin-Voigt rheological model), has been used. The structure is loaded by a uniform wind flow, acting normally to its plane. The nonlinear aeroelastic effects have been evaluated via the quasi-static theory, while structural nonlinearities have been ignored. The nonlinear equations ruling the dynamics of the aeroelastic system have been tackled by the Multiple Scale Method, directly applied to the partial differential equations, to get bifurcation equations for the Hopf bifurcation. The limit-cycle amplitude and its stability have been evaluated as functions of the mean wind velocity. The analytical solutions have been validated against numerical integration of the finite-difference discretized equations.
Numerical results have been obtained for a case study of suspension bridge, whose beam has a rectangular boxed cross-section, with an aspect ratio 2:1. It has been found that the bridge gallops in its first antisymmetric natural mode (having natural frequency ω = 2.12 rad/s). An explicit formula has been derived for the critical wind velocity, assuming the value U c = 34.13 m/s. The validity of the quasi-steady theory, when related to the vortex-induced vibration phenomenon, has been checked. The postcritical behavior has been investigated, and the bifurcation diagram provided. The well-known (and dangerous) phenomenon of "hard loss of stability" has been found to occur, according which, as far the critical wind velocity is passed, the motion amplitude reaches values equal to 1% of the bridge length (i.e., of the order of 2 m). Such a bad mechanical behavior is due to the unfavorable shape of the cross-section. All the asymptotic results, relevant to the transient as well to the steady motions, have been fully validated by the numerical analysis. This work is susceptible of farther insights, namely, (a) analysis of the influence of the flexibility of the pylons on the bridge dynamics; (b) study of the interaction between vortex-induced vibrations and galloping; (c) validation of the results by wind tunnel tests. Concerning item (a), elasticity of supports could cause (in principle, and to be investigated in the filed of the real bridges) an inversion between the now dominant antisymmetric and the first symmetric mode, as discussed in [6]. Such an inversion could also trigger (d) nonlinear interaction between the two modes. Concerning item (b) and (c), study of literature, as described in [55] could be exploited.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

MSM Multiple Scale Method NRT
Non Resonant Terms c.c. complex conjugate

Appendix A. Calibration of the Damping Coefficients
In aeroelastic stability analyses, the damping assumes a fundamental rule. In particular, for the j natural mode φ j (s), the modal damping ratio of the bridge is ξ j = −d 0 /ω j (according to Equation (25)), where the real coefficient d 0 depends by the internal and external damping parameters, η b , c eb (of the beam), and η c , c ec (of the cable) (see Equation (26)). These coefficients are chosen as follows: (a) the beam and the cable are considered independent each from the other; (b) a shape of vibration is enforced, equal to the modal shape of the bridge, and the modal quantities c, ω, m evaluated, parametrically dependent on internal and external damping coefficients, for the beam and the cable; (c) the beam and the cable modal damping ratios are computed and their value are assigned on experimental basis; (d) a system of equations is solved for the unknown damping coefficients; (e) the damping ratios relevant to the assembled bridge are evaluated for check purposes.
Referring to the j natural mode φ j (s) (either antisymmetric or symmetric) of the bridge, the modal quantities for the beam (index b) and cable (index c) are: