Efficient modal-based method for analyzing nonlinear aerostatic stability of long-span bridges

,


Introduction
The increasing main span lengths of recently projected long-span cable-stayed bridges demand a comprehensive analysis of their windinduced responses.Generally, the wind-related ultimate limit state in the design of long-span bridges has been flutter.According to Ge [1], flutter essentially controlled the design of four of the five bridges with the longest span in the world, including the Akashi, Xihoumen, Yi Sunsen, and Runyang bridges.This aeroelastic phenomenon can be generally defined, according to Simiu and Scanlan [2], as a "self-excited oscillation that a structural system by means of its deflections and their time derivatives taps off energy from the wind flow".When "this energy of motion extracted from the flow exceeds the energy dissipated by the system through mechanical damping", the flutter instability occurs, leading to an eventual structural failure.Further details can be found elsewhere, for instance, in Scanlan [3] and Larsen and Larose [4].Hence, specifications of bridge projects (see Stretto di Messina [5]) demand the bridge critical flutter velocity to be higher than a given threshold, in accordance with the wind characteristics of the project location.However, as main spans more frequently surpass 1 km in length, aerostatic stability arises as an important design limitation that can be in some cases even more governing than the dynamic flutter instability, noticeably for cable-stayed bridges.This wind-induced phenomenon, also known as static nonoscillatory divergence, and more traditionally, torsional divergence, is a stability problem that happens when "a velocity is reached at which the magnitude of the wind-induced moment, together with the tendency for twist to demand additional structural reaction, creates an unstable condition and the structure twists to destruction", as described for the single-degree-offreedom torsional divergence in Simiu and Scanlan [2].Aerostatic stability problems were noted in wind tunnel tests and reported in several studies, such as Hirai et al. [6], Cheng [7], and Cheng et al. [8].In the study by Boonyapinyo et al. [9], an aerostatic critical wind velocity of 76.5 m/s was analytically obtained for the Akashi Bridge, while the elastic flutter wind velocity obtained through wind tunnel tests was 92 m/s (Boonyapinyo et al. [10], Honshu-Shikoku Bridge Authority [11]).Another example can be found in the paper by Nagai et al. [12], where an aerostatic critical wind velocity of 80 m/s was reported for a singlebox cable-stayed bridge with a main span of 1400 m, while its flutter critical wind velocity was 120 m/s.It must be also noted that, in general, it is advised for the aerostatic critical wind velocity to be above the flutter critical wind velocity (Su et al. [13], Wind-resistant design specification for highway bridges [14]).These facts highlight the importance of considering simultaneously both the aerostatic and the aeroelastic ultimate limit states in the design of long-span cable-stayed bridges.
Aerostatic instability is a three-degree-of-freedom wind-induced phenomenon caused by the action of the steady aerodynamic forces on the flexible deck of a bridge.The implicit relationship between the three components of the displacements-dependent wind loads and the deck displacements leads to the instability of the structure at a particular wind velocity.The critical wind velocity associated with aerostatic instability U cr represents the critical bridge condition that must be considered in a bridge design framework in the same manner as it is the case with the flutter analysis.Hence, design-oriented analysis methods must focus on the efficient estimation of U cr , while tracing the deformation path as a function of the wind velocity.
Although aerostatic stability has been studied since the early years of aerospace engineering (see e.g.Bisplinghoff and Ashiley [15], Dowell [16]), the first methods used to analyze the aerostatic stability in longspan bridges were developed around the 1990's.The critical wind velocity was obtained using the so-called linear method (ASCE [17], Simiu and Scanlan [2]), which adopted a linear elastic stiffness matrix and modeled aerodynamic loads by a single component using a linearized derivative of the pitching moment.Later, the methodology was advanced by considering the three components of the wind displacement-dependent loads and the nonlinearities of the structure (Cheng et al. [8], Boonyapinyo et al. [9], Nagai et al. [12], Boonyapinyo et al. [18], Cheng et al. [19], Zhang [20]).The development of improved aerostatic instability analysis techniques is an active research topic and several recent contributions can be found, for instance, in Zhang et al. [21], Su et al. [13], Arena and Lacarbonara [22], Ge and Shao [23], Dong and Cheng [24], and Hu et al. [25], among several others.While some studies have focused more on the phenomenological analysis of the bridge displacements as a function of the wind velocity and the subsequent collapse mechanisms (e.g.Zhang et al. [26]), other studies have proposed alternative strategies for the fast identification of the critical wind velocity.In Boonyapinyo et al. [9,18], an outer loop combining the eigenvalue analysis under wind loads and the windvelocity bound algorithm is used to find the critical wind velocity.In Cheng et al. [27], an incremental-two-iterative procedure was proposed to speed up the convergence of the outer loop and increase the level of detail in tracing the wind-induced displacements when reaching U cr .This procedure was later enhanced in Zhang et al. [28].
The aforementioned methods address the aerostatic stability problem from the perspective of the stiffness, which is modeled considering nonlinear stiffness matrices or by introducing nonlinear finite element structural models.By contrast, the analysis of other aeroelastic phenomena, such as flutter, buffeting, or vortex-induced vibration, are often carried out adopting the dynamic equilibrium equation, and the common inputs are the modal properties of the bridge.Aiming at developing comprehensive design frameworks for preliminary design stages that may consider all the aeroelastic phenomena involved in the problem, the reformulation of the aerostatic stability analysis in terms of dynamic equations would be very useful.
This study proposes an alternative method for the assessment of the aerostatic instability wind velocity using the modal properties of the bridge as input parameters.A first contribution using dynamic properties to reformulate the linear method was proposed by Xiang [29], where the standard linear method (Simiu and Scanlan [2]) was modified by introducing the first symmetric torsional frequency (see Cheng et al. [19]).This simplifies the procedure by using as input values the same ones required in the analysis of other aeroelastic phenomena, but it keeps the intrinsic limitations of the linear method.The present paper further develops this approach by introducing a modal-based formulation that considers the three components of the wind displacementdependent loads, and the nonlinear aerodynamic and structural features of the problem.
A simplified version is first developed to solve the aerostatic analysis problem including the three nonlinear aerodynamic components but assuming a linear structural behavior of the bridge.An instabilityfinding outer loop is proposed using root-finding algorithms given the independence of the analyses at each wind velocity.The required inputs are the force coefficients as a function of the angle of attack, and the bridge's natural frequencies and modal matrix normalized to the mass.This information can be conveniently obtained from any commercial FEM software, which makes very simple the implementation of the method in design environments.Then, this algorithm is recast to carry out the analysis considering the nonlinear features of the structure, including: nonlinear behavior of the stays (sag effect), beam-column effect, and large displacements effect.These nonlinearities are modeled through a nonlinear FEM, and the sag effect is considered adopting the Ernst equation (Ernst [30]).Both the linear and nonlinear versions of the algorithm are developed adopting the instability-finding outer loop to identify the critical wind velocity, and an alternative incremental outer loop to trace the bridge displacements as a function of the wind velocity.The performance of the proposed methods are tested in three application examples, where the influence of the nonlinear features of the problem is analyzed.

