On the nonlinear effects of the mean wind force on the galloping onset in shallow cables

The critical aeroelastic behavior of horizontal, suspended, and shallow cables is analyzed via a continuous model accounting for both external and internal damping. Quasi-steady aerodynamic forces are considered, including their stationary contribution (mean wind force). This latter induces a rotation of the cable (steady swing) around the line connecting the suspension points, together with a deformation of the initial equilibrium profile under self-weight. First, by using perturbation methods, the nontrivial equilibrium configuration is analytically determined as a nonlinear function of the wind velocity. Then, the wind critical values at which bifurcations take place and the corresponding modal shapes are determined by solving a boundary value problem in the complex field. Numerical investigations are carried out to validate the perturbation solution. A preliminary nonlinear galloping analysis is also performed to verify the galloping onset in terms of non-trivial equilibrium path and critical modes. The nonlinear terms related to the fundamental path, from which bifurcations take place, play a key role revealing new insights. They are able to heavily influence the system bifurcation, making unstable configurations which instead would be aerodynamically stable without considering the effect of the mean wind force.


Introduction
Cables are widely used structures in civil and industrial applications. They assume a fundamental role in the carrying capacity of suspended and stayed bridges, as well as in the realization of cable-cars, electrical transmission lines and many other structures. A general characterization of statics and dynamics of cables can be found in the classic Irvine's book [1], where several applications are analyzed as well. Recent advances in the statics of cables under forces exclusively vertical but not necessary uniform are proposed in [2], then extended to the case of general 3D forces in [3], where a perturbation asymptotic expression for the elastic problem solution is proposed. At the same time, the dynamic characterization of suspension cable systems is very important (e.g., [4]). Concerning nonlinear vibrations of suspended cables, the quadratic terms play a fundamental role in presence of internal resonances (e.g., [5,6]), differently from taut strings for which only cubic terms are present (e.g., [7]). In the case of inclined cables, additional properties that cannot be obtained by using horizontal cable modeling can be observed (e.g., [8]).
The wind-induced vibration of cables is a classic topic in structural applications [9]. Various aeroelastic phenomena may arise including vortex-induced vibrations (e.g., [10,11]), rain-wind-induced vibrations (e.g., [12,13]), dry galloping (e.g., [14,15]), ice galloping (e.g., [16,17]), and wake galloping (e.g., [18,19]). When operating in cold regions, an important role is played by ice galloping. It occurs when there is an ice accretion on the cable surface that breaks the symmetry of the cross section. As a consequence, the lift force may induce aerodynamic instabilities even for moderate wind. This phenomenon is first described in [20] and is largely analyzed in the literature, in the framework of the quasi-steady theory [21][22][23]. Without wanting to be exhaustive, a few examples of iced cable galloping are given in [16,[24][25][26]. Effects of possible internal resonances, in discrete or continuous nonlinear models, are investigated by the authors in [17,27], whereas a stiff modeling, able to take into account in a consistent way both bending and torsional stiffness, are introduced in [28][29][30]. Concerning flexible cables, the cable rotation, and resulting variation of their wind attitude are however very important for changes in aerodynamic properties of the cross section (e.g., [27,31]). Finite-element solutions of the problem can be found, for example, in [32]. The role of the internal and external damping contributions is evaluated in [33], where in-plane galloping analysis of flexible cables is performed, pursuing a direct approach on the nonlinear pde's of motion.
Focusing the attention on galloping critical conditions, the 1:1 internal resonance effects are pointed out in [16,34,35] as concerns two degree-of-freedom models. A very recent paper [36] presents a perturbation galloping stability criterion for a three-degree-of-freedom coupled motion of an iced conductor. A study of multimodal galloping on an iced transmission line is carried out in [37] using a reduced model including the first four in-plane modes and the first torsional mode, which is obtained by using Galerkin spatial discretization. Critical conditions and nonlinear responses are determined in a purely numerical way; the actual role of the rotation caused by the static aerodynamic forces is not completely unveiled since it seems related only to the torsional mode considered.
In this paper, considering a perfectly flexible (i.e., with evanescent flexural stiffness) and torsionally rigid cable arranged horizontally, the influence of mean wind actions on ice galloping stability is carefully addressed in the context of a continuous approach. For this purpose, the structural model of cable is taken from [33], and here extended to possible occurrence of out-ofplane motions. Due to the acting forces, a configuration change of the static equilibrium of the cable occurs and bifurcations take place starting from a nonlinear, nontrivial fundamental path. After modeling the aeroelastic flexible cable (Sect. 2), the aim is, therefore, twofold: (1) to find the nontrivial path as a function of the mean wind velocity (Sect. 3), and, along this path, (2) to search for the critical wind velocities (Sect. 4). Both of these objectives will be pursued analytically with the use of appropriate perturbation techniques. Comparisons with a numerical solution are however performed (Sect. 5) in both linear and nonlinear field in order to verify the accuracy of the results obtained. Some final considerations are reported in the ending Sect. 6.

Equations of motion
A horizontal flexible cable, having small sag-to-span ratio and under the effect of gravity and uniform wind, is analyzed here. As long as the sole selfweight is considered, i.e., no wind blows, the cable hangs in the ā x ,ā y -plane (vertical), in the configuration referred to asC (see Fig. 1, solid thin line). A local triad is defined in such a configuration, constituted by the tangent, normal and binormal unit vectors, and denoted by (ā t ,ā n ,ā b ). When a uniform wind of velocity U = Uā z is considered as well, it induces a time-dependent displacement of the generic point of the center line of the cable from the configurationC, which is indicated as u = uā t + vā n + wā b (solid thick line in Fig. 1); the specific effect of the wind is described in detail in Sect.

2.3.
The equations of motion for the cable, which are well-known in the literature [1,5,27,38], are written here in terms of normal and binormal displacement only (v, w), i.e., with the tangent displacement (u) condensed, under the usual hypothesis that the celerity of the longitudinal waves is much greater than that asso- Fig. 1 Horizontal cable under self-weight and binormal wind flow ciated with the transversal ones. Moreover, bending and torsional stiffness is ignored, consistently with the model of flexible cables, and differently than what it is done in case of stiff cables [28][29][30]. When adapted to the problem at hand, the equations of motion read: Here, T 0 is the prestress, assumed uniform along the abscissa s ∈ [0, l]; l is the length of the cable, taken nearly equal to the chord (i.e., the distance between the supports); E A is the axial stiffness; m is the mass per unit length (including possible ice coating);k := mg T 0 is the prestress curvature, also assumed uniform; e(t) is the dynamic unit extension, uniform on s, from which the dynamic tension is evaluated asT = E Ae; f d n , f d b are damping forces and f a n , f a b aerodynamic forces, all per unit length, acting in the normal and binormal direction, respectively; the prime stands for s-derivative and the dot for time derivative.
Models for damping and aerodynamic forces are discussed below.

Internal and external damping forces
The damping scheme follows what it is proposed in [33], in the spirit of the Kelvin-Voigt rheological model (note that, in [33], only in-plane motions are involved while here the procedure is extended to possible outof-plane motions as well). In particular, both external and internal damping contributions are considered. The latter is essential to avoid, in a continuum approach, infinitely many coincident critical wind velocities, one for each mode and, therefore, to split the coalescence of bifurcation points.
Specifically, the external damping is assumed as proportional to velocities: f d e n = −c ev , f d e b = −c eẇ , with c e a damping coefficient. This leads to an external damping operator proportional to the mass operator. Conversely, the internal damping, which accounts for various dissipative phenomena occurring in the material, here is assumed as proportional to the stiffness operator, to make the analysis as simple as possible (see [33] for further comments on the internal damping contribution). Therefore, since from Eq. (1), the geometric and elastic linear force components are: then the damping force components become: with ζ an internal damping coefficient. By superimposing the internal and external mechanisms, the Rayleigh model of damping is obtained: The equations of motion (1) accordingly read:

Aerodynamic forces
Aerodynamic forces are formulated in the framework of the quasi-steady theory (e.g., [21,22]). They are influenced by the occurrence, at any point, of an ice coating which breaks the axis-symmetry of the circular cross section of the cable. The forces are assumed uniform along the centerline of the cable, due to its small sag-to-span ratio (ā n ā y ) and to the uniform shape of the ice coating.
As previously said, the gravity forceb = −mgā y (g is the gravity acceleration), which includes the own and the ice weights, lets the cable lie in the vertical plane (configurationC, light gray plane in Fig. 2a). The initial orientation of the iced cross section is described by the directionā s , which is inclined by the angle ϕ 0 to the binormal axisā b . In other words, ϕ 0 represents the initial orientation of the cross section when the cable hangs on the vertical plane, in the absence of wind (Fig. 2b).
The wind, which blows normally to the cable plane, produces aerodynamic forces which can be decomposed in steady plus unsteady parts. The steady part, i.e., the one not depending on the velocity of the body, can be combined to the gravity giving rise to a uniform and constant force fieldb, that results parallel to an inclined plane and modifies the equilibrium configuration. More specifically, a static swingθ(U ) of the cable is induced, i.e., the cable swings to this new plane, which still passes through the line connecting the supports but is inclined of the angleθ(U ) with respect to the vertical one. This new equilibrium configuration is referred to asC (dark gray plane in Fig. 2a). The rotated plane will be called swung plane in the following. As specified in Sect. 2.1, and contrary to what it is done in [30], the vertical planar configurationC is taken here as reference, and displacements are measured from it.
Under the aforementioned hypotheses, the aerodynamic force acting on the cross section, which is the sum of a drag (f d ) and a lift (f l ) component, is (see Fig. 2b): where ρ is the air density, b a characteristic radius of the cross section, × is the cross product, U r = U r a u r the relative velocity between the flow and the cable, defined as In Eq. (7) the contribution of the spinθ is ignored, as typically done for compact shapes, it being small as compared tou. As a consequence, it turns out that: In Eq. (6), C d (γ + ϕ 0 ) and C l (γ + ϕ 0 ) are the drag and lift coefficients, neglecting their intrinsic uncertainty [39]; note that, if the angle of attack γ = 0, i.e., if the cross section is in its fixed initial orientation, the face exposed to the wind would be described by the angle ϕ 0 . With reference to Fig. 2b, the expression for γ is evaluated from: where the dot indicates the scalar product, so that, from Eqs. (7) it turns out: From Eq. (10) (and from Fig. 2b), it is evident the contribution of the swing angleθ which, summed to the initial orientation of the cross section ϕ 0 , changes the face of the cross section to the wind. When the previous results are expanded for small velocity ratiosv/U,ẇ/U up to the first order (i.e., considering a purely linear contribution from aeroelastic actions), it is found: Consistently with the mechanical model of flexible cable, twist and torsional moment are not defined. From Eqs. (6) and (11)-(13), after series expansion with respect to the variablesθ,v/U,ẇ/U , the normal and binormal components of the vector of aerodynamic forces, referred to as f a n and f a b , respectively, become: where the terms uniquely depending onθ are collected inf α , while those which depend on the velocity componentsv,ẇ as well, are collected inf α , α = n, b.
The former are the steady components, depending on the swing angle (and, parametrically, on the wind velocity); the latter are the non-steady contributions, depending on structural velocities (and, parametrically, on the swing angle and wind velocity). In particular, the steady forces are: and the non-steady ones are: The coefficients in Eqs. (15) and (16) are given in Appendix A. The expression (16) of the aerodynamic forces, as a function of the swing angle, is asymptotically correct tillθ is small of the same order of the (nondimensionalized) velocitiesv/U,ẇ/U . Series expansions (15) and (16) are expected to be sufficiently accurate as long asθ is small. However, whenθ is large, a different approach must be followed, namely: (a) the exact (not expanded) expression of the steady forces are evaluated, with Eqs. (7)-(10), from Eq. (6), after neglecting all the structural velocity terms; (b) the non-steady forces are determined by a series expansion carried out on the structural velocities only, and not on the swing angle. In this case, the steady forces, which substitute Eqs. (15), become: They have the important drawback of not being explicit inθ and U , thus calling for a purely numeric approach. The unsteady forces, which substitute Eqs. (16), are:

The equilibrium equations
The nontrivial equilibrium configuration in the swung plane is described by the steady response, indicated as v(s),w(s),ȇ. This is the solution of the nonlinear equilibrium equations, which are obtained from Eqs. (1) when inertia, damping and non-steady components of the aerodynamic forces are neglected, namely: In order to evaluate the equilibrium path, that is the steady solutionv(s),w(s),ȇ as a function of the wind velocity U , Eqs. (19) are solved. Due to the dependence of the steady aerodynamic forces onθ, a further equation which states that the cable lies on the swung plane at equilibrium, i.e., relatingv(s) andw(s) toθ , should be considered.
However, as more convenient alternative, the swing angle can be first evaluated through a physical consideration on forces, and then, once the steady aerodynamic forces are known, Eqs. (19) are taken on, in order to obtain the displacement and strain. In particular, at the equilibrium, the total steady force, uniformly acting on the cable, isb = (f n − mg)ā y +f bāz (see Fig. 3), and this induces the cable to assume the planar parabolic configurationC in the swung plane, whose trigonometric tangent is the ratio between the two force components, namely: Equation (20), which is transcendent inθ , once solved in terms of the parameter U , allows one to evaluate the steady aerodynamic forces. Since it must bef n −mg < 0 as the cable remains in tension, it isθ < 0 in the velocity field of interest. Note that, in the swung plane, the intensity of the normal load per unit length changes from mg to The difference of load entails elastic strain in the cable, which modifies its sag (the profile still remaining parabolic) and its stress. Consequently, the dynamic characteristics of the swung cable modify, with respect to those originally owned in the vertical plane.
As a further comment on Eq. (20), it turns out that the value of the swing angle is independent of the initial curvature of the cable.

The linear incremental equations
The general dynamic response to Eqs. (1) can be expressed in incremental form from equilibrium, by splitting the solution as: whereṽ(s, t),w(s, t),ẽ (t) are incremental variables. By substituting Eqs. (21) in Eqs. (5), accounting for Eqs. (19) and linearizing in the increments, the linear incremental equations of motion are obtained (tilde omitted on the incremental variables): where the steady tensile force of the cable isT = E Aȇ. Note that, in writing Eq. (22-b), an integration by parts is performed, i.e., where the boundary terms go to zero. Equations (22) govern the small oscillations of the cable around the nontrivial equilibrium configuration, with the incremental displacement components referred to the original basis located in the vertical plane and the aerodynamic forces defined in Eq. (16).
Two main approaches can be possibly followed to get to the final goal, which is the evaluation of the bifurcation conditions of Eqs. (22). The first is valid for largeθ and it is purely numerical using the following steps: 1) pick of a wind velocity U ; 2) evaluation of the swing angle solving Eq. (20), by means of the Newton-Raphson method, using the exact and implicit definition of the forces (17); 3) evaluation of the equilibrium solution solving Eqs. (19) by means of the finite difference method; 4) solution of the infinite dimensional eigenvalue problem (22), with definition (18) for the forces and still using finite difference method. This approach requires a parametric sweep in terms of the wind velocity U , i.e., the procedure is entirely repeated for different values of U until critical conditions are found. The whole procedure is implemented in MATLAB c [40], and requires a very short computational time (a few seconds) for each wind velocity.
The second approach is analytical and makes use of perturbation methods; it is valid for moderately small θ and entails: 1) perturbation expansion and chainsolving of Eq. (20) to obtain the swing angle as a function of the wind velocity, using the approximated expression of forces in Eqs. (15); 2) perturbation solving of Eqs. (19) to obtain the equilibrium solution; 3) analytical seeking of the critical solution of the linear problem (22), with definition (16) for the aerodynamic forces.
In the following Sects. 3 and 4, the second, analytical procedure is described in detail, and finally, in Sect. 5, outcomes are compared to those given by the first, numeric, procedure.

