Flutter analysis methods for bridges stabilized with eccentric wings

Analysis methods for computing the flutter speed of bridges stabilized against flutter by stationary wings are presented. The wings are placed outboard the bridge deck to achieve a large lateral eccentricity, which enables them to produce enough aerodynamic damping to effectively raise the flutter speed. Given the focus on flutter, other wind effects are neglected. The analysis can thus be carried out in the frequency domain. The most so-phisticated method is based on a specially developed finite aeroelastic beam element, used for modelling a bridge-deck-plus-wings segment, leading to a multi-degree-of-freedom analysis. Such analysis is recommended if the wings do not extend over the full length of the bridge, a design choice that benefits cost efficiency. Second, a simplified two-degree-of-freedom flutter analysis method is described. Simplification is achieved by establishing the wind forces on the wings assuming quasi-steady, instead of unsteady, flow and taking them into account as additional damping and stiffness. Results of example calculations are compared to those of the multi-degree-of-freedom flutter analysis. Finally, it is shown how torsional flutter of a bridge equipped with such wings can be treated in a single-degree-of-freedom analysis. The method is applied to the first Tacoma Narrows Bridge.


Introduction
Bridges are exposed to various wind effects, in particular the static wind force, buffeting, vortex shedding, and flutter. In contrast to the first three effects, which correspond to mechanical response problems, flutter is a stability problemdue to the presence of motion-induced wind forcesand can lead to a sudden collapse. The longer the span of a bridge, the more flutter becomes a driving factor for design and cost. In the detailed design of a long-span bridge, all the effects listed above must be considered. This can be done by state-of-the-art methods for predicting the combined effect of buffeting, vortex shedding, and motion-induced forces and even for taking into account aerodynamic and structural non-linearities in time-domain approaches (Abbas et al., 2017). In the conceptual design of a very long-span bridge, however, it is reasonable to focus on flutter and leave other effects for later consideration. A pure flutter analysis is most conveniently performed in the frequency domain. The basis for this was created by Theodorsen (1934) with his airfoil flutter theory that includes complex functions describing the unsteady aerodynamic forces on a thin flat plate and leads to the critical wind speed for flutter onset (flutter speed). Scanlan and Tomko (1971) adapted this approach to bridge flutter. Instead of by closed-form complex functions, the unsteady aerodynamic forces are described by empirical real functions determined with sectional bridge deck models in the wind tunnel. In this approach, the spatial system is generalized to two degrees of freedom, heave and rotation, using the first vertical bending and torsional modes of vibration. If such a simplification is not justified, multi-mode and finite element approaches are available (Starossek, 1993;Vu et al., 2011;Tamura and Kareem, 2013). It is also possible to include the lateral motion of the bridge deck, but this normally does not much influence the computed flutter response of a bridge (Matsumoto et al., 2010;Vu et al., 2015).
The flutter speed of a bridge can be raised by lateral eccentric wings, running parallel to the bridge deck, that are rigidly attached to the deck by support structures longitudinally spaced at a certain distance (Fig. 1). The wings and their supports form a passive aerodynamic device called the eccentric-wing flutter stabilizer. The wings produce aerodynamic damping of the bridge deck motion that raises the flutter speed. The wing-produced damping of the bridge deck's rotational motion increases quadratically with the wing eccentricity, which explains the advantage of eccentricity. These statements were confirmed by wind tunnel tests and flutter analyses (Starossek et al., 2018). The eccentric-wing flutter stabilizer, in adapted form, can stabilize also other kinds of structures against flutter. The analytical methods presented here can possibly be used in such other applications.
For the conceptual design of the eccentric-wing flutter stabilizer in a project, its potential impact on flutter for conceivable wing configurations and bridge systems should be known. A parametric analytical study has been performed for this purpose (Starossek and Starossek, 2021). To work out the essential dependencies as clearly as possible, the computation focussed on flutter and neglected other effects. On the other hand, the number of different parameter sets to be considered as input parameters was high, which posed a special challenge to the computational tools in terms of ease of modelling and analysis. The computational tools developed and used for that parametric study are presented here.
The length of the wings is an essential design variable that should be accounted for in a flutter analysis. When this length is less than the length of the bridge, the aerodynamic properties are no longer constant along the bridge axis and, consequently, (mainly) the torsional component of the flutter mode shape differs from the respective vacuum mode shape. To accurately take this into account, a multi-degree-of-freedom (MDOF) analysis method is developed. It is based on a known finite element flutter analysis method (Starossek, 1992(Starossek, , 1993 that is extended to include the wings. This is achieved by a specially developed finite aeroelastic beam element capable of simultaneously modelling the bridge deck and the wings and hence facilitating the modelling of the many bridge-wing configurations to be considered in a parameter study. Moreover, a simplified two-degree-of-freedom (2-DOF) flutter analysis method is described. Simplification is achieved by establishing the wind forces on the wings on the assumption of quasi-steady flow past a thin flat plate and including them in a conventional 2-DOF flutter analysis as additional damping and stiffness. A reduction factor is suggested to approximately take into account the wing length in case it is less than the bridge length. Results of example calculations are compared to those of the MDOF analysis.
Finally, it is shown how torsional flutter of a bridge with wings can be treated as a single-degree-of-freedom problem leading to a simple flutter condition equation. The effect of the wings is reflected in that equation by an additional damping term. The method is applied to the first Tacoma Narrows Bridge.

Scope
The assumptions and methods of classical bridge flutter theory apply. That is, steady-state harmonic vibration is assumed and the analysis is performed in the frequency domain. Concerning the wind forces, only motion-induced lift forces and motion-induced aerodynamic moments are considered. They are linearly related to vertical displacements and rotations of the structure, and the respective velocities and accelerations, by complex aerodynamic coefficients that are functions of the degree of unsteadiness of the flow. In other words, the general framework of Theodorsen's (1934) flutter theory applies, even if the aerodynamic coefficient functions relating aerodynamic force quantities to displacement quantities are different from those originally given by Theodorsen (1934), but instead are determined specifically for the respective aerodynamic contour (Scanlan and Tomko, 1971). It is assumed that the oncoming wind is either non-turbulent or that any turbulence has little effect on the flutter speed or that it can be taken into account by the aerodynamic coefficient functions. Consequently, eigenvalue problems are established and solved.
With regard to the wings, it is assumed that there is no aerodynamic interference between the windward and leeward wings and the bridge deck so that classical bridge flutter theory can be applied separately to each of these three elements for determining the respective aerodynamic forces. In practice, interference can be prevented by positioning the wings above or below the bridge deck with sufficient vertical offset to the deck and between them. Furthermore, the wings are aerodynamically shaped so that drag and wake are small. Flutter analyses based on these assumptions have been validated by wind tunnel tests (Starossek et al., 2018). Turbulence of the oncoming wind acting on the wings and a possibly resulting effect on the buffeting-induced vibration of the bridge are neither considered here. Wind tunnel tests and analyses for Imja Bridge, Korea have shown that the response of the bridge to buffeting, and also to vortex shedding, is mostly reduced when wings are added (Meyer, 2019). This can be attributed to the aerodynamic damping produced by the wings that reduces vibrations due to any kind of excitation. While not within the scope of this paper, this means that the wings also improve serviceability and reduce fatigue damage.

General
The finite element developed here captures the structural forces and the motion-induced wind forces that arise in the bridge deck and the two wings. It represents a bridge-deck-plus-wings segment of the bridge extending in the longitudinal direction x and including rigid support structures at its ends (Fig. 2). A known finite aeroelastic beam element for modelling the bridge deck (Starossek, 1992(Starossek, , 1993 is extended to include the wings.
The element has seven independent degrees of freedom: the vertical displacements of the bridge deck end nodes, δ 1 and δ 2 , the angular displacements of the bridge deck end nodes about the transverse axis, δ 3 and δ 4 , and the angular displacements of the bridge deck end and centre nodes about the longitudinal axis, δ 5 , δ 6 , and δ 7 . For modelling the torsional deformation of the bridge deck, a minimum number of two degrees of freedom is required, which correspond to the angular displacements of the element end nodes about the longitudinal axis (a respective finite aeroelastic beam element without wings is described by Starossek (1992Starossek ( , 1993). The introduction of a third degree of freedom, that is, the angular displacement at the element centre, results in a more uniform accuracy in modelling the bending and torsional deformations and thus leads to a balanced relationship between numerical effort (which increases with the degree of segmentation of the bridge deck and hence the number of U. Starossek and R.T. Starossek beam elements) and accuracy (Starossek, 1992(Starossek, , 1993. Elements are connected to each other at the bridge deck end nodes only, that is, by the respective independent degrees of freedom δ 1 , δ 2 , δ 3 , δ 4 , δ 5 , and δ 7 .
Dependent degrees of freedom are the vertical displacements of the windward wing end nodes, δ W1 and δ W2 , the angular displacements of the windward wing end nodes about the longitudinal axis, δ W3 and δ W4 , the vertical displacements of the leeward wing end nodes, δ L1 and δ L2 , and the angular displacements of the leeward wing end nodes about the longitudinal axis, δ L3 and δ L4 .
This choice of independent and dependent degrees of freedom and the kinematic relationships between them are based on the following simplifying assumptions, which facilitate the representation of the combined bridge-deck-plus-wings segment by one finite element.
• The wings of adjacent segments are not connected to each other.
• The wings can freely rotate about their transverse axis at their end supports. • The wings are torsionally fixed at their ends to the support structures.
• The support structures are rigid in bending.
This has the following consequences.
• The deformation state is defined by only seven independent degrees of freedom. • The dependent degrees of freedom depend on δ 1 , δ 2 , δ 5 , and δ 7 .
• Node displacements do not cause bending deformations of the wings.
• The wings are twisted depending on δ 5 and δ 7 .
In view of the support conditions actually chosen in a project, the above simplifications are expected to be close to reality. In particular, structural continuity of the wings would involve them in the global load transfer of the bridge deck, which is not desirable. In any case, given the relatively small forces acting on the wings in comparison to the bridge deck, these simplifications only minimally affect accuracy, while they make additional independent degrees of freedom unnecessary. For further simplification, the torsional rigidity and the mass moment of inertia of the wings as well as the mass of the support structures are neglected and all cross sections are assumed to be doubly symmetric.
The independent and dependent displacements are summarized in respective vectors: The kinematic relationships of the latter to the former, expressed as linear transformations, are where a W c is the lateral distance of the windward wing centre, and a L c the lateral distance of the leeward wing centre, from the bridge deck axis (wing eccentricities). The transposed matrices t ⊤ W and t ⊤ L transform the forces that correspond to the respective dependent degrees of freedom to the forces that correspond to the independent degrees of freedom.

General
The 7 × 7 stiffness, mass, and aerodynamic matrices of the element, k, m, and a, provide linear transformations of the independent displacement vector δ into corresponding element force vectors f s , f i , and f a , which comprise the elastic restoring, inertia, and motioninduced wind forces, respectively, according to in which the factor − ω 2 results from the steady-state harmonic vibration approach and ω = circular frequency of motion (Starossek, 1992(Starossek, , 1993.
Each of the element matrices can be represented as the sum of two contributions, that is, where k d , m d , and a d are the contributions from the bridge deck, established in Section 3.2.2, and k c , m c , and a c are the contributions from the wings, established in Section 3.2.3.

Contributions of bridge deck
The contributions from the bridge deck were derived by Starossek (1992Starossek ( , 1993 using the principle of virtual displacements with cubic interpolation functions for vertical (bending) displacements and quadratic interpolation functions for angular (torsional) displacements (consistent matrices). They are as follows.
where the superscripts h and α of the submatrices refer to displacements and forces in the vertical or angular (torsional) direction, respectively, l = element length, E = Young's modulus, J = second moment of area, N = axial force (inducing geometric stiffness), G = shear modulus, J d = torsion constant, of the bridge deck, all of which are assumed to be constant over the element length, where m = mass per unit length, I = mass moment of inertia per unit length, of the bridge deck, both assumed to be constant over the element length, where ρ = air density, b = half chord of aerodynamic contour of bridge deck, c hh , c hα , c αh , c αα = complex unsteady aerodynamic coefficients of bridge deck contour that are functions of where k = reduced frequency and u = wind speed. Also b and u and hence k and the aerodynamic coefficients are assumed to be constant over the element length. They can vary between elements. In this case, each element p has its own reduced frequency k p and aerodynamic coefficients. The ratios of the k p to one another are defined by the spatial distribution of b and u and all k p and aerodynamic coefficients are determined by a single reference value k. Hence the solution method described in Section 4.1 below can still be used without major complication (Starossek, 1992(Starossek, , 1993. The complex unsteady aerodynamic coefficients used here are related to the real so-called flutter derivatives used by Scanlan (1992) by where i = imaginary unit. The merits of a complex versus real representation of the aerodynamic coefficients and forces are discussed by Starossek (1998). Closed-form expressions for the four complex unsteady aerodynamic coefficient functions were derived by Theodorsen (1934) and compiled by Starossek (1992). They are based on the assumption of potential flow past a thin flat plate and lead to fairly accurate results if the bridge deck contour is sufficiently streamlined. Otherwise, unsteady aerodynamic coefficient functions determined from wind or water tunnel measurements with sectional bridge deck models or from computational fluid dynamics calculations (Thiesemann, 2007;Starossek et al., 2009) are employed.

Contributions of wings
It follows from the assumptions listed in Section 3.1 that the stiffness contribution of the wings is zero: Before establishing the contributions m c and a c of the wings to the 7 × 7 mass and aerodynamic matrices of the element, the 4 × 4 matrices m W , a W , m L , and a L are established. They relate the dependent displacements of the windward and leeward wing ends, summarized as δ W and δ L , to the corresponding forces at the respective wing ends, summarized as f W,i , f W,a , f L,i , and f L,a , according to The task of establishing these matrices is identical for both wings given the choice of dependent degrees of freedom. Therefore, the subscripts W and L are temporarily omitted and the matrices m and a, valid for each of both wings, are derived.
Again referring to the assumptions listed in Section 3.1, it is found that the coefficients of m associated with angular (torsional) displacements and/or forces vanish. The remaining coefficients m hh ij , associated with vertical displacements and forces, are calculated using the expression obtained from the principle of virtual displacements, where m c = mass per unit length of wing on the respective side of the bridge deck (possibly including a contribution of the support structures). The interpolation functions ψ i and ψ j are linear functions of x, in accordance with the above assumptions. Assuming that m c is constant over the element length, the following result is obtained.
with the abbreviations Matrix a can be written Using the principle of virtual displacements, the following expressions for calculating the coefficients of the submatrices are obtained (Starossek, 1992(Starossek, , 1993.
where b c = half chord of a wing, c c,hh , c c,hα , c c,αh , c c,αα = complex unsteady aerodynamic coefficients of wing that are functions of which is the reduced frequency associated with the considered wing. The latter is related to the reduced frequency of the bridge deck, k, by Linear interpolation functions for both vertical and angular (torsional) displacements are used, in accordance with the assumptions of Section 3.1. When b c and u and hence k c and the aerodynamic coefficients are constant over the element length, the following result is obtained (again using the submatrix term z defined in Eq. (26)).
When the wings are sufficiently streamlined, Theodorsen's closedform expressions (Theodorsen, 1934;Starossek, 1992) can be used for the unsteady aerodynamic coefficients. Incidentally, the reduced frequency k c of a wing will often be on the order of 0.05 or less, in which case the unsteady aerodynamic forces considered here could be replaced, without much loss in accuracy, by quasi-steady aerodynamic forces (Försching, 1974) (see also Section 5).
The matrices m W and a W are obtained when substituting the corresponding windward wing properties m W c , b W c , c W c,hh , c W c,hα , c W c,αh , c W c,αα , the matrices m L and a L are obtained when substituting the corresponding leeward wing properties m L c , b L c , c L c,hh , c L c,hα , c L c,αh , c L c,αα , into the expressions for m and a given by Eqs. (25) and (34).
The contributions of both wings to the 7 × 7 mass and aerodynamic matrices of the element now follow from the transformations leading to where and in which the abbreviations are used, and where a hh While the flutter suppression effectiveness of the wings mainly stems from the aerodynamic damping they provide to the bridge deck's rotational motion, the wings also produce aerodynamic damping of the vertical motion and aerodynamic stiffness with respect to the vertical and rotational motion of the bridge deck. All these effects are reflected in a c .
When identical wings with the same mass per unit length m c and the same half chord b c are placed with the same lateral eccentricity a c on both sides of the bridge deck, the above submatrices take the following simpler forms.
and a hh

Modelling alternatives
The finite aeroelastic beam element developed here captures the forces that arise in the bridge deck and in the wings initially separately but then combines them by using kinematic relationships. Alternatively, the bridge deck and the two wings could each be separately modelled with an ordinary finite aeroelastic beam element, like the one described in Section 3.2.2, and then the full structural system could be assembled from these elements. The modelling and analysis effort increases accordingly, which would make this alternative unsuitable for a parametric flutter analysis study.
Another alternative is to model a combined bridge-deck-plus-wings segment by one ordinary finite aeroelastic beam element, like the one of Section 3.2.2, using unsteady aerodynamic coefficient functions that have been determined for the aerodynamic contour of the combined assemblage. Results reported by Meyer (2019) have been obtained in this way, using aerodynamic coefficient functions determined from wind tunnel measurements using a sectional bridge-deck-plus-wings model. This alternative can also be employed in a 2-DOF flutter analysis. Its flexibility with regard to the study of design variants is even smaller than that of the alternative described in the previous paragraph.

Method
The overall structural system is modelled with finite elements. The bridge deck with the wings attached to it is modelled using the finite aeroelastic beam element presented in Section 3. If wings are partly or completely absent, the same beam element is used, only that contributions k c , m c , and a c in Eq. (7) are zero for the elements that model the respective parts of the bridge deck. The stiffness, mass, and aerodynamic matrices of the system, K, M, and A, are assembled from the corresponding element matrices (Clough and Penzien, 1975), that is, as far as the bridge deck and the wings are concerned, from k, m, and a according to Eq. (7). The system matrices are n × n matrices, where n is the number of degrees of freedom of the system. An eigenvalue problem of the form { is established for solving this multi-degree-of-freedom (MDOF) problem, where Φ = system node displacement vector and g = structural damping parameter (Starossek, 1992(Starossek, , 1993. The structural damping represented by g is proportional to the elastic restoring forces and in counter-phase to velocity (Försching, 1974). Since the element aerodynamic matrices a are complex non-Hermitian matrices that depend on k, the same holds for the system aerodynamic matrix, A. When k ∈ R + is fixed, a linear complex non-Hermitian eigenvalue problem, with the eigenvalue ω 2 , is obtained. Its solution leads to n generally complex eigenfrequencies ω i = ω The minimum wind speed found in this manner is the flutter speed, which in the following is likewise called u. The corresponding ω ′ j is the circular flutter frequency, in the following simply called ω. For safely identifying the flutter speed, the search starts with large and then decreasing values of k or with small and then increasing values of u red = 1/k = reduced wind speed. The governing case is usually associated with one of the lowest eigenfrequencies and hence the eigenvalue analysis can be limited to a number of lowest, say two to eight, eigenmodes.
The above approach to account for structural damping leads to the simplest form of the MDOF flutter eigenvalue problem. Equivalent viscous modal damping ratios-to-critical, ξ i , for any eigenmode i can be determined as where ω 0 i = natural circular frequency of the undamped system in a vacuum (without aerodynamic forces) associated with eigenmode i (Starossek, 1992). If present, they would lead to the same flutter solution as g. When the structural damping is assumed to be viscous damping, the eigenvalue problem reads where C = system damping matrix. C can be constructed to possess special properties, assigning, for instance, a different damping ratio-tocritical ξ i to every eigenmode i (Clough and Penzien, 1975). A quadratic instead of a linear eigenvalue problem is now to be solved. Apart from this, the procedure outlined above for finding the flutter speed applies. A FORTRAN program implementing the finite aeroelastic beam element presented here and solving the MDOF flutter eigenvalue problem according to Eq. (56) using a vector iteration method has been developed. Because of the algebraic structure of the system matrices (non-normal matrix pair), iteration with right and left eigenvectors is performed (Zurmühl and Falk, 1986).

Example
The method is applied to an example problem. A simply supported girder with torsionally fixed ends and streamlined contour is considered. Streamlined wings are placed on both sides of the girder in a symmetrical configuration. The input data representing system properties and the results are presented in non-dimensional form. The respective framework of non-dimensional quantities is explained first.
The input data characterizing the girder are the frequency ratio, ε, defined as where ω α = natural circular frequency of torsional vibration, and ω h = natural circular frequency of vertical bending vibration, of the undamped system without wings in a vacuum (without aerodynamic forces), these frequencies being associated with the lowest modes of vibration, which, for the considered girder, are both symmetrical; the (structure-to-air) mass ratio, μ, defined as where m = mass per unit length of girder; ρ = air density, b = half chord of aerodynamic contour of girder; the reduced mass radius of gyration, r, defined as where I = mass moment of inertia per unit length of girder; and the structural damping parameter g defined in Section 4.1. Because the mass of the wings will be specified separately, parameters ε, μ, and r refer to the girder only, without considering the influence of the wing mass. The input data characterizing the wings are the relative wing eccentricity, ã c , defined as the relative wing width, b c , defined as U. Starossek and R.T. Starossek the relative wing length, L c , defined as where L c = length of wing on one side of girder, understood to be contiguous and centred in relation to midspan, L = length of girder; and the relative wing mass, m c , defined as where m c = mass per unit length of wing on one side of girder.
The results are presented in terms of the reduced frequency, k, defined by Eq. (19), the non-dimensional flutter frequency, ω, and the non-dimensional flutter speed, ζ. The latter two quantities relate the circular flutter frequency, ω, or the flutter speed, u, to system properties and are defined as A further non-dimensional result is the flutter speed increase ratio, R, which is defined as the flutter speed of the structure with wings to the flutter speed of the structure without wings, that is, and is an indicator of the flutter suppression effectiveness of the wings. Note that the non-dimensional results defined in this manner only depend on the non-dimensional input data defined above (Starossek and Starossek, 2021). Also note that the non-dimensional quantity defined by Eq. (68)  Properties μ, r, ã c , b c , and m c can in principle vary over the length of the girder, but are assumed to be constant here. The relative wing length, L c , is varied from 0 to 1 in 0.04-increments. Consequently, the structure is segmented into 50 identical elements leading to 199 degrees of freedom. Because the aerodynamic contours of both the girder and the wings are assumed to be streamlined, Theodorsen's unsteady aerodynamic coefficient functions (Theodorsen, 1934;Starossek, 1992) are used for both kinds of contours. They are computed based on a rational function approximation of the Theodorsen function (Bruno et al., 1987).
Properties L, b, E, G, and ω h are set to one and ρ to 1.225, which is immaterial to the non-dimensional results presented below because these only depend on non-dimensional properties. The axial force, N, is zero. The remaining input data are chosen such that the non-dimensional property values specified in Eqs. (70) and (71) are met.
For the representation of ω α and ω h , and thus ε, the respective analytical free-vibration frequency formulae are used for determining the input quantities J and J d . The ensuing error due to the interpolation functions of the finite element method being approximations of the true modes of vibration is negligible. 26 separate flutter analyses are performed, one for each value of L c . Selected results for three different wing lengths are shown in Table 1.
Thus, the maximum flutter speed increase ratio is found to be R = 8.5022 2.8348 = 3.00 which indicates a high flutter suppression effectiveness of the wings for the girder and wing properties chosen here. The further results are illustrated in Fig. 3, where the flutter speed increase ratio is plotted against the relative wing length. When the relative wing length is reduced, starting at 1, the flutter speed increase ratio initially decreases only slightly. This is because the reduction of the wing length begins at the girder ends, areas where both the vibration amplitudes, and thus the damping forces generated by the wings, and the effectiveness of these forces are small (peculiarities that are taken into account in the generalization formula Eq. (82)). It can also be seen in Fig. 3 that with a relative wing length of 0.5 the flutter speed increase ratio is about 2, that is, the flutter speed is still doubled when the wings are placed over only half the span.
Results obtained with the FORTRAN program, implementing the finite aeroelastic beam element and MDOF flutter analysis method described above, are compared to those obtained from 2-DOF flutter analyses using the method described in Section 5.1 below. The comparison is presented in Section 5.2.

Table 1
Results of MDOF flutter analysis for various relative wing lengths L c .

Method
The bridge system is generalized to a system with two degrees of freedom (2-DOF), heave and rotation, using a shape function that resembles either the lowest symmetric or the lowest antisymmetric modes of vertical and torsional vibration, whatever governs flutter. The corresponding natural frequencies of the bridge, ω h and ω α , are adopted for the generalized system. When the generalized system properties m and I are set equal to the prevailing mass and mass moment of inertia per unit length of the bridge, the generalized motion-induced wind forces are likewise to be established as forces per unit length as they act on bridge deck and wings. This assumes that every existing wing extends over the full length of the bridge.
The motion-induced wind forces on the wings could be taken into account in a similar manner as for the finite aeroelastic beam element presented above, that is, using unsteady aerodynamic coefficient functions established for the wings and a transformation from dependent to independent degrees of freedom. When a conventional 2-DOF flutter analysis program is at hand, it is more convenient to consider the motion-induced wind forces on the wings as being quasi-steady and to take these forces into account by simply modifying the input to such a program. (The motion-induced wind forces on the bridge deck are always considered as being unsteady.) The condition for this simplification is a small reduced frequency as explained further in Section 5.2 below. The procedure explained in the following assumes that the wings are sufficiently streamlined. According to potential theory, the lift force per unit length, L, acting on a thin flat plate with a half chord b c exposed to a steady flow with a speed u at a stationary angle of attack α is where ρ = air density. Its resultant acts at the leading quarter point of the plate (Försching, 1974). A similar force acts on a streamlined wing. An angular velocity of the bridge deck, α, leads to a vertical velocity of a wing, a cα , when the wing is placed at an eccentricity a c . Superposing this vertical motion onto the horizontal wind flow, an apparent stationary angle of attack is formed when the flow is thought of as quasi-steady. The resulting lift force per unit length on a wing and the respective moment per unit length on the bridge deck, M, both acting in opposite direction of the motion, are where the positive sign holds for a windward wing and the negative sign for a leeward wing. Hence a wing produces a viscous aerodynamic damping of the rotational motion of the structure that is given by the expression (related to unit length) Given that the wing is assumed to extend over the full length of the bridge, the corresponding damping ratio-to-critical is where ω α = natural circular frequency of rotational motion of an undamped system in a vacuum. The definitions in Section 4.2 of nondimensional input parameters and results are adopted. Hence the damping ratio-to-critical of the rotational motion contributed by one wing can also be expressed as where k = reduced frequency (defined by Eq. (19)), ω = non-dimensional flutter frequency. If the program at hand accepts as input damping parameters g as defined in Section 4.1, instead of damping ratios-tocritical ξ, the damping of the rotational motion contributed by one wing can be taken into account with the quantity g α,c computed from (compare Eq. (58)) The total aerodynamic damping contributed by the wings is the sum of the windward wing and the leeward wing contributions, which can both be determined from Eq. (78) or (79). This contribution is added to the respective structural damping of rotational motion quantity, that is, ξ α or g α , which is then used as modified program input. Further aerodynamic damping contributions of the wings exist but are small and neglected. This also holds for the damping of the vertical motion provided by the wings. Thus, ξ h,c and g h,c are set to zero and the input quantity for the structural damping of vertical motion, that is, ξ h or g h , remains unchanged.
The wings also produce aerodynamic stiffness that can be considered in the same quasi-steady fashion leading to the additive rotational stiffness term where, except for the expression in parentheses, the positive signs hold for the aerodynamic stiffness contribution of the leeward wing and the negative signs for that of the windward wing, and vice versa in the expression in parentheses. Both contributions are added to the structural rotational stiffness, k α . When non-dimensional input parameters are used, the expression is evaluated, with subscripts W and L referring to windward and leeward wing, respectively, and added to the frequency ratio, ε, to take into account the aerodynamic stiffness contributed by the wings. A correspondingly modified ε should also be used in Eqs. (78) and (79). When identical wings with the same eccentricity are provided on both sides of the bridge deck, the stiffness contributions of the wings cancel. (Without the approximation in Eq. (80), there remains a small, and usually negligible, residue.) Further aerodynamic stiffness contributions of the wings exist but are small and neglected. In such a generalized 2-DOF procedure, the input parameters refer to the entire system. When applied to a suspension bridge, for instance, the mass of the suspension cables should be included in m and I and in the ensuing quantities μ and r. (More generally, these quantities are generalized properties related to the distributed mass of the system by the respective mode of vibration.) Furthermore, the wing mass, m c , can be included in m, I, μ, r, and its effect on ω h , ω α , and ε can thus be taken into account, as further outlined by Starossek and Starossek (2021).
The above expressions for determining the aerodynamic damping and stiffness contributions of the wings contain the wind speed, u, or the reduced frequency, k, and the non-dimensional flutter frequency, ω. All refer to the flutter state and are unknown at the beginning of the calculation. Thus, the use of these expressions in a flutter analysis requires iteration.
So far it has been assumed that every existing wing extends over the full length of the bridge. The underlying generalization is revisited to expand it to cases where a wing is shorter than the bridge.
When the system properties are constant along the length of the bridge, the corresponding terms can be drawn in front of the generalization integrals. Furthermore, when the same shape function is used in all generalization integrals, these integrals cancel and thus the exact form of the shape function becomes irrelevant. This justifies the procedure carried out above, in which the bridge and wing properties are used directly as generalized system properties. It is accurate to the extent that 1) the system properties are constant along the bridge length, and 2) the vertical and torsional components of the flutter mode shape coincide, the latter being required to justify the use of the same shape function. Both conditions are mostly met when wings are absent or every existing wing extends over the full length of the bridge. Otherwise, both conditions are violated: the wing properties are not constant and, given that essentially only the torsional motion is affected by the wingproduced forces, the torsional component of the flutter mode shape will change when a wing is shorter than the bridge and differ from the vertical component that remains virtually unchanged (Starossek and Starossek, 2021). Therefore, generalisation becomes less accurate for such a case. If performed anyway, it leads to the reduction factor  82) can be expanded accordingly. The consideration of the reduction factor in the correspondingly modified input data is explained further by Starossek and Starossek (2021).

Example
The example problem of Section 4.2 is revisited. For the considered girder, flutter is governed by the lowest symmetric modes of vertical and torsional vibration, which are the lowest modes. The 2-DOF flutter analysis is done with a custom program implemented on a pocket calculator (Starossek, 1992), which uses the unsteady aerodynamic coefficient functions according to Theodorsen (1934) to take into account the motion-induced wind forces on the girder (as in Section 4.2). The program cannot take wings into account directly. However, it accepts as input quantities the structural damping parameter g separately for heave, g h , and rotation, g α . The damping parameter g α is zero if wings are absent (since the structural damping, represented by g, is set to zero). When wings are present, g α is determined from Eq. (79) multiplied by 2, since there are identical wings on both sides of the girder, and multiplied by F according to Eq. (82) for the partial-length wings condition (L c = 0.48). The damping parameter g h is zero in either case.
Because the wings are assumed to be massless, and neither add aerodynamic stiffness given their symmetrical layout, properties ε, μ, and r do not need to be modified and can directly be adopted from Eq.
(70) even if wings are present. These properties are accepted by the program as direct input quantities. Properties ã c and b c are given by Eq.
(71). The ensuing damping parameter g α and the results for three different wing lengths L c are shown in Table 2. The last column of Table 2 shows the relative deviation of the flutter speeds computed here from the values determined by MDOF flutter analyses reported in Table 1. For the system without wings (L c = 0), the results are virtually identical to those obtained from MDOF analysis. The small deviations are attributed to the MDOF discretization and the use of the analytical free-vibration frequency formulae in the representation of ε (see Section 4.2).
For the system with full-length wings (L c = 1), the results are in good agreement with those obtained from MDOF flutter analysis. The flutter speed differs by less than 1%. The good agreement is ascribed to a small reduced frequency of the wings, that is, k c =b c k = 0.013, which is much smaller than that of the bridge deck. In such cases, that is, more specifically, for k c < 0.05, the motion-induced wind forces can be determined on the basis of the quasi-steady flow assumption, as done here, without much loss of accuracy (Försching, 1974). The good agreement can also be considered as an element of validation of the further simplifications of the 2-DOF analysis method proposed here as well as of the FORTRAN program and the underlying MDOF flutter analyses method described above.
For the system with partial-length wings (L c = 0.48), the reduction factor according to Eq. (82) is F = 0.79768, assuming wings centred at midspan and choosing a sine half-wave as shape function ψ(x). The deviation of the values computed by 2-DOF flutter analysis from those determined by MDOF flutter analysis becomes noticeably larger, according to expectation (Section 5.1). The flutter speed is overestimated by 3.3%, which would be just about acceptable in practice.
The 2-DOF flutter analysis method proposed here was applied to a variety of other sets of girder properties and wing configurations and the results were compared to those obtained from MDOF flutter analysis. Again, excellent or good agreement was found for all cases with no-wing or full-length wing configurations. For certain partial-length wing configurations, deviations larger than the one reported in the previous paragraph were found, which in some cases would not be acceptable in practice (Starossek and Starossek, 2021).

Method
In bridges where the bridge deck has a bluff aerodynamic contour, the motion-induced wind forces can be such that vertical and torsional vibrations remain nearly uncoupled and flutter occurs in a purely or mainly torsional mode (torsional flutter). Due to the missing or weak aerodynamic coupling, generalization to a single-degree-of-freedom system is possible and only the solution of the equation for rotational motion is required (Starossek, 1994). In such cases, flutter is mostly governed by the aerodynamic moment induced by rotational (i.e., torsional) motion and in phase with angular velocity. This moment is determined by c ′′ αα , the imaginary part of the complex unsteady aerodynamic coefficient c αα of the bridge deck (related to the corresponding real-number flutter derivative (Scanlan, 1992) by c ′′ αα = 8/π⋅A * 2 ). When neglecting the real part of c αα , the flutter frequency becomes identical to the lowest natural frequency of torsional vibration and the flutter speed of the bridge without wings can be determined from the condition where ξ = structural damping ratio-to-critical, μ = (structure-to-air) mass ratio, r = reduced mass radius of gyration, as defined by Eqs. (61) and (62), and hence where μr 2 = (structure-to-air) mass moment of inertia ratio, I = mass moment of inertia per unit length of bridge deck (including suspension cables if present), ρ = air density, b = half chord of aerodynamic contour of bridge deck. For convenience, c ′′ αα is now considered a function of the reduced wind speed, u red , defined as rather than of the reduced frequency, k. The flutter condition expressed by Eq. (83) is met when the negative aerodynamic damping produced by the wind forces acting on the bridge deck (associated with positive values of c ′′ αα ) is as large as the inherent structural damping and thus the total damping is zero. The flutter speed, u, follows from solving Eq. (83) for u red and then Eq. (85) for u, setting the circular flutter frequency, ω, to the lowest natural circular frequency of torsional vibration, ω α .
When identical wings are added on both sides of the bridge deck along the full length of the bridge, Eq. (83) is replaced by the almost equally simple condition in which an additional term appears that represents the aerodynamic damping provided by the wings, with ã c = a c /b = relative wing eccentricity, b c = b c /b = relative wing width, a c = wing eccentricity, b c = wing width. It follows from Eq. (78) and hence is based on considering quasi-steady potential flow past a thin flat plate and assumes that the wings are sufficiently streamlined. The influence of the wing mass can be taken into account in ω α and μr 2 . The aerodynamic stiffness contributions of the wings and their influence on ω α could be considered using Eq. (80) or (81) (modifying ω α at the same ratio as ε). For the symmetrical wing configuration considered here (with wings on both sides of the bridge), however, these contributions cancel. If the wings are shorter than the bridge, this can approximately be taken into account by multiplying the last term in Eq. (86) by the reduction factor F according to Eq. (82).

Example
The collapse of the Tacoma Narrows Bridge in 1940 has been ascribed to torsional flutter. Thus, the above equations are tentatively applied for computing the flutter speed of that bridge. The structural data are b = 5.94m; I = 202, 400kgm 2 / m; f α = 0.233Hz; ξ = 0.54% where f α = lowest natural cyclic frequency of torsional vibration and thus ω α = 2πf α = 1.464 1/s (Farquharson et al., 1949). With ρ = 1.225 kg/m 3 follows μr 2 = 42.2. Some of the structural properties given in the literature differ from one another and are inconclusive. Therefore, the above mass moment of inertia I has been determined by the authors based on the design drawings reproduced by Farquharson et al. (1949). It includes the mass contribution of the suspension cables. The above torsional frequency (used here) was observed on the morning of the collapse when the bridge abruptly began to vibrate in the first antisymmetric mode. Immediately before collapse, this frequency changed to f α = 0.200 Hz (not used here, as presumably belonging to an already damaged system). For bluff bridge decks, the unsteady aerodynamic coefficient functions, including c ′′ αα (u red ), strongly differ from those determined by Theodorsen (1934). In such cases, they are established, for discrete values of u red , by measurements in a wind or water tunnel using a sectional bridge deck model or by computational fluid dynamics (CFD) analysis. Both methods were applied to a variety of bridge deck cross sections (Thiesemann, 2007;Starossek et al., 2009). For a simplified version of the first Tacoma Narrows Bridge deck (Fig. 4), the results shown in Table 3 were obtained from water tunnel testing for a Reynolds number of 200,000 and an angular displacement amplitude of 8 • .
When linear interpolation between the values given in Table 3 is used, solving Eq. (83) leads to u red = 1.468, corresponding to k = 0.681, and the flutter speed of the bridge without wings, computed from Eq. (85), becomes u = 12.8 m/s. This value is distinctly smaller than the wind speed of 18.8 m/s reported by Farquharson et al. (1949) to have caused the flutter-induced collapse of the Tacoma Narrows Bridge. However, it is in line with the analytical results of other authors (Larsen, 1998;Matsumoto et al., 2003). Reasons for the apparent discrepancy are suggested by Matsumoto et al. (2003). On this basis, it is concluded that the above equations are sufficiently accurate at least for estimating the relative flutter speed increase provided by the eccentric wings studied here.
The flutter speed of the Tacoma Narrows Bridge is now computed when identical wings are added on both sides of the bridge deck along the full length of the bridge. The same procedure applies except that Eq.
(83) is replaced by Eq. (86). The influence of the wing mass on ω α and μr 2 is neglected. Table 4 shows the results for the reduced wind speed, u red , the reduced frequency, k, the flutter speed, u, and the flutter speed increase ratio (relative to the flutter speed without wings), R, for various wing configurations defined by ã c and b c .
These results indicate that a substantial increase in flutter speed is achievable by adding wings provided certain minimum values are chosen for their eccentricity and width. If the wing configuration   defined by ã c = 2.00 and b c = 0.100 is selected, a flutter speed increase on the order of 280% is obtained. This configuration consists of wings that are 1.2 m wide with respective centres located about 6 m outside the edges of the bridge deck. A more general conclusion can be drawn from an inspection of Eq. (86). The terms on the right-hand side of that equation represent positive damping that raises the flutter speed. The first term on the righthand side stands for the inherent structural damping; the second term represents the aerodynamic damping provided by the wings. It is concluded that, in the case of uncoupled torsional flutter, the relative flutter speed increase due to the wings (i.e., the flutter speed increase ratio) is the larger, the smaller the structural properties ξ and μr 2 .
For comparison, analyses using the approaches of Sections 3 and 4 are performed for the above specified structural properties of the Tacoma Narrows Bridge, but assuming a streamlined (instead of bluff) bridge deck and using Theodorsen's unsteady aerodynamic coefficient functions. Since flutter now occurs as a coupled vertical-torsional vibration, the mass of the bridge deck and the vertical natural frequency are also of interest. They are m = 8,186 kg/m (determined by the authors based on the design drawings, including the suspension cables) and f h = 0.145 Hz (Farquharson et al., 1949). The flutter speed of the bridge without wings results in u = 37.5 m/s, almost three times the value u = 12.8 m/s determined above for uncoupled torsional flutter of the bluff bridge deck. Identical wings are now added on both sides of the deck, again along the full length of the bridge. The wing configuration is defined by ã c = 2.00 and b c = 0.100 and the wing mass is neglected. The flutter speed of the bridge becomes u = 45.2 m/s. The flutter speed increase ratio thus is R = 45.2/37.5 = 1.21, a value much smaller than the value 3.80 found for the same structural properties and wing configuration, but a bluff bridge deck.

Conclusions
Bridges provided with lateral eccentric wings rigidly attached to the bridge deck were considered and methods for computing the flutter speed of such systems were presented. When the wings do not extend over the full length of the bridge, a design choice that benefits cost efficiency, the flutter mode shape changes, which influences the flutter speed. This effect can accurately be taken into account in a multi-degreeof-freedom (MDOF) flutter analysis. Such analysis is conveniently performed using the specially developed finite aeroelastic beam element presented here. The element is used to simultaneously model the bridge deck and the wings, so that no additional independent degrees of freedom are required to represent the wings and modelling and analysis are simplified. Example calculations were performed for a girder-pluswings system, in which Theodorsen's unsteady aerodynamic coefficient functions were used for the girder as well as the wings and the length of the wings was varied between a no-wing and a full-length wing configuration. The results show that a substantial increase in flutter speed can be achieved by wings. They also confirm that the wing length is an important design variable.
Moreover, a two-degree-of-freedom (2-DOF) flutter analysis method was described in which the wind forces on the wings are established on the assumption of quasi-steady potential flow past a thin flat plate. The method is useful when a conventional 2-DOF flutter analysis program is at hand and the bridge system and the wing configuration allow for such a generalization and simplification. The motion-induced wind forces on the wings can then be taken into account by modifying the input to such a program. If the wings are shorter than the bridge, a reduction factor is used to approximately account for such a case. Results of example calculations were compared to those obtained from MDOF analysis. Good agreement was found for the no-wing and the full-length wing configurations. This confirms the validity of the quasi-steady flow assumption for the considered example, as was to be expected because a general limit condition for this validity was met. Reasonable agreement was found for a partial-length wing configuration confirming the suitability of the reduction factor.
Finally, uncoupled torsional flutter of bridges with bluff bridge deck was considered and it was shown how it can be treated in a singledegree-of-freedom analysis. The effect of the wings is reflected in the flutter condition equation by an additional damping term established on the assumption of quasi-steady potential flow past a thin flat plate. The method was applied to the first Tacoma Narrows Bridge for various wing configurations. The wind forces on the bridge deck were taken into account by unsteady aerodynamic coefficient values measured in a water tunnel using a sectional bridge deck model. It was found that the flutter speed is substantially raised by the wings. Compared to a bridge with a streamlined deck, but the same structural properties and wing configuration, a greater relative flutter speed increase was achieved. Based directly on the flutter condition equation, it was concluded that for bridges prone to torsional flutter the relative flutter speed increase is the larger, the smaller the structural damping and the (structure-to-air) mass moment of inertia ratio of the bridge.