Aerostatic equilibrium problem definition
The analysis of the aerostatic response of long-span bridges is a nonlinear problem caused by the dependency of the three components of the mean wind-induced loads on the angle of attack α.The equilibrium equation can be expressed as where K is the nonlinear stiffness matrix of the bridge, u stands for the vector of displacements, which has three components for each node j, as , and f s is the vector of nodal aerodynamic steady forces acting on the bridge deck, which is a function of the vector of angles of attack for each node α, and each component is defined as follows: where ρ is the air density, U is the wind velocity, l e is the span-wise length of the deck element, and C D , C L , and C M , are the time-averaged drag force, lift force, and pitching moment coefficients, respectively, which are a function of the wind angle of attack α.This definition can be made assuming the steady wind load model hypothesis, where the wind is considered to be constant without fluctuations, and the force coefficients are obtained assuming that no displacements are considered in the deck cross-section.Furthermore, the consideration of time-averaged aerodynamic forces allows the assessment of the bridge aerostatic response without adopting time-domain approaches, which are computationally much more demanding.

Aerostatic nonlinearities
There is a number of nonlinear features present in the aerostatic equilibrium of a bridge, which are sketched in Fig. 1 and can be summarized as: 1. Wind-induced displacement-dependent aerostatic loads (Fig. 1 (a)).
The relationship of the force coefficients with the wind angle of attack α represents a source of nonlinearities for the aerostatic analysis.While streamlined geometries usually present a linear relationship of C L and C M with α before the stall, and the drag response usually shows a U-shaped relationship with α, bluffer deck cross-sections can present aerodynamic responses far from being linear (see, for instance, Diana et al. [31]).2. Aerostatic equilibrium problem (Fig. 1 (b), Section 2.1).The assessment of the rotation of the deck at a wind velocity is a given nonlinear problem that must be solved iteratively since both the aerodynamic forces and the wind-induced displacements depend on each other.3. Nonlinear structural responses (Fig. 1 (c)).Focusing now on the left side of Eq. ( 1), the structural characterization of the bridge may comprise a number of nonlinearities that can affect the aerostatic response of the structure (Wang and Yang [32], Gimsing and Goergakis [33]).These nonlinearities can be listed as: 3.1 Interdependency of lateral and torsional deck displacements.
Cable-supported bridges usually present a torsional draginduced response which affects the aerostatic response of the bridge when the three components of the wind displacementdependent loads are considered.An example can be found in Boonyapinyo et al. [9], where a nose-down drag-induced rotation is described for the Akashi Kaikyo Bridge.3.2 Nonlinear behavior of the stays: The sag effect (Ernst [30], Wang and Yang [32], Karoumi [34], Thai and Kim [35]) modifies the Young's modulus E of each stay depending on its stress level.As the wind velocity increases, lift and moment components may induce upward displacements of the deck that can reduce the stays' stress level, leading to a drastic reduction of the global stiffness of the bridge.This is known as stiffness degradation (Zhang et al. [26]), and can drastically drop the aerostatic critical wind velocity.The sag effect can be modeled in different ways, as discussed in Freire et al. [36].One of the most popular is the Ernst equation (Ernst [30]), which can be written as: where E eq is the equivalent cable modulus of elasticity, E is the effective modulus, A is the stay cross-section area, q stands for the stay weight per unit of length, L is the horizontal projection of the stay length, and T is its tension force.3.3 Beam-column effect: The characteristic large compressions in the deck and towers of cable-stayed bridges explains the importance of considering this nonlinear effect (Fleming [37], Wang and Yang [32], Freire et al. [36]) that can be modeled via nonlinear FEM.3.4 Large displacements effect: Bridge responses can be sensitive to this effect in high wind load scenarios, which may increase the lateral deck displacement and subsequently the drag-induced rotation.It can be also modeled through nonlinear FEM.3.5 Nonlinear material behavior: Material nonlinearity can affect to the critical wind velocity, as discussed in the study by Ren [38].
All these sources of nonlinearities play a role in the aerostatic analysis of long-span bridges since large differences can be found when linear and nonlinear results are compared.For instance, a critical wind velocity of 135 m/s was identified in Boonyapinyo et al. [18] for a 1000m center span cable-stayed bridge using a nonlinear approach considering the three components of the displacement-dependent wind loads and the geometric nonlinearity in the analytical model, while the results obtained for the same bridge applying the linearized torsional divergence analysis was 519 m/s.This difference highlights the need for introducing nonlinear approaches in aerostatic-resistant design frameworks.

Aerostatic instability
The aerostatic instability problem consists in identifying the critical wind velocity at which the deck rotation caused by the wind displacement-dependent loads yields a positive feedback loop leading to the bridge collapse.Some of the most relevant methods are summarized below.

Linear formulation for a 1DoF-1node dynamical system
The critical wind speed U cr of a 1DoF-1node dynamical system can be evaluated following the linear method proposed by Simiu and Scanlan [2] as: Fig. 1.Schematic relationships between the terms of Eq. ( 1) that give place to the nonlinearities of the aerostatic analysis problem.
M. Cid Montoya et al.
where K α is the torsional stiffness, ρ is the density of the air, B is the width of the deck cross-section, and C ′ M,0 • is the slope of the moment coefficient with respect to the wind angle of attack.This expression was rewritten from a dynamical perspective (Xiang [29], Cheng et al. [8,19]) as: where f t is the first symmetric torsional frequency and where I m is the mass moment of inertia per unit of length, m is the mass per unit of length, r is the inertia radius of the section (I m = mr 2 ).

Linear formulation for full bridges: flexibility matrix and eigenvalue problem
The procedure for computing the aerostatic critical wind speed of full bridges was also introduced by Simiu and Scanlan [2], assuming linear structural behavior and 1 DoF linear aerodynamics.It consists in solving the following eigenvalue problem: where F α is the torsional flexibility matrix, I is the identity matrix, and ∂CM ∂α is the slope of the aerodynamic moment coefficient (Eq.( 2)).Then, the critical wind velocity is obtained as: where d t is the higher eigenvalue, and ΔL is the deck element discretization length.The full derivation of this expression is given in Appendix A. This procedure requires obtaining beforehand the flexibility matrix F α , which in the case of a full bridge usually entails carrying out a series of static analysis consisting of applying modal unitary moments along the deck and obtaining the modal rotation vectors to create the matrix.

