Time-domain simulations of turbulence effects on the aerodynamic flutter of long-span bridges

Aerodynamic flutter instability has been a major concern for long-span flexible bridges, such as suspension and cable-stayed bridges, subjected to wind actions that result in the so-called self-excited forces. Though turbulence effects on bridge flutter have been studied in the last few decades, its true effects remain a debate due to the limitation of previous wind tunnel facilities, such as using turbulence scales that are too small in these experiments. In this paper, the characterizations of self-excited forces are presented in both the frequency-domain and in the time-domain. Then, the flutter analysis is conducted under both smooth flow and turbulent flow in order to investigate the effect of wind turbulence on the flutter instability. The effect of wind turbulence is directly modeled in the time-domain in order to avoid the complicated random parametric excitation analysis of the equation of motion used in previous studies. By comparing the results of different turbulence intensities with that of the smooth flow, it is found that the turbulence has a stabilizing effect on bridge flutter. The turbulence can change the vibration patterns and weaken the spatial vibration correlation to some extent. As a result, the critical flutter velocity can be increased by 5% to 10% over that under smooth flow.


Introduction
The aerodynamic instability of long-span bridges under strong wind excitations has been a major concern in recent years, especially with the continuously increasing span length built around the world, especially in China. The well-known failure of the old Tacoma Narrows Bridge has heightened people's attention to wind-resistant design for long-span bridges and led to many research studies and investigations on bridge flutter performance (Zhang et al. 2017;Caracoglia 2018;Tang et al. 2019;Gao et al. 2020). At first, flutter analysis was based on the thin airfoil theory given by Theodorsen and Mutchler (1935). Afterwards, Scanlan and his co-authors (Scanlan and Tomo 1971;Scanlan and Jones 1990) developed the formulations of the lift, drag, and pitching moment motion-dependent forces, otherwise known as the self-excited forces, which are presently widely used. These equations involve the flutter derivatives obtained from experimental measurements on sectional bridge deck models. There are general two approaches for flutter analysis: the frequency-domain approach and the time-domain approach. In the frequency-domain method, the critical flutter velocity, the flutter mode shape, and the corresponding frequency can be obtained by conducting a complex eigenvalue analysis (Agar 1989;Miyata and Yamada 1990;Jain et al. 1996;Dung et al. 1998;Ge and Tanaka 2000;Liu et al. 2018;Han et al. 2015). In the time-domain method, the self-excited forces are represented in the form of the indicial functions (Scanlan et al. 1974;Zhang et al. 2010;Caracoglia and Jones 2013;Zhang et al. 2019) or rational functions (Chen et al. 2000;Chowdhury and Sarkar 2005), which can be identified through experimental tests or numerical approaches from available flutter derivatives.
While a few studies have investigated the effect of wind turbulence on bridge aerodynamic instability (Lin and Ariaratnam 1980;Bucher and Lin 1988a, 1988b, 1989Lin and Li 1993), for most of the previous work, the wind turbulence has been neglected in the flutter analysis. However, based on the experiments carried out in wind tunnels on full aeroelastic models of long-span suspended bridges, it is shown that the level of wind turbulence generated in wind tunnels influences the aerodynamic instability of the structure (Diana et al. 1993). Based on their observations on the measured flutter derivatives, Scanlan and Jones (Scanlan and Jones 1991;Scanlan 1997) found that the flutter instability performance can be enhanced in a turbulent wind field because the turbulence may weaken the inherent correlation of the self-excited forces along the bridge deck. At the same time, to theoretically investigate the wind turbulence effect on bridge flutter instability, complicated random parametric excitation analyses were conducted since the equation of motion becomes a randomly parametrically excited type of equation. The results showed that wind turbulence with high intensity might also have adverse effects that reduce bridge flutter instability Lin 1988a, 1988b;Cai et al. 1999). These results also raise the question: Can the influence of wind turbulence on flutter instability be fully reflected in the measured flutter derivatives? Huston (1986) has conducted a series of studies on the effects of large-scale turbulence on the flutter derivatives and found that the turbulence does not always stabilize bridge flutter. According to Li and Lin (1995), the presence of wind turbulence changes the combined structure-fluid critical mode and results in a new energy balance. They suggested it is the random deviation from the deterministic flutter mode that renders either the stabilizing or destabilizing effect of turbulence possible.
From these studies of turbulence effects on bridge flutter performance mentioned above, it can be found that the majority of studies rely on wind tunnel tests due to the inherent complexity of wind-bridge interactions, which makes their mathematical formulations extremely difficult. However, the accuracy of the results of wind tunnel tests largely depends on the matching degree of the turbulent atmospheric flows, which is influenced by many factors, such as the Reynolds numbers, the integral scales, the turbulence intensities, and the anisotropy, etc. Typically, only a fraction of the turbulence characteristics can be matched in a wind tunnel experiment (Haan 2000). Numerical analysis of bridge flutter performance in turbulent flow is relatively rare, and it has always been considered as a supplement of experimental simulation. These numerical simulation approaches, such as the random parametric excitation (RPE) analysis (Lin 1979;Cai et al. 1999), are usually too complex mathematically and/or too computationally consuming. In the present paper, a simplified numerical approach is proposed in which the influence of wind turbulence on bridge flutter instability is investigated numerically in the time-domain, which not only avoids the complicated stochastic solution process, but also can technically include any nonlinear effects (geometric and/or material nonlinear) in the analysis when deemed necessary. Taking the Karman spectrum as the target spectrum, turbulent wind fields with different turbulence intensities are simulated and utilized here. For comparison, three approaches: (i) the frequency-domain approach based on flutter derivatives, (ii) the time-domain approach based on rational functions under smooth flow, and (iii) the time-domain approach under turbulent flow, haven been applied to predict the critical flutter velocity. These results are compared and discussed based on the analysis of a prototype long-span bridge.

