Nonlinear interaction between self- and parametrically excited wind-induced vibrations

The aeroelastic behavior of a planar prismatic visco-elastic structure, subject to a turbulent wind, flowing orthogonally to its plane, is studied in the nonlinear field. The steady component of wind is responsible for a Hopf bifurcation occurring at a threshold critical value; the turbulent component, which is assumed to be a small harmonic perturbation of the former, is responsible for parametric excitation. The interaction between the two bifurcations is studied in a three-dimensional parameter space, made of the two wind amplitudes and the frequency of the turbulence. Aeroelastic forces are computed by the quasi-static theory. A one-D.O.F dynamical system, drawn by a Galerkin projection of the continuous model, is adopted. The multiple scale method is applied, to get a two-dimensional bifurcation equation. A linear stability analysis is carried out to determine the loci of periodic and quasi-periodic bifurcations. Limit cycles and tori are computed by exact, asymptotic, and numerical solutions of the bifurcation equations. Numerical results are obtained for a sample structure, and compared with finite-difference solutions of the original partial differential equation of motion.


Introduction
The study of the aeroelastic behavior of slender structures is a fascinating topic, of high scientific and technical value. The literature is rich of studies concerning specific flexible structures (cables [1][2][3][4][5][6][7], beams [8][9][10][11][12], plates [13][14][15][16]) and general aeroelastic phenomena (galloping [17][18][19][20][21], flutter [22][23][24][25], vortex-induced vibrations (VIV) [26][27][28][29][30]). The instability and bifurcation events can be related to different kinds of excitation. Self-excited autonomous systems, as structures subjected to steady wind, are prone to Hopf bifurcations, caused by zeroing of the total damping. Nonautonomous systems, as structures subjected to turbulent wind, are prone to parametric excitation phenomena, leading to divergence, flip and Neimark-Sacker bifurcations. Depending on the nature of the loads, the different kinds of excitation can interact. Some attention has been devoted in literature to interactive aeroelastic phenomena, as galloping parametric excitation [31][32][33][34][35][36][37][38][39][40] or galloping vortex-induced vibrations [41][42][43][44][45]. In particular, as regard the former class, the principal resonance of a single-degree-of-freedom system with two-frequency parametric and self-excitation is investigated in [31,34]; the method of multiple scales is used to determine the equations of modulation of amplitude and phase, and qualitative analyses are performed out to study steady state, limit cycle, and torus responses. In [32], analytical investigations of the system under parametric, self-excited, and inertial excitations were carried out. In [33], the galloping of tall prismatic cantilevered structures, due to unsteady wind, is analytically studied. The unsteady wind is considered by adding a time varying wind speed component to the mean wind speed. Consequently, the structure is subjected to multi-harmonic external and parametric excitation. The multiple scale method is used to study the effect of primary and secondary resonances on the galloping response of the structure. Starting from this last work, the parametric, external, and self-excitation of one/two towers under turbulent wind flow is studied in [35,36,39] and in [37,38], focusing the attention to the periodic and quasi-periodic galloping motions. It was shown that the unsteady component of wind can cause a significant decrease in the critical speed; however, the contribution of the unsteady component is less prominent at high wind velocities, at which large amplitude oscillations occur. It was detected, moreover, that periodic and quasi-periodic motions exist according to suitable combinations of the steady and unsteady wind parameters.
Here, the effect of a wind flow acting on a slender prismatic structure is investigated. A class of cross sections is considered, suffering sub-critical bifurcation followed by a regain of stability at high amplitudes, known in the literature as the hard-loss of stability [8,46,47]. When the turbulence frequency is nearly double a natural frequency of the structure, this latter is parametrically excited at the primary resonance. Therefore, self-excitation and parametric excitation interact, both in linear and nonlinear fields. The investigation, mainly carried out by perturbation methods, is aimed at evaluating: (a) as the turbulent component reduces the galloping critical velocity, and vice versa, (b) as the steady wind component modifies the parametric excitation instability domain. Moreover, (c) the analysis of the limit cycles generated by the flip bifurcations, and that of tori, generated by the Neimark-Sacker bifurcations, are also of interest.
The paper is organized as follows. In Sect. 2, a continuous aeroelastic model of a slender prismatic viscoelastic structure is reduced to a single-D.O.F system via the Galerkin method. In Sect. 3, the bifurcation equations, ruling the slow flow of the amplitude and phase of the involved mode, are derived via the multiple scale method. In Sect. 4, the linear stability analysis is carried out, and a three-dimensional stability domain is built-up. In Sect. 5, a nonlinear analysis of the bifurcation equation is performed; in particular, the bifurcation equations governing the slow-slow-flow on a torus, are derived. In Sect. 6, numerical results are reported for a sample structure; they are compared with direct numerical integration of the original partial integrodifferential equation of motion. In Sect. 7, some conclusions are drawn. In "Appendix A", the sample system is described. Two other "Appendixes", supplying details, close the paper.

