Preliminary Flutter Stability Assessment of the Double-Deck George Washington Bridge

Featured Application

This study provides a simplified approach to the flutter analysis of suspension bridges having two superposed decks. The George Washington Bridge engineering case is analyzed, in consideration of its historical relevance and age. The obtained results are compared with the predictions of simplified formulations available from other sources.

Abstract

We deal with the flutter analysis of the George Washington bridge, in both the single- and double-deck configurations of 1931 and 1962, respectively. The influence of the additional lower deck on the aerodynamic behavior is investigated. To overcome the lack of aerodynamic data, a simplified approach is followed based on Fung’s formulation, in which the flutter derivatives are expressed in terms of the real and imaginary parts of the Theodorsen function and of the steady-state aerodynamic coefficients of the deck cross-section. The latter are obtained by Computational Fluid Dynamics simulations conducted in ANSYS FLUENT, whereas the ANSYS Mechanical APDL finite element package is used to perform the flutter analyses. Two different methods for the application of the aeroelastic forces are employed for the double-deck configuration: (i) self-excited forces, based on flutter derivatives related to the whole cross-section, acting on the upper deck; and (ii) self-excited forces, based on flutter derivatives related to the single deck, simultaneously applied to the upper and lower decks. The obtained results are critically compared with theoretical predictions of simple formulas available from the literature; it is suggested that laboratory tests are needed since no experimental results seem to be available.

1. Introduction

Old bridges (conventionally, which are more than 50 years old), require special maintenance, accurate inspections, and a careful assessment of their safety conditions.

The conflict between economy and structural performance between the 19th and the 20th centuries led the design of long-span bridges to the development of very flexible and slender structures. The use of the elastic deflection theory allowed for very slender decks against static loads and shifted the design trend at that time from rigid trusses to slender edge girders. This evolution ended brutally with the Tacoma Narrows Bridge disaster due to the wind-induced flutter instability on 7 November 1940. From then on, any design of a flexible structure must assure the structure itself to be stable under the dynamic effects of wind loads. In fact, in addition to the already known static divergence due to the steady-state wind loads, flutter stability has become a governing criterion in the design of long-span suspension bridges. The objective of a (linear) flutter analysis is to predict the lowest critical flutter wind velocity and the corresponding flutter frequency.

A crucial point in flutter analysis is the definition of motion-dependent aerodynamic loads. The first analytical solutions were developed by Wagner [1] and Theodorsen [2] for thin airfoil. Theodorsen defined the self-excited forces as the superposition of circulatory and non-circulatory contributions, the former depending on the oscillation frequency and accounting for flow unsteady effects, and the latter independent of oscillation frequency and including inertial effects due to the moving fluid mass. Subsequently, Scanlan and Tomko exported some features of the Theodorsen’s results extending the formulation to bluff bodies, as bridge cross-sections. In this formulation, the wind loads induced by sectional harmonic motions are described by means of a linear format based on experimentally evaluated parameters, called “flutter derivatives”, that supply the lack of closed-form analytical formulations [3,4]. The relation between Wagner’s approach based on indicial functions, Theodorsen’s theory based on circulatory function and Scanlan’s flutter derivatives is well detailed in [5].

In the last decades, Scanlan’s formulation has been the most widely adopted, and the calculation of flutter derivatives has become an important step for any flutter analysis. Currently, the only method considered reliable for their calculation at the design stage is that of carrying out experimental tests on scale models in wind tunnels. Some notable examples are [6,7] for the Great Belt East Bridge and [8] for the Akashi Kaikyo Bridge; more recently [9] for the Jianghai Channel Bridge and [10] for the Hardanger Bridge.