Evaluation of the swing angle
The equilibrium equation (20), together with the definitions (15) of the steady forces, implicitly determines θ as function of U . It reads: Since the forces are expressed in polynomial form (Eqs. (15)), a perturbation method is conveniently applied. By rescaling the forces at order (i.e., by lettingf α → f α , so thatθ → 0 when → 0), expanding tanθ =θ +θ 3 3 + . . . and using the series expansion θ = θ 1 + 2θ 2 + 3θ 3 , the substitution in Eq. (24) and vanishing of the coefficients of different powers of produces the following perturbation equations: Chain solution leads to: from whichθ(U ) is reconstituted as: With this result, the steady forces can be expressed in terms of the wind velocity only. By substituting Eq. (27) in Eq. (15), one finally gets: where the positions (26) hold and h.o.t. stands for higher order terms.

Evaluation of the static displacements
To find the steady response of the cable, accounting for elasticity, the equilibrium equations (19) must be solved. The following trial solution is used: where y(s) =k 2 (s − l)s is the parabolic initial profile, withk = mg T 0 , and α n , α b are (small) nondimensional amplitudes, to be determined. Since y A = y B = 0, the boundary conditions in Eqs. (19) are satisfied.
By substituting Eq. (29) into Eqs. (19-a,b), the following nonlinear algebraic system in α n , α b ,ȇ is obtained: After eliminatingȇ, the system reduces to: where It is worth noticing that, as a check for the validity of the trial solution (29), since bothv(s),w(s) are proportional to y(s), this solution keeps the cable planar. Indeed, denoting by x(s) := (x(s) +ȗ(s))ā x + (y(s) + v(s))ā y +w(s)ā z the position of a generic point in the swung configuration, the vector n :=ā x × x(s) = y(s)[(1−α n )ā z +α bāy ] keeps its direction constant for any s, coincident with the normal to the swung plane. The rotation of the plane is the angle formed by n and a z , i.e.: Moreover, by combining Eqs. (31), it follows that , so that Eqs. (32) and (20) are consistent. Therefore, finding α b , α n and using Eq. (32) gives an alternative method (to Eq. (24)) to evaluateθ .
The solution of Eqs. (31) is now sought with a perturbation method. By rescaling the forces at orderand expanding the unknowns as α h = α h1 + 2 α h2 + 3 α h3 , substitution in Eqs. (31) and vanishing the coefficients of equal powers of gives the following perturbation equations: After chain-solving the linear systems (33), and reconstituting, the following expressions are obtained: Substitution of this results in Eq. (29) provides the static deflection of the cable, including rigid and elastic effects.
Concerning the increment of strainȇ, this can be evaluated by Eq. (30-c); the associated increment of tensionT = E Aȇ is: It is worth noticing from Eqs. (33)-(35) that, when Λ 2 → ∞ (inextensible cable) the linear part of the normal displacement ampliture α n1 disappears, as it happens in an inextensible pendulum performing small swing; furthermore the (nonlinear) elastic strain tends to zero, in a perturbation sense, when Λ 2 → ∞ (i.e., it contains terms of order 4 if the perturbation procedure is stopped at order 3 ), entailing no profile changes, whileT T 0 →f n mg , i.e., the tension ratio equates the load ratio; therefore a static effect only occurs. On the other hand, whenk → 0 (taut string), Λ 2 → 0 too, so that a deflection takes place (α n = 0), whileT → 0, i.e., kinematic effect only occurs.