Aeroelastic model
A continuous model of a slender prismatic visco-elastic structure, of length l, restrained to belong to the (i, j)plane and constrained at its ends A, B, subjected to a turbulent wind flow of velocity U (t)k, acting normally to the plane, is formulated (Fig. 1). The structural model is assumed to be linear (since, as it is well-known, elastic nonlinearities mainly affect the frequency, but not the amplitude of the aeroelastic response); the aerodynamic model is instead nonlinear. The generic equation of motion and the boundary conditions are: In these equations: v (s, t) is the deflection of the prism in the j-direction at the abscissa s ∈ (0, l) and time t ∈ (0, ∞); m is the mass per unit length; L e and L v are elastic and viscous linear operators, respectively; B eH and B v H are elastic and viscous linear operators acting at the boundaries; c e is an external damping coefficient; p a (s, t) are aerodynamic loads in the j-direction; a dot indicate partial differentiation with respect to the time. Equation (1) can be referred as a metamodel of prismatic structures, describing a class of real important structures, as columns, towers, chimneys, and bridges.

Aerodynamic model
The aerodynamic forces are determined according to the quasi-static theory [8]. This theory holds when the By assuming the wind velocity U (t) uniform in space, the quasi-static theory supplies the forces as nonlinearly depending on the structural velocityv (s, t) (which is responsible for modifications of the attack angle). By assuming thatv U 1, and the structure cross section is symmetric with respect to the horizontal axis, they are given by the odd-power series: where: a is the air mass density and A i are dimensionless aerodynamic coefficients, depending on the crosssection shape, of which D is a characteristic length. The aerodynamic forces include all the nonlinearities of the problem. By letting: in whichŪ is the (leading) steady-state wind velocity, and u(t) the (small) turbulent component, and linearizing Eq. (2) in the ratio ū U 1, it follows: where b i := 1 2 a D A i (for i = 1, 3, 5) has been introduced to simplify the notation.
In the following, it will be assumed that the excitation is harmonic, namely: withû Ū the amplitude and Ω the frequency of the turbulent component.

Single degree-of-freedom system
The equations of motion and boundary conditions Eq.
(1), by accounting for the aerodynamic forces (4) and the turbulence law (5), are recast in the following nondimensional form: where: and the following positions have been introduced: v * := v l , s * := s l , t * := tω r , are nondimensional operators; ω r is a reference frequency and T r a reference mechanical characteristic, having the dimension of a force; finally a dot denotes differentiation with respect t * .
The continuous system (6) is discretized as a single-D.O.F. system via the Galerkin method, by assuming: Here q (t * ) is a Lagrangian parameter and φ (s * ) is a natural mode of the undamped system, normalized according to max |φ| = 1, of frequency ω * , satisfying the boundary value problem: By weighting the residuals in the domain and at the boundaries, the following condition is enforced (Virtual Work Equation): from which a nondimensional ordinary differential equation is drawn: In this equation, the structural and wind coefficients are defined as: with I m := 1 0 φ 2 ds * . To simplify the notation, the star will be dropped ahead.

Bifurcation equation
The one-D.O.F. system (12) is prone to a Hopf bifurcation (galloping, in the technical literature) when b 1 < 0 (i.e., A 1 < 0) and the steady windŪ attains a critical valueŪ c . On the other hand, it is also susceptible of parametric excitation when the turbulence has a suitable frequency Ω and the intensityû overcomes a threshold, which depends on Ω. Here attention is devoted to the case of primary parametric resonance, for which: where ω is a natural frequency of the aeroelastic system, and σ a small detuning. The two types of bifurcation can interact in the nonlinear field, giving rise to a multiple bifurcation, described by the three bifurcation parametersŪ ,û, σ. It is worth noticing, however, that, while self-excitation can occur alone (when the turbulence is absent), the reverse does not hold, sincē U cannot be zero nor small in the model (due to the hypothesisv U 1). To analyze the nonlinear phenomenon, the Multiple Scale Method (MSM) is applied here [48] to capture the essential aspects of bifurcation (see, e.g., [49] for several applications of the MSM in bifurcation theory). Internal resonance cases, in principle possible, are excluded in this paper.

