A DYNAMIC PROBLEM INVOLVING A COUPLED SUSPENSION BRIDGE SYSTEM: NUMERICAL ANALYSIS AND COMPUTATIONAL EXPERIMENTS

. In this paper we study, from the numerical point of view, a dynamic problem which models a suspension bridge system. This problem is written as a nonlinear system of hyperbolic partial diﬀerential equations in terms of the displacements of the bridge and of the cable. By using the respec- tive velocities, its variational formulation leads to a coupled system of parabolic nonlinear variational equations. An existence and uniqueness result, and an exponential energy decay property, are recalled. Then, fully discrete approximations are introduced by using the classical ﬁnite element method and the implicit Euler scheme. A discrete stability property is shown and a priori error estimates are proved, from which the linear convergence of the algorithm is de- duced under suitable additional regularity conditions. Finally, some numerical results are shown to demonstrate the accuracy of the approximation and the behaviour of the solution.


(Communicated by Josef Malek )
Abstract. In this paper we study, from the numerical point of view, a dynamic problem which models a suspension bridge system. This problem is written as a nonlinear system of hyperbolic partial differential equations in terms of the displacements of the bridge and of the cable. By using the respective velocities, its variational formulation leads to a coupled system of parabolic nonlinear variational equations. An existence and uniqueness result, and an exponential energy decay property, are recalled. Then, fully discrete approximations are introduced by using the classical finite element method and the implicit Euler scheme. A discrete stability property is shown and a priori error estimates are proved, from which the linear convergence of the algorithm is deduced under suitable additional regularity conditions. Finally, some numerical results are shown to demonstrate the accuracy of the approximation and the behaviour of the solution.
1. Introduction. During the last decades the study of the so-called suspension bridges has received a large attention because this kind of bridges is a common type of civil engineering structure. It is well-known that these bridges may display certain oscillations under external aerodynamic forces like, for instance, it occurred in the famous Tacoma's bridge (see [2,5]), in which a strong wind caused the collapse of a narrow and very flexible suspension bridge.
Since the pioneering works by Lazer and McKenna (see, e.g., [14] or [18]), where mathematical models based on nonlinear partial differential equations were introduced to describe oscillations in suspension bridges including connection cables, a large number of papers have been published (see, for instance, [1,3,10,11,12,13,16,19,20,21,22,23,24,25] and the numerous references cited therein). Most of the above works deal with mathematical aspects as the application of analytical methods to scrutinize stationary solutions, periodic oscillations, longtime global dynamics, existence and uniqueness of solutions and so on. Moreover, we want to point out that in the literature other types of suspension bridges were studied, where the associated dynamic system was written in the vertical and torsional displacements of a cross-section of the bridge (see, e.g., [4,15]). However, to our knowledge the numerical analyses of the proposed models have not been performed yet.
In this paper, we revisit the problem considered in [7], where a model taking into account the coupling between the road bed and the suspension main cable was studied from the mathematical point of view, proving the existence and uniqueness of weak solutions by using the Faedo-Galerkin approximation procedure and Gronwall's lemma, and the exponential decay of the system energy. Moreover, the existence of a regular global attractor was also shown by using the semigroup theory and defining an adequate Lyapunov functional. Here, we aim to provide the numerical analysis of this dynamic problem, introducing fully discrete approximations, proving a discrete stability result and a priori error estimates, and performing numerical simulations which demonstrate the accuracy of the algorithm and the behaviour of the solution.
The paper is outlined as follows. The mathematical model is described in Section 2 following [7], deriving its variational formulation. An existence and uniqueness result, and an energy decay property, proved in [7] are also stated. Then, in Section 3 a numerical scheme is introduced, based on the finite element method to approximate the spatial domain and the forward Euler scheme to discretize the time derivatives. A discrete stability property is proved and a priori error estimates are deduced for the approximative solutions from which, under suitable regularity assumptions, the linear convergence of the algorithm is obtained. Finally, some numerical simulations are presented in Section 4.