The critical wind velocity
To compute the critical wind velocity and the corresponding critical mode, the incremental equations of motion (22) must be considered. By using Eq. (16) for the expression of the aerodynamic forces and Eqs. (29) for the equilibrium solution, they become: It is worth noticing that, in addition to the effects discussed in the previous section, the change in the equilibrium configuration (swing) induces a change also in the aerodynamic damping matrix, thus accounting for the modified exposure to wind of the cross section. Furthermore, it is important to point out that the analytical solving of the problem (36) is possible only as consequence of the evaluation of the non-trivial equilibrium path using the proposed perturbation method, as carried out in the previous Section. To solve the problem (36), variable separation is used, namely the substitution (v(s, t), w(s, t), e(t)) = (v(s),ŵ(s),ê) exp(λt), which produces: where: The eigenvalue problem (37)-(39) must be solved in the complex set; indeed, the galloping (critical) modes, which are the eigenfunctions associated to the eigenvalues with zero real part, are generally complex, due to the non-proportional nature of the damping operator (through its aerodynamic counterpart). Solution to the field equation assumes the form: where β 1 , β 2 , φ 1 , φ 2 are the eigenvalues and eigenvectors of the following 2 × 2 algebraic problem: Moreoverû is a particular solution to the non-homogeneous problem (consideringê as a known term), and C 1 , . . . , C 4 are arbitrary constants. By substituting the components ofû in the equation forê (38), this latter is drawn in terms of the arbitrary constants. Finally, from the boundary conditions (39), a 4 × 4 homogeneous algebraic problem follows. Explicit expressions of the solution are cumbersome, but help of an algebraic manipulator makes the problem easy to be solved. In the critical condition, it is λ = iω c , U = U c ; by splitting the characteristic equations in the real and imaginary parts, two transcendent real equations for ω c , U c , which are the critical frequency and wind velocity, respectively, are derived and solved by the Newton-Raphson method.