Cartesian form
To obtain the Cartesian form of the bifurcation equation, a new variable Z (t) is defined via: which transforms the non-autonomous Eq. (17) in an autonomous system. Then, is posed, and the real and imaginary parts of the equation are separated. The following two-dimensional dynamical system is obtained: where: is the Jacobian matrix at the origin, and n = n x , n y T is the vector of the nonlinear terms, defined as follows: Accordingly, the motion law (9) is expressed, at the leading order, by: where use has been made of Eqs. (16) and (14).

Polar form
To derive the polar form of the bifurcation equation (17), the complex amplitude is expressed as: with a (t) the modulus of the amplitude and ϕ (t) the phase. Here, according to Eq. (9), a(t) represents the motion of the material point at which φ attains its maximum modulus. Consequently, Eq. (17) splits into two real equations: (26) in which γ := 2ϕ (t) − σ t is the phase difference. The motion law (9), at the leading order, reads: Amplitudes and phases in the two representations are related by a = 2

Linear stability analysis
The rest position of the structure v (s, t) = 0 ∀ (s, t), is described by the trivial solution A = 0 of the bifurcation equation (17). Since the polar form (26) is singular at a = 0, the Cartesian form (21) is used. Stability of the origin is governed by the eigenvalues of the Jacobian matrix J 0 , which read: where the discriminant is: For later purposes, the associated (right) eigenvectors, solving (J 0 − λ ± I) x = 0, together with the left eigenvectors, solving the transpose conjugate eigenvalue problem J 0 −λ ± I y = 0, are evaluated as: In the absence of turbulence (û = 0 , σ = 0), the two eigenvalues Eq. (28) are real and coincident. A Hopf bifurcation occurs when λ + = λ − = 0, at the critical galloping velocityŪ =Ū c : The coefficients d 0 and d 1 (see Eq. (18)) depend on the the cross-section shape. Here, the interesting case d 0 < 0 (section prone to galloping) and d 1 > 0 (subcritical Hopf bifurcation) is considered.
The turbulent component of wind modifies the bifurcation condition, according to the nature of the eigenvalues (28). To investigate its effects, linear stability is studied in the σ,û,Ū -space, in which: (i)Ū is the distinguished parameter, and, (ii)û and σ are the splitting parameters (since they rule the coalescence, as it will appear clear soon). Results are qualitatively summarized in Fig. 2, where the linear stability domain is represented (a) by a 3D view, or (b,c,d) by cross sections parallel to the coordinate planes (here the eigenvalues are also sketched in the different regions).
In the û, σ -plane, the straight lines σ = ±d 1û are the locus at which Δ = 0. These lines separate the plane in an internal R P and an external R Q region, in which periodic and quasi-periodic motions, respectively, arise whenŪ takes a suitable value, as discussed ahead.
-In the internal region R P , it is û, σ > 0, so that the eigenvalues are real and distinct. Here, whenŪ is increased from zero, and it is such that λ ± = 0, two successive flip bifurcations occur, at which a periodic response, of period double of that of excitation, is triggered. This happen on the geometrical locus: which is a cone in the parameter space. The lowest bifurcation (λ + = 0) manifests at a steady wind velocityŪ û, σ <Ū c , so that turbulence reduces the galloping velocity. -In the external region R Q , it is û, σ < 0, so that the eigenvalues are complex conjugates. Here, whenŪ is increased from zero, and it is such that Re[λ + ] = Re[λ − ] = 0, a Neimark-Sacker bifurcation occurs, at which a quasi-periodic response is triggered. This happen on the geometrical locus: i.e., atŪ =Ū c ∀ û, σ . -At the interface between R P and R Q , it is û, σ = 0, so that the eigenvalues are real and coincident. They vanish atŪ =Ū c , where the Neimark-Sacker bifurcation degenerates in a flip bifurcation. The splitting parameters decide about the coalescence.
It is worth noticing that, referring to the slow flow ruled by the bifurcation equation (17), the flip bifurca-tion appears as a divergence (simple zero eigenvalue) and the Neimark-Sacker as a Hopf bifurcation. At the intersections of the relevant manifolds, a doublezero (Bogdanov-Takens) bifurcation occurs, where the Hopf bifurcation degenerates; here, the Hopf manifold dies (see, e.g., [51]). This mechanisms clearly appears in Fig. 2b, c. A mechanical interpretation of the instability domain can be given. WhenŪ =Ū c , the aerodynamic damping balances the structural damping, so that the monomodal response of the structure is undamped. Due to the parametric excitation, there exists a sector of theŪ =Ū c plane, in which the system is unstable, according to the well-known Mathieu equation (see Fig. 2a). When U <Ū c , the structural damping prevails over the aerodynamic damping, so that the unstable zone moves apart from the σ -axis. Consequently, a larger turbulence, or a smaller detuning, are needed to make the system unstable, namely: . When (ideal case, beyond the validity limits of the theory)Ū = 0, damping is only of structural type, so thatû takes its maximum valueû c =û c (σ ) (see Fig. 2b). Out of the parametric excitation zone, at U =Ū c , the phenomenon is lead by the self-excitation. However, the periodic motion experienced by the system under non-turbulent wind, is modified by the turbulence, which introduces a second frequency, Im[λ ± ], into the response, thus transforming the limit cycle in a torus.

