The role of aerodynamic forces in a mathematical model for suspension bridges

In a fish-bone model for suspension bridges studied by us in a previous paper we introduce linear aerodynamic forces. We numerically analyze the role of these forces and we theoretically show that they do not influence the onset of torsional oscillations. This suggests a new explanation for the origin of instability in suspension bridges: it is a combined interaction between structural nonlinearity and aerodynamics and it follows a precise pattern. This gives an answer to a long-standing question about the origin of torsional instability in suspension bridges.


Introduction
Since the Federal Report [1], it is known that the crucial event causing the collapse of the Tacoma Narrows Bridge was a sudden change from a vertical to a torsional mode of oscillation. Several studies were done on this topic, see [8,11,14] but a full explanation of the origin of torsional oscillations is nowadays still missing; see also the updated monograph [15] and references therein. In two recent papers the onset of torsional oscillations was attributed to a structural instability. In [2] a model of suspension bridge composed by several coupled (second order) nonlinear oscillators has been proposed. By using suitable Poincaré maps, it has been proved that when enough energy is present within the structure a resonance may occur, leading to an energy transfer between oscillators, from vertical to torsional. The results in [2] are purely numerical. We found a similar answer in [3] by analyzing a different mathematical model, named fish-bone. In this model, the main span of the bridge, which has a rectangular shape with two long edges and two shorter edges, is seen as a degenerate plate fixed and hinged between the towers. The midline of the roadway is seen as a beam with cross sections that are seen as rods free to rotate around their barycenters located on the beam. The degrees of freedom are the vertical displacement of the beam y, which is positive in the downwards direction, and the angle θ of rotation of the cross sections with respect to the horizontal position. The roadway is assumed to have length L and width 2 with 2 L. By considering the kinetic energy of a rotating object and the bending energy of a beam, the following system is obtained in [3]:    M y tt +EIy xxxx +f (y+ sin θ)+f (y− sin θ) = 0 0 < x < L, t > 0, where M is the mass of the rod, µ > 0 is a constant depending on the shear modulus and the moment of inertia of the pure torsion, EI > 0 is the flexural rigidity of the beam, f includes the restoring action of the prestressed hangers and the action of gravity. To (1) we associate the boundary-initial conditions For a linear force f the two equations in (1) decouple: this case was studied in [10]. In the nonlinear case, well-posedness of the problem was shown in [6]. For a suitable nonlinear f , in [3] we gave a detailed explanation of how internal resonances occur in (1), yielding instability. The aim of this analysis was purely qualitative and the bridge was seen as an isolated system with no dissipation and no interactions with the surrounding air. In particular, both theoretical and numerical results were given proving that the onset of large torsional oscillations is due to a resonance which generates an energy transfer between different oscillation modes. More precisely, when the bridge is oscillating vertically with sufficiently large amplitude, part of the energy is suddenly transferred to a torsional mode giving rise to wide torsional oscillations. Estimates of the energy threshold for stability were obtained both theoretically and numerically, see Section 2.
Our purpose in [3] was to emphasize the structural behavior of the bridge without inserting any interaction with the surrounding air. This procedure was also followed by Irvine [7, p.176] who comments his own approach by writing: In this formulation any damping of structural and aerodynamic origin has been ignored... We could include aerodynamic damping which is perhaps the most important of the omitted terms. However, this refinement, although frequently of significance, yields a messy flutter determinant that requires a numerical solution.
This comment says two things. First, that it was a good starting point to study (1) as an isolated system. Second, that, in order to have more accurate responses, the subsequent step should be to insert aerodynamic forces in the model. This refinement of the model (1) was also suggested to us by Paolo Mantegazza, a distinguished aerospace engineer at the Politecnico of Milan, and motivates the present paper. In order to better highlight the role of the aerodynamic forces, we do not insert in the model any other external action. We will show, both numerically and theoretically, that the threshold of instability of the system is independent of aerodynamic forces. This suggests a new pattern for the aerodynamic and structural mechanisms which give rise to oscillations in suspension bridges, see Section 6.