Description of the case studies
Geometrical and mechanical parameters of the cable assumed as case study are given here. The unstretched length is l = 267 m, the mass per unit length is m = 4.4 kg/m, and the axial stiffness is E A = 4.0×10 8 N. Two different values of initial sags d are considered, the first one giving rise to a very taut configuration and indicated as case S (smaller sag), the second one to a moderately taut configuration, referred to as case L (larger sag). The values for case S are: d = 3 m, which corresponds to an initial curvatureκ = 3.37 × 10 −4 m −1 , initial tensile force T 0 = 128.213 kN and the Irvine parameter /(2π) = 0.79, i.e., at the left of the first crossover point; the values for case L are: d = 26 m, i.e., initial curvatureκ = 2.92 × 10 −3 m −1 , initial tensile force T 0 = 14.734 kN and the Irvine parameter /(2π) = 20.39, i.e., far to the right of the first cross-over point. In the two cases, different values of external and internal damping coefficients are assumed, respectively, so as to have equal modal structural damping ratio for the first two natural modes (without wind), in the amount of 0.5%: they are c e = 0.056 kg×s/m and ζ = 239.4 kg×m×s for case S, c e = 0.020 kg×s/m and ζ = 73.0 kg×m×s for case L. The aerodynamic parameters are: air mass per unit volume ρ = 1.25 kg/m 3 , radius b = 0.102 m; the drag C d (ϕ) and lift C l (ϕ) coefficients refer to the NDT cross section in [16,41], and they are shown in Fig. 4a. Moreover, for both the cases S and L, the analyses are carried out for two different values of initial angle of the ice coating, which correspond to configurations not prone to galloping in the Den Hartog meaning [20], i.e., where C d (ϕ) + C l (ϕ) > 0: they are (1) ϕ 0 = π/8 rad = 22.5 • and (2) ϕ 0 = 11π/48 rad = 41.25 • (red dashed lines in Fig. 4b). Summarizing, a total number of four cases are analyzed, namely S.1, S.2, L.1, and L.2, where 1 or 2 indicates the initial angle of the ice coating. In all the cases, the outcomes of the perturbation procedures are compared to those of the numeric approach, as performed through a finite difference code implemented in MATLAB c [40].

