On the discontinuous dynamics of a class of 2-DOF frictional vibration systems with asymmetric elastic constraints

1.   Introduction

Among bridges, suspension bridges are the only usable to cover long distances (above 1 km) and they turn out to be quite unstable. In Figure 1 we sketch the components of a suspension bridge. The deck (or roadway) is supported from below by a girder that may be composed by stiffening trusses. Four towers (or piers) sustain two (or more) parallel cables which, in turn, sustain the hangers. At their lower endpoints the hangers are linked to the deck and sustain it from above.

The idea of suspension bridge appears in 1595, in pictures by the Italian engineer Fausto Veranzio ^[103], although these bridges were never built. The Jacob Creek Bridge, erected in Pennsylvania in 1801 following a project of the Irish engineer James Finley ^[50], is the first suspension bridge ever built. Since the very beginning, many suspension bridges had serious problems under the action of the wind or of traffic loads, see ^[1] and [59,Chapter 1]. The most impressive failure of history is certainly the Tacoma Narrows Bridge (TNB) collapse in November 1940, see ^[101]. But other bridges collapsed in a similar way due to hurricanes. The Brighton Chain Pier (UK) collapsed in 1836 ^[89], the Menai Straits Bridge (UK) collapsed in 1839 ^[88], the Wheeling Suspension Bridge (West Virginia) collapsed in 1854 ^[39], the Matukituki Suspension Footbridge (New Zealand) collapsed in 1977 [71,p.180]. We refer again to [59,Chapter 1] for more details on these events. Prior to all these collapses, unexpected and destructive torsional oscillations suddenly appeared. This naturally raises the question: is there some way to prevent or reduce them?

To answer this question, we need first to understand the phenomena which lead to failures. Nowadays the scientific community agrees that the main culprit for the TNB collapse were the torsional oscillations that suddenly appeared, see [94,p.63]. Just after the collapse, a Board of Experts was appointed for a Federal Report ^[2] that, however, was unable to explain the sudden appearance of torsional oscillations which remained a mystery for long time.

There have been many attempts to explain this phenomenon. The first explanation is through some external resonance but this was immediately discarded with simple arguments, see Section 2.1. Subsequently, Billah-Scanlan [27,p.121] claimed that Scanlan-Tomko ^[93]

...demonstrated conclusively that the catastrophic mode of the old Tacoma Narrows Bridge was a case of what they termed single-degree-of-freedom torsional flutter due to complex, separated flow.

But McKenna ^[83] wrote that Billah-Scanlan ^[27]

...offered a mathematical model which is only valid for very small displacements and can only be verified in ideal wind tunnel experiments of "in torsion $0\leq\alpha\leq \pm 3^o$". We are asked to believe that these "penetrating insights" explain the Tacoma Narrows oscillation. To us, the case is less than convincing.

With some sarcasm McKenna also writes apparently the authors were not familiar with the concept of absolute value and he concludes by saying that the explanation in ^[27]

...is a perfectly good explanation of something that was never observed, namely small torsional oscillations, and no explanation of what did occur, namely a large vertical oscillation with a double amplitude of five ft. and a frequency of 38 per min. followed by a change to the torsional.

Indeed, the Report [2,p.31] shows that torsional oscillations were never recorded prior to the day of the collapse. Scanlan ^[91,92] also wrote that

...the original Tacoma Narrows Bridge withstood random buffeting for some hours with relatively little harm until some fortuitous condition "broke" the bridge action over into its low antisymmetrical torsion flutter mode.

The words "fortuitous condition" denote a lack of rigorous scientific explanation. As a conclusion to this diatribe, the civil engineer Scott ^[94] and McKenna ^[82] recently wrote that

...opinion on the exact cause of the TNB collapse is even today not unanimously shared...

...there is no consensus on what caused the sudden change to torsional motion.

Whence, the mystery of the TNB collapse remained unsolved for about 80 years. All the suggested explanations were of purely aerodynamic type and all of them were discarded either because the quantitative parameters did not fit the theoretical explanations or because the experiments in wind tunnels did not confirm the underlying theory. Irvine [71,p.176] suggested that a full explanation needs also to take into account the nonlinear behavior of structures; see also ^[76]. To this end, a first step is to isolate the bridge from forcing and damping. Following this suggestion, intensive structural studies on different suspension bridge models were performed in the last decade. Let us mention studies through the Melan equation ^[84] (both analytically ^{[38,61,65,68]} and numerically ^{[42,96,97,98]}), through coupled oscillators ^{[7,8,26,41,63]}, beam equations ^{[12,14,22,24,25,26,53,54,56,57,62]}, plate equations with two degrees of freedom ^{[9,23,28,29,37,46,49,64]}; see also the monographs ^[55,59]. These different models and the overall mathematical description of certain phenomena found consensus among engineers ^{[6,32,33,41,42,44,73,78,106,107]}.

The first purpose of this paper is to combine together all these phenomena and to give a new fully detailed and unpublished explanation of the TNB collapse. This requires the definition of three different stability thresholds for each of the vertical modes of oscillation. This is done in Section 2 where we also show that our explanation matches the description of the collapse given in ^[2,94]. Hence, our (qualitative and quantitative) mathematical explanation appears sound and reliable, thereby solving the long standing mystery of the TNB collapse.

But understanding a problem is only the first part of the work, the second part consists in solving the problem. In this respect, a possibility is to introduce a damping parameter. In Section 3 we consider a degenerate plate equation introduced in ^[24,55] and its damped version ^[64]. We set up a finite horizon optimal control problem in which the damping term plays the role of the control. Theorem 3.6 shows that the optimal control does not exist; this result is complemented with some numerics that highlight the role played by the structural parameters. A further possibility to improve the stability of bridges is to find "optimal decks" able to convey the energy from the wind to suitable oscillation modes. This problem is addressed in Section 4 where different shape optimization problems are suggested and partially solved. We believe that a full solution of these problems, both theoretically and practically, would increase the stability of future suspension bridges; in particular, in windy regions such as the Arabian Peninsula and the Sahara Desert.

2.   The conclusive explanation of the TNB collapse

2.1.   Why not a simple resonance

Two days after the TNB collapse, the New York Times published an article ^[10] attempting to describe the causes of the collapse. It says:

Like all suspension bridges, that at Tacoma both heaved and swayed with a high wind. It takes only a tap to start a pendulum swinging. Time successive taps correctly and soon the pendulum swings with its maximum amplitude. So with this bridge. What physicists call resonance was established, with the result that the swaying and heaving exceeded the limits of safety.

This comment is quite puzzling. First, the simplified comparison between a pendulum and a bridge. Second, the pendulum has to be tapped "correctly", did the wind tap correctly the TNB that day? Third, what is the maximum amplitude of swing of a pendulum? Fourth, resonance is introduced as if it was a whim of physicists. Serious doubts about the above explanation were also raised by Lazer-McKenna [77,$\S$ 1] who observe that

the phenomenon of linear resonance is very precise. Could it really be that such precise conditions existed in the middle of the Tacoma Narrows, in an extremely powerful storm?

The engineers Billah-Scanlan ^[27] make a fool of physics textbooks who attempt to explain the TNB collapse by means of a resonance, while also the physicists Green-Unruh ^[69] wrote that

making the comparison to a forced harmonic oscillator requires that the wind generates a periodic force tuned to the natural frequency of the bridge.

The erroneous explanation through resonance has its origin from different accidents occurred in other bridges. The Broughton Suspension Bridge (Manchester) collapsed in 1831 due to an external resonance caused by a troop marching over the bridge in step ^[4]. As a consequence of this incident, it was decided that troops should "break step" when crossing a bridge. But this did not solve the problem. The Pont Suspendu de la Baisse-Chaîne (France) collapsed in 1850 while a battalion of soldiers was marching across it, killing 226 of them. From ^[3] we learn that, even if the soldiers were told to break the step, in order to withstand they went back to step by following the structural frequency of the bridge: this created a (real!) resonance and amplified the oscillations. The brains of the soldiers induced them to do so, but the wind has no brain!

