Lyapunov stability of suspension bridges in turbulent �ow

In the era of sleek, super slender suspension bridges, facing the issue of stability against dynamic wind actions represents an increasingly complex challenge. Despite the signiﬁcant progress over the last decades, the impact of atmospheric turbulence on bridge stability remains partially not understood, evoking the need for innovative research approaches. This study aims to address a gap in current research by investigating the random ﬂutter stability associated with variations in the angle of attack due to turbulence, which has not formally been addressed yet. The present investigation employs the 2D rational function approximation (2D RFA) model to express self-excited forces in a turbulent ﬂow. The application of this type of models to bridge dynamics yields a viscoelastic coupled dynamic system characterized by memory eﬀects and driven by broad-band long-time-scale noise, described here by a linear homogeneous time-variant diﬀerential equation, which shows apparent nonlinear features, and which has rarely been matter of research. Utilizing a Monte Carlo methodology, this work innovates in applying the largest Lyapunov exponent (LE) and the moment Lyapunov exponents (MLE) to study bridge random ﬂutter stability. The calculation of LE and MLE under diverse turbulent wind conditions uncovers lower ﬂutter stability than without turbulence eﬀects. In most cases, sample and low-order p -th moment stability thresholds closely align with the bridge dynamic response pattern; therefore, the ﬂutter critical wind speed is unequivocal. However, under certain turbulence scenarios, it is necessary to resort to MLE for a complete description of stability, evoking some additional consideration of which statistical moments should be considered for the engineering assessment of the ﬂutter limit. Finally, this work provides a qualitative insight into the instability mechanisms