Equilibrium configuration
As a first result, the swing angleθ evaluated from the analytical solution, Eq. (27), is shown in Fig. 5, as a function of the wind velocity U . Note thatθ is independent of the sag-to-span ratio, therefore Fig. 5a holds for both cases S.1 and L.1, while Fig. 5b for both cases S.2 and L.2. It turns out that the perturbation solution (in blue solid line) is in good agreement with the numeric one (in red dashed line) up to wind velocities which do not overtake about U = 8 m/s in Fig. 5a and U = 10 m/s in Fig. 5b, corresponding to angles of magnitude 0.1 rad and 0.2 rad, respectively. These angles represent the upper limit where the cubic series expansion is valid; moreover, different limit values are obtained for the two initial angles as a consequence of the different magnitudes of the derivatives of the aerodynamic coefficients, which lead to a better correspondence to the assumed ordering for the case (2). As a further comment, the swing angle is negative as expected for the convention assumed.
The amplitudes of the parabolic shape at equilibrium, α n and α b , as given by Eq. (34), are shown in Figs. 6 and 7, as function of the wind velocity, for the four analyzed cases. As reasonable, the binormal component is always larger than the normal one. Moreover, among the four cases, the one in Fig. 6b, i.e., S.2 appears to have larger amplitudes in both the components. The perturbation outcomes are in good agreement with those relevant to the numeric procedure, even if some discrepancies are obtained for the α n component in case S.2 for wind velocities larger than U = 8 m/s, due to the loss of accuracy in the evaluation of the swing angle.
It is worth noticing how, in cases S.1 and S.2 (Fig. 6), the binormal displacements are much smaller than in cases L.1 and L.2 (Fig. 7), as a consequence of the different sag-to-span ratios. Anyway, no significant differences occur keeping fixed the sag-to-span ratio and changing the initial angle of the ice coating, i.e., between case S.1 and S.2, or L.1 and L.2.
It is interesting to evaluate the total tensile force of the cable at equilibrium in the swung plane, namely T 0 +T , as a function of the wind velocity. This is shown in Fig. 8a for cases S.1 and L.1, and in Fig. 8b for cases S.2 and L.2. It can be observed that, in most of the range of the considered velocities (U ≤ 8 m/s), the agreement between analytical and numeric results is very good, as well as modifications in the tensile force are of a few percent. On the other hand, bad agreement is found for larger velocities, particularly for cases S.2 and L.2 (Fig. 8b). Actually, the perturbation solution gives a totally wrong outcome in case L.2, and this is due to the evaluation ofT as a term proportional to the difference of α n and 1 2 (α 2 n + α 2 b ) (see Eq. (35)), which both assume similar values, very close to 0.025 at U = 10 m/s (α b about 0.25), in violation of the perturbation order. In these cases, further perturbation orders should be considered, so as to give consistency to the analytical solution even for larger wind velocities.