Summarizing, a simple resonance as described in the New York Times, is not the reason of the TNB collapse. In order to understand the real phenomena which lead to collapses, we need to discuss some models.

2.2.   Plate models for bridges

A rectangular plate $\Omega = (0, L)\times(-\ell, \ell)\subset \mathbb{R}^2$ resists to vertical loads exclusively by means of bending. The flexural properties of a plate depend on its thickness, which we denote by $d$, compared with its width $2\ell$ and its length $L$. Since for bridges modeled through plates one has that $2\ell\ll L$, the plate should be classified according to the ratio $2\ell/d$, see [102,$\S$ 1.1]. Modern suspension bridges have quite small $d$: for instance, the Messina Bridge (see Figure 2) was supposed to have $d = 4.68$m, $2\ell = 52$m, $L = 3300$m, $2\ell/L\approx0.016$, $2\ell/d\approx11.11$ so that it may be classified as a thin plate, see again [102,$\S$ 1.1]. On the other hand, the main span of the original TNB had the measures $L\approx853.44$m, $2\ell\approx11.89$m, $d\approx1.33$m, see p.11 and Drawings 2 and 3 in ^[2]. Therefore, $2\ell/d\approx8.94$ and also the TNB may be considered as a thin plate.

We may then view the deck of the bridge as a thin rectangular plate and, for simplicity, we normalize its length to be $L = \pi$:

The boundary constraints (on $\partial\Omega$) of its out-of-plane displacement $u$ should model a plate hinged on its short edges, where the bridge is supported by the ground, and free on its long edges. If an external vertical load $f\in L^2(\Omega)$ acts on the plate, then, according to the Kirchhoff-Love theory of elasticity ^[74,80] (see also ^[52,60]), up to some multiplicative constants the total energy of the plate is given by

Here, $\sigma\in(-1, 1)$ is the Poisson ratio; hence, the quadratic part of the functional $\mathbb E(u)$ is positive. For a mixture of steel and concrete, such as the deck of a bridge, one has $\sigma\approx0.2$.

For the stability analysis of the deck, a crucial role is played by the vibrating modes of the plate and we are thus led to analyze the eigenvalue problem

The following result holds:

Theorem 2.1. (^[49]) The set of eigenvalues of $(2.2)$ may be ordered in an increasing sequence of strictly positive numbers diverging to $+\infty$ and any eigenfunction belongs to $C^\infty(\overline\Omega)$; the set of eigenfunctions is a complete system in a suitable closed subspace of $H^2(\Omega)$.

Moreover, all the eigenvalues and eigenfunctions belong to one of the following families:

where $m\ge1$ are integers, $A_m^k$ and $B_m^k$ ($k\ge1$ being an integer) are explicit real constants, while $\mu_{m, k}$ (resp. $\nu_{m, k}$) are the eigenvalues associated to the eigenfunctions $(2.3)$–$(2.4)$ (resp. $(2.5)$–$(2.6)$).

Hence, all the eigenfunctions of $(2.2)$ have the form

where the $\varphi_{m, i}$'s are linear combinations of hyperbolic and trigonometric sines and cosines, either even or odd with respect to $y$. Since $\ell\ll\pi$, even functions of $y$ are almost constant, whereas odd functions of $y$ have a linear approximation. The two families $(2.3)$–$(2.4)$ (of the kind $\approx c_m^k\sin(mx)$) represent the vertical modes while two families $(2.5)$–$(2.6)$ (of the kind $\approx c_m^ky\sin(mx)$) represent the torsional modes, see the qualitative pictures in Figure 3.

The eigenvalues of $(2.2)$ solve explicit equations and, by taking the parameters of the collapsed TNB, it is possible to compute them numerically ^[49]. It turns out that the $10$ lowest eigenvalues are all vertical, the least vertical mode is positive and corresponds to the eigenvalue $\mu_{1, 1}$ and to the eigenfunction $\approx C\sin x$, the least torsional mode is the 11th and corresponds to the eigenvalue $\nu_{1, 2}$ and to the eigenfunction $\approx Cy\sin x$.

With this static plate model at hand, the subsequent step is to study the dynamics of the plate. This is why, in the next subsection, we describe in detail how the wind inserts energy in the bridge.

2.3.   From drag to lift

When the wind hits the deck of a bridge the flow is modified and goes around the deck. In the (downwind) leeward, the flow creates vortices which are, in general, asymmetric. The same phenomena may be highlighted in lab experiments, see Figure 4.

In most cases, the wind hits the deck transversally with a force named drag $\mathcal{D}$ and the asymmetry of the vortices behind the deck generates an orthogonal force named lift $\mathcal{L}$, which triggers vertical oscillations, see the sketch of a flow around a disk and against a stadium-like cross section of the deck of a bridge in Figure 5.

In the experiments sketched in Figure 5, one can vary the value of the wind speed $W$ hitting the obstacle and observe that the vortex pattern depends on the Reynolds number $\text{Re} = C\, W/ \eta$, in which $\eta > 0$ is the viscosity of the fluid and $C > 0$ is a geometric constant depending on the shape (or aspect ratio) of the obstacle. The pictures in Figure 6 summarize the observations and the symmetry breaking according to the Reynolds number when the obstacle is a disk. In the left picture (${\text{Re}} = 0$) the flow has a double symmetry, both vertical and horizontal, and the regime is that of unseparated flow. In the middle picture ($0 < {\text{Re}} < 40$) there is symmetry breaking in the horizontal direction, with the appearance of a pair of vortices in the wake. In the right picture ($40 < {\text{Re}} < 150$) there is symmetry breaking in the vertical direction leading to a laminar vortex street. For larger ${\text{Re}}$ the picture becomes even more disordered with the appearance of turbulence and chaos, see ^[79]. Whence, the lift appears for ${\text{Re}}\approx40$, see ^[66] for a detailed explanation and ^[67] for a theoretical quantitative analysis of symmetry breaking.

In bridges, the induced vortices generate a lift that triggers the vertical oscillations of the deck. After a transition time, called the Wagner effect ^[105], the vertical oscillations reach a periodic-like motion which is maintained in amplitude by a somehow perfect equilibrium between the input of energy from the wind and the internal dissipation. If the resulting vertical oscillation is sufficiently large then a structural instability appears: there is a sudden transition from vertical to torsional oscillations. This purely qualitative description requires a refinement through a quantitative analysis of the characteristic frequencies, both of the wind and the structure. This may be achieved through the Floquet theory applied to coupled oscillators, see ^{[7,8,26,41,63]}. With this theory one finds the thresholds of instability in dependence of the structural parameters. This will be used in connection with Proposition 2.4 below.

As sketched in Figures 5 and 6 the wind hitting the deck gives rise to periodic vortices whose frequency increasingly depends on the wind velocity. This law also depends on the so-called aspect ratio of the cross section, namely a number that can be computed from the shape of the section. Once the bridge is given (and, therefore, also the aspect ratio), from the European Eurocode1 ^[40] we learn that there are reliable ways to determine the "explicit" law. Summarizing,

More precisely, the displacements $u = u(x, , y, t)$ of the rectangular deck $\Omega$ are governed by an evolution plate equation such as

see ^{[9,23,28,29,37,46,49,64]} (and Section 2.2), complemented with initial-boundary conditions. In ($2.8$), $\gamma > 0$ is a geometric parameter depending on the shape of the cross section and on the Poisson constant of the material composing the deck, $\delta > 0$ is a damping parameter representing the structural friction (dissipation), while $f(u) = f(u, u_{xx})$ is a (possibly nonlocal) nonlinearity that takes into account the behavior of structures ^[76]. Finally, the right hand side in ($2.8$) is the lift force $\mathcal{L} = \mathcal{L}(t)$ due to the vortices with $W$ and $\omega$ being related by ($2.7$): in Eurocode1 ^[40] the amplitude of the oscillation is considered to be the squared wind velocity whereas the engineering literature (see e.g., [99,Chapter 6]), assumes that the lift force varies periodically in time with the same frequency governing the vortex shedding, this justifies the choice $\mathcal{L}(t) = W^2 \sin(\omega t)$.