Description of bridge and dynamic characteristics
The Taihong Bridge analyzed in the present study is a single span suspension bridge that has a span length of 808 m, a streamlined steel box girder with a width of 37.5 m, and a height of 3 m. The two main cable planes are 33.6 m apart, and the bridge deck is suspended by hangers at intervals of 12 m. The two bridge towers are reinforced concrete structures with a height of 112.7 m and 107.6 m, respectively. A sketch of the bridge is shown below in Fig. 1.
The commercial finite element software ANSYS is used here to establish the 3D finite element (FE) model of the Taihong Bridge. In the finite element model, the main girder and the towers are simulated by Beam4 elements, and the main cables and suspension cables are simulated by Link10 elements. Figure 2 shows the finite element model of the Taihong Bridge.
The dynamic properties of the Taihong Bridge, including its natural vibration frequencies and mode shapes are analyzed based on the dead load deformed configuration. The results are shown in Table 1. The vertical bending and torsional vibrations are usually the critical modes of flutter for suspension bridges. It can be found that the

Identification of flutter derivatives using forced vibration method in time domain
The flutter derivatives represent the characteristics of the self-excited aerodynamic forces of the bridge section and are the most important parameters in the flutter instability analysis of long-span bridges. Many studies have focused on the identification of the flutter derivatives (Yamada and Ichikawa 1992; Sarkar et al. 1994;Gu et al. 2000;   Chen et al. 2002). Self-excited lift force L se and pitching moment M se per unit length are defined as (Scanlan and Tomo 1971) in which h is the vertical or heaving displacement; α is the torsional or pitching displacement;ḣ is the vertical or heaving velocity;α is the torsional or pitching velocity; ρ is the air density; U is the mean wind velocity; B is the bridge deck width; is the reduced circular frequency; ω h and ω α are the circular frequencies of the heaving and pitching motions, respectively; and H * i and A * i (i = 1,2, …,4) are the flutter derivatives (H * 1 , H * 4 , A * 2 and A * 3 are called direct flutter derivatives while H * 2 , H * 3 , A * 1 and A * 4 are called cross flutter derivatives). By combining the coefficients of flutter derivatives, Eqs. (1a) and (1b) can be expressed as: in which H i and A i (i = 1,2, …,4) are equivalent flutter derivatives whose relationships with the flutter derivatives are given in matrix forms as: According to Eq. (2), H i and A i (i = 1,2, …,4) can be obtained from the self-excited force time-histories and vibration time-histories by the principle of least square fitting. Then H * i and A * i (i = 1,2, …,4) can be obtained according to Eq. (3) based on the identified H i and A i (i = 1,2, …,4). The time-domain forced vibration method can be divided into single DOF and coupled vibration identification methods according to the forced vibration mode of the model. In this paper, 1-DOF harmonic oscillation of the tested model is simulated and the single DOF identification method is used. More specific details can be found in previous studies (Han et al. 2014).
The results of the 8 identified flutter derivatives under different wind attack angles are shown in Fig. 4, though 18 flutter derivatives are used in the general formulation described later.