Nonlinear approaches
The nonlinear features of the problem described in Section 2.2 limit the applicability of the aforementioned classical linear approaches.In order to take into account all these nonlinearities, several methods have been proposed adopting a double loop scheme.The inner loop is dedicated to solving the equilibrium problem (Eq.( 1)) for a given wind velocity, while the outer loop updates the wind velocity value for the inner loop and seeks to find the critical wind velocity U cr .
On the one hand, the solution of the equilibrium problem in the inner loop is usually addressed using the linearized incremental equilibrium equation: where Δu is the displacement increment vector, and Δf s is the nodal force increment vector induced by the increment in the displacements.
The nonlinear features described in Section 2.2 can be introduced here by defining f s considering the aerodynamic nonlinearities and introducing in K(u) the full range of structural nonlinearities.The solution of Eq. ( 9) can be achieved by the nonlinear finite element method, as the one used in Cheng et al. [8], Boonyapinyo et al. [9,18], and Ge and Shao [23].Nonlinear finite element methods are computationally very intensive and rely on the definition of stiffness or flexibility matrices, or on directly iterating over the FEM of the bridge.
On the other hand, the outer loop has been traditionally managed by adopting an incremental strategy by updating the wind velocity adding a wind increment step as U iout+1 = U iout + ΔU, where i out is the iteration of the outer loop and ΔU the wind increment step.Examples adopting this sequential strategy can be found, for instance, in Cheng et al. [8,19].Later contributions improved the efficiency of the sequential approach by adopting a variable wind step ΔU to speed up the identification of the critical wind velocity, such as Cheng et al. [27] and Zhang et al. [28].Alternatively, the outer loop strategy proposed by Boonyapinyo et al. [18,9] permits a faster convergence to U cr by combining an eigenvalue Fig. 2. Standard response of torsional divergence (in red) and convergence (in green) of the nonlinear calculation of static rotation in a linear structure: (a) deck rotation as a function of the angle of attack; (b) convergence of the inner loop iterative nonlinear calculation of the static rotation (Eq.( 10)) when U < U cr ; and (c) divergence of the inner loop nonlinear calculation when U > U cr .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)analysis under wind loads and the updated wind-velocity bound algorithm to identify the next velocity increment, improving the performance of the sequential method.This approach effectively reduced the number of outer iterations but at the cost of a high computational burden in each iteration to define the next value of wind velocity to be analyzed.
This paper seeks to develop an efficient method that can be used in design applications.The goal is to obtain U cr in an efficient way using information that can be readily obtained from a FEM, in the same way as in the analysis of other wind-induced phenomena, such as flutter or buffeting.Hence, the proposed method combines modal formulation with a root-finding algorithm that speeds up the convergence of the outer loop.This avoids the definition of the stiffness or flexibility matrices and permits the identification of the critical wind velocity from the result of dynamic eigenvalue analyses of the bridge nonlinear structural FEM.

Proposed modal-based approach for aerostatic stability analysis
The purpose of the proposed method to study the aerostatic stability of long-span bridges is twofold: (1) to unify the structural inputs required to analyze aerostatic stability and other aeroelastic phenomena, such as flutter or buffeting: the modal matrix normalized to the mass and the vector of natural frequencies, and (2) to obtain the critical wind velocity in an efficient way in order to assess the performance of the structure in design frameworks.In the following, the fundamentals of the method are described in Section 3.1; then, the solution procedure assuming linear structural behavior is presented in Section 3.2, while an advanced nonlinear approach is developed in Section 3.3.

Basic features of the proposed modal-based approach
The vector of displacements and rotations of a given node j along the deck u j = [ u j , w j , α j ] under the action of the steady wind displacement- ) ] (Eq.( 2)) can be ob- where N is the number of modes, ϕ n is the modal vector of the n-th mode normalized to mass, and ω n stands for natural frequency (see Appendix B).These input values can be conveniently obtained from any FEM commercial software, making straightforward the implementation of the method.Since the aerostatic loads are a function of the angle of attack, this implicit expression must be solved iteratively to obtain the nonlinear nodal displacements of the deck.On the other hand, an efficient algorithm to find the critical wind velocity U cr can be formulated by analyzing the deck rotation path as a function of the wind velocity and the convergence of the inner rotationload loop.The standard response of a deck section under the action of aerodynamic loads for growing values of the wind speed U is shown in Fig. 2 in terms of deck rotation.It can be observed in Fig. 2 (a), that the response can be obtained for wind velocity values lower than U cr , where the nonlinear problem in Eq. ( 10) convergences to a real value, as shown in Fig. 2 (b).On the contrary, Fig. 2 (c) shows a classical divergence pattern for wind velocities higher than U cr .The differences in the convergence/divergence patterns can be used to identify when a given wind velocity U is higher or lower than U cr .Therefore, the wind velocity domain can be divided into two different regions: where the nonlinear problem converges, and where it does not.The instability bifurcation point can be found using any mathematical root-finding method such as the bisection method, the Newton-Raphson method, or optimization methods.This would represent an outer loop in the nonlinear problem described in Eq. (10).

Sequential increment of the wind velocity
A sequential method adopting an increasing wind velocity U iout+1 = U iout +ΔU until reaching the instability (divergence) condition can be implemented to find U cr .This procedure also provides the evolution of the nodal rotations and displacements of the bridge as a function of the wind velocity, in exchange of a higher computational burden.This is the approach adopted in most of the published studies addressing the aerostatic stability of long-span bridges, as discussed in Section 2.3.3.
Fig. 3 shows a flowchart of this procedure, where two loops can be identified.The inputs of the problem are the modal matrix ϕ n , the natural frequencies ω n , the force coefficients that are dependent of the angle of attack C D (α), C L (α), and C M (α), and the initial angle of attack α 0 .The inner loop is an implicit problem which seeks to obtain the deck displacements, particularly the deck rotation α, at a given wind velocity U.This is carried out iterating over Eq. (10), until the steady wind loads f s and the deck displacements u converge.Then, the outer loop moves to the next outer iteration i out by relaunching the inner loop with a higher wind velocity: U iout+1 = U iout + ΔU.This process is repeated until Fig. 3. Flowchart of the sequential increment solution procedure to solve the aerostatic stability problem using the proposed modal-based approach.
M. Cid Montoya et al. reaching a wind velocity for which the internal loop diverges, identifying the aerostatic instability (U cr ).