2.4.   Prevailing vertical mode

From [2,p.20,p.28] we learn that, usually, in an oscillating suspension bridge, there is a prevailing mode of oscillation. In other words, the solution of ($2.8$) is similar to a multiple of one of the fundamental (linear) modes of the deck; see Section 2.2. This is why, in ^[54,55,56], the prevailing mode was characterized as the mode holding the largest part of the energy in the Fourier decomposition of the oscillation, resulting after solving ($2.8$). As shown in ^[23,49], the lowest frequencies (corresponding to physical wind velocities), are all associated to vertical modes and this is why, originally, a vertical mode is prevailing. So, let us introduce a criterion able to predict which will be the prevailing mode of oscillation.

Initially, the deck is still and the wind starts its transversal action on the deck. For some time, the oscillation of the deck remains small. This suggests to linearize $(2.8)$. By taking into account a possible prestressing $P\ge0$ in the $x$-direction, this leads to consider the linear problem

From [46,Theorem 7], we learn that both the torsional and the vertical skew-symmetric components of the solution are zero. Therefore, we may write the solution of (2.9) as

where $w_k$ are the symmetric vertical eigenfunctions, namely ($2.3$) and ($2.4$) in Theorem 2.1 with $m$ odd. Denote by $\lambda_k$ the eigenvalue of $(2.2)$ associated with $w_k$. Let

then the Fourier coefficients $S_k(t)$ satisfy the ODE

The explicit solutions of (2.10) are given by

and are composed by a damped part (multiplying the negative exponential) and a linear combination of trigonometric functions. In fact,

Hence, the parameter measuring the amplitude of each of the $S_k$'s is

But we also need to take into account the size of the $w_k$'s, originally normalised in $L^2$: therefore, the amplitude of oscillation of each mode is

It is readily seen that

In Figure 7 we represent the functions $A_1$, $A_3$, $A_5$, $A_7$, as in (2.11), for $\omega\in(0, 60)$ and for $P = 0$.

It turns out that these functions all have a steep spike close to their maximum and are fairly small elsewhere, of several orders of magnitude less. The peak (maximum) represents a kind of "structural resonance" called lock-in in the engineering literature, see e.g., ^[33]. The height of the spikes is decreasing with respect to $k$ and the maximum is moving to the right (larger $\omega$).

We are now in position to explain how the prevailing mode is determined by the wind: for a given $\omega > 0$, the prevailing mode $w_k$ is the one maximizing $A_k(\omega)$. For each $\omega > 0$, numerical experiments in ^[29] determine which $k$ maximises $A_k(\omega)$ and show that the prestressing constant $P$ in (2.9) does not influence significantly the prevailing mode. Combined with ($2.7$), this shows that the wind velocity determines the prevailing vertical mode. In fact, as shown in ^[29], there is a whole interval of wind velocities $W$ for each vertical mode $v$:

Proposition 2.2. Let $v$ be a vertical mode, corresponding to some eigenfunction of $(2.2)$ of the kind $(2.3)$ or $(2.4)$.

$(i)$ There exist $0 < L_v < H_v$ (depending on $v$) such that $W\in(L_v, H_v)\Longleftrightarrow\mathit{\mbox{the prevailing mode is }}v$.

$(ii)$ There exists $A_v > 0$ such that, as $W$ varies in the interval $(L_v, H_v)$, the amplitude of oscillation of $v$ varies in $(0, A_v]$.

In particular, Proposition 2.2 states that the mode $v$ cannot oscillate with an amplitude larger than $A_v$.

2.5.   Thresholds and frequencies leading to torsional instability

The classical elastic/plastic principle states that when an elastic body becomes plastic it does not come back to its original shape (or position) and it may crash. Although a suspension bridge is quite flexible, due to this principle it cannot withstand too large energy inputs. For each mode there exists a maximal amplitude of oscillation below which the bridge maintains its elastic behavior and above which the bridge collapses. For a given vertical mode $v$ let us call this maximum amplitude the plastic threshold and let us denote by $P_v$ this threshold. Therefore,

Proposition 2.3. Each vertical mode $v$ has its own plastic threshold $P_v > 0$.

An isolated dynamical system has a constant energy split on its vibrating modes (Fourier coefficients). However, the energy itself may move from one mode to another. The energy transfer may be studied by using the (linear) Floquet theory for differential equations with periodic coefficients and the Hill equations ^[36,104]. For some bridge models this phenomenon has been studied in a large number of papers ^{[7,8,9,12,23,24,25,26,29,37,41,46,49,53,54,56,57,64]}. It turns out that the transfer of energy from vertical to torsional modes occurs when the former reach a given amplitude of oscillations. This happens for each vertical mode $v$: we call torsional instability threshold this critical amplitude and we denote it by $T_v$. Whence,

Proposition 2.4. Each vertical mode $v$ has its own torsional instability threshold $T_v > 0$.

The next result explains why some oscillations are not visiblle in suspension bridges and why not all the vertical oscillations may switch to torsional oscillations.

Proposition 2.5. Let $v$ be a vertical mode, corresponding to some eigenfunction of $(2.2)$ of the kind $(2.3)$ or $(2.4)$ and assume that $W\in(L_v, H_v)$ so that $v$ is prevailing, see Proposition 2.2. Let $A_v$, $P_v$, $T_v$ be, respectively, as in Propositions 2.2, 2.3, 2.4. Then:

$(i)$ if $A_v < \min\{P_v, T_v\}$, then the vertical mode $v$ can neither lead to a collapse, nor switch to a torsional oscillation;

$(ii)$ if $T_v < \min\{A_v, P_v\}$, then the vertical mode $v$ cannot lead to a collapse and switches to a torsional oscillation as soon as $W$ generates a vertical oscillation with amplitude larger than $T_v$;

$(iii)$ if $P_v < \min\{A_v, T_v\}$, then the vertical mode $v$ leads to a collapse as soon as $W$ generates a vertical oscillation with amplitude larger than $P_v$.

Remark 2.6. The minimal request for the design of a bridge is that $A_v < P_v$ for any vertical mode $v$, which prevents the bridge to collapse due to a "simple" vertical oscillation; this fact is well emphasized in the engineering literature ^[1,6,76,87]. Hence, in most bridges, case $(iii)$ does not occur for any vertical mode $v$. We are so left with cases $(i)$ and $(ii)$ and which of these two cases occurs usually depends on the vertical mode $v$ considered.

Remark 2.7. The Report [2,p.31] shows that torsional oscillations were never recorded at the TNB prior to the day of the collapse. In all the bridge collapses due to a torsional oscillation, the involved torsional mode was the second one, see the left picture in Figure 3. According to the detailed analysis on the TNB by Smith-Vincent [100,p.21], this is also the only torsional oscillation ever seen in bridges or on their models ("The only torsional mode which developed under wind action on the bridge or on the model is that with a single node at the center of the main span"). The reason for this is explained in [59,Section 5.3.3]: due to the nonlocal effect of the sustaining cables, it is the least energy torsional configuration.

In the next subsection we use Proposition 2.5 and these remarks in order to interpret all the events seen during the TNB collapse.

2.6.   Why is this explanation consistent

Thanks to the videos available on the web ^[101], most people have seen the spectacular collapse of the TNB and the torsional oscillations were considered the main cause of this dramatic event ^[2,94]. With a bullet $\bullet$ we recall here the phenomena observed, whereas $©$ contains our comments showing that they can be explained by means of the discussion in the previous subsections.