Flutter instability analysis
The equation for motion of a bridge under smooth flow can be expressed as: where M, C, and K are the global mass, damping, and stiffness matrices, respectively; q, q, and € q represent the nodal displacement, velocity, and acceleration vectors, respectively; and F se denotes the vector of the nodal aero-elastic forces. Self-excited lift force, L se , drag force, D se , and pitching moment, M se , per unit length are defined in a general form with 18 flutter derivatives as (Scanlan 1978;Jain et al. 1996) in which H Ã i , P Ã i , and A Ã i (i =1 to 6) are the aerodynamic flutter derivatives related to the vertical, lateral, and torsional directions, respectively; h, p, and α are the vertical, lateral, and torsional displacements of the bridge deck, respectively; and the dot on the cap denotes the derivative with respect to time. An iterative process is usually used to define the oscillation frequency present in K and the aerodynamic flutter derivatives H Ã i , P Ã i , and A Ã i Approach I-Frequency-Domain Approach using Smooth Flow Flutter Derivatives.
The 3D finite element flutter analysis is performed to determine the critical flutter wind speed and flutter modes of the Taihong Bridge under different attack angles of wind flows. In the analysis, a pair of user-defined Matrix27 elements are used to simulate the aerodynamic forces acting on each element of the main girder. One Matrix27 element is used to simulate the aerodynamic stiffness, and the other one is used to simulate the aerodynamic damping (Han et al. 2015). By solving the equation of motion after assembling the aerodynamic stiffness and damping matrices, the flutter instability and critical wind velocity of the system can be determined from the values of the real part of the complex eigenvalues. The system is dynamically stable if the real part of all eigenvalues is negative and dynamically unstable if the real part of one or more eigenvalues is positive. At a certain wind velocity, the real part becomes zero, which means the system is on the critical state and the corresponding wind velocity is defined as the critical flutter wind velocity.
In this study, damped complex eigenvalue analysis is carried out under wind velocities ranging from 0 to 150 m/s of different wind attack angles by assuming that the structural damping ratio is 0.5%. The complex eigenvalues of multiple modes under a wind attack angle of − 3°are shown in Fig. 5 as an example. It is found that the critical flutter mode is mode 16, i.e., the first symmetric torsion mode. The deformation shape of mode 16 under the critical flutter velocity is shown in Fig. 6, which can be found to be the coupling of the symmetric vertical bending and the symmetric torsion. The critical flutter wind velocity and corresponding frequencies under different wind attack angles are summarized in Table 2. To compare the frequency-domain and time-domain approaches, the self-excited forces per unit span can also be expressed in terms of impulse response function as follows (Lin and Yang 1983): where I denotes the impulse function of the self-excited forces, and the subscripts  According to the equivalency of the spectral characteristics between the aerodynamic self-excited forces expressed by the impulse response function and the self-excited forces defined with the flutter derivatives by Scanlan (1978) as shown in Eq. (6), the relationship between the impulse response function and the flutter derivatives can be obtained as follows: where the I is the Fourier transform of I, and the subscripts denote the corresponding force components.
The Roger's approximation, a kind of rational function approximation approach, is often utilized to express I fx ð f ¼ L; D; M; x ¼ h; p; αÞ (Chen et al. 2000). Taking the impulse function of the lift force induced by the vertical motion as an example, I L seh can be expressed as: where d l and A n (l = 1 to m; n = 1 to m + 3) are frequency independent coefficients, which can be determined through parameter fitting of the flutter derivatives obtained by the wind tunnel test. The value of m is user-defined, which determines the approximation accuracy. Among the rational function, the third term of (i.e., the A 3 term) represents the additional aerodynamic mass and is normally negligible. It should be noted that for each force component, the coefficients of d l and A n (l = 1 to m; n = 1 to m + 3) are different, and m = 2 is used in the present study. The nonlinear leastsquares method is used here to determine these coefficients in Eq. (8). Then, to validate the procedure, the flutter derivatives are back calculated based on Eq. (7) and called simulation values here. Figure 7 shows the comparison of experimental values and simulation values of the flutter derivatives of the Taihong Bridge under a wind attack angle of − 3°.
By taking the Fourier transform of (8) after obtaining the corresponding parameters, the impulse response function can be expressed as: By substituting Eq. (9) into Eq. (6a), the aerodynamic lift force induced by the vertical motion can be obtained as follows: Similarly, the expressions for the other self-excited force components can be obtained. After obtaining the expressions of all the aerodynamic self-excited forces, flutter where F(t i ) represents the equivalent nodal forces induced by external forces at time t i , and R(t i ) denotes the equivalent nodal resistance of the structure. At time t i + Δt, the equation of motion becomes as follows: The Newmark-β method is used to solve this equation of motion considering the geometry nonlinearity of the structure and the nonlinearity of aerodynamic loads. The nodal acceleration and velocity can be expressed as follows: According to Eq. (13) and Eq. (14), the acceleration and velocity of the structure at a given time can be obtained through iterations. By increasing the wind speed gradually and conducting transient dynamic analysis of the structure, the flutter critical wind velocity can be determined. Figures 8, 9, and 10 show the time histories of the displacements of the mid-span under a wind attack angle of 0°at pre-critical, critical, and post-critical stages, respectively. The vertical and torsional displacements are defined at the center of the bridge deck. It can be found that when the wind velocity is 77 m/s, 79 m/s, and 80 m/s, the displacement amplitudes are convergent, becoming constant, and changing to divergent, respectively. Therefore, 79 m/s is identified as the onset wind velocity of bridge flutter. Figure 11 shows the spectrum of displacements at the critical wind velocity, and it is found that the vibration is a vertical-torsional coupled one with the dominant frequency being around 0.448 Hz for the torsional vibration (Fig. 11a) and 0.10 Hz and 0.446 Hz for the vertical vibration (Fig. 11b). It is noted that the first symmetric torsion mode (Mode 16) has a natural frequency of 0.451 Hz and this frequency reduces slightly when it approaches the flutter stage as shown in Fig. 5. Meanwhile, the 1st anti-symmetric vertical bending mode (Mode 3) has a frequency of 0.173 Hz and it approaches 0.1 Hz at the occurrence of flutter. The flutter in the time-domain analysis shows a coupled vibration involving a few modes similar to that in the frequencydomain. A few other modes might be also involved since they have a frequency near the neighborhood of 0.1 Hz (see Table 1). The critical flutter velocities under different wind attack angles and corresponding dominant frequencies are summarized in Table 3. It can be found that the results of the time-domain analysis are consistently slightly lower than that from the frequency-domain analysis.