Direct search of U cr : Instability-finding algorithm
An alternative strategy can be adopted when the displacements trace is of no interest by taking advantage of the independence of the calculations for each wind speed.This is carried out herein by substituting the outer loop in Fig. 3 by a root-finding algorithm able to identify the instability bifurcation point, namely the wind velocity for which the convergence of the inner loop changes from convergence to divergence, as depicted in Fig. 2.
Algorithm 1 scheme shows how this method can be implemented adopting the bisection method.The inputs are again the modal matrix ϕ n , the natural frequencies ω n , the force coefficients C D (α), C L (α), and C M (α), and the initial angle of attack α 0 of each deck node j.The lower U a and upper U b bound of the domain are the initial wind velocity U 0 , which is a low value close to 0 m/s, and a very high value where the instability is guaranteed U ∞ .Then, a middle value is obtained as U m = (U a + U b )/2, and the bridge response at U m is calculated using Eq.(10).
A convergence condition C for the f s -u inner loop is defined to evaluate if the deck rotation at midspan α mid converges or diverges for each iteration of the outer loop i out .To do so, a maximum number of iterations of the inner loop is established I in .Then, the rotation of the deck at midspan at that number of iterations α mid (I in ), which is shown as a black diamond in Fig. 4 (a), is compared with the rotation at some point of the inner loop convergence history (α mid (γ⋅I in )), represented as blue dots.If α mid (γ⋅I in ) is higher than α mid (I in ) multiplied by γ (e.g.γ = 0.5), the inner loop is converging, while otherwise is diverging.The condition is: It can be seen in Fig. 4 (a), that for the iterations of the outer loop (i out ) where the inner loop converges (C > 0, green cases, as previously shown in Fig. 2 (b)), the normalized rotation (blue dot) will be higher than half the value of the normalized rotation at I in , represented as a red star.This leads to a positive value of C. On the other hand, the value of C is negative in the divergent cases (C < 0, red lines, as in Fig. 2 (c)).This difference can be used to divide the domain into two regions for the inner loop iterations as shown in Fig. 4 (b): a convergence region (green), and a divergence region (red).With this criterion established, the root-finding algorithm can easily identify the wind velocity at which the instability occurs.) ) ) Fig. 4 (b) shows how the bisection method is able to iterate using the condition defined in Eq. ( 11) until reaching convergence in a reduced number of outer iterations.The outer loop convergence is achieved when the difference between U b and U a is lower than an imposed convergence tolerance ∊ U or when the number of outer iterations i out reach the maximum number of outer iterations allowed I out , as described in Algorithm 1. Several root-finding algorithms can be used for this task; however, as illustrated in Fig. 4 (b), the condition is flat in two wide regions of the domain: for relatively low wind velocities and relatively high wind velocities.In the low velocity range, this is because the value of α mid (γ⋅I in ) is equal to α mid (I in ) since the convergence of the inner loop is achieved in a reduced number of inner iterations lower than γ⋅I in .Hence, C takes the value of γ.In the high velocity range, this occurs when the value of γ⋅α mid (I in ) is much higher than α mid (γ⋅I in ) (see i out #2 in Fig. 4 (a)) Fig. 4. Representation of the convergence criterion of the root-finding algorithm (Eq.( 11)).(a) shows the convergence of the f s -u inner loop, in terms of α, for different wind velocities (outer loop iterations i out ).It must be noted that they are normalized to the deck rotation value at midspan (α mid (I in ), black diamond) in order to compare the convergence histories (see Fig. 2).Blue dots identify the value of α mid (γ⋅I in ) used in Eq. ( 11 and therefore C→ − γ (see Fig. 4 (b)).In both cases, this leads to a constant value of C (see Eq. ( 11)) for that range of wind velocities, as shown in Fig. 4 (b).This encourages the use of gradient-free algorithms like the bisection method.

Proposed modal-based iterative approach for aerostatic stability analysis considering nonlinear structural behavior
An advanced version of the method for the assessment of U cr in nonlinear structures is presented here.This approach is sketched in Fig. 5 and consists of two stages.The first phase seeks to obtain the estimated critical wind velocity taking into account the structural nonlinear properties of the bridge at the no-wind or initial wind scenario U cr,0 .Hence, before conducting the eigenvalue analysis in the linear modal-based approach (see Algorithm 1), the prestressing forces and the Young's moduli of each stay must be updated considering the cable-sag effect and other nonlinear structural features such as beam-column effect and large displacements through the nonlinear structural FEM.This is depicted in Fig. 5 in the "Deck zero-displacement N-E eq loop", where the prestressing forces N and the equivalent Young's moduli E eq are iteratively updated.The prestressing forces are obtained by solving the system of equations that impose zero nodal displacements in the deck under self weight load adopting the zero displacement method (Wang et al. [39]), and their values are kept constant along the aerostatic stability analysis since they are independent of the wind scenario.The equivalent moduli are then obtained using Eq. ( 3) for each stay, depending on its cross-section area, its Young's modulus, its length, and stress level.These values will be modified as the wind loads change for each wind velocity.Once this loop converges, the eigenvalue analysis of the bridge FEM provides the modal matrix Φ normalized to the mass and the natural frequencies ω, which are different from the ones obtained using the linear analysis due to slight changes in the stays' stiffness.Then, the root-finding modal-based approach (Algorithm 1) finds the estimated critical wind velocity considering the bridge properties at the initial wind scenario U cr,0 .
In the second phase, a root-finding method based on the bisection scheme seeks the identification of the critical wind velocity updating the structural nonlinearities for each wind scenario.The lower and upper wind velocities are set as the initial wind velocity U 0 and the critical wind velocity obtained in the previous phase U cr,0 .This allows the definition of a mean wind velocity U m,0 = ( U 0 + U cr,0 )/ 2, which is the wind scenario to be evaluated in the first iteration of the outer loop i out .Then, the next task is to complete the nonlinear rotation f s -u-E eq loop, which finds the deck rotation equilibrium considering the changes in the stays' stiffness caused by the wind-induced displacements.The first step in this process is the completion of the f s -u loop, which can be carried out using the nonlinear FEM or Eq. ( 10) with the dynamic information taken from the nonlinear FEM.The resulting displacements require updating the equivalent stiffness of the stays, which is identified in Fig. 5 as "Update E eq ".The most influential effect on the cable stiffness is the vertical displacements of the deck anchorages, which is caused by the lift-induced vertical displacement of the deck, and the rotation caused by the moment component.If upward vertical displacements occur, the tension of the stays will decrease according to Eq. ( 3) and subsequently, its stiffness may be reduced.On the other hand, when the wind loads cause a downward vertical displacement, the stiffness of the stay is expected to increase according to Eq. ( 3), unless the maximum stress level is surpassed leading to the collapse of the stay.Furthermore, horizontal displacements in the tower may cause similar effects, which are particularly important for the backstays' stiffness.Once the stiffness of the entire cable-supporting system is updated for the current wind velocity U m , the eigenvalue analysis is conducted to obtain the bridge modal properties.Then, the linear modal-based approach (Section 3.2.2) is used to obtain an estimated critical wind velocity U L cr for the current properties of the bridge at that wind velocity.
As the wind velocity grows, the stiffness of the bridge is expected to decrease (see Zhang et al. [26]).Therefore, the estimated critical wind velocity U L cr for the updated bridge dynamic properties will decline.Hence, there will be a wind velocity value where the estimated critical wind velocity U L cr is equal or lower than the current wind velocity U, which will correspond to the instability point.This situation is modeled by the following nonlinear convergence criterion C NL : ) This process allows the estimation of the critical wind velocity without defining any nonlinear stiffness matrix, only using nonlinear FEM evaluations.The linear approach contributes to speeding up the nonlinear calculation of the critical wind velocity.