$\bullet$ From [2,p.20,p.28] we learn that, in the months prior to the collapse,

...one principal mode of oscillation prevailed... the modes of oscillation frequently changed... seven different motions have definitely been identified on the main span... from the simplest, that of no nodes, to the most complex, that of seven nodes .

$©$ The first part of this description shows that all the vertical oscillations of the TNB were close to a pure mode, which justifies the introduction of a rigorous mathematical definition of prevailing mode ^[54,55,56]; roughly speaking, it is the mode containing the largest part of the energy, see Subsection 2.4. Besides the obvious remark that from no nodes to seven nodes there are eight (and not seven!) different motions, the other parts of this description tell us that vertical oscillations on higher modes were never observed before the day of the collapse and, hence, only the vertical eigenfunctions of the kind ($2.3$) and ($2.4$) with $m\le8$ played the role of prevailing mode. This means that $A_v < T_v$ for these modes, corresponding to case $(i)$ in Proposition 2.5, recall that Remark 2.6 rules out case $(iii)$.

$\bullet$ In view of the above observation, the transition from vertical to torsional oscillations was unexpected and [2,p.31] reports that

Prior to 10:00 A.M. on the day of the failure, there were no recorded instances of the oscillations being otherwise than the two cables in phase and with no torsional motions.

$©$ This completes the previous description and shows that the vertical modes ($2.3$) and ($2.4$) of the TNB with $m\le8$ satisfied $A_v < T_v$ and could not reach an amplitude of oscillation yielding the transition to a torsional mode. But Propositions 2.2 and 2.4 state that other vertical modes $v$ (with $m\ge9$) may have thresholds satifying $A_v > T_v$.

$\bullet$ On the day of the collapse, the bridge was oscillating vertically as many other times before. According to Eldridge [2,V-3], a witness of the TNB collapse,

the bridge appeared to be behaving in the customary manner... these motions, however, were considerably less than had occurred many times before.

$©$ We believe that "customary manner" means that only vertical oscillations were involved. Moreover, "considerably less" certainly refers to the amplitude of oscillations since the energy is not visible: recall that also the mode (and not only the amplitude) determine the energy. Following Proposition 2.2, that day the wind velocity $W$ fell inside some interval $(L_v, H_v)$ and generated a vertical prevailing mode $v$ with an amplitude of oscillation smaller than $A_v$ (and smaller than the amplitude of any vertical oscillation previously seen at the TNB).

$\bullet$ Then a sudden change in the motion was alarming, a violent destructive torsional movement started. A witness to the collapse was Farquharson, the man escaping in the video ^[101]. According to his detailed testimony in ^[45],

a violent change in the motion was noted. This change appeared to take place without any intermediate stages and with such extreme violence that the span appeared to be about to roll completely over.... The motion, which a moment before had involved a number of waves (nine or ten) had shifted almost instantly to two.

$©$ The "violent change" has to be attributed to a loss of stability, namely the exact lock-in instant: the amplitude of the ninth (or tenth) vertical mode generated a frequency identical to the second torsional mode. That day the wind velocity $W$ fell into an interval $(L_v, H_v)$ associated with a vertical mode $v$ of the kind of ($2.3$) and ($2.4$) with $m = 9$ or $m = 10$. For these modes one has $T_v < A_v$ and, hence, they belong to case $(ii)$ in Proposition 2.5 (recall Remark 2.6). At some time, the wind velocity $W$ generated a vertical oscillation of amplitude larger than $T_v$ and, without any intermediate stages, the oscillation became torsional, see again Proposition 2.5. Why this new born oscillation was on the second torsional mode is explained in Remark 2.7. We also point out that the results obtained in ^[23,29] show that the 9th and 10th mode are more prone to switch to torsional oscillations because their thresholds $T_v$ are smaller than those of lower modes.

$\bullet$ McKenna ^[83] wrote that

... if the explanation in ^[27] has any validity, why were small torsional oscillations never observed? After all, the bridge was known to have oscillated vertically in winds of 3 m.p.h., and remained motionless in winds of 35 m.p.h., (when, according to ^[27], "divergent amplitudes" are reached). It is also worth noting that the bridge had survived winds of 48 m.p.h. without undergoing torsional oscillations, ^[2], page 28.

$©$ Small torsional oscillations were never observed because when the amplitude of the vertical oscillation mode $v$ reaches the torsional instability threshold $T_v > 0$ (see Proposition 2.4) a parametric resonance occurs and almost all the energy suddenly transfers to torsional. Winds of 3 m.p.h. triggered visible vertical oscillations because the resulting frequency $\omega$ of the vortices fell into intervals containing the spikes in Figure 7. Winds of 35 m.p.h. did not trigger visible vertical oscillations because the resulting frequency $\omega$ of the vortices fell into the intervals far away from the spikes in Figure 7. Winds of 48 m.p.h. did not create torsional oscillations because the resulting frequency $\omega$ of the vortices triggered a vertical mode $v$ for which $A_v < \min\{P_v, T_v\}$, see case $(i)$ in Proposition 2.5.

The detailed description of the TNB collapse given in this section shows the necessity of finding devices in order to improve the stability of bridges. This leads to some natural shape optimization problems that we analyze in the next sections.

3.   An optimal control problem for suspension bridges

3.1.   A degenerate plate model with intermediate piers

The plate model introduced in Section 2.2 fails to take into account two essential features of suspension bridges. On the one hand, a suspension bridge is usually divided in three adjacent spans, separated by two piers (or towers) where the deck is hinged, see Figure 1. On the other hand, the rectangular plate $\Omega$ in ($2.1$) has edges with huge discrepancies ($\ell\ll\pi$), of about two orders of magnitude: this means that, in fact, the plate is quite similar to a beam. But since beams do not display torsional oscillations, they also fail to model a suspension bridge, see the left picture in Figure 8 in which a torsional oscillation of the central span is sketched. A possible compromise is then to consider a degenerate rectangular plate $\Omega = (0, \pi)\times(-\ell, \ell)$ composed by a central beam in position $(0, \pi)\times\{0\}$ and by a continuum of cross sections of length $2\ell$ free to rotate around their center placed on the beam. Assuming that the piers are in the symmetric positions $x = \pm a\pi$ (for some $a\in(0, 1)$), the view from above of this degenerate plate and the position of the piers are drawn in the right picture of Figure 8 in which the beam is the white midline and some cross sections are represented by the orthogonal black thin segments. This model was called a fish-bone in ^[24], whereas the piers were introduced in ^[55,56].

While the vertical displacement $u$ can be damped by stiffening the beam or with stronger elastic connections to the ground at $x = \pm\pi$ and to the piers at $x = \pm a\pi$, there is no simple way to damp the torsional displacement $\theta$ since the endpoints of the cross sections at $y = \pm\ell$ are free to move, see the boundary conditions in ($2.2$). This leads to an initial value problem for the following system of PDE's:

where $I = (-\pi, \pi)$ and $\delta > 0$, complemented with the boundary-interior conditions

Compared with ($2.8$), the system $(3.1)$ is degenerate since the vertical displacement $u$ is governed by the beam-type equation $(3.1)$$_1$ while the torsional angle $\theta$ is governed by the wave-type equation $(3.1)$$_2$. A full physical justification of $(3.1)$–$(3.2)$ is given in ^[55,56]. Due to the presence of the constraints at $x = \pm a\pi$, the solutions of $(3.1)$ fail to be smooth and a suitable notion of weak solution is needed.

Definition 3.1. Endowed with the scalar products $(u, v)_V = \int_Iu"v"$ and $(u, v)_W = \int_Iu'v'$, consider the Hilbert spaces

and let $\langle\cdot, \cdot\rangle_V$ and $\langle\cdot, \cdot\rangle_W$ be the duality pairings in $V(I)$ and $W(I)$. We say that the functions

are weak solutions of $(3.1)$–$(3.2)$ if