One mode approximation of the fish-bone model
We first introduce some simplifications of the model which, however, maintain its original essence and its main structural features. First, up to scaling we may assume that L = π and M = 1. Then, since we are willing to describe how small torsional oscillations may suddenly become larger ones, we use the following approximations: cos θ ∼ = 1 and sin θ ∼ = θ; see [3] for a rigorous justification of this choice. Since our purpose is merely to describe the qualitative phenomenon, we may take EI π L 4 = 3µ π L 2 = 1 although these parameters may be fairly different in actual bridges. For the same reason, the choice of the nonlinearity is not of fundamental importance; it is shown in [2] that several different nonlinearities yield the same qualitative behavior for the solutions. Whence, as suggested by Plaut-Davis [12, Section 3.5], we take f (s) = s + s 3 . Finally, we set z := θ and the system (1) becomes    y tt + y xxxx + 2y(1 + y 2 + 3z 2 ) = 0 0 < x < π, t > 0, To (4) we associate some initial conditions which determine the conserved energy E of the system, that is, Existence and uniqueness of solutions were proved in [6] by performing a suitable Galerkin procedure, see also [3] where more regularity was obtained. The proof is constructive: to (4) we associate the functions and the approximated m-mode system where j = 1, ..., m. Then, suitable a priori estimates allow to prove that ) as m → +∞ for all T > 0. Hence, the functions in (5) approximate the solutions of (4). The error committed when replacing y with y m and z with z m can be rigorously estimated, see [3,Theorem 2].
In what follows we focus our attention to the simplest case m = 1. Then, system (6) reads with some initial conditions If we take ζ 0 = ζ 1 = 0, then the unique solution of (7)-(8) is (y 1 , z 1 ) = (ȳ, 0) withȳ =ȳ(η 0 , η 1 ) being the unique (periodic) solution of We callȳ the first vertical mode with associated energy Since we are interested in the stability of the solution z 1 ≡ 0 corresponding to ζ 0 = ζ 1 = 0, we linearize system (7) around (ȳ, 0). The torsional component of the linearized system is the following Hill equation [5]:ξ We say that the first vertical modeȳ at energy E(η 0 , η 1 ) is torsionally stable if the trivial solution of (11) is stable. By exploiting a stability criterion by Zhukovskii [16], in [3] we obtained the following theoretical estimates. Proposition 1. The first vertical mode y at energy E(η 0 , η 1 ) (that is, the solution of (9)) is torsionally stable provided that Proposition 1 gives a sufficient condition for the torsional stability. The numerical results obtained in [3] show that the threshold of instability could be larger. We quote a couple of them in Figure 1. We plot the solution of (7) with initial conditions for different values of y 1 ∞ . The green plot is y 1 and the black plot is z 1 . For y 1 ∞ = 1.45 no 3 How to introduce the aerodynamic forces into the model?
Even in absence of wind, an aerodynamic force is exerted on the bridge by the surrounding air in which the structure is immersed, and is due to the relative motion between the bridge and the air.
Pugsley [13, § 12.7] assumes that the aerodynamic forces depend linearly on the "cross" derivatives and functions. Similarly, Scanlan-Tomko [14] obtain the following equations satisfied by the torsional angle θ: where I, ζ θ , ω θ are, respectively, associated inertia, damping ratio, and natural frequency. The r.h.s. of (13) represent the aerodynamic force which is postulated to depend linearly on bothθ and θ with A, B > 0 depending on the structural parameters of the bridge. Let us mention that the arguments used in [14] to reach the l.h.s. of (13) have been the object of severe criticisms (see [9]), due to some rough approximations and questionable arguments. Nevertheless, the r.h.s. of (13) is nowadays recognized as a satisfactory description of aerodynamic forces. Following these suggestions, we insert the aerodynamic forces in the 1-mode system (7). We first consider the case where only the cross-derivatives are involved. This leads to the following modified system:   ÿ 1 + 3y 1 + 3 2 y 3 1 + 9 2 y 1 z 2 1 + δż 1 = 0 z 1 + 7z 1 + 9 2 z 3 1 + 27 2 z 1 y 2 1 + δẏ 1 = 0 with δ > 0. As in (12), we take the initial conditions for different values of σ and we wish to highlight the differences, if any, between (7) and (14). For (14) we have no energy conservation; however, let us consider the (variable) energy function Let us now consider the case where also the cross-terms of order 0 are involved. Then, instead of (14) we obtain the system   ÿ 1 + 3y 1 + 3 2 y 3 1 + 9 2 y 1 z 2 1 + δ(ż 1 + z 1 ) = 0 z 1 + 7z 1 + 9 2 z 3 1 + 27 2 z 1 y 2 1 + 3δ(ẏ 1 + y 1 ) = 0 where the coefficient 3 in the second equation comes from the variation of the energy Also for (17) we do not have energy conservation but the function E in (18) better approximates the internal energy. It may be questionable whether to include the last term δy 1 z 1 into E since this term depends on the aerodynamic forces. However, the behavior of E, which we analyze in the next section, does not depend on the presence of this term.

Numerical results
For (14) we first take σ = 1.47 and we modify the aerodynamic parameter δ. To motivate this choice we note that no energy transfer seems to occur for σ below this threshold, furthermore σ = 1.47 is also the numerical threshold found when no aereodynamic force is inserted in the model, see Section 2. In Figures 2 and 3 we plot both the behavior of the solutions (first line) and the behavior of the energy E(t) (second line), for increasing values of δ. The first lines in Figures 2 and 3 should be compared with the second picture in Figure 1 (case δ = 0). We note that, as the aerodynamic parameter increases, the transfer of energy is anticipated but it is not amplified. Quite surprisingly, on the second line we see that the energy E(t) remains almost constant except in the interval of time where the transfer of energy occurs: for increasing aerodynamic parameters δ we observe increasing variations in the energy behavior.
Then we maintain fixed δ = 0.01 and we increase the initial energy, that is, the initial amplitude of oscillation. In Figures 4 and 5 we plot both the behavior of the solutions (first line) and the behavior of the energy E(t) (second line), for increasing values of σ.
It turns out that all the phenomena are anticipated (in time) and amplified (in width) and reach a quite chaotic behavior for σ = 3 where we had to stop the numerical integration at t = 90.  For (17) we take again as initial conditions (15) but with σ ≥ 1.47 so that we are above the energy threshold of instability, see Section 2. In Figure 6 we plot both the behavior of the solution (first line) and the behavior of the energy E(t) (second line) of (17)- (15).
It is quite visible that the instability is further anticipated but now also the amplitude is enlarged. Moreover, the energy increases also in absence of torsional instability: this variation is due to the cross-derivatives since all the other terms appear in the energy (18). The very same behavior is obtained for the internal energy, namely the energy (18) without the last term δy 1 z 1 . We also remark that the energy E fails to follow a regular pattern only in presence of instability. We only quote these numerical results because all the other experiments gave completely similar responses.