Critical conditions
The critical modes and relevant critical conditions in terms of wind velocity and modal frequency are now analyzed, as solutions of the eigenvalue problem (22). In particular, in Table 1, the critical conditions are reported for both analytical and numeric solutions, showing good agreement, specifically in the modal frequency where the percentage error is practically zero.
For case S.1, the corresponding modal shape is shown in Fig. 9. It is anti-symmetric and real, involv-  ing both normal (v) and binormal (ŵ) components with the same order of amplitude, in good agreement with numeric outcomes. Moreover, in Fig. 10 the evolution of all the eigenvalues which become unstable are shown in the complex plane, as the parameter U is increased and assumes the values in Table 2; concurrently, the corresponding modal shapes are shown at the right side of the same Figure. From this Figure, which is realized with data coming from the numeric procedure, it turns out that the eigenvalues which become unstable are four; furthermore, a regain of stability is achieved as the wind velocity is increased, due to the progressive change of orientation of the cross section with the increasing of the swing angle and, as a consequence, to the modification of the intensity of the aerodynamic forces. It is worth noticing that the order in which the eigenvalues regain stability is different than that of their loss of stability. Moreover, the regain of stability is not (a) (b) Fig. 9 Critical mode for case S.1: a real and imaginary parts ofv; b real and imaginary parts ofŵ necessarily associated to a comeback of the evolution of the system on the trivial dynamic solution. The system might evolve on other possible, stable and coexisting attractors, deriving from secondary bifurcations here not investigated. For case S.2, the critical mode is shown in Fig. 11: it turns out to be symmetric and complex, with only one semi-wave along the span, coherently with a configuration slightly at the left of the first crossover point. In this case, the agreement with numeric outcomes is good, in particular for the normal component (Fig. 11a). The eigenvalues which lose stability are shown in the complex plane in Fig. 12, at the wind velocities reported in Table 3; relevant eigenfunctions are shown as well at the right part of the figure. In this case, the eigenvalues which lose stability are six, and the order in which they regain stability is the opposite of that in which they lose it.
For case L.1, in Fig. 13, the critical mode is again anti-symmetric, while it turns out to be symmetric for    Fig. 14, having a three semi-waves shape, which is coherent with the configuration at the right of the cross-over. In the latter case, as for case S.2, still the normal component has better agreement with the numeric outcomes than the binormal component, which presents some quantitative differences, even if the qualitative agreement is always very satisfactory. As a summary, for both cases S.1 and L.1, i.e., for the first value of ice-coating initial angle, and independently of the sag-to-span ratio, it turns out that the critical mode is anti-symmetric and real; on the other hand, for the second value of ice-coating initial angle, the critical mode is symmetric and complex. Therefore, it appears that the modal shape of the cable at bifurcation is strictly related to its initial attitude to wind and, then, to the value of the aerodynamic forces, which have the ability to significantly change the structure of the Hamiltonian dynamical system and, consequently, the nature of the eigenfunctions, triggering symmetric or anti-symmetric modal shapes. Consistently, critical modes might be strongly different from natural ones, in absence of wind.
In this respect, it is interesting to look at the value of the steady angle under which the ice-coating faces the wind at bifurcation, i.e., −θ + ϕ 0 , checking the corresponding value of the Den Hartog's coefficient C d +C l , valid for 1-d.o.f. vertical galloping (see Fig. 15): it is evident how the effect ofθ (and, then, of the mean wind force) is to reduce the Den Hartog coefficient, so as to take the cross section closer to the 1-d.o.f. instability region, even if it can still be non-negative (e.g., for cases S.2 and L.2). Specifically, it is evident that, due to the generally concurrent presence of both the modal componentsv,ŵ, the mechanism of instabilization is more complicated than that which only relies on the Den Hartog's coefficient (see, e.g., [35] for 2-d.o.f. systems).