Application examples
Three application examples are used in this paper to demonstrate the capabilities of the proposed method to directly find the aerostatic critical wind velocity, which are summarized in Table 1.The first case is a 1node system representing a bridge deck restrained by a torsional spring under the action of the wind-induced torsional moment, resulting in a linear one-degree-of-freedom problem.For this application case, the Fig. 5. Flowchart of the proposed iterative nonlinear root-finding modal-based approach.
M. Cid Montoya et al. performance of the proposed algorithm is compared with the linear method for 1 DoF (Eq.( 4)).The second example shows the performance of the proposed approach considering a full bridge model, hence an Nnodes dynamical system.The results are compared with the ones obtained solving the eigenvalues problem based on the flexibility matrix.In this case, the structural performance is assumed to be linear, while the nonlinear aerodynamics effects and the influence of the three components of the wind-induced static loads are discussed.In the last example, all nonlinearities are considered.The nonlinear behaviour of the bridge is modeled by conducting nonlinear FEM calculations, including beamcolumn effects, large displacements, and the stays' sag effect using Ernst equation (Eq.( 3)).
The deck cross-section adopted for the three application examples is the Scanlan's G1 section ( [40]) with a deck width B = 40 m.The aerodynamic coefficients are the ones reported in Cid Montoya et al. [41] obtained by sectional model wind tunnel tests at a Reynolds number Re = 4⋅10 5 , which are shown in Fig. 6, following the definition in Eq. ( 2).It must be noted that the drag and moment coefficients are multiplied by 2 in this figure.The sign convention adopted is that the positive values are down-wind for the drag, upward for the vertical force, and nose-up for the moment.The slopes of the experimental force It should be noted that any complex aerodynamic behavior that may be caused by complex bluff deck geometries, ancillary facilities or appendages, such as crash-barriers, handrails, wind barriers, or guide vanes, among others, will be reflected in the force coefficients (see, for instance, Buljac et al. [42]).Hence, the applicability of the method is independent of the bridge deck aerodynamics.On the other hand, the aerodynamic contribution of moving elements such as road and railway traffic to the bridge aerodynamics (see Cai et al. [43]) may be omitted since aerostatic stability is an ultimate limit state that may only occur at high wind velocities when the bridge traffic is typically closed.

Example #1: 1DoF-1node system
This application example simulates the response of a cable-stayed bridge at midspan by adopting a 1DoF-1node dynamical system, as the one shown in Fig. 7, with the mechanical properties reported in Table 2.This means that only the torsional component of the aerodynamic load is considered in this application case, and it is modeled by the moment coefficient at α = 0 • (C M,0 • ) and its slope computed with the values at α = 0 • and 2 resulting in a linear relationship for the moment with the angle of attack that can be written as: The value of the aerostatic critical velocity obtained for this application case using the linear method (Eq.( 4)) is U cr = 165.115m/s.This problem can be also solved by using the algorithms proposed in Section 3.2 based on the formulation developed to obtain the responses of the bridge using the modal analysis (Eq. ( 10)).Fig. 8 shows the results obtained using both the sequential algorithm, which allows the estimation of the entire rotation curve as a function of the wind velocity (see Section 3.2.1),and the root-finding algorithm, which permits a fast convergence to identify the critical wind speed in a reduced number of iterations (see Section 3.2.2).The sequential approach calculates the deck rotation with a wind increment of ΔU = 1 m/s.The root-finding algorithm uses the bisection method within an interval of wind velocities between U a = 1 m/s and U b = 250 m/s, the number of iterations of the M-α inner loop has been I in = 100, and the factor γ is 0.5 (Eq.( 11)).It can be observed in Fig. 8 that the three methods converge to the same value, and the root-finding approach is able to reach the vicinity of the critical velocity in just a few iterations.As it can be seen in Table 3, a difference lower than 1 m/s in the estimation of the critical velocity is obtained in the eighth iteration.This ability will be later of utmost importance in the nonlinear application cases, where the computational Fig. 7. Representation of the 1DoF-1node dynamical system.

Table 2
Structural and aerodynamic properties used in the first application example.
M. Cid Montoya et al.
burden of the FEM analysis will highlight the advantages of this method.
It must be noted that the value of the inner loop convergence condition (Eq.( 11)) when approaching the solution may be affected by the maximum number of iterations used in the inner loop I in and the value of the parameter γ.However, since the convergence criterion is based on the convergence or divergence of the inner loop, the solution will not be affected.

Example #2: 3DoF-Nnode cable-stayed bridge with linear structural behavior
The proposed approach is also tested on a long-span cable-stayed bridge with a main span of 1316 m, in the order of magnitude of the Russky, the Sutong, or the Stonecutters bridges, or the bridge studied in the paper by Nagai et al. [12].The deck cross-section is the same as the one described in Section 4 (see Fig. 6).
The layout of the bridge is shown in Fig. 9, and the sizing of the cablesupporting system and the deck plate thickness are the ones obtained after completing a structural optimization in Cid Montoya et al. [44].The main values are the following: deck plate thickness t = 2.27 cm, backstay cross-section area A B = 0.543 m 2 , and average cross-section area of the stays A s = 0.04 m 2 .A deck non-structural mass of 8 T/m and a non-structural mass moment of 2300 Tm 2 /m are adopted.The Young's modulus for the steel of the stays is E stays = 190 GPa.A summary of the main natural frequencies and mode shapes is provided in Fig. 10.
The structural response of the bridge is assumed to be linear in this application example.Hence, the only input data required to carry out all the calculations reported in this section is a single FEM dynamic eigenvalue analysis with the aforementioned characteristics.

Table 3
Convergence evolution of the root-finding algorithm in terms of the interval middle velocity U m (see Algorithm 1), the interval of each outer iteration U b − U a , the angle of attack at the midspan α mid and the value of the convergence condition C according to Eq. (11).Fig. 9. Layout of the cable-stayed bridge used as application example with the Scanlan's G1 deck cross-section.Dimensions in m.The numeration adopted for the stays is provided along the deck.
Fig. 8. Aerostatic response of the 1DoF-1node dynamical system.Comparison of the results obtained using the linear method, and the sequential and root-finding modal-based approaches.

Results considering linear 1DoF aerodynamics
The evaluation of the aerostatic stability of the bridge using the flexibility matrix method described in Section 2.3.2 is used as a reference to assess the capabilities of the proposed modal-based approach.The torsional flexibility matrix obtained from the FEM of the full bridge is shown in Fig. 11, where the torsional rotation along the deck for unitary torsional moments (M t = 1 Nm) applied along the deck central X axis (see Fig. 9) at the location of the stays anchorages is shown.Therefore, this flexibility matrix is defined using 40 control points.As an illustrative example, the torsional rotation obtained along the deck when the unitary moment is applied to the deck longitudinal axis at point #16 (stay #16) is depicted in red color in Fig. 11.The problem described in Eq. ( 7) and Eq. ( 8) is solved first assuming linear aerodynamics and 1DoF; hence the moment component of the steady aerodynamic action is modeled by the moment coefficient and its slope at α = 0 • (see Table 2).The critical wind velocity obtained adopting the above assumption is U cr = 198.91m/s.
Alternatively, the modal-based root-finding approach for the full bridge model is applied in the following.The input data are the dynamic properties of the full bridge considering the first 22 modes, the most relevant are provided in Fig. 10, and the aerodynamic properties given in Table 2.The mode shapes are defined by 171 nodes uniformly distributed along the deck.In order to compare these results with the ones obtained using the flexibility matrix method, only the moment coefficient and its slope have been adopted to model the aerodynamic loads.The performance of the modal-based root-finding algorithm is reported in Fig. 12.The maximum number of iteration in the rotation inner loop is I in = 100, and the interval of wind velocities ranges from It can be noted in the convergence history shown in Fig. 12 (a), which is similar to the conceptual convergence diagram sketched in Fig. 4, that the critical wind velocity is identified in a small number of iterations.In Fig. 12, the convergence criterion is calculated considering the deck rotation at midspan α mid .The result obtained is U cr = 198.65 m/s, very close to the results obtained applying the flexibility matrix method, with a difference in absolute value of Δ Ucr = 0.26 m/s.This validates the proposed modal-based procedure for the full bridge model.Fig. 12 (b) shows the torsional rotation along the deck evaluated for the first two iterations of the root-finding algorithm.It can be observed that in the first iteration of the outer loop (U = 125.5 m/s), the rotation of the deck is relatively low, while for the second iteration, which is close to the instability point, the rotations are very high.It must be remarked that the goal of the proposed approach is the identification of the instability point, and for wind velocities close to the critical wind velocity, the structural response could be unrealistic.This will be further addressed when the application of the nonlinear version of this method Fig. 11.Torsional flexibility matrix F α (Eq.( 7)) represented as the rotation caused by each unitary moment over all sections.The rotation along the deck caused by the unitary torsional moment at point # 16 is depicted in red as an illustrative example. is discussed.Fig. 12 (c) shows the convergence history of the rotation inner loop for these two first iterations, where it can be noted that, for structures assuming linear behavior, the convergence of the inner loop is slower as the outer loop algorithm gets closer to the solution, as anticipated in Fig. 4 (a).