Nonlinear analysis
Bifurcation analysis calls for determining: (a) the limit cycles (periodic solutions) generated by the flip bifurcation; (b) the tori (quasi-periodic solutions) generated by the Neimark-Sacker bifurcations. The periodic solutions, which are equilibrium points for the bifurcation equations, are found by analytically or numerically solving nonlinear algebraic equations. The quasiperiodic solutions, which are periodic orbits for the bifurcation equations, are found (i) via numerical integration of the ordinary differential equations, and (ii) asymptotically, by carrying out a new perturbation analysis (as, e.g., done in [52], and, then, in [37,38] for a different engineering application).

Limit cycles
Periodic motions are the equilibrium points a = const, γ = const of the polar Eq. (26). By requiringȧ = 0 andγ = 0 and by zeroing the right-hand side, sin γ and cos γ are obtained and the variable γ condensed using the relation sin 2 γ + cos 2 γ = 1. The resulting equation for the unknown a is the following: In general, the roots of the Eq. (34) cannot be determined analytically. However, this is viable in the particular case of perfect resonance (i.e., σ = 0), for which the solutions read: together with γ 1 = 0 and γ 2 = π , respectively. The Eq. (35) consist in multi-valued functions, whose existence depends on the wind parameters and on the nature of the aerodynamics coefficients (i.e., by the cross-section geometry).
To detect stability of the periodic solution, the dynamical system Eq. (26) is perturbed by letting a (t) = a e + δa (t) and γ (t) = γ e + δγ (t), where δa, δγ are small deviations from the equilibrium solution (a e , γ e ). Substituting in the equation and linearizing in the perturbation, leads to the following variational equation: where J e := J i j (a e , γ e ) is the 2 × 2 Jacobian matrix at the equilibrium, whose coefficients are: By denoting with μ the eigenvalues of J e , the periodic solution is stable if Re (μ) < 0 for all μ, and unstable if Re (μ) > 0 for at least one μ. In the resonant case, the Jacobian matrix is diagonal, with eigenvalues J 11 (a i± , γ i ) and J 22 (a i± , γ i ).

Tori
Quasi-periodic motions are found as limit cycles for the bifurcation equations, via a perturbation analysis. Use is made of the Cartesian form (21) of the bifurcation equation and the MSM is newly applied to build-up periodic solutions arising at Neimark-Sacker points.
Here, the eigenvalues of the Jacobian matrix J 0 are purely imaginary, i.e., λ ± = ±i , with := 1 2 and whereŪ =Ū c . At the bifurcation, the motion is: where R (t) = 1 2 r (t) e iθ(t) is a complex amplitude, and where use has been made of Eq. (30). By lettinḡ U =Ū c + ΔU to explore the neighborhood of the bifurcation point, and newly applying the MSM (see "Appendix B"), the following polar form of the bifurcation equation (sometimes referred to as ruling the slow-slow-motion), is drawn: Quasi-periodic solutions are found as equilibrium points r = const of Eq. (39), i.e., where θ 0 is an initial phase and: is a frequency correction.
The motion has therefore two frequencies, the driving one Ω 2 , and the modulating one * . At any given s =s, the modulating amplitude: exists only when σ > σ M û and it spans the range: whose extreme values represent the lengths of the axes of an elliptical trajectory traveled in the (X ,Y ) statespace (see the later Fig. 13). It should be noticed, however, that the perturbation solution carried out here (local bifurcation analysis) is unable to capture possible interactions (global bifurcations analysis) between tori and limit cycles. Therefore, it is expected that it loses validity at high amplitudes of motion. The question will be investigated via numerical analyses.