for all $(\varphi, \psi)\in V(I)\times W(I)$ and all $t > 0$.

We complement ($3.4$)–($3.5$) with some initial conditions

assuming that $u_0 \in V(I)$, $\theta_0 \in W(I)$, $u_1, \theta_1 \in L^2(I)$. The well-posedness of this problem, may be obtained through the standard Galerkin method.

Proposition 3.2. ([55,Proposition 4.1]) For all $u_0\in V(I)$, $\theta_0\in W(I)$, $u_1, \theta_1\in L^2(I)$, there exists a unique weak solution $(u, \theta)$ of $(3.4)$-$(3.5)$-$(3.6)$. Moreover, $u\in C^2(\overline{I}\times \mathbb{R}_+)$ and $u_{xx}(-\pi, t) = u_{xx}(\pi, t) = 0$ for all $t > 0$.

Associated to weak solutions, we have the natural energy

The "partially damped" feature of $(3.1)$ raises several natural questions on the long-time behavior:

These questions have the following answers.

Theorem 3.3. (^[64]) For all $u_0\in V(I)$, $\theta_0\in W(I)$, $u_1, \theta_1\in L^2(I)$, the vertical component of the corresponding solution $(u, \theta)$ of $(3.4)$-$(3.5)$-$(3.6)$ satisfies

Moreover, there exists $\beta > 0$ such that if $u_0\in V(I)$, $\theta_0\in W(I)$, $u_1, \theta_1\in L^2(I)$ are small enough in such a way that

then there exist $C, \eta > 0$ such that the vertical component of the corresponding solution $(u, \theta)$ of $(3.4)$-$(3.5)$-$(3.6)$ satisfies

If $(\theta_0, \theta_1)\neq(0, 0)$, then the torsional component $\theta = \theta(t)$ of the corresponding solution $(u, \theta)$ of $(3.4)$-$(3.5)$-$(3.6)$ satisfies

Finally, if $(u_0, u_1)\neq(0, 0)$, then the energy $E(t)$ in $(3.7)$ is strictly decreasing, while if $(\theta_0, \theta_1)\neq(0, 0)$, then $E(t)\to E_\infty > 0$ as $t\to\infty$.

3.2.   The role of the damping term as a control parameter

By combining Theorem 3.3 with some numerical results and by comparing $(3.1)$ with its undamped version ($\delta = 0$), the conclusion in ^[64] was that partial damping can lead to disasters, namely that partial damping may enlarge the torsional oscillations even in absence of external forces. In this subsections we go a step further with the analysis of $(3.4)$-$(3.5)$-$(3.6)$, by considering $\delta$ as a control parameter and by wondering if it can be optimized in order to reduce the torsional oscillations.

To set up rigorously this optimal control problem, let us first rewrite the energy $E$ in ($3.7$) as

where $E_u$, $E_\theta$, $E_c$ denote, respectively, the vertical energy, the torsional energy, the coupling energy.

A straightforward consequence of Theorem 3.3 reads

Corollary 3.4. Let $u_0\in V(I)$, $\theta_0\in W(I)$, $u_1, \theta_1\in L^2(I)$, and let $(u, \theta)$ be the corresponding solution of $(3.4)$-$(3.5)$-$(3.6)$. Then the following limits exist and

Moreover, if $(u_0, u_1) = (0, 0)$ then $E_\infty = E(0)$, while if $(\theta_0, \theta_1) = (0, 0)$ then $E_\infty = 0$.

The last statement may be obtained after noticing that, by Proposition 3.2 (uniqueness), $(u_0, u_1) = (0, 0)$ implies $u(t)\equiv0$ while $(\theta_0, \theta_1) = (0, 0)$ implies $\theta(t)\equiv0$.

From a physical point of view, one is mainly interested in a fixed horizon optimal control problem, see e.g., ^[31]. More precisely, one fixes some finite $T > 0$ and considers $(3.1)$ only on the interval $[0, T]$. Then the optimal control problem reads

This means that if the wind inserts an energy $E(0)$ within the bridge, we aim to determine an optimal partial damping $\delta\ge0$ such that, regardless of the initial data $(u_0, u_1, \theta_0, \theta_1)$ whose global energy equals $E(0)$, the fixed horizon energy $E(T)$ is minimal. We set up this problem by introducing the closed bounded manifold ${\mathcal M} = {\mathcal M}\big(E(0), T\big)$ defined by

We then define

and address the following

Problem 3.5. Is the infimum in $(3.12)$ attained? That is, does there exist initial data $(u_0, \theta_0, u_1, \theta_1)\in {\mathcal M}$ such that the corresponding energy $(3.7)$ satisfies $E(T) = \Gamma_T$?

The numerical results in Section 3.3 suggest that the answer to Problem 3.5 might be positive, see Remark 3.8. In Section 3.4 we prove

Theorem 3.6. Let $u_0\in V(I)$, $\theta_0\in W(I)$, $u_1, \theta_1\in L^2(I)$ be such that $E_u(0)+E_\theta(0)+E_c(0) = E(0)$ for some $E(0) > 0$. Let $\big(u(t), \theta(t)\big)$ be the corresponding weak solution of $(3.4)$-$(3.5)$-$(3.6)$, with the associated energies as in $(3.7)$ and $(3.10)$. Let $T > 0$ be fixed.

If $\delta = 0$, then $E(T) = E(0)$.

If $\delta > 0$, then:

$\bullet$ for any $\alpha\in(\Gamma_T, E(0)]$ there exist initial data $(u_0, \theta_0, u_1, \theta_1)\in {\mathcal M}$ such that $E(T) = \alpha$;

$\bullet$ $\Gamma_T\to0$ as $T\to\infty$.

Theorem 3.6 shows that the optimal control problem ($3.11$) has no solution, there exist initial data with nontrivial torsional component for which the energy $E(T)$ is close to $E(0)$. Hence, as soon as the bridge has some torsion (even very small!), it may collapse.

3.3.   Numerics for bimodal solutions

Consider the two eigenvalue problems

The sequences of the associated eigenfunctions form an Hilbertian basis of, respectively, $V(I)$ and $W(I)$. Let $e_\lambda$ be an $L^2$-normalized eigenfunction of $(3.13)$$_1$ related to the eigenvalue $\lambda$ and let $\eta_\kappa$ be an $L^2$-normalized eigenfunction of $(3.13)$$_2$ related to the eigenvalue $\kappa$. There is a natural coupling between these modes and, for some of them, the presence of the piers yields an additional coupling measured by the coefficient

There exist couples $(\lambda, \kappa)$ for which $A_{\lambda, \kappa} = 0$ and couples for which $A_{\lambda, \kappa}\neq0$, see ^[56]. If $A_{\lambda, \kappa} = 0$, then the space $\langle e_\lambda\rangle\times\langle\eta_\kappa\rangle$ is invariant. This means that if, for some real numbers $c_1, c_2, c_3, c_4$, we take initial data such as

then the solution of $(3.4)$-$(3.5)$-$(3.6)$ has the form

We call solutions such as ($3.15$) bimodal solutions.

For the collapsed TNB it was $a = 14/25$ so that the ninth vertical eigenvalue was $\lambda\approx633$ whereas the second torsional eigenvalue was $\kappa\approx3.189$, see ^[56]. These eigenvalues correspond to the oscillations seen the day of the collapse ^[2], as recalled in Section 2.6 and in Remark 2.7; in this case, we have $A_{\lambda, \kappa} = 0$, see again ^[56]. Then we take initial data as ($3.14$) with $(c_3, c_4)\neq(0, 0)$ so that the solution of $(3.4)$-$(3.5)$-$(3.6)$ has the form ($3.15$) with $(w, z)$ solving the system (for $t\ge0$)