Introduction
The quest to span ever greater distances through the construction of super-long suspension bridges has led civil engineering into previously uncharted territories.Notably, the evolution from the Storebaelt Bridge (Denmark), once representing a benchmark for aerodynamically streamlined long-span bridges, to the recently completed C ¸anakkale Bridge (Turkey) illustrates a significant progression in engineering ambition and capability [1].The C ¸anakkale Bridge, with its main span exceeding that of Storebaelt by about 25%, set a new standard for what is achievable with a new concept of streamlined bridge cross sections.However, the proposed spans for projects like the Sulafjørden Bridge (Norway; for information about the Coastal Highway Route E39 project, see [2]) and Messina Strait Bridge (Italy; see, e.g., [3]), potentially surpassing 3000 m, are poised to redefine even further these standards.This jump (Messina Bridge would be 63% longer than the C ¸anakkale Bridge) signifies a previously unimaginable change of perspective concerning the complex challenges of wind-structure interaction.Indeed, as the span lengths of bridges venture into new magnitudes, the implications of nonlinear aeroelastic effects may intensify, transitioning from mere considerations to potentially dominant factors that govern structural behavior, evoking a deeper understanding and innovative approaches to ensure their safe design.
Particular attention is given to wind-induced catastrophic instabilities, such as classical flutter for long-span suspension bridges equipped with modern aerodynamic deck cross sections.Determining the flutter stability threshold is often crucial before evaluating the dynamic response to turbulent wind (buffeting) of long-span suspension bridges.Indeed, it stands as a primary concern for designers, who need to make decisions regarding the aerodynamic feasibility of the bridge.Moreover, the importance of an accurate estimate of flutter critical wind velocity is underscored by codes and design specifications, which establish different return periods to check bridge buffeting response and flutter stability (see, e.g., [4,5]).
Traditionally, the analysis of bridge structures exposed to turbulent wind has been hinged on the linearization of self-excited forces (i.e., motion-induced aerodynamic forces), assuming small oscillations of the bridge around a steady state, and often disregarding the effects of turbulence on the bridge aeroelastic behavior.These force models can be either steady or unsteady and are responsible for aeroelastic instabilities, in particular flutter.From a mathematical point of view, flutter is a Hopf bifurcation, marked by a pair of complex conjugate eigenvalues, the real part of which becomes positive when reaching the critical wind velocity.Based on the linear assumption, past research has primarily concentrated on incipient and deterministic flutter of bridges ( [6][7][8][9][10] among others), in the wake of the pioneering works in aeronautics, where this problem has become well-known for airplane wings since the 1920s [11,12].However, experimental observations have shown that the post-critical flutter response tends to evolve towards limit cycles, after either subcritical (see, e.g., [13,14] for a flat plate and a bridge, respectively) or supercritical bifurcations (see, e.g., [15,16] for some bridges).Furthermore, the onset of flutter instability is often significantly influenced by fluctuations in wind velocity due to turbulence, as observed by Diana [17].These insights have led to the emergence of two distinct branches of aerodynamic force models that aim at including such nonlinearities.This work specifically focuses on the nonlinearity resulting from turbulence.Synoptic turbulent wind is typically represented as a superposition of a mean velocity and fluctuating components across three directions, often modeled as broad-band multivariate ergodic random processes.Nakamura and Ozono [18] made a distinction between the effects of small-scale turbulence, whose size is comparable to the thickness of the shear layers, and large-scale turbulence, whose size is larger than the body dimensions.Small scales may have an impact on the general bridge section aerodynamics, affecting shear layer separation (see also [19,20]), while large scales can be seen by the body as an instantaneous change in kinetic pressure and wind slope.The latter may produce a random variation over time in the fluid-structure interaction parameters, particularly when self-excited forces exhibit a strong sensitivity to the mean angle of attack (e.g.[21][22][23][24]), resulting in time-varying dynamic systems, which can be susceptible to instabilities such as random flutter (e.g.[25,26]).
Since the early 1980s, there has been a focused investigation into the parametric effects caused by kinetic pressure variations due to longitudinal turbulence, starting with the seminal work of Lin and Ariaratnam in 1980 [27].Their work delved into stability against torsional random flutter of bridges, leading to the assessment of statistical moment boundaries.In their model, turbulence was represented as a white-noise randomness in wind velocity magnitude, assuming the bridge response as a Markov vector.This simplification aids in the stability analysis of a diffusion process described by an Itô stochastic differential equation (through Stratonovich or Wong and Zakai conversion rules [28]).Later, expanding the analysis to aerodynamically coupled two-degree-of-freedom systems (including vertical bending and torsion) and modeling the turbulence-induced noise as an Ornstein-Uhlenbeck process, moment stability was extensively studied.In particular, Bucher and Lin [29,30] and Lin and Li [31] highlighted the energy transfer between modes due to noise in the self-excited forces, which, in the considered cases, was beneficial for system stability.However, these studies also demonstrated that turbulence can both stabilize and destabilize structural modes, depending on the aerodynamics of the bridge section (see also [32]).In the early 2000s, Poirel and Price [25,33] assessed the random flutter stability of a twodegree-of-freedom airfoil by numerically calculating the largest Lyapunov exponent, which is a measure of sample stability.Later, they explored how stochastic stability, in terms of phenomenological Hopf bifurcation (P-bifurcation, [34,35]), is affected by turbulence characteristics such as intensity and integral length scale, along with structural parameters and an additional stiffness nonlinearity.They found that an increase in turbulence noise in the self-excited forces tends to bring forward the onset of sample stability and modify the shape of the joint probability density function (PDF) of the airfoil pitching response at the limit cycle oscillation.
Subsequent research has emphasized the crucial role of the time-varying angle of attack induced by atmospheric turbulence.Indeed, large-scale turbulence can slowly (relatively to the structural time scale) alter the wind incidence, thereby affecting the aerodynamics of the bridge section and, consequently, the self-excited forces.Examples of studies exploring such nonlinear effects are the works by Diana et al. [21], Chen and Kareem [36], and Barni et al. [22,37].These studies differ in the model used to translate the slowly-varying turbulence-induced modulation of unsteady self-excited forces from the frequency domain, where experimental aerodynamic derivatives are usually available, to the time domain.Specifically, rheological models are employed in [21,23,38] and rational functions in [22,36,37].The parametric effects associated with the angle of attack were found to be more significant for bridge sections than those mentioned earlier related to the fluctuations in the kinetic pressure [39].Notably, the model presented and experimentally validated by Barni et al. [22,37] relies on a linear time-variant framework.Indeed, it accounts for the variation in the angle of attack due to turbulence through the coefficients of the rational approximation of the selfexcited force transfer function.In this case, the angle of attack also alters the terms related to fluid memory in the convolution integrals, introducing an additional layer of complexity compared to the aerodynamic force model used by Poirel and Price in their studies on airfoils [25,[33][34][35].Building on this model, Barni and Mannini [39] provided additional insight into the stabilizing and destabilizing roles of large-scale turbulence, comprehensively addressing the parametric excitation related to time-variations in kinetic pressure, angle of attack and reduced velocity.They found that the turbulenceinduced stabilizing or destabilizing effects primarily arise from the so-called average parametric effect associated with the variations in the angle of attack.This effect introduces time-invariant aerodynamic damping and stiffness, which are not affected by the spatial correlation of turbulence.The study also discusses the role played by parametric resonances in the first torsional mode and by resonances of combination type between vertical bending and torsional modes, which are however generally minor in most turbulence scenarios, especially up to a certain value of turbulence intensity.
Finally, this analysis clearly highlights how such time-variant linear systems markedly exhibit some classical features typical of nonlinear systems.
Despite the availability for over 20 years of models that account for the effects of turbulence-induced angle of attack on bridge aeroelasticity, methods to formally evaluate random flutter stability have not been explored in this context yet.This gap likely arises from the additional complexity in self-excited force models introduced by a time-variant angle of attack, which makes the determination of statistical moment stability an intricate task, in contrast to the relatively straightforward cases involving only the variations in the wind kinetic pressure.This work represents a first endeavor of an extensive Lyapunov stability analysis, incorporating both the largest Lyapunov exponent (LE) and moment Lyapunov exponents (MLE), through Monte Carlo methods for a suspension bridge exhibiting a nonlinear aeroelastic behavior due to changes in the angle of attack induced by turbulence.MLE comprehensively characterize the stability of a noisy system, as they can detect situations where, though almost surely stable, an unacceptable transient response can occur due to the asymptotic divergence of some statistical moments.In different contexts, these mathematical tools have already been applied to study the stability of coupled viscoelastic systems under parametric noise excitation, but examples of such MLE applications remain scarce in the literature.Among these, remarkable contributions are those by Doyle and Sri Namachchivaya [40] and Sri Namachchivaya et al. [41], who examined the p-th moment stability through moment Lyapunov exponents (MLE) in a two-degree-of-freedom system with stiffness coupling, parametrically excited by a small noise characterized by a realistic spectrum.Liu et al. [42] employed MLE to study the stochastic stability of coupled elastic systems with non-viscous damping, driven by a white noise.Additionally, MLE were calculated by Li and Liu [43] for the airfoil model previously considered in [34].However, the stochastic stability of coupled viscoelastic systems with timevariant memory, such as a suspension bridge under turbulent wind, has never been rigorously studied with these methods.
In this work, a simplified equivalent 2D model of the Hardanger Bridge, in Norway, considering its three sectional degrees of freedom (lateral and vertical displacements, and torsional rotation), is assumed as case study.The variations in reduced velocity and kinetic pressure due to turbulence are neglected, as it was shown in [39] that they have a minimal impact on the Hardanger Bridge stability.Sample and p-th moment stability are discussed under various turbulent wind scenarios, modeling the flow velocity fluctuations as Ornstein-Uhlenbeck processes.Finally, for the sake of physical understanding of stochastic stability results, an analysis of Floquet exponents is conducted, focusing on a simplified model of turbulence parametric excitation.

Background: stochastic state-space model
The self-excited forces acting on a bridge deck are typically expressed using the mixed time-frequency formulation, which employs aerodynamic derivatives as suggested by Scanlan in [6].This linearized model is valid for small-amplitude harmonic motion and incorporates the unsteadiness resulting from fluid memory.Therefore, by referencing the force system illustrated in Fig. 1, the self-excited force vector can be expressed as follows: Fig. 1: Sketch of the Hardanger Bridge section with the reference system for displacements, forces and wind velocities where q se = q y q z q θ T , r = y z θ Here V m is the mean wind velocity, ρ is the air density, B is the bridge cross-section width, and r = [y z θ] T represents the vector of deck displacements (see Fig. 1).P j , H j and A j with j ∈ {1, 2, . . ., 6} are the dimensionless aerodynamic derivatives, which depend on the reduced frequency of oscillation K = ωB/V m , where ω is the circular frequency of oscillation, and on the mean angle of attack.These coefficients are usually measured through wind tunnel tests.q y , q z and q θ are the self-excited drag, lift and moment.Alternatively, for any motion, Eq. ( 1) can be transferred in the frequency domain and expressed solely as a function of K (e.g., [44]).
Aiming at a time-domain formulation of self-excited forces that accounts for the effect of large-scale turbulence, the 2D rational function approximation (RFA) model, proposed by Barni et al. [22], postulates that the transfer function between self-excited forces and bridge motion is modulated over time in a quasi-steady fashion by the angle of attack resulting from low-frequency wind velocity fluctuations.Although these forces are assumed to be linear with respect to the bridge motion, time-varying parameters are introduced as a function of the angle of attack.This marks a significant shift from conventional models and implies nonlinear features in the bridge dynamic behavior and stability, as will be highlighted later in this work.It is important to note that the turbulent wind field exhibits not only temporal but also spatial variability (multivariate process), leading to random fluctuations of the angle of attack along the length of the bridge girder.This spatial variation is not accounted for in the simplified 2D bridge model considered here (refer to Fig. 1); thus, compared to [37], the degrees of freedom in modeling both the structural behavior and the self-excited forces are significantly reduced, along with the computational effort of the analysis.However, this choice is justified by the findings in [39], which demonstrate that the partial correlation of turbulent fluctuations has a negligible impact on self-excited forces and bridge stability.Conversely, the assumed perfect correlation of the buffeting load (external forces) leads to a strong overestimation of the dynamic response to turbulent wind.
However, as will be discussed later, this does not pose any issues in terms of stability.
According to the 2D RFA model, the time-domain self-excited forces read: where In this model, the angle of attack is exclusively associated with turbulence, and therefore given by where u(t) and w(t) are the longitudinal and vertical turbulent wind velocity components, respectively.In contrast, the contributions from bridge torsional motion are considered negligible under stable conditions [36,37,45].
For a three-degree-of-freedom dynamic system, the equations of motion can be written in a matrix form as follows: where M, C, and K denote the structural mass, damping, and stiffness matrices, respectively (see Section 4.1).In contrast, Ĉae and Kae represent the memoryless components of aerodynamic damping and stiffness that directly stem from the RFA approximation.q b is an external load vector.Therefore, applying the state-space transformation γ 1 = r, γ 2 = ṙ and γ l+2 = Ψ l , l ∈ {1, 2, . . ., N − 2} to Eqs. ( 3) and ( 4), the following system is obtained: where I ∈ R 3×3 is the identity matrix.Eq. ( 5) can be expressed in compact form as follows: where is the time-variant state matrix of the system, and the input matrix.The system described by Eq. ( 6) is linear time-variant, implying that the stability of the bridge does not depend on initial conditions and external forces (e.g., buffeting forces q b ).Consequently, the mean wind velocity and the parametric effect of turbulence associated with the angle of attack emerges as the sole external factor impacting system stability.Notably, as state at the beginning of this section, this aspect is traditionally viewed as a nonlinear feature, despite the system linearity in terms of state variables.

Lyapunov exponents
The flutter stability threshold of a time-invariant system can be determined by the sign of the real part of the eigenvalues of the state matrix.Nevertheless, for a timevariant stochastic system like the one of Eq. ( 6), these eigenvalues fluctuate over time, requiring more complex stability considerations, such as sample or p-th moment stability (see [46]).Pivotal in the analysis of sample stability of stochastic systems perturbed by ergodic processes is the concept of Lyapunov exponents, which serve as a measure of the system response sensitivity to a perturbation.Indeed, the Multiplicative Ergodic Theorem [47] establishes the existence of the Lyapunov exponents, which are deterministic numbers even if the system is stochastic.For an n-dimensional system, there exist n Lyapunov exponents, represented as λ 1 , . . ., λ n ∈ R, which delineate the average exponential rates at which the n random subspaces asymptotically either expand or contract.It is a well-established fact that after a long-enough time almost all sample trajectories will expand or contract in the direction indicated by the largest LE, the sign of which determines the sample or almost sure stability of the system.
The largest LE associated with the homogeneous part of Eq. ( 6) is defines as: A negative largest Lyapunov exponent indicates that the system is stable with probability one, while a positive exponent means that sooner or later the system response will diverge.
LE can be derived using analytical methods, such as the well-known Khasminskii procedure [48], applicable to simple dynamic systems described by linear Itô stochastic differential equations.However, the problem specified in Eq. ( 6) presents considerable complexity for analytical approaches, involving multiple degrees of freedom, aerodynamic coupling, and a complex nonlinear relationship between the state matrix and the stochastic angle of attack.Consequently, a well-established numerical method is employed to calculate the largest LE, following the algorithm described in [49].In this case, the numerical estimator for the largest LE is defined as: where E[•] denotes the expected value, K∆t = T is the observation time, which can be extended as necessary to ensure the convergence of the estimator.The term x k denotes the solution at the k-th timestep, normalized with its Euclidean norm, so that x k resides on the unit sphere S 3[2+(N −2)]−1 .Eq. ( 8) is based on the fact that the norm of the state vector at the K-th step is calculated as the product of the norms of all preceding steps: Such a normalization is instrumental to avoid data overflow or underflow when the system is strongly unstable or stable, respectively.

Moment Lyapunov exponents
Even if the solution of Eq. ( 6) is almost surely stable when λ γ < 0, i.e., ∥γ(t)∥ → 0 as t → ∞ with probability one at the exponential rate λ γ , it is still possible for the system to exhibit other kinds of instabilities, such as those related to the statistical moments of the state vector norm.To visualize this behavior, Fig. 2 provides a preview of a case study that will be elaborated upon in Section 5. Fig. 2a shows the torsional response of the bridge that decays from unit-sphere initial conditions with λ γ < 0 for V m = 52 m/s and grows with λ γ > 0 for 64 m/s (λ γ = 0 for V m = 61 m/s).In contrast, Fig. 2b reports three samples of the bridge response, all of which decay to zero as t approaches infinity, but two of them exhibit a large transient response.Therefore, since almostsure convergence does not necessarily imply moment convergence, moment Lyapunov exponents (MLE) are essential for a comprehensive description of the stability of the stochastic system in Eq. ( 6).These exponents are defined as follows: It is demonstrated in [50] that the slope of the MLE function Λ γ (p) at p = 0, Λ ′ γ (0), is equal to the largest LE λ γ .
To the best of the authors' knowledge, there are only two numerical methods for the evaluation of the MLE using a Monte Carlo approach, both proposed by Xie [51,52].Although some studies in the literature have applied the first method [51] ( to determine MLE in the stability analysis of nonlinear aeroelastic systems described by stochastic differential equations (see, e.g., [53,54]), Xie himself in [52] pointed out some issues in estimating MLE with this approach.First, estimating the p-th statistical moment of the norm of the state vector as a sample mean, the variance of the MLE estimator is equal to the variance of the population divided by the sample size.
Consequently, when the abovementioned statistical moments are unstable, the variance of ∥γ∥ p significantly increases with time, and an acceptable reduction of the variance of the MLE estimates may require an unreasonably large number of samples.Moreover, due to the finite length of floating-point representation of numbers in computers, once the solution of the equations attains very large values, the contribution of all other samples disappears.Therefore, a correct estimation of MLE becomes impossible when there is a tendency to rare, extremely large, system response (i.e., when the statistical moments of the norm of the state vector tend to get unstable).These issues cannot partially be circumvented by reducing the observation time (like, for instance, the example in Fig. 2 of [52]) because atmospheric turbulence noise presents time scales much longer than those characterizing the associated dynamic system.
To circumvent these issues, this study utilizes the method proposed by Xie and Huang [52] for numerically estimating the MLE in linear homogeneous stochastic systems.Since the abovementioned problems arise from large values of ∥γ∥ p , the statistics of log∥γ∥ are leveraged to estimate the MLE.If log∥γ∥ is normally distributed for t → ∞, a good estimate of the MLE can then be determined from the relation below: Here, once again T = K∆t is the observation time.For a linear homogeneous system parametrically excited by stationary ergodic diffusion processes, the normal distribution of log∥γ∥ and the validity of Eq. ( 10) is guaranteed by some conditions, but these are difficult to verify for a complex problem such as that governed by Eq. ( 6).However, as noted by Arnold [46], such conditions are commonly satisfied in engineering systems, and Xie and Huang [52] suggest that the asymptotic normality of logarithm of norm can also numerically verified a posteriori (this will be done for our system in Section 5.2).The solution γ(T ) of Eq. ( 6) is determined through the normalization procedure outlined in [52] (Section 4, Step 3) to prevent numerical underflow or overflow.The MATLAB➤ codes employed for calculating both LE and MLE for a linear stochastic differential equation subject to random noise are included as supplementary material.
4 Case study

Hardanger Bridge
The previously outlined model is applied to evaluate the aeroelastic stability of the The finite element model of the bridge, reported in [37], was used to perform a modal analysis after the application of dead loads to include the nonlinear stiffness effect of the cables.To simplify the problem and the computational cost of the Monte Carlo analysis, the bridge is then modeled as an equivalent two-dimensional dynamic system with three degrees of freedom.The mass matrix M is composed by the equivalent masses, obtained by dividing the modal masses by the integral of the squared modal displacements along the bridge deck.The three components of the vector r signify the mode contributions to the deck response at the points of maximum modal displacement.The modes selected for this analysis are the first lateral, vertical and torsional symmetric modes (modes 1, 4 and 15 reported in [37]).The stiffness and damping matrices K and C in Eq. ( 5) are deduced from the bridge modal properties.
Fig. 3: Lift and moment aerodynamic derivatives associated with the torsional motion for the Hardanger Bridge (see also [22]).Results of wind tunnel measurements are reported as a function of the reduced wind velocity V r for various mean angles of attack.2D RFA approximations are also shown.
A structural modal damping ratio of 0.5% is assumed for all modes.

Turbulent wind model
The characteristics of wind turbulence at the Hardanger Bridge site and deck level are provided in terms of Kaimal spectra: where fi = f Li Vm , i = u, w is the nondimensional frequency, based on the mean wind velocity V m , the integral length scale of turbulence L i and the generic frequency f (in Hz).σ 2 i is the variance of the associated process, in wind engineering often expressed in terms of turbulence intensity Finally, according to the wind field measurements in [55], an upward mean wind velocity inclination of 2 to 2.5 deg relative to the horizontal plane is typically observed at the Hardanger Bridge site.Therefore, α m = 2.5 deg is assumed as the mean value of α(t).
In the present work, in order to be in the conditions to calculate the MLE according to Eq. ( 10), the longitudinal and vertical turbulent wind velocity components, which define the angle of attack α(t), are modeled as Ornstein-Uhlenbeck (OU) processes.
Then, the time evolution of the turbulent velocity components is governed by the following equations: W u and W w represent two scalar Wiener noises with unit variance.G i,1 ≥ 0 (inverse of a time-scale) and G i,2 ≥ 0 (diffusion), where i ∈ {u, w}, are the coefficients of the linear filter.When the processes get stationary (theoretically for t → ∞, in fact G i,1 t ≫ 1 is enough), the associated one-sided power spectral densities read: where ) is again the process variance (for t → ∞).
The coefficients G i,1 and G i,2 are determined by pursuing a match between the OU spectra and the given Kaimal spectra.First, the same total variance σ 2 i is imposed by setting G 2 i,2 = 2 G i,1 σ 2 i .Then, the parameters G i,1 are obtained through a nonlinear least-squares fitting procedure, as in several previous works (see, e.g., [32,56,57]).It is worth noting that the spectrum of an OU process mainly differs from the Kaimal spectrum in the fact that the latter decays for high frequencies with a power of −5/3 while the former with a power of −2, thus introducing non-negligible differences in the energy distribution among the frequencies, as shown in Fig. 4a.Nevertheless, this still represents a reasonable model of atmospheric turbulence (see Dryden turbulence model [58]).In this context, however, a characteristic nondimensional frequency, Gθ w,1 = G w,1 / (2πf θ ), is preferred to the integral length scale to characterize the various turbulent wind scenarios.This nondimensional frequency gives a measure of the energy distribution of the noise across the frequencies in relation to the torsional mode frequency f θ , which is a pivotal parameter for system stability, as will be shown later.
On the other hand, the vertical turbulence intensity I w is a key indicator of the total variance of noise fluctuations (considering that generally α ≃ w/V m ).
Turbulent wind velocity time histories are generated by numerically integrating the stochastic differential equations, Eqs. ( 12) and ( 13), using the Euler scheme and a time step of 0.001 s.Fig. 4b showcases the simulated power spectral densities of the vertical wind velocity component compared to the corresponding target OU spectra.
Clearly, decreasing the nondimensional frequency Gθ Fig. 5: RMS of the buffeting torsional response of the Hardanger Bridge, obtained using either an LTI approach (gray circles) or the formulation of Eq. ( 2) for the selfexcited forces (red crosses).The solid and dash-dot red lines represent the sample and 2nd-moment stability thresholds, respectively.Finally, the blue triangles indicate the LTI eq stability thresholds [39].
Finally, one may remark that, without any additional low-pass filter, the angle of attack partially contravenes the slowly-varying assumption of the model leading to Eq. ( 2) [22].Indeed, the 2D RFA model has experimentally been validated for variations in the angle of attack up to a frequency equal to two-thirds of the bridge motion frequency (in the experiments, a limitation arose because the setup did not allow the investigation of higher ratios; see [59]).Nonetheless, this work assumes that the model validity can be extended to higher ratios, considering that the majority of turbulence energy is at low frequency, where the assumption of slow variation in the angle of attack holds.Indeed, the contribution of turbulent fluctuations at frequencies higher than the first bridge mode frequencies was found not to significantly influence the flutter stability and the dynamic response of the Hardanger Bridge [39].