Numerical results
Here, by referring to a sample structure described in "Appendix A", the nondimensional coefficients appearing in Eq. (18) are taken to assume the following values: d 0 = −0.14, d 1 = 0.19, d 3 = 4717.47, d 5 = −9.19 × 10 7 . These numerical values are consistent with the ordering performed in the perturbation analyses.
The asymptotic solutions, previously determined, are validated ahead against numerical results. These latter consist of: (i) numerical integrations of the nonlinear partial differential equation of motion Eq. (6), discretized by finite-differences (see "Appendix C"), and, (ii) numerical integrations of the bifurcation equations Eq. (26) and (21), for given initial conditions.

Linear analysis
In the case study, the first natural frequency is ω = 2.12 rad/s (ω * = 8.77 in dimensionless form) and the corresponding galloping velocity (in the absence of turbulence) isŪ c = 34.13 m/s (U * c = 0.7 in dimensionless form). The influence of the turbulent component of the wind is illustrated in the (quantitative) linear instability domain in Fig. 3, where (a) a 3D-plot and (b-d) planar contour plots are shown. Plots are extended beyond the limits of validity of the theory, which requiresû Ū .
It is seen that the critical wind velocity decreases when the turbulence increases and/or the detuning decreases. For example, whenû = 0.15 and σ = 0.02, it is

Limit cycle analysis in the resonant case
The nonlinear behavior of the system at the perfect resonance (σ = 0) is first analyzed in the (Ū ,û)-plane. The relevant bifurcation analysis calls for determining limit cycles and analyzing their stability. Amplitude and phase of the limit cycles are defined by the explicit expressions a i± , γ i (i = 1, 2) in Eq. (35). First, the existence domain of each solution branch in Eq. (35) is studied in the parameter plane (see Fig.  4). In the same figure, the bifurcation locus, made by two generatrices of the cone, are marked in red. The plot displays the existence of a complex scenario, in which from zero to four limit cycles coexist (the cir-cled numbers denote the number of solutions in each region). To detect the stability of the periodic solutions, the sign of the corresponding eigenvalues J 11 and J 22 of the Jacobian matrix Eq. (36) is evaluated: They are found to be both negative only when evaluated at the (a 2+ , γ 2 ) limit cycle, which is therefore the unique stable bifurcated branch.
Planar bifurcation diagrams are built-up and shown in Figs. 5 and 6, namely: (i) a versusŪ , for selected values ofû (ticks I, II in Fig. 4), and, (ii) a versusû, for selected values ofŪ (ticks III to VI in Fig. 4).
The a,Ū diagrams are commented first. In case (I) there exist three families of limit cycles, all branching from the trivial path. Two of them are close to the galloping curve (û = 0) (one on the left, the other on the right), generated, whenŪ is increased, at entering and coming out the cone of Fig. 3. A third family bifurcates fromŪ = 0, at which the theory loses validity; how- ever, the existence of this branch has been confirmed by numerical results, as it will be discussed later. The dangerous phenomenon of hard loss of stability, already manifesting in absence of turbulence, is observed. It is generated by an (unstable) subcritical bifurcation, followed by a regain of stability, which entails a finite jump of the amplitude, just passed the critical velocity. In the numerical example, this amplitude is about 1% of the length (i.e., double of the cross-section height). Moreover, the existence of a stable portion of the lowest branch, creates a potentially dangerous competitive attractor at high amplitude, although existing in a very narrow range of the velocity (0.25, 0.26). In case (II), in which the turbulence is higher, the scenario changes, since the stable and unstable portions of the two lowest curves weld together, the highest remaining undisturbed. Here, once the trivial path loses stability, the system jumps to an even higher-amplitude limit cycle, equal to about 1.5% of the length. The dangerous effect of the turbulence is thus highlighted.
When the bifurcation diagrams a,û are plotted, Fig. 6 is obtained. In all cases, just one family of limit cycles exists, generated, forû increasing, at the (unique) crossing of the cone. WhenŪ <Ū c (cases (III) and (IV)), a hard loss of stability is observed. However, in case (III), in whichŪ is far from the galloping value, the stable branch fails to overlap the sub-critical region. This aspect could be interpreted as due to the prevailing influence of the parametric excitation with respect the galloping phenomenon. In contrast, in case (IV), the stable branch extends on a small portion of the super-critical region. WhenŪ ≥Ū c (cases (V) e (VI)), the trivial solution is unstable, so that motion stabilizes on limit cycles of great amplitude, which increases with the turbulence.
To sketch the evolution of the system in presence of one or more attractors, the transients asymptotic solutions, obtained by direct numerical integrations of the bifurcation equation (26), are investigated in the statespace at varying of the initial conditions. By referring to the bifurcation diagram in Fig. 5I, the (a, γ ) trajectories are shown in Fig. 7 for the steady wind velocities: (a)Ū = 0.5, at which two stable steady solutions coexist; (b)Ū = 0.84, at which only one stable steady solution exists. It is seen, in Fig. 7a, that trajectories originating at low amplitudes (i.e., below that of the unstable solution) are attracted by the trivial solution, irrespectively of γ ; in contrast, trajectories originating at higher amplitudes can be attracted by one of the two stable solutions, according to the initial gamma. In Fig.  7b it appears that there exist a trajectory leading to the stable solution which attracts all the surrounding orbits. Finally, the steady-state and transients asymptotic solutions, are validated against exact finite-difference solutions of the partial differential equations (6). In particular, for given initial conditions (a = 0.01, γ = 0.05 in Fig. 7b), the exact response and the asymptotic a (t) time-histories are compared in Fig. 7c. After a transient has been exhausted, the motion stabilizes on a limit cycle of amplitude a, depending on the wind parameters. When the limit cycle amplitudes are extracted by the recorded responses, relevant to different wind parameters, the (green) bullets reported in Figs. 5I and 6IV are found. In all cases, the agreement between exact and asymptotic solutions is excellent.