Discussion on the influence of 3DoF nonlinear aerodynamics in the aerostatic stability
The modal-based approach based on Eq. ( 10) permits obtaining the bridge response including the three components of the aerodynamic steady loads.Therefore, the critical velocity can be estimated considering the 3DoF nonlinear aerodynamic response using the root-finding algorithm, and its influence in the bridge response can be studied taking advantage of the sequential method (see Section 3.2.1).
Five cases are studied herein using the modal-based approach to study the influence of nonlinear aerodynamics, which are summarized in Table 4. Case #1 is the one reported in Section 6.1, which models the moment wind load component using the moment coefficient at α = 0 • and its slope.Case #2 improves the previous one by adopting the experimental nonlinear curve of the moment coefficient as a function of the angle of attack, as depicted in Fig. 6.The third case introduces the three components of the aerodynamic steady wind load, considering the linear approximations based on the force coefficients at α = 0 • and their slopes.
Therefore, the lift and drag coefficients are estimated from their slopes ) as: Case #4 advances the previous one by estimating the drag force coefficient with a quadratic approximation fitting the values of the moment coefficient at α = 0 • , 2 • , and 4 • , whose values are 0.057, 0.059, and 0.067, respectively, resulting in: a "L" = linear approximations provided in Eq. ( 13)-( 15); "Q" = quadratic aproximation given in Eq. ( 16); "Exp" = experimental curves reported in Fig. 6.The approximations for the lift and moment coefficients are kept linear since this assumption is sufficient to accurately model the aerodynamic response of the streamlined deck used as application example.The last case (Case #5) improves the previous ones based on linear and quadratic aerodynamics by adopting the experimental nonlinear aerodynamics for the three force components.The result obtained using the flexibility matrix is identified as Case #0 and it is used as reference.
It can be seen that when the experimental curve of the moment coefficient (case 2) is adopted instead of the approximation based on the moment and the slope at α = 0 • (case 1 in Table 4, see Eq. ( 13)), the critical velocity grows.This is so because the linear approximation of the moment coefficient based on the slope (Eq.( 13)) overestimates the value of C M at large angles of attach, as shown in Fig. 6.On the other hand, when the nonlinear formulation for the three components is adopted, the critical velocity decreases due to the contribution of the draginduced nose-up rotation α D , as shown in Fig. 13.This effect, which is controlled by the structural configuration of the bridge, can be also identified in mode #1 as shown in Fig. 10, where positive lateral displacements provoke negative rotation according to the global coordinate system (Fig. 9), which result in nose up rotations of the deck (see Fig. 13).The opposite effect was found in Boonyapinyo et al. [9] for the Akashi Kaikyo Bridge, where the drag induced rotation delayed the critical speed.
It must be noted that the lateral stiffness of the bridge is much lower than the torsional and vertical ones, as it can be deduced from the values of the natural frequencies shown in Fig. 10.This remarks the noticeable influence of the drag load on the rotation of the deck, and the necessity of considering the 3 DoF wind load model even when adopting linear structural analyses.This effect is highlighted in Cases #4 and #5, where the drag influence is even more important as the angle of attack grows.
The aforementioned cases can be further analyzed through the study of the evolution of the deck rotation as a function of the wind velocity applying the sequential method.In Fig. 14, the large difference between the critical speeds of cases 1 and 2 is clearly caused by the change in the model relating the moment coefficient and the angle of attack for α⩾6 • , as shown in Fig. 6.On the other hand, the influence of the drag-induced Fig. 14.Deck rotation at midspan obtained for the five cases described in Table 4. Fig. 15.Effects of the initial wind angle of attack on the critical wind velocity (a) and the torsional behavior at midspan (b).The sign criterion adopted is shown in (a).rotation is noteworthy when comparing cases 1 and 3, where the larger rotation in case 3 is due to the drag contribution.However, this effect is more important when the complete curve of the drag coefficient is considered, as shown in case 5.This is so because the slope of the drag increases noticeably for α⩾8 • (see Fig. 6), and the torsional rotation remarkably grows, leading to an important reduction of the critical wind speed.It must be noted that the quadratic approximation for the drag coefficient curve combined with the linear approximations of the lift and moment coefficients (case 4) provides a very accurate estimation of the critical wind velocity (see Table 4).However, the rotation versus the wind velocity curve presents some discrepancies for large angles of attack.The applicability of these approximations will be further discussed in Section 7, where the nonlinear structural features of the bridge will also be considered to show how important is the structural response of the bridge for large angles of attack.

Discussion on the influence of the initial angle of attack
The torsional behavior of the deck at midspan considering several initial angles of attack has been analyzed and the results are reported in Fig. 15.The initial angle of attack α 0 makes an impact in the aerodynamic nonlinearities since the nonlinear effects produced by the change in the force coefficients as a function of the angle of attack will be anticipated or delayed depending on the initial value of the wind angle of attack.The direction of the deck rotation depends on the curve of the moment coefficient, which changes its sign around α = − 2 • , and the drag-induced rotation α D of the deck.The critical initial wind angle of attack is around α 0 = − 3.66 • , as shown in Fig. 15.When the initial angle of attack reaches lower values, the torsional divergence is delayed since the sign of the moment coefficient turns negative and the draginduced rotation now hinders the deck rotation.Similar behavior was found in Arena and Lacarbonara [22] for the Hu Men suspension bridge.On the other hand, the critical wind velocity obtained for positive initial angles of attack is very similar, around U cr = 174 m/s, but the rotation of the deck is larger as the initial angle of attack grows.