Nonlinear galloping
In order to have a first insight into the cable nonlinear behavior close to the first bifurcation point, a spatial finite-difference approximation of the nonlinear equations of motion (5) is carried out and a numerical integration of the resulting equations in time is performed (details on the numerical procedure are given in Appendix 2). In this way, it is also possible to fully validate, in a completely independent way, the results obtained with respect to the equilibrium configurations (Sect. 5.2) and to the first critical modes (Sect. 5.3). Results are limited to the case S for the sake of brevity; anyway, outcomes show that aerodynamic forces greatly influence the critical response despite the feature of the Hamiltonian system.
Concerning the case study S.1, Fig. 16 shows the time histories for steady oscillations at a mean wind velocity equal to 5 m/s, slightly higher than the first critical speed. Cable oscillations occur around the equilibrium configuration (marked by a red dashed line in the figure), with values fully consistent with those of Fig. 6a and Eq. (29), and are characterized by a single frequency corresponding to the first critical one (see Table 1). From the direct integration of the nonlinear motion equations, the cable dynamics during a period of oscillation, subtracted to the equilibrium configuration determined by the mean wind, is presented in Fig. 17: the shapes are entirely in accordance with the first (normalized) critical modes of the system (Fig. 9), pointing out a not negligible coupling between the two planes, already highlighted by the critical modes. It should be noted that the vibration shapes remain almost constant during the oscillation confirming the real nature of the corresponding eigenfunctions (the imaginary part of which is almost zero, see Fig. 9).
Moving on to the case study S.2, Fig. 18 shows the steady oscillations at a mean wind velocity equal to 6.4 m/s, immediately after the first bifurcation point, and just prior to the second one. Also in this case the average value of the oscillations is completely consistent with the values obtained for the equilibrium configuration ( Fig. 6b and Eq. (29)), with a clear prevalence of the horizontal component. The motion is monomodal periodic, with the presence of a single vibration frequency (i.e., the first critical frequency in Table 2). The corresponding cable deformations during a vibration period Table 3 Numerical values of the velocities for the positions in Fig. 12 Fig. 19, still subtracted to the equilibrium configuration associated to the mean wind, allowing one to highlight the dynamics of the motion. It occurs almost completely in the vertical plane with a modest coupling in the horizontal direction. Moreover, the vertical oscillations maintain almost the same shape over time, while the horizontal oscillations present remarkable variations in shape during a vibration period. This result is in excellent agreement with what was obtained in the critical condition analysis (Fig. 11 -normalized critical mode), in which the complex nature of the horizontal critical mode was highlighted, and just an order of magnitude lower than the vertical mode (which instead was almost real, as confirmed by this nonlinear galloping analysis).