To overcome costs and difficulties arising from wind tunnel testing, several efforts were made with the aim of developing some simplified methods for the calculation of flutter derivatives. Although these methods are not appropriate for the final design stage, they allow useful and versatile studies in early design phases. Some authors utilized the derivatives of the thin airfoil, as [11] for the Humber bridge. This simplification leads to relatively small errors when the bridges are characterized by streamlined cross-sections. Al-Assaf [12] adopted an alternative approach for open-truss stiffened suspension bridges: he estimated the aerodynamic derivatives based on the correlation between the thin airfoil derivatives and those of other bridges having a similar deck configuration. The method was applied to evaluate the flutter stability of the second Tacoma Narrows Bridge, with special focus on the effects of the side grates. Another simplified approach is given by the quasi-steady theory, in which the aerodynamic loads depend on the instantaneous relative velocity between the flow and the cross-section. In this framework the flow-induced forces are described by means of non-linear static relationships involving the wind angle of attack and the displacements of the structure. A linearization of this model allows a comparison with Scanlan’s semi-empirical approach and furnishes an expression of the flutter derivatives as functions of the steady-state aerodynamic coefficients of the cross-section [13]. These simplified expressions of the flutter derivatives are widespread for the streamlined deck cross-sections. Some works on the validity of this approach are [14,15,16] and some applications for the determination of the critical state can be found in [17,18,19]. In a recent work [20], a linear superposition of flat plate aerodynamics was adopted to estimate the flutter derivatives of streamlined multi-box deck sections. In that case, correction factors were introduced into the proposed analytical formulas to better fit the available experimental data. In the present paper, a simplified approach is used based on Fung’s formulation [21], in which the flutter derivatives are expressed in terms of the steady-state aerodynamic coefficients and of the real and imaginary parts of the Theodorsen function.

Two-dimensional Computational Fluid Dynamics (CFD) simulations were conducted in ANSYS FLUENT for the determination of the lift and moment coefficients, varying the wind attack angle. In recent years, the CFD approach gained importance with respect to the traditional experimental investigation. Some applications related to the bridge cross-sections are [22,23]. More recently, Brusiani et al. [24] pointed out that Reynolds-averaged Navier-Stokes (RANS) approach coupled with the k−ω SST turbulence model, identify the best compromise between accuracy and computational cost. The CFD framework also provides different approximate techniques for the direct calculation of flutter derivatives, some applications are: [25,26,27,28]. Nevertheless, the calculation of the steady-state aerodynamic coefficients requires simpler simulations compared to those required to compute the flutter derivatives, especially in terms of computational effort. In addition, the approximations introduced by a 2D geometrical representation, neglecting several appendages that could affect the aerodynamic behavior, do not justify the use of more refined methods.

Concerning flutter analyses, several methods can be found in literature [29,30], a method developed by Hua and Chen [31], which allows the analysis to be carried out with the FE package ANSYS Mechanical APDL, is adopted in this paper. The software, in fact, provides specific user-defined elements, through which it is possible to implement the motion-dependent aeroelastic forces as expressed in Scanlan’s formulation.

The George Washington suspension bridge was chosen as a case study, partly because of its historical significance. In fact, it was the first bridge whose main span exceeded one kilometer. The bridge was opened to traffic in 1931. It is characterized by a total length of 1450 m and a mid-span of 1067 m. Initially it was composed only of the upper level, whereas the lower deck was constructed from 1958 to 1962 because of the increasing traffic flow. Both configurations will be analyzed in the following: the original one with the single deck and the stiffened one with the two decks (Figure 1). The choice is also motivated by the peculiarity of the current bridge having two superposed decks. A similar configuration can be detected in the Verrazano-Narrows Bridge (USA), also having a great historical relevance, and in the Yangsigang Yangtze River bridge (China), which currently has the third longest span in the world after the Akashi Kaikyo bridge (Japan) and the recent 1915 Çanakkale bridge (Turkey), the current World record.

As a matter of fact, several investigations on the aerodynamic performances of trussed girders have been made, but only a few for double-deck bridges [32]. Due to the lack of aerodynamic data and studies on the subject, two simplified ways for the application of the aeroelastic forces are followed in the present study with the aim of obtaining approximate predictions. The first way considers the self-excited forces, based on the flutter derivatives related to the whole cross-section, acting on the upper deck. In the second case, the self-excited forces, based on the flutter derivatives related to the single deck, are simultaneously applied to the upper and lower decks. Both approaches produce plausible results, which can serve as a reference for future comparisons with alternative investigation methods. In principle, the proposed approach can be used for analyzing other double-deck bridges, and deserves further attention for proper validation.

2. Motion Related Wind Load

The equation describing the motion of the bridge in the smooth flow can be expressed as:

$$M\ddot{X}+C\dot{X}+KX=F_{ae}.$$

(1)

where, in a Finite Element framework, $M$, $K$, and $C$ represent the global mass, stiffness, and damping matrices, respectively; $\ddot{X}$, $\dot{X}$, and $X$ are the nodal acceleration, velocity, and displacement vectors, respectively; and $F_{ae}$ denotes the vector of nodal self-excited forces. The three degrees of freedom of the bridge deck cross-section, namely the vertical displacement $h$, the torsional rotation $\alpha$, and the horizontal displacement $p$, to which correspond the aeroelastic forces of lift $L_{SE}$, moment $M_{SE}$, and drag $D_{SE}$, are reported in Figure 2, where $U$ indicates the undisturbed mean wind velocity.