2.
The model and its variational formulation. In this section, we present briefly the model, the required assumptions and the variational formulation of the mechanical problem, and we state an existence and uniqueness result. We refer the reader to [7] for details.
Let [0, ], > 0, be the one-dimensional beam (bridge) or rod (cable) of length and denote by [0, T ], T > 0, the time interval of interest. Moreover, let x ∈ [0, ] and t ∈ [0, T ] be the spatial and time variables, respectively. In order to simplify the writing, we do not indicate the dependence of the functions on x and t, and a subscript under a variable represents its derivative with respect to the prescribed variable.
Let u denote the downward deflection of the deck midline (of unitary length) in the vertical plane with respect to his reference configuration (that is, it represents the bending displacement of the road bed), and let w be the vertical displacement of the vibrating spring standing for the main cable. Thus, the mechanical problem which simulates the coupled system of a vibrating suspension bridge and a cable is written in the following form (see [7]).
Here, the notation (F ) + stands for the positive part of a function F , i.e. (F ) + = max{0, F }. Moreover, the term −k 2 * (u − w) + models a restoring force due to the cables, f is the given vertical dead load distribution on the deck, the constant p represents the axial force acting at the ends of the road bed in the reference configuration (being negative when the bridge is stretched and positive when compressed), and g denotes an external source applied in the cable. Finally, following [7] we assumed that all the constants in the model were equal to 1 for the sake of simplicity in the writing. However, in order to simplify the calculations in this work we note that we have modified slightly the boundary conditions employed in [7], where the ends of the vibrating beam were considered pinned, assuming instead that they are rigidly fixed now. The analysis performed there could be adapted with some minor changes.
In order to obtain the variational formulation of Problem P, let Y = L 2 (0, ) and denote by (·, ·) the scalar product in this space, with corresponding norm · . Moreover, let us define the variational spaces V and E as follows, with scalar product (·, ·) V (resp. (·, ·) E ) and norm · V (resp. · E ) defined in V (resp. E).
By using the integration by parts and the Dirichlet boundary conditions at x = 0, , we write the variational formulation of Problem P in terms of the bending velocity v = u t and the vertical velocity e = w t .
where the bending displacement and the vertical displacement are then recovered from the relations: The following theorem has been proved in [7], providing the existence of a unique solution to Problem VP as well as an energy decay property.
and assume that the initial data have the regularity Then, there exists a unique solution to Problem VP with the following regularity: Moreover, this solution is continuously dependent on the given data.
In addition, if we assume that f = g = 0 (i.e. there are not volume forces) and that the axial force p is small enough, denoting by E(t) the energy of the system given by 2 , then it decays exponentially; i.e. there exist two positive constants Q and ω such that 3. Numerical analysis: Fully discrete approximations and a priori error estimates. In this section we consider a fully discrete approximation of Problem VP. This is done in two steps. First, we assume that the interval [0, ] is divided into M subintervals a 0 = 0 < a 1 < . . . < a M = of length h = a i+1 − a i = /M and so, to approximate the variational spaces V and E, we construct the finite dimensional spaces V h ⊂ V and E h ⊂ E given by where P r ([a i , a i+1 ]) represents the space of polynomials of degree less or equal to r in the subinterval [a i , a i+1 ]; i.e. V h is composed of C 1 and piecewise cubic functions and E h is composed of continuous and piecewise affine functions. Here, h > 0 denotes the spatial discretization parameter. Moreover, we assume that the discrete initial conditions, denoted by u h 0 , v h 0 , w h 0 and e h 0 , are given by P h 1 and P h 2 being the classical finite element interpolation operators over V h and E h , respectively (see [9]).
Secondly, we consider a partition of the time interval [0, T ], denoted by 0 = t 0 < t 1 < · · · < t N = T . In this case, we use a uniform partition with step size k = T /N and nodes t n = n k for n = 0, 1, . . . , N . For a continuous function z(t), we use the notation z n = z(t n ) and, for a sequence {w n } N n=0 , we denote by δw n = (w n − w n−1 )/k its divided differences.
Therefore, using the implicit Euler scheme, assuming that f, g ∈ C([0, T ]; Y ) the fully discrete approximation of Problem VP is as follows.
where the discrete bending displacement and the discrete vertical displacement are then recovered from the relations: The existence of a unique solution to Problem VP hk can be obtained proceeding as in the continuous case, applying a standard Faedo-Galerkin approximation procedure (see, for instance, [6,17]).
The aim of this section is to obtain some a priori error estimates on the numerical errors u n − u hk n , v n − v hk n , w n − w hk n and e n − e hk n . Now, we have the following discrete stability property. Proof. In order to simplify the writing of this proof, we remove the superscripts h and k in all the variables and we assume that f = g = 0.
Taking as a test function ξ h = v n in equation (14) we have where we used Cauchy-Schwarz and Young inequalities, it follows that Now, taking ψ h = e n as a test function in equation (15) we find that (δe n , e n ) + ((w n ) x , (e n ) x ) + (e n , e n ) − k 2 * ((u n − w n ) + , e n ) = 0, and so, taking into account that (δe n , e n ) ≥ 1 2k e n 2 − e n−1 2 , Combining now (17) and (18) we can conclude that Thus, by induction and using Poincare's inequality we find that and, applying a discrete version of Gronwall's inequality (see, for instance, [8]) we find the desired stability estimates.
Remark 1. We note that we could modify Problem VP hk by using a semi-implicit scheme in the following form: Find the discrete bending velocity v hk = {v hk n } N n=0 ⊂ V h and the discrete vertical velocity e hk = {e hk n } N n=0 ⊂ E h such that v hk 0 = v h 0 , e hk 0 = e h 0 and, for n = 1, . . . , N , where the discrete bending displacement and the discrete vertical displacement are again recovered from relations (16).
This new problem is now uncoupled and linear and so, the numerical resolution is easier. However, Lemma 3.1 should be modified and the resulting estimates will vary accordingly. Now, we are able to prove the following a priori error estimates result.