Results
Prior to examining the system stability using Lyapunov exponents and moment Lyapunov exponents, the bridge dynamic response to the external buffeting load vector q b is calculated for various turbulent wind scenarios according to the buffeting load model in [37].Fig. 5 presents the results in terms of root mean square (RMS).For the sake of brevity, only the torsional response is reported here, as it is generally the most representative for flutter analyses (e.g., for the Hardanger Bridge, see [37]).This figure aims to provide immediate and practical insight into the correlation between the subsequent stability analyses and the actual dynamic response of the system under external load excitation.RMS values are obtained averaging the results for 20 onehour-long time histories calculated at a sampling rate of 12 Hz.Each frame of Fig. 5 includes the bridge dynamic response calculated according to Eq. ( 6) and the associated non-parametrically excited case (linear time-invariant, LTI), where the angle of attack due to turbulence is considered constant at 2.5 deg (as per Eq. ( 2  w,1 = 1.28 to Gθ w,1 = 0.64.Barni and Mannini [39] showed that this destabilizing effect of turbulence is mostly ascribed to the so-called "average parametric effect", with parametric resonances playing only a secondary role.This effect arises from the fact that the aerodynamic derivatives averaged over the time history of the angle of attack α(t) typically differ from those corresponding to the mean angle of attack α m .This discrepancy results in average aerodynamic stiffness and damping that significantly diverge from the values used in the classical time-invariant approach.This phenomenon becomes pronounced if the aerodynamic derivatives exhibit substantial nonlinear and non-symmetric variations around α m , and it intensifies for increasing turbulence intensity.However, this effect is independent of the spectral distribution of parametric excitation energy.In this regard, the blue triangular markers in Fig. 5 indicate the deterministic flutter stability threshold when using averaged aerodynamic derivatives (LTI eq approach, see [39]), as opposed to the values at α m = 2.5 deg (LTI case).The remarkable contribution of the average parametric effect in varying bridge stability is apparent.The gap between the LTI eq thresholds and the vertical asymptotes of the nonlinear buffeting response (red crosses) hints at additional effects, likely driven by parametric resonances.Their impact can be twofold; indeed, they can both stabilize the bridge, as noted for Gθ w,1 = 1.28 (turbulent energy at higher frequencies), and destabilize it, as for Gθ w,1 = 0.43 (turbulent energy at lower frequencies).

Sample stability
The largest Lyapunov exponents, as well as the MLE in the next section, are determined for a system undergoing free-decay vibrations, starting from unit-norm initial conditions (∥γ∥= 1).To ascertain the sample or almost sure stability threshold of the bridge in turbulent flow, the LE is calculated using Eq. ( 8).Fig. 6 shows the convergence of the numerical estimator as a function of the observation time T .For the sake of conciseness, the figure is based on a representative set of turbulence parameters and mean wind velocity, but similar trends were observed in other cases.Given that the LE is a numerical estimator and infinite time histories are impractical as per the LE definition, we consider 1000 samples of the LE and examine the estimated PDF every T G w,1 /(2π) = 500.As the system behaves linearly with respect to the state variables, the magnitude of the initial conditions does not impact the stability of the system.Fig. 6 illustrates how the dispersion of the estimator decreases with the extension of the observation time and the mean becomes stable after approximately T G w,1 /(2π) = 4500.Consequently, a time window of about 21600 s corresponding to T G w,1 /(2π) = 5000 is selected as the standard for subsequent LE calculations.Fig. 7 shows how the largest LE varies with mean wind velocity under various turbulent wind conditions.Additionally, the figure considers the associated LTI system, which is useful for interpreting the results.In this case, the Lyapunov exponents coincide with the real part of the eigenvalues of the state matrix and, up to a certain wind Fig. 7: Largest LE as a function of the mean wind velocity V m for various turbulent wind conditions (colored thick lines with markers).The Lyapunov exponents associated with the LTI system are also displayed (thick gray line) together with those calculated for the system using average aerodynamic derivatives (LTI eq approach; dashed lines for the largest LE, dotted lines for the other Lyapunov exponents; the colors indicate again the associated turbulent wind scenario).velocity, the largest LE is that associated with the horizontal mode, which exhibits lower damping.However, for high wind velocity, due to aerodynamic interaction with the vertical mode, the damping in the torsional mode progressively reduces, so that the associated LE becomes dominant at around 63 m/s and flutter instability occurs at 71.2 m/s.The average parametric effect on the Lyapunov exponents is isolated in the figure reporting the pattern of the real part of the state-matrix eigenvalues associated with lateral and torsional modes, calculated according to the LTI eq approach [39].While the impact on the lateral mode branch is negligible (due to the essential insensitivity of the aerodynamic derivative P * 1 to the average parametric effect), the torsional LE branch rises up at significantly lower wind velocity compared to the classical LTI system, leading to an earlier bridge instability.This can mainly be ascribed to the increase in turbulent flow in the important aerodynamic derivative A * 2 , which can even become positive (which means a negative aerodynamic damping contribution in torsion) for high turbulence intensity (see Fig. 10 in [39]).
When the parametric excitation induced by the fluctuating angle of attack α(t) is fully accounted for in the calculation of the largest LE (Eq.( 8)), the largest LE pattern exhibits a significantly lower gradient with the mean wind velocity and a smoother transition from the horizontal branch to the torsional branch.Importantly, the results demonstrate a sensitivity to the nondimensional characteristic frequency of turbulence, especially between Gθ w,1 = 1.28 and Gθ w,1 = 0.64, which cannot be appreciated by considering the average parametric effect only (the latter is clearly independent of the frequency distribution of turbulent fluctuations).Such phenomenon, likely related to parametric resonances (Section 5.3 will further elucidate this aspect), postpones the change of sign of the largest LE, thereby enhancing bridge stability for Gθ w,1 = 1.28.
In this case, although not as crucial for bridge stability, a damping reduction in the lateral mode is also promoted by turbulence (upward shift of the largest LE pattern for V m ≲ 40 m/s).As will be better explained in Section 5.3, this is attributed to the energy transferred from the torsional mode to the lateral mode promoted by the parametric excitation induced by α(t) [30,39].
The sample stability thresholds are then marked with vertical red solid lines in Fig. 5.These thresholds comply with the pattern of the torsional dynamic response of the bridge, representing the ideal vertical asymptotes of the curves.However, in some turbulence scenarios, particularly for Gθ w,1 = 1.28 and I w > 7.5%, the increase in the RMS of the torsional response is less sudden, and very large values of the response are attained before the sample stability threshold.In such cases, the latter does not seem conservative for design purposes.

p-th moment stability
As explained in Section 3.2, the algorithm proposed by Xie and Huang [52] is utilized in this study for the calculation of MLE.This method requires that log∥γ∥ is normally distributed, a condition that is satisfied for the current problem, as exemplified by the histograms in Fig. 8a and Fig. 8b.Moreover, Fig. 9 shows the results of a study to establish the number of samples and the observation time T (simulation length) needed for reaching convergence in the calculation of MLE.5000 samples and a time T = 7200 s were judged a conservative choice and used in all of the following analyses.
For a specific turbulence scenario and two wind velocities, Fig. 8 compares the MLE curves obtained with the two algorithms discussed in Section 3.2 and here denoted as Xie (2005) [51] and Xie (2009) [52].While in both cases the slope in the origin of the MLE curve, Λ ′ γ (0), coincides with the values of LE reported in Fig. 7, the two MLE estimators already diverge for small values of p. Notably, the pattern obtained with Xie (2005) algorithm is qualitatively very similar to the example reported in [52] for a completely different problem and large values of T (Section 2, Fig. 2 thereinto).The reasons for this discrepancy and the clearly wrong pattern obtained with Xie (2005) algorithm have already been explained in Section 3.2 and, as expected, the more unstable is the system and the higher is the order p of the statistical moment considered, the larger is the difference between the two results.It is worth adding here that no improvement is obtained by significantly increasing the number of samples; moreover, unlike different physical systems, here the long time scale of the noise does not allow fixing the problem by reducing the observation time T (see Fig. 9), thereby making absolutely necessary the use of Xie (2009) algorithm.Fig. 8: Comparison between the MLE estimators according to the two algorithm described in [51] (black line) and [52] (red line) for 5000 samples and an observation time T = 7200 s.The turbulence scenario considered is that with I w = 7.5% and Gθ w,1 = 0.64.The PDF of log∥γ(T )∥ is also reported with its fitted normal distribution.

Further insight via Floquet exponents
Previous analysis underscored specific results that cannot fully be explained by the average parametric effect and that, in fact, depend on the frequency distribution of noise variance.To shed some light on these parametric effects, a simplified system parametrically excited by a sinusoidal gust is studied through Floquet theory, in a similar way to [39].Specifically, a sinusoidal variation in the angle of attack α(t) is assumed in Eq. ( 6), parametrized by frequency f * (the so-called pumping frequency), amplitude α 0 , and mean value α m = 2.5 deg.The resulting system is then time-periodic.
Obviously, the impact on bridge stability of broad-band turbulence parametric excitation cannot be simply analyzed as a superposition of effects of individual harmonics composing the random process, akin to Fourier analysis.In fact, the effects of these harmonics combine in a nonlinear manner when forming the eigenvalues of the system transition matrix (see, e.g., [60]).However, the simplified time-periodic system is expected to provide some qualitative insight into the behavior of its noisy counterpart.
Floquet approach requires the analysis of one period of bridge response, 1/f * (for more details, see, e.g., [60]).The asymptotic behavior of the system is ruled by Floquet multipliers ζ j (for j = 1, . . ., n, being n = 3[2 + (N − 2)] the size of the state vector  modes.This phenomenon, called parametric coupling resonance by Barni and Mannini [39], is responsible for the bowl-shaped unstable regions, which expands with an increase in mean wind velocity.
At high mean wind velocity, the stability maps in Fig. 12 are significantly influenced by the average parametric effect, making it challenging to extract insight into the contribution of parametric resonances.For this reason, Fig. 13 presents stability maps for the same bridge structure but excluding the average parametric effect on self-excited forces.This is done by linearly approximating with the 2D RFA the experimental aerodynamic derivatives with respect to the angle of attack.This new rational approximation also slightly alters the system behavior at the mean angle of attack α m = 2.5 deg and consequently the system LTI stability threshold (V LTI cr = 66.8 m/s instead of 71.2 m/s).However, due to the nonlinearity of the problem, it is important to note that a rigorous distinction between parametric resonances and average parametric effect is impractical with numerical approaches.Indeed, this analysis only seeks to provide a qualitative understanding of the transition in the LE pattern from that dominated by average aerodynamic derivatives, as per the LTI eq approach, to that only ruled by parametric resonances.In this case, the parametric resonance of difference type associated with the nondimensional frequency fθ − fz becomes the most significant one for system destabilization, growing wider and the attendant LE getting larger when the mean wind velocity increases.However, for wind speeds exceeding approximately 55 m/s, a broad and stabilizing parametric resonance, centered around the sum fθ + fz , also starts to be visible.It can be recognized in the map as a dark blue region, which becomes a real stability island once the mean wind velocity exceeds the LTI flutter stability threshold (V LTI cr = 66.8 m/s).This result also represents a clear evidence of the marked nonlinear behavior of the system, as the same range of pumping frequencies has a destabilizing impact in presence of the average parametric effect (see the panel for V m = 50 m/s in Fig. 12).
These stability maps enable us to speculate about some unexplained features of the noisy system.To aid this analysis, Fig. 14 presents the dimensional spectra for the vertical turbulent wind velocity, which is useful for localizing the energy of the parametric excitation across the nondimensional frequencies f = f /f θ .For a mean wind velocity of 60 m/s, about one third of the total turbulent energy is found at f < 1 for Gθ w,1 = 1.28, compared to slightly less and more that two thirds for Gθ w,1 = 0.64 and 0.43, respectively.The energy distribution slightly changes with the mean wind velocity, but this does not significantly alter the substance of the discussion for V m in the range 40 to 75 m/s.Referring to the map presented in Fig. 13 and trying to qualitatively extend the results to broad-band turbulence, one may observe that for a nondimensional characteristic spectral frequency Gθ w,1 = 1.28, average parametric effect aside (which is independent of the spectral distribution of fluctuation energy), the parametric excitation mainly contributes to stabilization.In contrast, for Gθ w,1 = 0.64 and 0.43, a significant portion of the energy is located in the zone of the destabilizing parametric resonance fθ − fz .Although this is only a conjecture, it qualitatively explains the behavior of the dynamic response in Fig. 5 and that of the largest Lyapunov exponent in Fig. 7, which show that the system parametrically excited by broad-band random noise is generally more stable for Gθ w,1 = 1.28 than for Gθ w,1 = 0.64 and 0.43.

Conclusions
This work represents the first application of largest Lyapunov exponents and moment Lyapunov exponents to assess the random flutter stability of a suspension bridge subjected to the parametric excitation due to turbulence-induced angle of attack.The following main conclusions can be highlighted.
❼ For the Hardanger Bridge, the turbulence-induced parametric excitation anticipates flutter instability compared to the time-invariant condition, confirming previous findings for this case study.
❼ The random flutter critical velocity determined by LE closely aligns with the bridge dynamic response pattern and with low-order p-th moment stability thresholds for most turbulent wind conditions considered.In these cases, the flutter stability threshold is unequivocal.However, certain turbulence scenarios reveal that sample stability is not conservative, and it is necessary to resort to MLE for a complete description of stability.In these circumstances, the stability index shows a reduced slope with the mean wind velocity, evoking some additional consideration of which statistical moments should be considered for the engineering assessment of flutter stability.
❼ Bridge random flutter stability is often predominantly influenced by the average parametric effect (always destabilizing for the current case study), but additive and subtractive parametric resonances between torsional and vertical modes can also play an important role for some turbulent wind conditions.❼ Floquet analysis, though it is not able to fully address the problem of random flutter induced by parametric variations in the angle of attack, aids in providing a qualitative understanding of the parametric effects governing the stability of the bridge system.The behavior of LE and MLE in relation to the frequency-distribution of parametric excitation energy is consistent with the findings of Floquet stability maps, which revealed the contribution of additive and subtractive parametric resonances to system stabilization or destabilization.
❼ The behavior of Floquet exponents also underscores the complexity of the dynamic system considered.Specifically, the nature of the parametric resonance of additive kind between torsional and vertical modes hinges on system average damping, exhibiting a destabilizing role in the presence of an average parametric effect and a stabilizing role when it is absent.
This work shows a viable strategy to determine the flutter stability threshold for long-span suspension bridges, considering turbulence-induced parametric excitation.This approach could also be applied to other time-varying self-excited force formulations, as long as the dynamic system without external loads can be modeled as linear.Furthermore, the nonlinear features exhibited by these time-variant systems are intriguing, and their parametric resonance mechanisms deserve further investigation in the future, possibly employing simplified analytical methods to enhance understanding of such complex behaviors.

Declarations
where N − 2 is the number of additional aeroelastic states, and d l , l ∈ {1, 2, . . ., N − 2}, represent the rational function coefficients, which are polynomials (ζ and ι are the polynomial degrees) of the time-variant angle of attack α = α(t).Ψ l ∈ R 3 , l ∈ {1, 2, . . ., N − 2}, represent the additional aeroelastic states used to account for fluid memory effects, which are defined through a time-varying convolution:

Fig. 2 :
Fig.2: Torsional dynamic bridge response due to unit-sphere initial conditions.The response refers to a turbulent characterized by I w = 12.5% and a non-dimensional frequency Gθ w,1 = 1.28

Hardanger
Bridge, subjected to turbulent wind.This suspension bridge, crossing the Norwegian Hardanger Fjord, is characterized by a main span of 1310 m and towers reaching a height of 186 m, establishing it as Norway's longest and Europe's thirdlongest bridge.It also features a fairly streamlined single-box steel deck with a height of 3.2 m and a width of 18.3 m.


The aerodynamic derivatives are necessary to identify the 2D RFA model.As detailed in[22], these aerodynamic coefficients were measured in the wind tunnel through forced-vibration tests for 11 different mean angles of attack, ranging from −8 to +8 deg, and for reduced velocities (V r = V m /(f B), where f denotes the vibration frequency) up to 55.The rational functions associated with the self-excited force model include N = 4 terms, which comprise two additional aeroelastic states, and adopt 5th-order polynomials for all A i (i = 1, . . ., N ) and d l (l = 1, . . ., N − 2) parameters appearing in Eq. (2).Consequently, a total of 228 coefficients are determined, ensuring a high-fidelity fit to the experimental data, as exemplified in Fig.3.If the angle of attack exceeds the available experimental interval [−8, +8], the 2D RFA coefficients are assumed to remain constant, equal to the values at either −8 or +8 deg.Two critical observations emerge from Fig. 3. Firstly, the aerodynamic coefficients exhibit a pronounced dependence on the mean angle of attack, indicating the potential for a significant time-variant behavior of the dynamic system.Secondly, the crucial coefficient A * 2 , intimately connected to the aerodynamic damping in torsion, registers positive values (indicative of negative aerodynamic damping contributions) for mean angles of attack approximately exceeding 5 deg.This fact underscores the importance of incorporating such nonlinearity in the self-excited force model.
8 and δ w = 9.4 are the spectral parameters suggested for the specific bridge site[55].Longitudinal and vertical wind velocity fluctuations, u and w, respectively, are assumed as statistically independent.To examine the bridge stability behavior under different turbulent wind conditions, a range of turbulence intensities and integral length scales are considered.Specifically, the longitudinal turbulence intensity I u and integral length scale L u are let vary between 10% and 25%, and between 100 m and 300 m, respectively.These intervals reflect the common atmospheric turbulence parameters for long-span suspension bridges like the Hardanger Bridge.The vertical turbulence intensity and integral length scale are set to be half and 10% of the longitudinal values, respectively ([55]).

w, 1 Fig. 4 :
Fig. 4: Nondimensional power spectral density of turbulent vertical wind velocity fluctuations.(a) Reference turbulence spectrum (thick gray line) compared with the OU spectrum (blu line).(b) Simulated spectra (colored lines), averaged over 5000 time series, compared with the corresponding target OU spectra (thick gray lines) for various nondimensional frequencies Gθ w,1 .f y , f z and f θ indicate the lateral, vertical and torsional frequencies of bridge vibration modes.
) with α ≡ 2.5 deg).In the LTI scenario, the deterministic flutter bifurcation occurs at a critical wind velocity V LTI cr = 71.2m/s, and the corresponding vibration frequency is f cr = 0.289 Hz.It is important to remark once again that these RMS values are somewhat exaggerated and not realistic due to the perfect correlation of the external load, a limitation inherent to the simplification of a 2D bridge model.Nevertheless, they are valuable for visualizing the instability onset.

Fig. 5
Fig.5clearly shows that the flutter stability of the Hardanger Bridge is reduced due to the random parametric effects of turbulence.The destabilizing impact increases with the energy of the parametric excitation, namely the turbulence intensity.Additionally, an effect of the energy distribution across the frequencies is observable, particularly passing from Gθ w,1 = 1.28 to Gθ w,1 = 0.64.Barni and Mannini[39] showed that this

Fig. 6 :
Fig. 6: PDF of the LE estimator, calculated for 1000 samples and reported for various nondimensional observation times spaced T G w,1 /(2π) = 500 apart.The turbulence parameters are I w = 7.5% and Gθ w,1 = 0.64, and the mean wind velocity is 54 m/s.µ and σ denote the mean and the standard deviation of the estimator.

Fig. 10 illustrates
Fig.10illustrates the patterns of MLE for all the previously considered turbulent wind scenarios and mean wind velocities close to, and slightly below, the sample stability threshold.p-values in the range −2 ≤ p ≤ 5 are considered.MLE exhibit fast variations with the mean wind velocity under most turbulent conditions, making low-order p-th moments unstable for just a slight decrease in the mean wind velocity with respect to the critical sample value.To visualize the 2-nd moment stability with respect to the dynamic response curves reported in Fig.5, a vertical dash-dotted line is added at the mean wind velocity for which the associated MLE becomes positive.This threshold refers to the asymptotic unbounded growth of the norm of the state vector, which, in the case of bridge flutter, is clearly tightly connected to the variance of vertical and above all torsional response of the bridge.It is apparent that 2-nd moment stability complies with the evolution of the RMS of the bridge torsional response better than sample stability in those cases in which in the bridge dynamic response increases in a smoother way at high wind velocity.The calculated MLE are also employed to determine the stability index, which corresponds to the non-trivial zero of the MLE curve and which is shown in Fig.11 as a

(Fig. 9 :
Fig. 9: Numerical convergence of MLE with respect to the observation time T and number of samples S. The MLE in (a) refers to 5000 samples, while (b) considers a simulation length of 7200 s.In both cases turbulence is charcterized by I w = 7.5% and Gθ w,1 = 0.64.

Fig. 10 :
Fig. 10: MLE patterns close to the sample stability threshold for various turbulent wind conditions.The labels next to the curves represent the mean wind velocity V m (in m/s).

Fig. 11 :
Fig. 11: Stability index as a function of the mean wind velocity in the proximity of the sample stability threshold, for various turbulent wind conditions.

Fig. 12 :Fig. 13 :
Fig. 12: Stability chart of the bridge dynamic system (normalized largest LE) for different amplitudes and nondimensional frequencies f * = f * /f θ of the sinusoidal variation in the angle of attack modulating the self-excited forces.The gray line highlights the stability boundary, while V LTI cr represents the time-invariant flutter threshold.

Fig. 14 :
Fig. 14: Power spectral density of the vertical wind velocity for different nondimensional characteristic frequencies of turbulence, Gθ w,1 (I w = 7.5% and V m = 60 m/s).The abscissa is the normalized frequency f = f /f θ .