Limit cycle analysis in the quasi-resonant case
The nonlinear periodic motions in the quasi-resonant case are analyzed. Differently from the resonant case, explicit forms of amplitude and phase of the limit cycles are not available, so that numerical solutions are obtained. The same values of the turbulence,û = 0.1, 0.2, already considered in Fig. 5a, b, are considered in Figs. 8 and 9, respectively. In each of them, three different values of detuning σ have been considered, such that the û, σ point in the parameter space is (a) inside the region R P , (b) at the interface between R P and R Q , (c) inside the region R Q . For comparison, in addition to the galloping curve (thin green line), the curves relevant to σ = 0 (thin red lines) have also been reported. At the bottom of the figures, the phases are plotted. Referring to Fig. 8a, it is seen that, in the periodic region R P , the limit cycle are only slightly affected by the small detuning. In contrast, when the detuning is large (Fig. 8b, c), due to the disappearance of the flip bifurcation points B 2 , B 3 from the trivial path, limit cycles not crossing this path arise, generated by welding of the former curves, which, however, leaves almost unaltered their stable parts. Concerning the phase, it appears that the detuning alters the constant values γ 1,2 = 0, π of the perfect resonance, so that the phase changes with the parameter. The phases, as the amplitudes, show a welding mechanism. A very similar behavior is exhibited in Fig. 9.
Some bifurcation diagrams, plotting the steady amplitude and phase of the limit cycle versus the turbulence amplitudeû are shown in Fig. 10 for an assigned detuning σ = 0.05 and steady wind velocity: (a) A comparison between quasi-resonant (black lines) and resonant (thin red lines) cases is made. It is observed, with respect to the resonant case, that the detuning rounds the cusps of the amplitude bifurcation diagrams, but does not involve qualitative changes. It mainly modifies the unstable branches of the bifurcation diagrams, pushing them forward.
The plots of the steady amplitude and phase of the limit cycle versus the detuning parameter σ are shown in Fig. 11 for an assigned turbulence amplitudê u = 0.1 and steady wind velocity: (a)Ū = 0.85Ū c ; (b)Ū =Ū c ; (c)Ū = 1.15Ū c . In these bifurcation diagrams it emerges that the detuning strongly influences the amplitude and phase of the limit cycles. In particular, there exist two closed curves (for amplitudes and phases) in the subcritical range and just one curve in the supercritical range. This circumstance proves the existence of a tubular surface in the Ū , σ, a -space, emerging from the plane of the trivial solution (a = 0), whose transverse and longitudinal sections are the bifurcation diagrams in Figs. 8, 9 and 11. The asymptotic solutions are validated against numerical results: green bullets are superimposed to the asymptotic bifurcation diagrams in Figs. 8c, 9c, 10c and 11c. The exact results are found to be in excellent accordance with the analytical predictions.