Theorem 3.2. Let the assumptions of Theorem 2.1 still hold and denote by (v, e)
and (v hk , e hk ) the solutions to problems VP and V P hk , respectively. Assume that the following additional regularity is satisfied: Then we have the following a priori error estimates, for all Proof. Taking as a test function ξ = ξ h ∈ V h ⊂ V in equation (8) at time t = t n and subtracting it to equation (14) we have * ((u n − w n ) + − (u hk n − w hk n ) + , ξ h ) = 0, and therefore, where we used Lemma 3.1 and the notation δu n = (u n − u n−1 )/k, and the following estimates have been used, we find that Now, taking as a test function ψ = ψ h ∈ E h ⊂ E in equation (9) at time t = t n and subtracting it from equation (15) we have ((e t ) n − δe hk n , ψ h ) + ((w n − w hk n ) x , ψ h x ) − k 2 * ((u n − w n ) + − (u hk n − w hk n ) + , ψ h ) +(e n − e hk n , ψ h ) = 0 ∀ψ h ∈ E h , and so, we obtain, for all ψ h ∈ E h , ((e t ) n − δe hk n , e n − e hk n ) + ((w n − w hk n ) x , (e n − e hk n ) x ) + (e n − e hk n , e n − e hk n ) −k 2 * ((u n − w n ) + − (u hk n − w hk n ) + , e n − e hk n ) = ((e t ) n − δe hk n , e n − ψ h ) + ((w n − w hk n ) x , (e n − ψ h ) x ) + (e n − e hk n , e n − ψ h ) −k 2 * ((u n − w n ) + − (u hk n − w hk n ) + , e n − ψ h ). Keeping in mind that (δe n − δe hk n , e n − e hk n ) ≥ where δw n = (w n − w n−1 )/k, using several times the above Cauchy's inequality it follows that 1 2k e n − e hk n 2 − e n−1 − e hk n−1 ≤ C (e t ) n − δe n 2 + e n − e hk n 2 + w n − w hk n 2 + u n − u hk Combining estimates (20) and (21) we have, for all Thus, by induction we obtain, for all Finally, taking into account that = (e n − e hk n , e n − ψ h n ) + (e h 0 − e 0 , e 1 − ψ h 1 ) applying a discrete version of Gronwall's inequality (see again [8]) and Poincare's inequality, we conclude estimates (19).
We note that estimates (19) can be used to perform an error analysis. Thus, as an example, assume the following additional regularity on the continuous solution: We note that with this regularity we immediately find that (see [9]): Therefore, we have the following. Corollary 1. Let the assumptions of Theorem 3.2 and the additional regularity (22) hold. Then, the linear convergence of the algorithm is deduced; i.e. there exists a positive constant C, independent of the discretization parameters h and k, such that The proof of the above result is done in several steps. First, we have the following property of approximation by finite elements (see, e.g., [9]): Moreover, using again regularity conditions (22), we have . Finally, the remaining terms in estimates (19) can be bounded as follows (see [8] for details), . Combining all these estimates, the linear convergence is deduced.
4. Numerical results. In order to verify the behaviour of the numerical method described in the previous section, some numerical experiments have been performed. Moreover, for the sake of generality in the simulations, we have included coefficients in the definition of the mechanical problem P. 4.1. Numerical scheme. Given the solution u hk n−1 , v hk n−1 , w hk n−1 and e hk n−1 at time t n−1 , the discrete bending velocity is obtained from the discrete nonlinear variational equation . Later, we get the vertical displacement from the discrete nonlinear variational equation ρ 2 (e hk n , ψ h ) + k 2 a 2 ((e hk n ) x , ψ h x ) + c 2 k(e hk n , ψ h ) − k k 2 * ((u hk n − w hk n ) + , ψ h ) = ρ 2 (e hk n−1 , ψ h ) + k(g n , ψ h ) − k a 2 ((w hk n−1 ) x , ψ h x ), where the discrete bending displacement and the discrete vertical displacement are then recovered from the relations We note that both numerical problems consist of nonlinear symmetric systems, and so a fixed point method was applied for their solution.

4.2.
A first example: Numerical convergence. Our aim with this first example is to verify the numerical convergence of the numerical scheme. In this sense, the following problem is considered.
Using the discretization parameters k = h = 10 −3 , in Fig. 2 the evolution in time of the vertical displacement of the bridge center is shown for several values of the axial force p. As can be observed, the increment in the compression forces makes the oscillating response to be softer when the bridge is released.

4.4.
A third example: The coupling term. As a final test we analyze the dependence of the model with respect to the parameter k 2 * . Taking both the bridge and the cable in their reference configuration as initial conditions, a positive linearly increasing force acts over the bridge for one second. Taking the discretization parameters k = h = 10 −3 , running several simulations for different values of k 2 * the effect of the coupling term is easily noticed: when this value raises, the cable restricts the deformation of the bridge and so the displacement decreases, while at the same time the cable deformation increases, as can be seen in Fig. 3.