Approach III -time-domain approach under turbulent flow
Since flutter derivatives are a direct description of the self-excited forces, a number of studies have been conducted to investigate the influence of turbulence on flutter derivatives. Diana et al. (1992) found that turbulence had little effect on flutter derivatives by measuring the flutter derivatives in a turbulent wind field from the full-scale test. Sarkar et al. (1994) also concluded that turbulence did not significantly affect the flutter derivatives under smooth flow. Therefore, it seems like turbulence does not significantly influence the self-excited forces through the flutter derivatives directly. Due to the lack of experimentally obtained flutter derivatives under fully "matched" turbulent flow, the effect of wind turbulence on the flutter derivatives might not fully reflect its effects on the bridge flutter performance. Therefore, in the following analysis, the flutter derivatives measured in the smooth flow instead of in turbulent flow are used in the formulations of the self-excited forces, i.e., they remain the same as in Approach II. In Fig. 8 Time-histories of motion at mid-span (pre-critical, U = 77 m/s) order to investigate the effect of wind turbulence on the critical flutter velocity, a turbulent wind velocity component is added on the mean wind velocity of Approach II in the formulations of the self-excited forces. This is a similar approach to that used by Lin (1988a, 1989) in treating the self-excited forces in turbulent flow. However, Lin (1988a, 1989) used a complicated stochastic solution process for the parametrically excited differential equations of motion. In comparison, the present approach is to numerically solve the nonlinear equations of motion in the timedomain.
Taking the aerodynamic lift force induced by vertical motion as an example, it can be expressed as: in which u(t) denotes the turbulent component of wind velocity, which can be obtained by field measurements or numerical simulations. Due to the lack of field measured wind data, the harmonic synthesis method is used here to generate a sample of the turbulent wind velocity time-history of the bridge main beam. Four turbulent wind fields, with turbulence intensities ranging from 5% to 20%, are simulated here to investigate their effects on flutter instability. The main parameters of the turbulent wind field are listed in Table 4. Figure 12 shows the turbulent  wind velocity time-history (I u = 5%) at the mid-span and the 1/4-span of the bridge. Figure 13 shows the simulated turbulent wind spectrum, and it can be found that the simulated turbulent wind spectrum agrees well with the target spectrum.
The results of flutter analysis under different turbulent flows, compared with the results under smooth flow, are shown in Fig. 14. Specific critical flutter velocities under different wind attack angles are summarized in Table 5. It can be found that wind turbulences do have stabilizing effects on bridge flutter. In general, the critical flutter velocity increases monotonically with the increase of turbulence intensity from 5% to 20%. It should also be noted that with relatively low turbulence intensities (such as typical 5% and 10%), the critical flutter velocities are almost the same with that under the smooth flow, and the stabilizing effect is not obvious. With relatively high turbulence intensities (such as 15% and 20%), the critical flutter velocity can be increased by 5% to 10% due to the turbulence effect. The stabilizing effect is mainly from the decorrelation effect of the turbulence on the aerodynamic forces as discussed below. Figure 15 shows the post-critical time histories of the displacements at the mid-span under turbulent flow (I u = 15%). It can be found that the motion of the bridge deck becomes much more irregular and random in the turbulent flow compared to that in the smooth flow. Hence, turbulence not only increases the critical flutter velocity, but also  changes the vibration patterns to some extent. Figure 16 shows the comparison of the correlation of the motion at mid-span and 1/4-span under the smooth flow and turbulent flow. The correlation coefficient of the vertical displacements at these two points under both the smooth flow and the turbulent flow are calculated in Table 6. It can be found that, under the smooth flow, the vibration is highly synchronized at each location of the bridge deck, which is more likely to generate stronger wind-bridge interactions. However, this consistent pace in vibration has been broken in the turbulent flow. This    kind of inconsistency in motion affects the motion-dependent flutter forces and contributes to the increase of the critical flutter velocity.