Torus analysis
The nonlinear quasi-periodic motions, occurring on tori of the state-space, are studied. These solutions are gen-  (21)), and asymptotically, as described in Sect. 5.2. Since the bifurcation is of codimension-1, the bifurcation diagram in Fig. 12a is sufficient to explain the phenomenon. Here, in addition to the quasiperiodic solution (exact and asymptotic), represented via the interval [a min , a max ], also the periodic solutions are plotted. AtŪ =Ū c , a torus bifurcates from the trivial solution. It is of subcritical type, and manifests again the hard loss of stability phenomenon. After regaining stability, a collision of the torus with a limit cycle occurs, leading to the disappearance of the torus. Exact results (green segments and green bullets) clearly show the stable part of both torus and limit cycle branches. The asymptotic solution displays the qualitative mechanism. The torus branch plotted in Fig. 12a is just that corresponding to a û, σ pair; if this latter is changed, other branches of similar type are found. If a crosssection atŪ =Ū c of all these curves is made, together with σ = const orû = const, the plots in Fig. 12b, c are obtained. They, in addition to the collision already described, show the existence of stable tori external to the limit cycle tube (see Fig. 12c).
The transients asymptotic solutions, drawn by the bifurcation equation (39) through direct numerical integration, are investigated in the (X, Y ) state space, at varying of the initial conditions. Although the motion essentially occurs in the one-dimensional space r (t), as the first of numbered Eq. (39) proves, the twodimensional representation adopted here is more clear, since it allows distinguishing closed orbits from equilibria. In particular, the periodic (stable and unstable) solutions are represented by inclined ellipsis, whose axis lengths are the maximum and minimum value of the amplitude, according Eq. (45). By referring to the bifurcation diagram in Fig. 12a, the (X, Y ) trajectories are shown in Fig. 13 for the steady wind velocities: (a)Ū = 0.59, at which two periodic stable solutions coexist; (b)Ū = 0.75, at which only one periodic stable solution exists. Figure 13a shows that, when the initial conditions are inside the unstable closed orbit, the motion is attracted by the trivial solution. In contrast, when they are external to the unstable orbit, the trajectories are attracted by the stable orbit. In Fig. 13b, it appears that all the trajectories dye on the stable orbit.
Moreover, concerning the asymptotic solution, when compared with the exact one, it is observed that: (i) it captures the qualitative trend, but it is not very accurate for steady wind velocities far from the galloping value; (ii) it is unable to describe the interaction with limit cycles (as a matter of fact, it exists after the collision). Item (i) clearly appears in Fig. 14, where a comparison between the time-histories of the exact (green lines) and analytical (black lines) amplitudes X, Y, a are shown for given X = 0, Y=0.0025 initial conditions, represented in Fig. 13a. The rough approximation of the asymptotic solution is mainly due to the fact it has been limited to the first-order. It is guessed that pushing forward the perturbation expansion one step further would give better results. With reference to item (ii), the exact and asymptotic (Eq. (43)) slow frequency * is plotted in Fig. 15 versus the wind parameters.

(a) (b) (c)
It appears that, although an ordering violation occurs in the perturbation scheme, which leads the frequency * to vanish, the asymptotic solutions in Fig. 15b, c, capture with reasonably accuracy the zero-frequency occurrence, denoting the collision of torus and limitcycle. Of course, negative frequencies, beyond the collision, have no meaning. In conclusion, even a loworder perturbation solution, as expected, gives useful information.