From Theorem 3.3 we know that $w(t)\to0$ as $t\to\infty$ (vanishing vertical component) and that $z(t)\not\to0$ (nonvanishing torsional component). We numerically studied system ($3.14$)–($3.16$) by varying the involved parameters. In all the experiments below we started with null kinetic energy, namely, $\dot{w}(0) = \dot{z}(0) = 0$.

We first fixed $S = 50$, $\delta = 1$, $z(0) = 0.1$ and we varied the initial longitudinal amplitude of oscillation $w(0)$, obtaining the plot in the left picture of Figure 9. The results were identical by taking $\delta = 0.1$, the only change being the speed of convergence of $E(t)$ towards $E_\infty$. The plot shows that the residual energy $E_\infty$ increases as the initial energy increases, but $E_\infty$ always remains smaller than some threshold ($\approx0.285$) even for huge initial energies $E(0)$. This leads to the following

Conjecture 3.7. The residual energy $E_\infty$ of $(3.16)$, which is pure torsional energy in view of Theorem 3.3, is always less than some $\overline{E} > 0$, regardless of $w(0)$.

If this conjecture were true, this would add a further conclusion to Section 2, that is, although the switch from vertical to torsional oscillations asymptotically breaks down the energy of several orders of magnitude, the remaining torsional energy is enough to make the bridge collapse.

Then we fixed $w(0) = 5$, $\delta = 1$, $z(0) = 0.1$ and let $S$ vary. As expected, the residual energy increases with $S$ since this increases both the nonlinearity and the interaction between the two equations in ($3.16$). Starting from $S = 0$ (linear decoupled system), we increased the strength of the interaction until $S = 100$ and we obtained the plot in the right picture of Figure 9. It turns out that, again, the residual energy remains bounded ($E_\infty\le0.26$), much smaller than the initial energy $E(0) = \frac{633}{2}w(0)^2+\frac{3.189}{2}z(0)^2+\frac{25}{4}\big(w(0)^4+6w(0)^2z(0)^2+z(0)^4\big)\approx11821$.

We also fixed $S = 50$, $w(0) = 5$, $z(0) = 0.1$, and we varied $\delta$, obtaining the results in Table 1 where we see that, for increasing $\delta$, there is no clear trend. Sometimes the initial torsional energy increases of a factor 10 or 20, or even 30. For $\delta = 50$ (huge) it even increases of a factor 100. Our feeling is that larger $\delta$ transfer more quickly the energy from vertical to torsional, but the amount also depends on the possible phase displacement between the two oscillations.

Remark 3.8. Theorem 3.3 suggests that a possible minimizer for Problem 3.5 has to be sought among purely vertical solutions, namely for $(u_0, \theta_0, u_1, \theta_1)\in {\mathcal M}$ with $\theta_0 = \theta_1 = 0$. Moreover, the above numerical results suggest to focus the attention on unimodal solutions, namely solutions of $(3.4)$-$(3.5)$-$(3.6)$ in the form $\big(u(x, t), \theta(x, t)\big) = \big(w(t)e_\lambda(x), 0\big)$ with $\lambda$ being an eigenvalue of $(3.13)$$_1$. These solutions appear when both $u_0$ and $u_1$ are proportional to $e_\lambda$, which leads to the following damped Duffing ODE:

for some $c_1, c_2\in \mathbb{R}$. We numerically studied this problem and found that $E(T)$ was increasing with $\lambda$. This leads to the conjecture that the infimum in $(3.12)$ is attained by initial data $(c_1e_{\lambda_1}, 0, c_2e_{\lambda_1}, 0)\in {\mathcal M}$, where $e_{\lambda_1}$ is a normalized eigenfunction corresponding to the least eigenvalue of $(3.13)$$_1$.

3.4.   Proof of Theorem 3.6

By using the Galerkin approximation, one sees that $t\mapsto E(t)$ is differentiable and

We point out that ($3.17$) is not straightforward due to the fact that weak solutions of $(3.4)$–$(3.5)$ are nonsmooth and one cannot use as test functions $u_t$ and $\theta_t$; this is why one needs to use the finite-dimensional Galerkin approximation.

By ($3.17$), the energy $E(t)$ is constant if $\delta = 0$, which proves the first statement in Theorem 3.6.

The case $\delta > 0$ is much more complicated. First, we need a continuous dependence result.

Lemma 3.9. Let $(u_0, \theta_0, u_1, \theta_1)\in V(I)\times W(I)\times L^2(I)\times L^2(I)$ and let $\{(u_0^n, \theta_0^n, u_1^n, \theta_1^n)\}_n\subset V(I)\times W(I)\times L^2(I)\times L^2(I)$ be a sequence such that

Let $(u, \theta)$ and $(u^n, \theta^n)$ be the corresponding weak solutions, with associated energies $(3.7)$ $E = E(t)$ and $E^n = E^n(t)$ such that $E^n(0) = E(0)$. Then $E^n\to E$ uniformly in $[0, T]$ for all $T > 0$.

Proof. We subtract from both ($3.4$) and ($3.5$) the corresponding equations satisfied by $(u^n, \theta^n)$ and we obtain

for all $(\varphi, \psi)\in V(I)\times W(I)$ and all $t > 0$. We set $w^n = u-u^n$ and $z^n = \theta-\theta^n$; then we choose $\varphi = w^n_t(t)$ in (3.19) and $\psi = z^n_t(t)$ in (3.20). In fact, this choice is not legitimate since $w^n_t$ and $z^n_t$ do not have the needed regularity; however, this may be justified by using (again!) the finite dimensional Galerkin approximation and by letting the dimension tend to infinity. In order to avoid heavy notations (we would need four approximating sequences for $u$, $u^n$, $\theta$, $\theta^n$), we proceed formally knowing that each step can be rigorously justified.

With this choice of $\varphi$, and with iterated use of the Hölder inequality, (3.19) yields

where we omitted the dependence on $t$ in the norms involved in the upper bounds. By ($3.17$), we know that both $E(t)\le E(0)$ and $E^n(t)\le E^n(0)$ for all $t\ge0$. Therefore, the four quantities

are uniformly bounded for $t\ge0$. Hence, there exists $K > 0$ such that the previous inequality implies that

With the above choice for $\psi$ and by arguing in a completely similar way, from (3.20) we infer that (for a possibly different $K > 0$)

By setting $F_n(t): = \|w^n_t(t)\|_2^2+\|w^n_{xx}(t)\|_2^2+\|z^n_t(t)\|_2^2+\|z^n_x(t)\|_2^2$, by adding the two above inequalities, and by using the Young and the Poincaré inequalities, we infer that (for a possibly different $K > 0$)

Upon integration over $(0, t)$, this yields $F_n(t)\le F_n(0)e^{Kt}$ for all $t\ge0$. Since ($3.18$) implies that $F_n(0)\to0$ as $n\to\infty$, this shows that $F_n\to0$ in $L^\infty(0, T)$ for any $T > 0$. By recalling the expression of $F_n$, this implies that

In turn, by the Poincaré inequality, we also have that

All together, these convergence prove that $E^n\to E$ uniformly in $[0, T]$ for all $T > 0$.

The next step is to combine Corollary 3.4 with ($3.17$) to infer that the initial data $(0, \theta_0, 0, \theta_1)\in {\mathcal M}$ yield a constant energy ($3.7$), that is, $E(T) = E(0)$. Let $\alpha\in(\Gamma_T, E(0)]$; then there exists $(u_0^\alpha, \theta_0^\alpha, u_1^\alpha, \theta_1^\alpha)\in {\mathcal M}$ such that the corresponding energy $E^\alpha$ satisfies $E^\alpha(T) = \overline{\alpha}\in(\Gamma_T, \alpha)$. Consider a continuous curve

and, for any $s\in[0, 1]$, consider the weak solution $(u_s, \theta_s)$ of ($3.4$)-($3.5$) with initial data $r(s)$. This solution has an energy $E^s = E^s(t)$, as in ($3.7$). By Lemma 3.9, the map $s\mapsto E^s(T)$ is continuous over $[0, 1]$ so that it takes any intermediate value $\alpha\in[\overline{\alpha}, E(0)]$ for some $s$. This proves the first item for $\delta > 0$ in Theorem 3.6.