Theoretical results
As pointed out by Irvine [7, p.176], the numerical approach is probably the most appropriate to analyze a model which also involves aerodynamic forces. The reason is that a satisfactory stability theory for systems such as (14) and (17) is not available. Nevertheless, some theoretical conclusions can be drawn also for these systems, in the spirit of the results obtained in [3] (see also Section 2), where the main idea is to study the solutions of the systems near the pure vertical modeȳ. Here "pure" means that no interactions with the surrounding are admitted and only the structural behavior of the bridge is considered. Let us explain how some of the results for the isolated system (7) may be extended to (14); one can then proceed similarly for (17).
For system (7), two steps are necessary to define the torsional stability of the unique (periodic) solutionȳ of (9) (the pure mode): (i) we linearize the torsional equation of the system (7) around (ȳ, 0), see (11); (ii) we say that the pure vertical modeȳ at energy E(η 0 , η 1 ) is torsionally stable if the trivial solution of (11) is stable.
We point out that the system (7) is isolated and that (11) is unforced. In this situation, the above steps (i)-(ii) are equivalent to: (I) in the torsional equation of the system (7) we drop all the z 1 -terms of order greater than or equal to one and we replace y 1 withȳ; (II) we say that the pure vertical modeȳ at energy E(η 0 , η 1 ) is torsionally stable if all the solutions of (11) are globally bounded.
If we replace (7) with the system (14), then (i)-(ii) make no sense while (I)-(II) do. A linearization as in (i) would exclude the aerodynamic forces, while acting as in (I) preserves them and gives rise to whereȳ is the unique periodic solution of (9). Needless to say, f and a have the same period. The definition (ii) is inapplicable to (19) since ξ ≡ 0 is not a solution, while (II) is a verifiable property for (19). This is the definition of stability that we adopt for the vertical modeȳ in presence of aerodynamic forces. With this definition we can prove the following statement.
Proposition 2. Let y be the pure vertical mode at energy E(η 0 , η 1 ), that is, the solution of (9). Then y is torsionally stable for (7) if and only if it is torsionally stable for (14).
Assume now thatȳ is torsionally stable for (14). Then ξ is bounded for any choice of ξ h , that is, any choice of the constants A and B in (20). This shows that also ξ h is bounded and proves the stability ofȳ for (7). 2 This statement deserves a couple of straightforward comments.
• In agreement with the numerical results described in Section 4, Proposition 2 shows that the energy threshold for stability does not depend on the strength of the aerodynamic forces; in particular, an isolated system has the same energy threshold.
• By applying Propositions 1 and 2 we infer that the energy threshold for the stability of (14) is at least 235 294 .

Conclusions
It is clear that in absence of wind or external sources a bridge remains still. A vertical load, such as a vehicle, bends the bridge and creates a bending energy. Less obvious is the way the wind inserts energy into the bridge: let us outline how this happens. When a fluid hits a bluff body its flow is modified and goes around the body. Behind the body, or a "hidden part" of the body, the flow creates vortices which are, in general, asymmetric. This asymmetry generates a forcing lift which starts the vertical oscillations of the body. Up to some minor details, this explanation is shared by the whole community and it has been studied with great precision in wind tunnel tests, see e.g. [8,15]. The vortices induced by the wind increase the internal energy of the structure by generating wide vertical oscillations. When the amount of energy reaches a critical threshold our results in [3] show that a structural instability appears: this is the onset of torsional oscillations. The results in the present paper show that, at this stage, the aerodynamic forces excite the internal energy irregularly giving rise to further self-excited oscillations.
The whole energy-oscillation mechanism is here described through a very simplified model which certainly needs to be significantly improved. But, at least qualitatively, we believe that the "true" mechanism in a suspension bridge will follow this pattern: 1) the interaction of the wind with the structure creates vortices; 2) vortices create a lift which starts vertical oscillations of the bridge; 3) when vertical oscillations are sufficiently large, torsional oscillations may appear; 4) the onset of torsional instability is of structural nature; 5) the aerodynamic forces excite the energy only when the structural torsional instability appears; 6) the energy threshold of stability is independent of the strength of aerodynamic forces.