Conclusions
Although the frequency-domain approach and the time-domain approach are theoretically equivalent since the rational functions are extracted from the experimentally obtained flutter derivatives, it can be found that there are still some differences in the numerical values of the critical flutter velocities by comparing the results of Approach I and Approach II. The differences may be caused by the following reasons. Firstly, the numerical identification of the parameters in the rational functions introduces a degree of approximation and results in numerical errors due to the highly irregular behavior of several flutter derivatives. Secondly, the frequency-domain approach is based on the linear-elasticity theory and decomposes the structural response into multiple main participation modes by using the linear modal decomposition technique, while the timedomain approach considers the nonlinearities of the geometry and aerodynamic loads. The frequency-domain method is much more straightforward since it directly relies on experimental data and is less expensive from a computational point of view. However, the time-domain method is more flexible and powerful in some bridge analysis, such as the nonlinearity analysis, coupled flutter, and buffeting analysis, etc. By comparing the critical flutter velocities of Approach II and Approach III, it can be found that the turbulence in the cross wind can raise the critical flutter velocity. The increase in amplitude is dependent on the turbulence intensity. It is also found that the turbulence influences the vibration patterns and the spatial vibration correlation. However, it should be noted that the effect of turbulence on flutter derivatives is not considered in this study. Many studies have been conducted to investigate the flutter derivatives under the turbulent flow, and different research conclusions have been made. It is still a debate whether the effect of wind turbulence on flutter instability can be fully reflected by the measured flutter derivatives. Therefore, presently, the mechanism of the influence of turbulence on flutter derivatives is still not clear, and further study is needed.
The effect of the turbulent wind field on the aerodynamic instability of bridges is a complex process, which needs to take into account its effect on the flutter derivatives, the effect on the spatial correlation of self-excited forces, the effect of turbulenceinduced buffeting on the aerodynamic instability, and other factors comprehensively. In the present study, the turbulence effects on flutter instability are numerically simulated in the time domain, which considers the nonlinear effects and spatial correlations. The effect of buffeting vibrations on flutter instability is under study in this developed numerical framework and will be reported later. Future studies, involving both wind tunnel experiments and numerical analysis, are needed to further advance the understanding of the effects of turbulence on the aerodynamic instability of bridges.