Conclusions and final remarks
This paper deals with the research of critical conditions of ice galloping in suspended cables. The key point of the work is the nonlinear contribution of the non-trivial equilibrium path that influences the terms governing galloping onset conditions. It reveals new insights on the bifurcation phenomenon, being able to make unstable cable configurations that instead would be aerodynamically stable without the effects of the mean wind force.
Concerning the modeling aspects, the main novelty of the paper is the use of a continuous model, unlike what happens in current literature. In this way it is not necessary to select in advance a certain number of modal shapes (as usually happens in discrete Galerkin models, e.g., [37]). The proposed method is able to naturally follow the evolution of critical conditions as the wind speed increases. The problem is for now restricted to the case of horizontal cables to have an exact plane description of the non-trivial fundamental path, even if the proposed approach can be easily used in the presence of non-planarity of the static configuration of the cable.
A second strength of the work is that the search for critical conditions is proposed analytically through the use of perturbation techniques. In particular, the nonlinear, non-trivial fundamental path from which the bifurcation takes place is first determined in two steps: (a) determination of the swing angle that identifies the plane on which the cable lies due to the effect of the mean wind force (swung plane), and (b) evaluation of the corresponding static displacements of the cable, together with the associated increment of tension. An eigenvalue problem is then deduced starting from the continuous formulation of the problem, which leads to transcendental characteristic equations in analytical form.
Concerning numerical results, some points are worth noting: the vibration modes obtained are naturally complex due to the non-classical nature of aerodynamic damping; the agreement of the analytical solution of the problem with the numerical one is excellent on critical instability conditions (i.e., when the eigenvalues become unstable), while a greater error could be present for higher mean wind speeds due to the analytical calculation of the swing angle (truncated to the third order); thanks to the continuous formulation, the influence of steady swing is naturally present on all modes of the system (unlike what happens in the literature for discrete multimodal galloping, e.g., [37]); therefore, all critical eigenvalues return to the stability domain for sufficiently high mean wind velocities; the numerical integration of the nonlinear equations of motion (in both mechanical and aerodynamic terms) has allowed a further validation of the goodness of the results obtained from the critical conditions, at least limited to the first bifurcation point. Indeed, through a completely independent numerical procedure, full correspondences on both equilibrium path values and critical modes of the system are determined, even in cases where the oscillation modes assume significant complex representations. In fact, the steady swing strongly influences the aerodynamic forces, since it changes the exposition to wind of the body. Therefore, its consistent inclusion in the model is fundamental, and the nonlinear terms arising from it turn out to strongly affect the onset of galloping. A special comment deserves the fact that, in the presented examples, the unstable mode can be symmetric or anti-symmetric despite the lower vibration frequency of the cable always corresponds to a symmetric out-of-plane mode for any configuration (smaller or larger sag) considered. This can be attributable to the presence of the aerodynamic damping operator, which in turn depends on the non-trivial equilibrium path through the aerodynamic coefficients, capable of significantly modifying the critical behavior of the system. This result, which may significantly affect design choices for suspended cables in cold regions, is worthy of future investigations, also using particular perturbation techniques (e.g., [42]). In this context it is also worth noting that, as a consequence of the proposed approach, the motion must be described in terms of incremental variables, which are superimposed to the static configuration in order to find the critical conditions and carry out possible nonlinear bifurcation analyses. The incremental variables are here expressed in terms of the original extrinsic basis, in the vertical plane where the cable lies in absence of wind. Therefore, the locutions 'in-plane' and 'out-of-plane', relating to cable dynamic displacements, lose their meaning since they should be referred to the vertical, not to the currently rotated plane when the wind is blowing. Then, the cable natural modes will appear coupled in the two displacement components. But this approach, apparently more complicated than the intrinsic one (e.g., [29]), has the great advantage not to require any change of basis, nor projection of forces, while keeping the original meaning of the variables.
The first results of nonlinear galloping oscillations, limited to the first bifurcation, has to be extended to the analysis of the limit-cycle arising at successive Hopf bifurcations. However, this is a problem which appears far from being trivial, since a cluster of several modes (from four to six in the examined examples) could in principle interact in the postcritical field. Moreover, the regain of stability of the fundamental nontrivial path does not exclude the existence of supercritical competitive attractors, produced by the previous bifurcations.
Funding Open access funding provided by Universitá degli Studi di Genova within the CRUI-CARE Agreement.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.