Example #3: 3DoF-Nnode cable-stayed bridge with nonlinear structural behavior
The last example is presented aiming at testing the proposed approach in a problem considering simultaneously aerodynamic and structural nonlinearities.The nonlinear behavior of the stays, particularly the sag effect, is modeled by the Ernst equation (Ernst [30]) as previously described in Eq. ( 3).Gravity loads, stays prestressing forces, and other structural nonlinearities are considered in the FEM, such as the beam-column effect and large deformations.The material behavior is assumed to be linear.Deck aerodynamic characteristics are the ones described in the previous example.The FEM is the same as in the example used in Section 6.However, as described in Fig. 5, the sag effect should be evaluated at the initial stage of the process, so slight variations in the stiffness of the stays may be found, as it will be shown in Fig.  First, the evolution of the bridge deformation with the increase of the wind velocity until reaching the collapse situation is analyzed through nonlinear FEM analyses in Section 7.1.Then, the nonlinear version of the root-finding algorithm is applied to obtain the critical wind velocity.Finally, results are presented and compared in Section 7.2.

Responses obtained using sequential analysis
Fig. 16 shows the evolution of the deck displacements as the value of the wind velocity increases using the nonlinear method reported in Section 3.3 but adopting a sequential wind velocity increase to trace the bridge responses versus the wind velocity.The results at midspan are compared in Fig. 16 (a, c, and e) with the linear response obtained in Section 6.On the other hand, Fig. 17 shows the evolution of several important bridge structural properties as a function of the wind velocity, such as the main vertical and torsional natural frequencies, the variation of the stiffness of the stays due to the sag effect, and the loss of stress as the wind speed increases.
Focusing first on the deck displacements, it is interesting to note that the response for linear and nonlinear cases are very similar for wind velocities lower than 100 m/s, which shows that nonlinear effects do not play a relevant role for those wind velocities.Differences in the vertical response shown in Fig. 16 (c) for the range of studied wind velocities are due to the lower value of Young's modulus of the backstays caused by the sag effect, as reported in Fig. 17 (c).For wind velocities higher than 100 m/s, some differences are apparent in Fig. 16, which are caused by the nonlinear structural behavior of the bridge, particularly for the vertical response.It may be worth mentioning that a growing separation is found between the linear and nonlinear vertical responses between U = 100 m/s and 125 m/s.This is caused by the continuous reduction in the stiffness of the backstays (see Fig. 17 (c)), which leads to a reduction in the vertical stiffness of the bridge, as described in Fig. 17 (a) for the first vertical natural frequency (mode 5, see Fig. 10).This phenomenon is explained by the abrupt loss of stiffness of the backstays and the continuous loss of stress of the central span stays shown in Fig. 17  (d), leading to the eventual reduction of the stiffness of those stays.The whole mechanism is explained in Fig. 18, where it can be observed that the upward displacement of the deck leads to the loss of stress at the stays, horizontal displacements at the top of the towers also take place, producing the subsequent reduction of the stresses in the backstays.As the lift grows with higher speeds, the load generates a reduction in the stays stress which slows the growing vertical displacement, as shown in Fig. 16 (c).
It must be highlighted a behavioral change that is found for U = 130 m/s and U = 135 m/s.It can be seen in Fig. 16 (d) that the vertical displacement of the side span changes from downwards to upwards due to the change in the sign of the lift load that happens around α = 0.25 • (see Fig. 6).Given the non-symmetrical distribution of stays' crosssection area along the deck (see Cid Montoya et al. [44]), the change   The torsional response is also affected by the loss of stress and stiffness in the windward stays, as shown in Fig. 18 (b) and Fig. 17.First, the loss of stress in the windward stays (Fig. 17 (d)) slows the increasing torsional response with the wind velocity increase.However, as the backstays lose their stiffness, the overall torsional stiffness of the bridge decreases, as it can be observed through the natural frequency of mode 12 in Fig. 17 (b), leading to the abrupt loss of stiffness in other modes, such as mode 11, and the bridge instability.
It must be noted that the instability occurs with reduced values of displacements and rotations, as it was also found in the cases studied in Boonyapinyo et al. [18] and Nagai et al. [12].The phenomena is controlled by the stiffness degradation shown in Fig. 17, as it was also described in the paper by Zhang et al. [26].

Root-finding approach
The information provided in Fig. 16 and Fig. 17 shows that the critical wind velocity should be in the range U = [135, 140] m/s.However, the precise value should be found by an efficient automatic procedure without requiring the study of the entire collapse sequence, as analyzed in Section 7.1.As described in Fig. 5, the nonlinear rootfinding approach presented in this study is based on the estimation of the critical wind velocity of the structure for a given wind velocity.Fig. 19 shows the value of the estimated critical wind velocity U L cr when the bridge properties at different wind velocities are considered.It may be worth mentioning that before the nonlinear structural effects become relevant, i.e.U < 100 m/s, the estimated critical velocity is U L cr = 175.59m/s, similar to the linear value reported in Table 4 for case #5, as expected.However, as the wind velocity U increases and consequently the nonlinear effects diminish the stiffness of the bridge (see Fig. 17), U L cr drops.When the estimated critical wind velocity U L cr is equal or lower than the current wind velocity U, the instability is expected to occur, which is reflected in the convergence criterion C NL (Eq.( 12)).
Fig. 20 reports the convergence history of the nonlinear root-finding algorithm.The values of the convergence condition C NL are obtained following Eq.( 12) and can be inferred from those previously reported in Fig. 19.The bisection method is applied in the range U = [1, 175.59] m/s, and the value obtained as critical wind velocity is U cr = 138.5 m/s.It can be observed that a fast convergence is obtained with a relative error lower than 1 m/s in only 6 iterations, highlighting the efficiency of the approach.
It must be noted the large difference between the result obtained using the linear method (U cr = 176.04m/s, see Section 6) and the one reported here considering the structural nonlinearities (U cr = 138.5 m/s), which represents a reduction of 21.32%.This emphasizes the importance of considering all the nonlinearities involved in the phenomenon.Furthermore, the relevance of the aerostatic stability analysis in the bridge design must be highlighted, since the flutter critical velocity reported in Cid Montoya et al. [44] for this application case is U Exp f,cr = 127.36m/s.The similarity between these values of critical wind velocities emphasizes that in the design of ultra long-span cable-stayed bridges the phenomenon of aerostatic stability can be as relevant as flutter.