Conclusions
The aeroelastic stability of a slender prismatic viscoelastic structure under turbulent wind flow has been analyzed by a one-D.O.F. dynamical system, drawn by a Galerkin projection of the continuous model. Aerodynamic forces have been evaluated via the quasi-static theory. The unsteady wind has been split in a steady part and in a fluctuating part, assumed harmonic, of frequency close to the double of the fundamental frequency of the structure. The steady part is responsible for galloping; the harmonic part for principal parametric excitation. The two bifurcation mechanisms interact in the nonlinear field. The multiple scale method has been applied to find bifurcation equations. These have been solved for analyzing linear stability, limit cycles (periodic motions) and tori (quasi-periodic motions), bifurcating from the trivial path. The analytical solu-tions have been validated against numerical integration of the finite-difference discretized equations. The following main results have been found.
1. The parameter space is made of two regions: In each of them, periodic and quasi-periodic solutions bifurcate from the equilibrium. In the periodic region, the parametric excitation mechanism prevails over galloping. Here, however, the steady wind influences the parametric excitation instability domain, since it modifies the resultant (structural plus aerodynamic) damping. The higher the wind velocity, the wider the instability domain.
On the other side, turbulence reduces the galloping velocity of the steady wind case. In the quasiperiodic region, the galloping mechanism prevails over parametric excitation. However, the turbulent component changes the Hopf bifurcation of the autonomous system into a Neimark-Sacker bifurcation for the non-autonomous system, but it does not affect the critical velocity. 2. In the perfect resonant case, there are several limit cycle branches. They are multi-valued, and exhibit the hard-loss of stability phenomenon. Just one branch is stable. When a detuning is introduced, either it weakly affects the resonant solutions (small detuning), or it dramatically changes them (sufficiently large detuning), explaining as unstable bifurcated paths merge and move away from the trivial path. The agreement between exact and asymptotic periodic solutions is excellent. 3. Torus branches bifurcate from the trivial path. They also manifest the hard loss of stability phenomenon. After having regained stability, they collide with a limit cycle and disappear.

Appendix A: Illustrative example
A narrow pipeline, box-girder suspension bridge, is considered as an illustrative example of the theory developed. Such kind of industrial bridges is commonly adopted in the world to overcome ground depressions (or rivers) in constructing lines for liquid and gas transportation. Narrowness of bridge assures validity of the the aerodynamic quasi-static theory, since the time a fluid particle takes to cross the bridge is small with respect the natural period of the structure. The boxgirder hypothesis also entails that, for sufficiently long bridges, the torsional natural frequency is much higher than the flexural frequency, this circumstance preventing flutter (for which a coalescence between the two frequencies is requested), and allowing considerations of transverse vibrations only (galloping).

Background
A standard, single-span, pipeline suspension bridge is considered, consisting of a pair of identical parallel cables, a stiffening box-girder, equispaced hangers and two supporting towers (Fig. 16). Here, the classical continuum model of suspension bridge is adopted, originally derived by Melan in 1906 [53], revised by Bleich et al. in 1950 [54] for in-plane linear vibrations analysis, and recently re-examined by Luco and Turmo in 2010 [55]. Moreover, the Irvine's theory [56] for torsional vibrations of box-girder suspension bridges is used, aimed to investigate possible dynamic coupling, triggered by aeroelastic forces, between flexural and torsional behavior. According to these works, the stiffening girder is modeled as an Euler-Bernoulli beam; the vertical hangers are assumed to be inextensible and massless and are uniformly distributed along the span (i.e., they work as an inextensible and shear-flexible curtain); consequently, cable and beam undergo the same vertical displacement field v (s, t); the horizontal displacements of the cable are free, not being restrained by the curtain; the towers are assumed of negligible flexibility, so that the cables are fixed at the ends; out-ofplane displacements are ignored, but twist is accounted for; any structural (geometrical and material) nonlinearities are ignored. According to the Irvine's theory of shallow cables [57], the equations of motion of a single cable read:s Here: a dash denotes spatial-differentiation, a dot timedifferentiation and an overbar a quantity related to a single cable of the pair;T 0 is the pretension, k := w T 0 the curvature under self-weightw (including half the weight of the deck),T d the dynamic tension (all assumed constant in space) andm c is the cable mass The equations of motion of the box girder, ignoring non-uniform torsion effects, are: where: E b I b is the bending stiffness, G b J b the Bredt torsional stiffness, m b the mass per unit length of the beam, I Gb the centroidal mass inertia moment of the box, p b the external vertical force, c b the external torsional couple per unit length. By eliminating the reactive forcer between the first of numbered Eqs. (46) and (47), and appending boundary conditions, the following equations for the planar motion are obtained: where m := m b +2m c , p := p b +2p c are the total mass and forces, respectively, and where unbarred quantities refer to the pair of cables (e.g., T 0 = 2T 0 and E c A c = 2E cĀc ). By proceeding in the same way for the second of numbered Eqs. (46) and (47), and enforc-ing compatibility between the cross section and the cables (i.e., by letting v = − b 2 θ at the left cable), the equations ruling the twist motion are found: where G J := G b J b + b 2 4 T 0 and I G := I Gb + b 2 4 m c are the total torsional stiffness and inertia moment of the bridge, respectively.