To prove the second item, take initial data $(u_0, 0, u_1, 0)\in {\mathcal M}$ so that, by Corollary 3.4, we have $E_\infty = 0$ and there exists $\tau > 0$ such that $E(\tau)\le\beta$, see ($3.8$). By Theorem 3.3 we know that there exists $C, \eta > 0$ such that

yielding, in particular, $\Gamma_T\le E(T)\le C e^{-\eta T}$ whenever $T\ge\tau$; therefore $\Gamma_T\to0$ as $T\to\infty$.

4.   Shape optimization problems for future designs

4.1.   Minimization of the drag and lift

We introduce here a 3D model able to describe the wind forces acting on the bridge. As in Section 2.2 we view the deck of the bridge as a thin plate defined by

Compared with the 2D rectangular plate $\Omega$ in ($2.1$), $D$ has a thickness $d > 0$. We consider the approximation of an unbounded region surrounding the deck

The domains $D$ and $\Omega$, as well as their intersections $K = (-\ell, \ell) \times (-d, d)$ and $\Sigma = (-L, L)^2 \setminus \overline{K}$ with the plane $x = 0$, are represented in Figure 10 (not in scale).

This particular rectangular shapes for $D$ and $K$ have been chosen to explain the model. It is now our purpose to vary these shapes and to set up some shape optimization problems. We assume that the motion of the wind is governed by the 3D Navier-Stokes equations

where $\eta > 0$ is the kinematic viscosity and $W$ is the transversal velocity of the wind (in the $y$-direction).

The rate of strain tensor $\sigma$ and the stress tensor $\mathbb{T}$ are given by

where $\mathbb{I}$ is the identity $3 \times 3$ matrix. In a viscous fluid, the total force exerted by the fluid over the obstacle $D$ is given by the vector field

where the minus sign is due to the fact that the outward unit normal $n$ to $\Omega$ is directed towards the interior of $D$. Since $\partial D$ is nonsmooth, a generalization of ($4.4$) is needed to define $F_D$, see ^[51,66].

Aiming to prevent collapses, the first target is to minimize the drag in order to lower the energy input within the structure; from the engineering point of view, this amounts to seek optimal aerodynamic shapes. If the inflow is transversal, as in $(4.3)$, the drag force is the horizontal component of $F_D$. By ^[13], for smooth obstacles $D \subset \mathbb{R}^3$ the drag force may also be written as

The explicit form of $F_D$ can be used to seek the optimal shape of the body $D$ (with given volume) that minimizes the drag. The following necessary condition for the obstacle $D$ to minimize the drag holds.

Theorem 4.1. (^[85,86]) Assume that the boundary datum $W$ is sufficiently small, so that there exists a unique solution $u$ of $(4.3)$, and consider the linear problem

Then the body $D_0$ of given volume minimizing the drag, if any, satisfies

The results in ^[13] prove the differentiability of the drag with respect to the variations of a Lipschitz domain and this is certainly a strong hint on how to set up the full drag minimization problem. Nevertheless, they leave open the following question.

Problem 4.2. Find a class of domains $D$ in which the minimization of the drag has to be sought.

While the drag force acts in the direction of the flow and, hence, in a uniquely determined one-dimensional direction, the lift force is orthogonal to the drag and has two degrees of freedom in a 3D setting. This is why it is more convenient to focus on the 2D cross section $K$ of the obstacle, whenever the obstacle itself is the right cylinder $D = (-\Lambda, \Lambda)\times K$ as in ($4.1$). As explained in Section 2.3, due to the vortex shedding the drag generates the lift and, then, one also needs to find shapes of the cross section $K$ that minimize the lift generating vertical oscillations, see Figure 11. Then, the total force ($4.4$) is obtained by computing the integral over the boundary $\partial K$ of the cross section, see again Figure 10 (right picture). Related to $(4.3)$, we assume that the 2D motion of the wind is governed by the (planar) Navier-Stokes equations

Even in this simplified form, the lift force (in the $z$-direction) is quite difficult to compute and one then proceeds numerically, for instance through Computational Fluid Dynamics (CFD). But the cross section of the deck of a suspension bridge can have very unpleasant shapes due to side-walks, lamps, guard-rails, etc... In such case, one obtains many "bad" vortices. In Figure 12 we quote some of the results obtained in ^[30]. This suggests the following

Problem 4.3. Find the optimal shape for the cross section that maximizes the wind velocity threshold generating asymmetric vortices.

As sketched in Figure 6, the threshold is determined in case of a circular cross section. But solving Problem 4.3 needs a rigorous characterization of the class of cross section where to seek the optimal one. A possible strategy is to reduce the admissible shapes to rectangular cross sections that are completed on the short edges by two caps of given measure (or mass) but with different shapes. In Table 2 we quote and complement some numerical results from ^[66] obtained for plates having thickness $d = 0.25$m and for $\eta = 1.5 \times 10^{-5}$m$^2$/s, the kinematic viscosity of air. The first column contains the shapes of the obstacles, all having two caps of total area $\pi d^2/8$ of equal mass. The second column contains the range of numerically found critical velocities $W_{*}$ yielding symmetry breaking of the flow.

But this should only be seen as the very beginning of the story. These results are just for few shapes and much more has to be done.

Problem 4.4. Through CFD, complement Table 2 with different shapes for the cross section and derive a rule in order to maximize $W_*$.

4.2.   Maximizing the stability by varying the shape of the deck

As we have seen in Section 2.5, vertical oscillations may suddenly transform into torsional oscillations, see the sketch of the cross section in Figure 13, and the law governing this phenomenon is related to the frequencies of the fundamental vertical and torsional modes.

Since the squared frequencies of the oscillation modes are equal to the eigenvalues of ($2.2$), see Theorem 2.1, one is so led to modify the shape of the plate in order to modify also the eigenvalues of ($2.2$). In fact, also engineers are nowadays strongly interested in the shape optimization for the design of bridges and, in particular, on the sensitivity analysis of certain eigenvalue problems ^[72]. In this respect, Banerjee ^[11] goes even further: he claims that the free vibration analysis is a fundamental pre-requisite before carrying out a flutter analysis.

Whence, in order to improve the torsional stability of the plate one needs to analyze the behavior of the eigenvalues of ($2.2$) with respect to perturbations of the rectangular plate $\Omega = (0, \pi)\times(-\ell, \ell)$. The simplest possible perturbation consists in varying the ratio between width and length, namely $\ell/\pi$. For this problem, we have

Theorem 4.5. (^[16]) Let $\mu_{m, k}(\ell)$ and $\nu_{n, j}(\ell)$ be the eigenvalues of $(2.2)$. Then there exist explicit formulas to compute their derivatives with respect to $\ell$.

Theorem 4.5 solves the technical problem of computing the derivatives of the eigenvalues with respect to the variations of $\ell$ but fails to give ideas on how to improve the torsional stability. Some hints come from the physical/engineering literature. Rocard ^[90] considers the natural vertical ($\omega_v$) and torsional ($\omega_t$) frequencies of the bridge. He claims that for common bridges one has $\omega_t > \omega_v$ and that the critical velocity $W_c$ of the wind for which the bridge undergoes to flutter may be computed through the formula