Discussion on the influence of the aerodynamic nonlinearities in structures showing nonlinear behavior
When structural nonlinearities are considered, the aerostatic instability occurs due to stiffness degradation at low values of displacements, particularly for the deck rotation, whose value is around 2.5 • (see Fig. 16).Similar behavior can be found in the investigations by Cheng et al. [8], Boonyapinyo et al. [9], Cheng et al. [19].This fact can affect the influence of the nonlinear features of the deck aerodynamics on the instability mechanism.
Table 5 reports the critical wind velocity obtained considering four different sets of curves for the force coefficients.In Case I, the experimental force coefficients reported in Fig. 6 are used, which is the same case reported above in Section 7.2, and it is used here as a reference.Case II reports the result obtained when linear approximations are adopted for the three force coefficients (Eq.( 13)-( 15)).It must be noted that there is an important difference in the estimated critical wind velocity, mistakenly increasing the critical wind velocity and the safety of the bridge design.Therefore, the deck aerodynamics cannot be linearly approximated.Case III seeks to analyze the accuracy of the linear approximation in the lift and moment coefficients by adopting these approximations along with the drag experimental curve.The result is exactly the same as the one obtained in Case I, which reveals that the differences between Case I and II are caused by the drag coefficient.Case IV substitutes the drag experimental curve by the quadratic approximation given by Eq. ( 16), and the result is very accurate.This shows that the linear approximations for the lift and moment coefficients and the quadric approximation for the drag coefficient guarantee an accurate approximation of the experimental coefficients in the range α = [ − 4 • , 6 • ], which is sufficient for the accurate identification of the aerostatic critical wind velocity the bridge.The last conclusion can be considered for design purposes.

Concluding remarks
This study proposes an alternative approach for the estimation of the aerostatic stability critical wind velocity of long-span bridges.The method consists of the assessment of the bridge static deflection using modal formulation instead of the classical stiffness or flexibility matrices.The input parameters related to the structural properties of the bridge are the modal matrix normalized to the mass and the natural frequencies, which facilitates the aerostatic stability analysis and helps to simplify an eventual multi-response aeroelastic analysis of a longspan bridge, required in comprehensive design frameworks for preliminary design stages.A nonlinear version of the method is proposed to take into account the structural nonlinearities of the bridge.The proposed method is tested for three application examples considering an increasing number of nonlinearities and its accuracy is assessed by comparison with the results obtained using traditional methods.These comparisons have shown that the proposed method is simple, accurate, efficient, and reliable, and it is very suitable for being used within comprehensive design frameworks taking advantage of commercial FEM software.The method enables the identification of the aerostatic critical wind velocity requiring about 6 iterations of the instability-finding algorithm with an expected error lower than 1 m/s for the critical wind velocity.Forthcoming studies may address the use of the proposed algorithm in design frameworks to carry out the aerostatic shape a "L" = linear approximations provided in Eq. ( 13)- (15); "Q" = quadratic aproximation given in Eq. ( 16); "Exp" = experimental curves reported in Fig. 6.
M. Cid Montoya et al. optimization of the bridge deck.Furthermore, the tests and discussions presented in this study provided some interesting findings related to the aerostatic stability phenomenon and the analysis of long-span cable-stayed bridges.First, it is important to mention the important role of the drag coefficient in the aerostatic stability when the three components of the wind displacement-dependent loads are considered.The sign of the draginduced rotation determines if the structure will increase or decrease the critical wind velocity.On the other hand, the lift coefficient, which does not show any influence on the results when assuming linear structural behavior, is very important when structural nonlinearities are considered.This is so because when the lift component induces deck vertical upward displacements, the stays of the bridge reduce their level of stress, leading to a subsequent loss of stiffness due to the sag effect.This phenomenon is also influenced by the deck rotation, caused by both the pitching moment and drag-induced rotation, which increases the vertical displacement of the windward deck anchorages reducing, even more, the stiffness of the stays.This eventually leads to the stiffness degradation of the deck, which can be clearly identified by the drop in the main natural frequencies of the bridge, speeding-up the collapse of the structure.The aerostatic instability can take place even with a reduced value of the rotation of the deck, and showing lower displacements than in the linear analysis.This reduces the influence of the aerodynamic nonlinearities and consequently the amount of required data to characterize the aerodynamic response.In the application case adopted in this study, linear approximations for the moment and lift coefficients, and the quadratic approximation for the drag coefficient were sufficient for accurate identification of the critical wind velocity of the bridge.
The methodology proposed in this study was developed considering steady wind loads, which is the general approach to estimate the aerostatic stability critical wind velocity of long-span bridges.Future work would advance the proposed framework to consider the analysis under unsteady wind loads, turbulent fluctuations, non-stationarity, and transient winds, among others.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Algorithm 1 .
Instability-finding modal-based algorithm procedure IDENTIFICATION OF CRITICAL WIND SPEED INPUT: ϕ n , ωn, CD(α), CL(α), CM(α) Ua = U0; U b = U∞; α = α0; iout = 1 while iout < Iout do [Outer loop] Um = (Ua + U b )/2 for iin = 1 : Iin [Inner loop] for j = 1 : ) to calculate the value of C for each iteration of the outer loop.(b) shows how the root-finding algorithm (outer loop) finds U cr using the value of C. Inner loop convergence is represented in green, while divergence is shown in red.The outer loop iterations (Bisection convergence) are plotted in blue.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)M. Cid Montoya et al.

Fig. 6 .
Fig. 6.Force coefficients as a function of the angle of attack α according to Eq.

Fig. 10 .
Fig. 10.Mode shapes and natural frequencies of the bridge at wind velocity U = 0 m/s.Natural frequencies in Hz.The vertical continuous gray lines show the position of the towers, while the gray dashed lines indicate the midspan location.

Fig. 13 .
Fig. 13.Drag-induced deck rotation α D caused by the standalone steady drag force D s .The red lines represent the transversal bars used in the FEM to link the stays with the longitudinal beam that models the bridge deck, which are used here to visualize the deck rotation.Figure (a) shows the FEM; (b) shows a detail of the rotation at midspan; and (c) shows a rendered visualization of the deck rotation.

Fig. 12 .
Fig. 12. Results obtained for the full bridge model with linear structural behavior and deck aerodynamics using the modal-based root-finding algorithm.(a) shows the convergence of the root-finding algorithm (outer loop); (b) shows the rotation along the deck for the two first iterations; and (c) shows their convergence in the rotation M(α)-α loop (inner loop).

Fig. 17 .
Fig. 17.Evolution of bridge properties under increasing wind velocity.(a) shows the first vertical natural frequency; (b) describes the evolution of relevant torsional natural frequencies; (c) shows the equivalent Young's modulus E eq of some stays; and (d) shows the stress level of some windward stays.Stays numeration according to Fig. 9.

Fig. 18 .
Fig. 18.Nonlinear deformation of the bridge at wind velocity U = 135 m/s.(a) shows the deformation obtained with the FEM with a deformation scale factor 50 and the mechanism that leads the stays to a loss of tension and ultimately of stiffness is shown.(b) provides a schematic detail of the deck deformation at the deck location of stay #20.

Fig. 19 .Fig. 20 .
Fig.19.Evolution of the estimated critical wind velocity U L cr for each value of U compared with the linear (L α mid ) and nonlinear (NL α mid ) deck rotations at midspan as a function of the wind velocity U.The stability limit based on the convergence condition C NL = 0 (Eq.(12)), this is U L cr = U, is shown as a reference. .

Table 1
Summary of the application examples used in this study.

Table 4
Critical wind velocities obtained by the classical flexibility matrix based approach (Section 2.3.2) and the proposed modal-based approach (Section 3.2.2) considering different aerodynamic models.

Table 5
Critical wind velocities considering different configurations of aerodynamic forces.