According to Scanlan’s formulation [33], the self-excited aeroelastic forces per unit deck length can be expressed as a linear function of the displacement and velocity parameters of the deck as follows:

$$L_{SE}=\frac{1}{2}\rho U^{2}B\left[K{H}_1\frac{\dot{ h}}{U}+K{H}_2\frac{B\dot{ \alpha }}{U}+K^2{H}_3\alpha +K^2{H}_4\frac{h}{B}+K{H}_5\frac{\dot{ p}}{U}+K^2{H}_6\frac{p}{B} \right],$$

(2a)

$$D_{SE}=\frac{1}{2}\rho U^{2}B\left[K{P}_1\frac{\dot{ p}}{U}+K{P}_2\frac{B\dot{ \alpha }}{U}+K^2{P}_3\alpha +K^2{P}_4\frac{p}{B}+K{P}_5\frac{\dot{ h}}{U}+K^2{P}_6\frac{h}{B} \right],$$

(2b)

$$M_{SE}=\frac{1}{2}\rho U^2{B}^2\left[K{A}_1\frac{\dot{ h}}{U}+K{A}_2\frac{B\dot{ \alpha }}{U}+K^2{A}_3\alpha +K^2{A}_4\frac{h}{B}+K{A}_5\frac{\dot{ p}}{U}+K^2{A}_6\frac{p}{B} \right],$$

(2c)

where $\rho$ is the air density, $K=\frac{\omega B}{U}$ is the reduced circular frequency, and $H_i\left(K\right), P_i\left(K\right), A_i\left(K\right)$ are the flutter derivatives. In most practical applications, $p$,$\dot{ p}$, and $P_i$ ($D_{SE}$) can be ignored.

As mentioned in the Introduction, the most reliable way to calculate flutter derivatives is based on experimental tests made in the wind tunnel. For the bridge chosen as a case study, no data regarding flutter derivatives were found in the literature, therefore a simplified approach was followed for their calculation. It was decided to utilize a formulation developed by Fung [21] for the thin airfoil because it allows one to input the slopes of the steady-state aerodynamic coefficients, and, unlike the quasi-steady approach, it provides a unique formulation for the determination of the derivatives $H_2, A_2$. The expressions of the flutter derivatives are the following:

$$K{H}_1^{\#}\left(K\right)=-C_l^{\prime } \tilde{F}\left(K\right),$$

(3a)

$$K{H}_2^{\#}\left(K\right)=C_l^{\prime } \left[\left(\frac{1}{4}-a\right)\tilde{F}\left(K\right)+\frac{\tilde{G}\left(K\right)}{K}\right]+\frac{C_l^{\prime } }{4},$$

(3b)

$$K^2{H}_3^{\#}\left(K\right)=C_l^{\prime } \left[\tilde{F}\left(K\right)-K\tilde{G}\left(K\right)\left(\frac{1}{4}-a\right)\right]+C_l^{\prime } K^2\frac{a}{4},$$

(3c)

$$K^2{H}_4^{\#}\left(K\right)=C_l^{\prime } K\tilde{G}\left(K\right)+\frac{C_l^{\prime } K^2}{4},$$

(3d)

$$K{A}_1^{\#}\left(K\right)=C_m^{\prime }\tilde{F}\left(K\right),$$

(3e)

$$K{A}_2^{\#}\left(K\right)=-C_m^{\prime }\left[\left(\frac{1}{4}-a\right)\tilde{F}\left(K\right)+\frac{\tilde{G}\left(K\right)}{K}\right]+\frac{C_l^{\prime } }{4}\left(a-\frac{1}{4}\right),$$

(3f)

$$K{A}_3^{\#}\left(K\right)=-C_m^{\prime }\left[\tilde{F}\left(K\right)-K\tilde{G}\left(K\right)\left(\frac{1}{4}-a\right)\right]+C_l^{\prime } K^2\frac{a}{4}+\frac{C_l^{\prime } K^2}{128},$$

(3g)

$$K^2{A}_4^{\#}\left(K\right)=C_m^{\prime } K\tilde{G}\left(K\right)-C_l^{\prime } K^2\frac{a}{4},$$

(3h)

where $A_i^{\#}\left(K\right)$ and $H_i^{\#}\left(K\right)$ are the flutter derivatives for the thin airfoil as expressed by Fung ($K=\omega B/U)$, $\tilde{F}\left(K\right)$ and $\tilde{G}\left(K\right)$ are respectively the real and the imaginary parts of the Theodorsen function, $C_l^{\prime }$ and $C_m^{\prime }$ are the derivatives of the lift and moment coefficients with respect to the wind angle of attack and $a$ is the distance between the shear center and the centroid of the airfoil, normalized with respect to the chord $B$. For low values of reduced frequency, $\tilde{F}\left(K\right)$ tends to one, $\tilde{G}\left(K\right)$ tends to zero and the Fung formulation tends to the quasi-steady one.

These formulas were tested for two bridges of which both steady-state aerodynamic coefficients and flutter derivatives are available: the Great Belt East in Denmark, with a streamlined deck [6,34], and the Akashi Kaikyo in Japan, with a truss girder [35,36]. The derivatives of the aerodynamic coefficients, expressed in $\left[ra{d}^{-1}\right]$, are: $C_l^{\prime }=-4.37$ and $C_m^{\prime }=1.17$ for the Great Belt Bridge; and $C_l^{\prime }=-1.19$ and $C_m^{\prime }=0.3$ for the Akashi Kaikyo Bridge. The curves in Figure 3 and Figure 4 show a comparison between the flutter derivatives: obtained experimentally (solid crossed line), calculated by the Fung formulas (solid line), calculated by the quasi-steady formulation (dashed line) [26], and predicted by the Theodorsen theory for the thin airfoil (dash dotted line). From the comparisons in Figure 3 and Figure 4, it is possible to notice that the Fung formulation leads generally to a good alignment with the experimental curves: the average trend is always respected except for $H_4$, that in some cases can be neglected [12]. A remarkable superposition is obtained for $A_1$ and $A_3$ of both bridges, and also for $H_3$ in the case of the Great Belt Bridge. Good predictions are also obtained for $A_2$ derivative.

As expected, for a truss-stiffened girder, Theodorsen formulation for thin airfoil becomes unable to reliably predict the aerodynamic behavior. Even though the quasi-steady formulation also shows a good agreement, the Fung formulation was chosen over the quasi-steady approach because this latter is better suited for low values of reduced frequency and high values of reduced wind velocities [16]. In fact, the bridge chosen as a case study is relatively stiff and characterized by high values of torsional and vertical eigenfrequencies, therefore flutter instability is expected to occur at a relatively low reduced velocity.

3. Full-Order Flutter Analysis Using ANSYS Mechanical APDL

The main lack of almost all the FE packages commonly used in civil engineering is the impossibility of including motion-dependent aeroelastic loads into the model and then of performing a complex eigenvalue analysis. Hua and Chen [31] proposed a method that allows one to perform a flutter analysis using the commercial finite element package ANSYS Mechanical APDL. The method is based on the definition of the aeroelastic loads by means of a particular user-defined element as shortly illustrated below. According to the method developed by Hua and Chen [31], the motion-dependent effect due to the wind load is considered for each deck element by the element aeroelastic stiffness and damping matrices modeled by the user defined MATRIX27 elements. Following Scanlan’s formulation, the self-excited forces are expressed as functions of the flutter derivatives, as recalled in Section 2.

It is possible to manipulate Equations (2a)–(2c) in order to obtain the equivalent nodal forces acting on the ends of the generic deck element, hence a lumped formulation can be used to derive the element aeroelastic stiffness and damping matrices [30]. Assembling the element matrices into global aeroelastic stiffness ($K_{ae}$) and damping ($C_{ae}$) matrices, the mathematical model of the system integrated with the effect of aeroelasticity is obtained [31]:

$$M\ddot{X}+\left(C-C_{ae}\right)\dot{X}+\left(K-K_{ae}\right)X=0$$

(4)

The global aeroelastic stiffness and damping matrices contain flutter derivatives, which depend on the reduced wind velocity $U_R=\frac{U}{fB}$ ($f=\frac{\omega }{2\pi }$, oscillation frequency). Therefore, the system expressed by Equation (4) is parameterized by wind velocity and vibration frequency. By Equation (4), a complex eigenvalue analysis can be carried out to determine the critical wind velocity and vibration frequency. Being the conjugate pairs of complex eigenvalues ${\lambda }_j={\sigma }_j\pm i{\omega }_j$, and the conjugate pairs of complex eigenvectors ${\varphi }_j=p_j\pm i{q}_j$, the system will be dynamically unstable if the real part of any eigenvalue is positive. Hence, the condition for the onset of flutter instability is stated as follows: at a certain critical wind velocity $U_f$, the system has only one eigenvalue ${\lambda }_f$ with zero real part, and the imaginary part ${\omega }_j$ of the complex eigenvalue ${\lambda }_f$ becomes the flutter frequency. It is necessary to provide the variation of both wind velocity and vibration frequency in the complex eigenvalue analysis, so that a mode-by-mode tracking method can be employed to iteratively search the flutter frequency and the flutter velocity [31].

The motion-dependent effect due to the wind is taken into account for each deck element by the element aeroelastic stiffness and damping matrices. These matrices should be properly compiled to implement the aeroelastic motion-dependent load by means of the flutter derivatives. User defined MATRIX27 element in ANSYS Mechanical APDL can only model either an aeroelastic stiffness matrix or an aeroelastic damping matrix instead of both of them simultaneously, so a pair of MATRIX27 elements must be attached to each node of a generic bridge deck element as illustrated in Figure 5. The MATRIX27 elements $e1$ and $e3$ represent respectively the aeroelastic stiffness and damping of the node $i$, as the MATRIX27 elements $e2$ and $e4$ represent respectively the aeroelastic stiffness and damping of the node $j$.

4. The George Washington Bridge

As mentioned in the Introduction, the two configurations of the George Washington bridge were chosen as case studies (Figure 1). The deck cross-sections of both configurations were modeled in ANSYS FLUENT, where the steady-state aerodynamic coefficients were calculated via CFD analyses. Finally, full bridge models were developed in ANSYS Mechanical APDL to perform structural and flutter analyses. Flutter derivatives characterizing the aeroelastic load of full bridge models were calculated by means of Equation (3), with $C_l^{\prime }$ and $C_m^{\prime }$ resulting from the CFD analyses.

4.1. Steady-State Aerodynamic Coefficients and Flutter Derivatives

Once the Fung formulation was validated, the steady-state aerodynamic coefficients of both the George Washington deck cross-sections were calculated in ANSYS FLUENT. For each section model and different wind attack angles, RANS simulations were performed modeling the turbulence by the $k-\omega SST$ method, following the indications in [24,37]. In order not to interrupt the logical flow of the speech, more information is given in Appendix A. The resulting steady-state aerodynamic coefficients of lift, $C_l$, and moment, $C_m$, are shown in Figure 6, where dashed lines represent the linear regressions. In the case of the double deck (Figure 6b), static coefficients deviate from linearity for higher angles of attack, so the linear regressions are restricted to the intervals $\left[-2°;2°\right]$ and $\left[-1°;1°\right]$ for lift and moment coefficients, respectively. The first derivatives of static coefficients at the origin, calculated as the slopes of linear regressions, are: $C_l^{\prime }=-6.291$, $C_m^{\prime }=1.902$ for the single-deck and $C_l^{\prime }=-11.024$, $C_m^{\prime +}=1.954$ for the double-deck configurations.

Flutter derivatives, calculated accordingly by means of Equations (3a)–(3h), are shown in Figure 7 and Figure 8. These were used to calculate the aerodynamic lift and moment actions to be applied to the bridge section.

For the double-deck section, the lift and moment actions, based on static coefficients of the whole section, are applied to the top deck only, as sketched in Figure 9a. In addition, a second alternative method was adopted, which allows avoiding the use of data in Figure 6b and provides an additional model for comparison and validation of results. This latter method permits the same aeroelastic forces calculated for the single deck to be simultaneously applied to both upper and lower decks, as illustrated in Figure 9b. This procedure assumes that there is no aerodynamic interference between the two decks, whereas the mechanical interference is correctly modeled by the elements composing the truss. Furthermore, the pressure coefficients of the lower deck are assumed to be equal to those of the upper one and the aerodynamic influence of the secondary truss elements is also neglected. Despite the previous simplifications, the advantage of this procedure lies in the greater reliability of the results of CFD analysis for the single deck. Moreover, it can be conjectured that the simplified formulation previously introduced for the calculation of flutter derivatives (see Equations (3a)–(3h)), being based on the aerodynamics of the thin airfoil, is more suitable for the description of the behavior of a single deck rather than of a double deck.

4.2. Finite Element Models

The ANSYS Mechanical APDL numerical models were realized with the geometrical properties found in literature [38,39,40]. The suspension cables were modeled as beam elements (BEAM188), with a circular cross-section of radius 0.573 m. An equivalent elastic modulus of 1.07 × 10¹¹ Pa was adopted in order to take in to account the influence of the side spans. The hangers were modeled using LINK180 elements, that transmit only tensile forces between the deck and the main cables. The floor system was modeled by means of equivalent beams having the properties listed in

Table 1. Geometrical properties of the truss girder of the double-deck George Washington Bridge.

Elements	Area [m2]	Bending Moments of Inertia [m4]
		Horizontal Axis	Vertical Axis
Equivalent longitudinal upper beam	0.4465	Iy = 0.1867	Iz = 38.6601
Upper chord	0.0549	Iy = 0.0049	Iz = 0.0060
Equivalent transversal upper beam	0.237	Ix = 0.9498	Iz = 3.8102
V-shaped up-winds	0.0374	I1 = 0.0016	I2 = 0.0006
Equivalent longitudinal lower beam	0.3439	Iy = 0.1108	Iz = 29.3241
Lower chord	0.0484	Iy = 0.0029	Iz = 0.0051
Equivalent transversal lower beam	0.0816	Ix = 0.1559	Iz = 1.3645

, wherein the V-shaped up-winds and the lower elements refer to the double-deck configuration. The following material properties were assigned to each element except the suspension cables: Young’s modulus = 2.1 × 10¹¹ Pa, Poisson’s ratio = 0.3 and mass density = 7860 kg/m³. Point mass elements (MASS21) were attached to the upper chords and to the equivalent longitudinal beams, both lower and upper, in order to model the inertial contribution of the sidewalk slabs, roadway slabs and other secondary elements. According to Dana et al. [39], a total dead load of 417 kN/m and 569 kN/m were obtained for the single-deck and the double-deck configurations, respectively. Lastly, the MATRIX27 elements were attached to each node of the equivalent longitudinal beams to model the aeroelastic forces in terms of aeroelastic stiffness and damping.

According to the two different methods introduced in the previous section for the double-deck configuration, two different models were realized. In the first model, the pairs of MATRIX27 elements, compiled with the flutter derivatives in Figure 8, were attached only to the upper deck. Whereas, according to the second method, MATRIX27 elements, compiled with the flutter derivatives plotted in Figure 7, were attached to both upper and lower decks. The finite element models are shown in Figure 10.

5. Flutter Analyses

Before introducing the aeroelastic nodal forces by MATRIX27 elements, a modal analysis with no wind was performed: the results are summarized in Figure 11 and Figure 12, where the labels V, T, L stand for vertical, torsional and lateral, respectively; and the labels A and S stand for antisymmetric and symmetric, respectively. Since it is well known that flutter instability usually involves mainly the first modes at lower frequencies [31], only the first six modes were extracted for both configurations. The geometric nonlinearity of the structure was taken into account running the modal analysis in a preloaded configuration accounting for the gravity loads. Once eigenfrequencies and eigenmodes have been evaluated, the flutter analysis was performed by the method described in [31]. The aeroelastic loads were modeled as nodal forces affecting the stiffness and damping matrices of the element to which the nodes belong, so MATRIX27 elements were incorporated into the structural model to perform the damped eigenvalue analyses.

Finally, the iterative procedure for the determination of flutter speed and frequency was carried out. The previous computational steps were followed for the eigenmodes collected in Figure 11 and Figure 12. The damped complex eigenvalue analyses were conducted for the model under different wind velocities. The step increment of wind velocity was set variable, from a maximum value of 10 m/s in the ranges far from the flutter instability to a minimum of 1 m/s close to the instability value. Accordingly, the resulting flutter wind velocity has the accuracy of one m/s. The variation of the complex eigenvalues vs. wind velocity is plotted in Figure 13, Figure 14 and Figure 15, where the labels within parentheses refer to the mode shapes under wind action: for instance, the label VT indicates a mode characterized by both vertical and torsional components, in which the vertical component is predominant. The real and imaginary parts of the eigenvalues represent the logarithm decay rates and damped vibration frequencies of these modes, respectively. As stated above, the flutter condition occurs when the real part of any eigenvalue becomes positive. The results obtained for the 1962 configuration refer to the two methods previously introduced and recalled hereafter:

• The first method consists of the application of the aeroelastic forces to the upper deck, where the latter are expressed via the flutter derivatives evaluated using the steady-state aerodynamic coefficients obtained for the whole two-deck cross-section (see Figure 6b, Figure 8 and Figure 9a).
• The second method consists of the application of the aeroelastic forces to both decks, where the latter are expressed via the flutter derivatives evaluated using the steady-state aerodynamic coefficients obtained for the upper deck (see Figure 6a, Figure 7 and Figure 9b).

The analyses furnished a critical flutter speed of 40 m/s (144 km/h) and a critical flutter frequency of 0.130 Hz for the single-deck configuration. As regards the double-deck configuration, flutter speeds of 107 m/s (385.2 km/h) and 79 m/s (284.4 km/h), and flutter frequencies of 0.203 Hz and 0.195 Hz, were provided by the first and the second method, respectively.

As already said, the variations of the imaginary part of the complex eigenvalues in Figure 13, Figure 14 and Figure 15 describe the change in the oscillation frequencies with the increasing wind speed, whereas the variations of the real part tell us about stability of each single vibration mode, and thus of the bridge deck. Accordingly, the first torsional-vertical anti-symmetric mode (TVA) is responsible for flutter instability in all the analyzed cases. For higher wind speeds, a second flutter mode is attained in the torsional-vertical symmetric mode (TVS); see Figure 13, Figure 14 and Figure 15. For the double-deck configuration, the first and the second flutter velocities are rather close, being close also the corresponding eigenfrequencies; see Figure 14 and Figure 15. In relation to that, we can say that interaction between modes having similar frequencies could be a triggering factor for instability, thus lowering the critical wind speed.

The two different methods adopted for the double-deck configuration furnish a difference of about 26% and 4% in the estimated critical wind speed and flutter frequency, respectively. Moreover, they provide similar frequency trends as wind speed increases and the same prediction of the flutter mode. Both introduce simplifications from the aerodynamic standpoint: in addition to the approximate evaluation of the flutter derivatives, the former method imposes the torsion rotation axis of the cross-section to coincide with the upper deck axis; the latter assumes that there is no aerodynamic interference between the two decks and neglects the influence of secondary truss-elements. Lastly, it should be noted that the use of flutter derivatives itself represents a strong simplification of the complex aerodynamic behavior of bluff bodies, as bridge decks.

Simplified Formulations and Comparisons

Given the absence of reference data for a direct validation, different simplified formulas available in the literature for predicting the flutter wind velocity were employed for a comparison. The results are collected in

Table 2. Flutter velocity according to different simplified formulations.

Bridge Configuration	Flutter Velocity [m/s]
	Frandsen [41]	Selberg [42]	Rocard [42]	Matsumoto [43]	Van der Put [44]	Fu and Wang [45]	Present Work
Single-deck	39.56	36.08	36.21	36.09	38.70	46.45	40
Double-deck	121	113.86	114.44	113.87	89.23	97.76	107–79

. The simplified formulas adopted are in Equations (5a)–(5f). The first one is based on the reduction of the divergence wind velocity ($U_d)$ and is the only one involving the aerodynamics of the deck by means of the derivative of the moment coefficient ($C_m^{\prime }$). Equations (5b)–(5d) are empirical formulas fitted for thin airfoils, while Equations (5e) and (5f) contain a coefficient fitted for suspension bridge decks.

Frandsen [41]:

$$U_f=\left(B{\omega }_t\sqrt{\frac{2I}{\rho B^4{C}_m^{\prime }}}\right)\sqrt{1-{\left( \frac{{\omega }_v}{{\omega }_t} \right)}^2}$$

(5a)

Selberg [42]:

$$U_f=2.623f_{t}B\sqrt{\left(1-\frac{1}{{\gamma }_{\omega }^2}\right)r_{\alpha }{\mu }_1}.$$

(5b)

Rocard [42]:

$$U_f=6.282f_{t}B\sqrt{\left(1-\frac{1}{{\gamma }_{\omega }^2}\right)\frac{r_{\alpha }^2{\mu }_1}{1+8r_{\alpha }^2}}.$$

(5c)

Matsumoto [43]:

$$U_f=3.71B{f}_t\sqrt{\frac{\sqrt{Im}}{\rho B^3}\left[1-{\left( \frac{f_v}{f_t} \right)}^2\right]}.$$

(5d)

Van der Put [44]

$$U_f=\eta \left[ 1+\left( \frac{{\omega }_t}{{\omega }_v}-0.5\right)\sqrt{\frac{2r}{B}0.72{\mu }_2} \right]\frac{B}{2}{\omega }_v.$$

(5e)

Fu and Wang [45]:

$$U_f=\eta \left(2.5\sqrt{\frac{2r}{B}{\mu }_2}\right)B{f}_t.$$

(5f)

where ${\gamma }_{\omega }=\frac{{\omega }_t}{{\omega }_v}$, $r_{\alpha }=\sqrt{\frac{I}{m{B}^2}}$, ${\mu }_1=\frac{2m}{\rho B^2}$, ${\mu }_2=\frac{4m}{\pi \rho B^2}$, $r=\sqrt{\frac{I}{m}}$. I is the polar mass moment of inertia per unit span length, m is the mass per unit span length, ${\omega }_v, {\omega }_t$ are the angular frequencies of the first vertical mode and the first torsional mode, $f_v, f_t$ are the respective frequencies and $\eta$ is an empirical coefficient representing the difference between the flat plate and a generic profile (equal to 0.7 for suspension bridges [44]). Other simplified formulations based on flutter derivatives can be found in [46].

Results provided by the different formulations for the single-deck configuration are similar to the ones obtained in the present work within the range $\left[-10\%;15\%\right]$. Larger variations are found for the double-deck configuration. In particular, formulas fitted for the thin airfoil provide a considerable overestimation of the critical speed with respect to the formulas fitted for the suspension bridges. The flutter limit provided by the first method is close to the prediction obtained by Equations (5b)–(5d), whereas the flutter velocity provided by the second method is closer to the predictions of Equations (5e) and (5f), and is the most conservative prediction.

Hence, although for a single-deck configuration having the aerodynamic coefficients similar to those of the airfoil it is possible to obtain reasonable results using simplified formulations, for a double-deck configuration it is necessary to perform more detailed analyses.

6. Discussion

The flutter stability of the George Washington bridge was investigated for the two different configurations, i.e., the original single-deck and the current double-deck. Due to the lack of literature data on the aerodynamic parameters for the bridge under consideration, two simplifications were introduced for the description of the motion-related wind loads: one for the calculation of flutter derivatives and the other for the application of aerodynamic loads on the two decks of the current configuration.

The method adopted for the calculation of flutter derivatives allows them to be obtained by analytical formulations based on the steady-state aerodynamic coefficients. The latter parameters were calculated by CFD simulations using the finite element software ANSYS FLUENT. The application of the method provided good results on the bridges chosen for validation, namely the Great Belt and the Akashi Kaykio. The main advantage lies in the relative simplicity of the calculation of the steady-state aerodynamic coefficients, on which the formulation is based. In fact, the computational cost required is much lower than that which characterizes the CFD analyses necessary for a direct numerical calculation of flutter derivatives. Of course, this approach furnishes approximate predictions. To further validate the method, comparisons should be extended to several bridge decks for which both steady-state aerodynamic coefficients and flutter derivatives are available. Once the flutter derivatives were defined, the flutter analysis was performed by a finite element ANSYS Mechanical APDL model endowed with MATRIX27 elements [31].

With regards to the double-deck configuration, two alternative methods were adopted for the definition and application of aeroelastic loads: a first one where lift and moment actions, based on static coefficients of the whole section, are applied to the top deck only; and a second one where lift and moment actions, based on static coefficients of the single-deck section, are applied to both upper and lower decks. Both methods introduce simplifications from the aerodynamic viewpoint. The former provides predictions similar to those of the simplified formulas calibrated for the airfoil. The latter is more aligned with the results provided by simplified formulations accounting for the difference between airfoil and deck cross-sections. With respect to the modal frequencies, the two methods provide similar trends for increasing wind speed, a 4% difference in the critical frequency estimation, and the same prediction of the flutter mode. According to the results of this study, the presence of the lower deck has raised the flutter wind velocity by more than 200%.

Although a definitive validation of the numerical results was not possible, due to the lack of experimental data, the authors believe that the results obtained here can be used for future comparisons with others obtained by different methods. In addition, the proposed approach provides an approximate way for estimating the flutter velocity of double-deck long-span bridges in the absence of more detailed analyses.

7. Concluding Remarks

Two simplified methods were adopted to achieve a preliminary flutter stability assessment of the double-deck George Washington Bridge. Both methods provided a realistic outcome, the second one seeming the most reliable and being the most conservative. The predicted critical flutter wind velocity is about 79 m/s (284 km/h). According to [47], for the bridge site, the wind speed corresponding to approximately a 1.6% probability of exceedance in 50 years is 129 mph (about 57.7 m/s). Given the lack of experimental data, it is suggested that wind tunnel tests are needed to obtain the flutter derivatives of the bridge deck.