where $C > 0$ is a physical constant depending on the shape and the material composing the deck. This formula was later modified ^[81,95] by putting a different multiplicative constant $C$. Moreover, Holmes [70,p.293] shows that the Rocard-Selberg formula does not always agree with experimental measurements, which makes ($4.6$) unreliable from a theoretical point of view. Indeed, the values of the "characteristic parameters" of a bridge are not found theoretically but with wind tunnel tests, that allow to compute the so-called flutter derivatives, also known as aeroelastic (or aerodynamic) derivatives. These are the coefficients to be inserted in suitable linear ODE's used to find the vertical and torsional frequencies, as first suggested by Scanlan-Tomko ^[93]. The flutter speed may then be computed through closed formulas, which, however, are very different from ($4.6$); see ^[76]. Since several vertical and torsional frequencies ($\omega_v$ and $\omega_t$) are involved in the oscillations of a bridge, see Theorem 2.1, this suggests the following problem whose solution needs the explicit formulas of the derivatives from Theorem 4.5.

Problem 4.6. For any vertical mode $v$, let $T_v > 0$ be its torsional instability threshold, see Proposition 2.4. Find the optimal width of a bridge maximizing $T_v$ for the relevant vertical modes $v$; "relevant" refers here to the modes which are more prone to switch to torsional oscillations, see Proposition 2.5.

So far, we have only considered rectangular plates $\Omega$ but some recent glass suspension bridges have been built with different shapes, in particular in China, see the two pictures in Figure 14. Therefore, we are also interested in modifying the shape of the free edges of the plate without altering the width of the hinged ones. Physically meaningful shapes should be symmetric with respect to the midline of the deck, as in the right picture in Figure 14, namely

for some $\phi\in C^0[0, \pi]$ such that $\phi(0) = \phi(\pi) = 0$. Note that $\Omega_{0} = \Omega$.

The first step consists in solving the following problem.

Problem 4.7. Find a domain functional able to quantify how much a plate is prone to transform vertical oscillations into torsional ones. What properties are to be expected from this functional?

Such a functional is expected to depend on the particular couple of (vertical, torsional) modes considered. First attempts to solve Problem 4.7 were made in ^[15] by requiring the functional to satisfy several properties commonly accepted in literature ^{[11,70,72,81,90,93,95]} and in practice ^[40]. Although [16,Theorem 3.3] provides explicit formulas for the computations of the derivatives of the eigenvalues with respect to shape variations, the final conclusion in ^[16] was that there exists no domain functional satisfying all the expected properties. Therefore, Problem 4.7 appears quite challenging.

4.3.   Maximizing the stability through an optimal stiffening of the deck

For the plate $\Omega$ in ($2.1$), consider the functional space

We also denote by $H^{-2}_*(\Omega)$ the dual space of $H^2_*(\Omega)$ and by $\langle\cdot, \cdot\rangle$ the corresponding duality. Let $f\in L^2(\Omega)$ and consider the problem

which contains the same boundary conditions as ($2.2$). For $f\in H^{-2}_*(\Omega)$, this problem may be written in weak form as

This formulation is obtained by minimizing the energy of the deformation $u\in H^2_*(\Omega)$ subject to a force $f\in H^{-2}_*(\Omega)$, namely

It was shown in ^[49] that ($4.8$) admits a unique solution $u\in H^2_*(\Omega)$ and, since the boundary conditions in ($4.7$) satisfy the complementing conditions, elliptic regularity applies. This means that for $f\in L^2(\Omega)$ the solution of ($4.7$) satisfies $u\in H^4\cap H^2_*(\Omega)$.

Since $\Omega$ is a planar domain, we have the embedding $H^2_*(\Omega)\subset C^0(\overline{\Omega})$ so that weak solutions of ($4.7$) are continuous. Therefore, for any weak solution $u$ we can compute the gap function

that was introduced in ^[15] in order to measure the torsional performances of bridges. Clearly, to compare the performances for different forces $f$, one should normalize them in a suitable space.

A challenging shape optimization problem is then

Problem 4.8. Is it possible to reduce the torsional oscillations in bridges by inserting a strong material in an open subset $D\subset\Omega$ of given measure $|D| = \alpha\in(0, 2\pi\ell)$?

Problem 4.8 needs to be set up with extreme precision. Let $\chi_D$ be the characteristic function of $D$ and $d > 0$ be the stiffness of the strong material. One has the choice of stiffening the plate, in which case the energy becomes

or weakening the action of the wind with some device, which leads to the energy

There are several essential differences between these two energies. For $f\in H^{-2}_*(\Omega)$, the minimizer $u_{f, D}\in H^2_*(\Omega)$ of $\mathbb{E}_1$ satisfies

but, even if $f$ is smooth, this weak formulation admits no strong form. On the other hand, it appears difficult to define $\mathbb{E}_2$ if $f\not\in L^1(\Omega)$, while if $f\in L^1(\Omega)$ the minimization of $\mathbb{E}_2$ leads to the (strong) equation

complemented with the boundary conditions in ($4.7$). For both the energies $\mathbb{E}_1$ and $\mathbb{E}_2$, after finding the minimizer $u_{f, D}\in H^2_*(\Omega)\subset C^0(\overline{\Omega})$ for given $D$ and $f$, we compute its gap function:

It appears of great engineering interest to solve the following problem.

Problem 4.9. Find two classes ${\mathcal D}$ (of domains of given measure $\alpha\in(0, 2\pi\ell)$) and ${\mathcal F}$ (of normalized functions) such that

Then study the minimization problem $\min\{ \mathcal{G}^\infty_D:\, D\in {\mathcal D}, \, |D| = \alpha\}$.

In other words, we need to solve the minimaxmax problem

that was introduced in ^[17] and that has common points with some torsion problems considered in ^[34,35,75]. From a practical (engineering) point of view, homogenization should be avoided: this means that

With this additional constraint, some partial answers to Problem 4.9 are known.

Theorem 4.10. (^[17]) For a given open set $D\subset\Omega$, the problem

has a solution (case $\mathbb{E}_1$). For a given open set $D\subset\Omega$ and $p\in (1, +\infty]$, the problem

has a solution (cases $\mathbb{E}_1$ and $\mathbb{E}_2$).

Note that, in both cases, there is no uniqueness since $\pm f$ generate opposite minimizers and opposite gap functions. Moreover, the class of open sets $D$ appears too large for applications and some geometric constraints should be added. Some particular classes with geometric constraints yielding a solution of ($4.9$) were found in ^[17] and the role of symmetry was emphasized. Further classes needing a strong numerical support were introduced in ^[5,20] but a fully satisfactory answer to Problem 4.9 is still missing. We feel that polygonal trusses might play an essential role, both because hexagons have several geometric and energetic performances ^[58] and because common suspension bridges mostly use triangles, see Figure 15.

But adding a stiffening truss is not the only way to improve the stability of the deck. For the next shape optimization problem, we need to analyze how a suspension bridge is erected. The deck segments are put in position one aside to the other and have the shape of (large) rectangles, see Figure 16.

The purpose is to decide which of these segments (rectangles) should be made of a more stiff material. To this end, the rectangular deck $\Omega$ is divided into some smaller rectangles (the segments) as in Figure 17.

If $\Omega = (0, \pi)\times(-\frac{\pi}{150}, \frac{\pi}{150})$, then the rectangles are given by

Problem 4.11. Find the optimal position of $m$ stiff rectangles ($1 < m < n$) in order to minimize the torsion of the deck.

For this problem, we have

and the optimal shape will have a form such as in Figure 18.

For ${\mathcal F} = \{f\in L^2(\Omega); \, \|f\|_2 = 1\}$, the gap function for these composite rectangles was computed numerically in ^[5] for different $m$ and $n$. The results suggested the conjecture that with no lower bound on the minimum size of the deck segments, the optimal shape does not exist and a minimizing sequence leads to homogenization, precisely what has to be avoided for practical reasons. We also refer to ^{[18,19,21,43]} for more shape optimization problems and to ^[47,48] for further recent developments.

Acknowledgments

This research was funded by the National Plan for Science, Technology and Innovation (MAARIFAH), King Abdulaziz City for Science and Technology, Kingdom of Saudi Arabia, Award Number (13-MAT887-02).

Conflict of interest

The authors declare no conflict of interest.