Aeroelastic Response from Indicial Functions with a Finite Element Model of a Suspension Bridge

The present paper describes a comprehensive analysis of the aeroelastic bridge response in time-domain, with a finite element model of the structure. The main focus is on the analysis of flutter instability, accounting for the wind forces generated by the bridge motion, including twisting as well as vertical and horizontal translation, i.e. all three global degrees of freedom. The solution is obtained by direct integration of the equations of motion for the bridge-wind system, with motion-dependent forces approximated from flutter derivatives in terms of rational functions. For the streamlined bridge box-girder investigated, the motion dependent wind forces related to the along-wind response are found to have a limited influence on the flutter velocity. The flutter mode shapes in the time-domain and the frequency domain are consistent, and composed of the three lowest symmetrical vertical modes coupled with the first torsional symmetric mode. The method applied in this study provides detailed response estimates and contributes to an increased understanding of the complex aeroelastic behaviour of long-span bridges.


Introduction
The verification of wind-resistant design of a long-span bridge relies on experimentally obtained aerodynamic data from wind tunnel studies, applied in computational models of different levels of complexity. Established wind load models are based on the principle of superposition of wind loads separated into mean wind loads, buffeting loads and motion-induced loads [1]. While the 10 minutes mean wind velocity with an annual probability of exceedance of 0.02 is typically considered when studying the buffeting load effects [2], the aeroelastic stability has to be ensured at more extreme wind speeds.
In the present work, the motion-induced wind loads on the bridge deck are based on the formulation by Scanlan [3]. The load model defines the self-excited forces as linear functions of the displacements and velocities based on the assumption of small oscillations. In the time-domain, the aeroleastic loads can be represented by indicial functions in convolutions with the time-history response of the bridge [4].
Strong winds generate large displacements of the bridge girder and the cables. The time-domain finite element analysis (FEA) performed, include modifications of the structural stiffness due to large displacements, i.e. geometric non-linearity, through iterative solutions of the system equations. While this involves higher computational effort in comparison to modal based analyses in the frequencydomain, time-domain analyses enables a direct inspection of wind-induced stresses in different parts of

The Hardanger Bridge girder section
The computational method developed is applied in a case study considering the Hardanger Bridge, in service since 2013. The suspension bridge has nearly 200 m tall towers and a main span of 1310 m, connected to the 45 m and 25 m long viaducts. The girder is a wedge-shaped 18.3 m wide steel box with a width to depth ratio of 5.5. The mass per unit length is 12800 kg/m and the mass moment of inertia is 428400 kg·m 2 /m. The bridge girder carries two lanes of the road traffic, one in each direction, and a separate lane for bicycles and pedestrians. The bridge girder cross-section is illustrated in figure 1.
A series of wind tunnel tests [17] was dedicated to investigate the box-girder's aerodynamic behaviour, including the vortex-shedding response, the buffeting response and the flutter instability. The aerodynamic characteristics of an updated bridge girder design were investigated by Svend Ole Hansen ApS [16], in collaboration with the Norwegian Public Road Administration (NPRA) [18]. The wind tunnel data was utilized in various studies to estimate the full-scale aerodynamic bridge response [13,14,15,19]. In parallel with extensive analyses based on modal coordinates, a detailed verification of the bridge design was launched by NPRA, using a finite element model for the response analysis in timedomain. The present work was initiated to implement a comprehensive representation of the motiondependent loads, through the use of indicial functions, for the time-domain FEA. The original windresistant design of the Hardanger Bridge was based on flutter derivatives estimated from ambient vibration data from 2 DOF section model tests [18] in which the lowest vertical and torsional modes were represented. The objective of the present work was to implement all self-excited load components based on experimental data. An additional and complete set of flutter derivatives was provided as part of a cooperation project between Ruhr-University of Bochum and the University of Stavanger (UiS). The forced vibrations tests, performed with a section model in nominally smooth flow, and the identification of a complete set of 18 flutter derivatives were carried out by the UiS partners [20].

Finite element model
The finite element model, illustrated in figure 2, is defined in the general purpose finite element software Abaqus [21]. The model was made available for the present investigation by NPRA. The cables, hangers and the bridge girder are modeled with beam elements, of type B31 in the Abaqus notation. The elements are defined by two nodes, with 6 DOF per node, providing linear interpolation of displacements within the element. The bridge girder and the main cables are discretized in 64 sections corresponding to 20 m long elements, following the spacing between the hangers. The cables of the side spans are discretized with 20 elements. In total, 64 pairs of hangers are suspended from the main cables. Each hanger is modeled by one beam element. The cables and the hangers are assigned a negligible bending stiffness, thus effectively acting as truss elements.
The model devoted to the aeroelastic bridge analyses includes a simplified representation of the bridge towers, by linear springs and point masses. The structural damping of the bridge girder is differentiated for the lateral, the vertical and the twisting response. A damping ratio, relative to critical, was assumed 0.35% for both the fundamental horizontal and the lowest vertical mode and 0.32% for the first torsional mode. The structural damping is implemented with dashpot elements added to the girder and the cable elements in combination with stiffness proportional damping distributed to the beam elements along the girder. Input data i.e. dashpot coefficients and Rayleigh parameters corresponding to the assumed damping ratios is determined by iterative procedures analyzing the forced vibration response of the FE-model.

Time-domain solutions
The equation of motion of the bridge-wind system driven by the mean wind loads and the aeroelastic loads, Fq and Fae respectively, is given in equation (1): Here, M, C, and K are the global mass, damping and stiffness matrices and R(t) is the nodal displacement vector of the structure. The non-linear solver option in Abaqus [21] relevant for direct implicit integration of the equation of motion is an implementation of the Hilber-Hughes-Taylor method [22]. The algorithm has options to use a numerical damping and an adaptive time stepping scheme. In the present analyses, a fixed time step of 0.008 s and zero numerical damping were applied, and the numerical integration carried out using the trapezoidal rule. Wind load on the bridge girder and the cable elements was implemented by developing special purpose elements in Fortran. The uel subroutine is called at each time-step during iterations and the update of the self excited loads. The aeroelastic loads, Fae(t), are defined by convolutions implementing indicial functions and the motion of the girder. Selfexcited drag loads acting on the cables and the hangers are defined by the quasi-steady theory.

Motion-induced forces
The self-excited lift force Lae, per unit length of the bridge girder, is given in equation (2) following the notation by Scanlan [3]: where dynamic pressure, mean wind velocity and the deck width are denoted q0, U and B respectively. Flutter derivatives Hi * are functions of the reduced frequency k=ωB/U.
The motion-induced twisting moment Mae is expressed analogously to equation (2) by replacing Hi * with Ai * and the drag Dae by replacing Hi * with Pi * .

Formulation in terms of rational functions
By defining a non-dimensional displacement amplitude vector r0={z0/B 0 x0/B} T , assuming harmonic oscillations at an arbitrary frequency , the motion-induced loads can be expressed with (3) [23]: Equation (3) defines nine complex transfer functions Qmn for the motion-induced loads. Roger [24] proposed a rational function approximation of the aeoroelastic loads, given in equation (4) as presented by Abel [25]: The aeroelastic loads determined by convolutions with indicial functions and the structural response are given in equation (6) [26].
When transferred to physical time t, the loads take the form of equation (7) where n specifies nodal location along the bridge girder.
The transfer functions in equation (4) can be separated into real and imaginary parts as shown in equation (8) and (9): The transfer functions Qmn are fitted to the experimentally obtained flutter derivatives along the imaginary axis k in the complex plane. The real parts of the functions characterize aeroedynamic stiffness, whereas the imaginary parts characterize aeroedynamic damping. The parameter matrix A2 is interpreted as aerodynamic mass and the corresponding load term is assumed negligible in the context of the bridge response analysis. The coefficient matrices Ai and the parameters l are approximated with the extended minimum state method [23]. The A0 matrix incorporate the steady values of the flutter derivatives. To ensure decaying impulse functions, i. e. stability of the aeroelastic system, l is restricted to have positive values. Mathematically, the l parameters define the pole locations of Qmn(p) along the negative real axis. The transfer functions are constrained through a set of common poles l. The summation terms of the functions, called lag terms, considerably extend the computational effort with respect to the response analyses in time-domain. Furthermore, in a state-space format, the number of lag terms contribute to the dimension of the aeroelastic system. In this context Karpel [27] proposed a minimum-state method to reduce the dimension of the state-space equations, which introduced a coupling between the coefficients in the lag terms. Tiffany and Adams [28] suggested to use the extended minimum state method, splitting up the minimization problem into a linear process to determine the coefficient matrices and a non-linear process for finding the location of the poles [23].

Previous studies using rational function approximations for aeroelastic response analyses of bridges
Wilde and Fujino [29] and Chen et al. [8] utilized Roger's approximation [24] in the application of windloads in the bridge response analyses. Wilde and Fujino based their work on the flutter derivatives corresponding to the theoretical load formulation by Theodorsen [30]. Chen et al. utilized both the experimentally obtained flutter derivatives for a twin-box girder and the Theodorsens functions and presented the selected approximations for H3 * , A1 * , A2 * and A3 * .
Peil and Kirch [31] investigated the suitability of the rational function approximations to properly describe the motion-induced wind loads on bridge girder sections as well as the flat plate. Based on the use of experimental quantified flutter derivatives, they demonstrated that the approximation error could be reduced by increasing the number of lag terms, but at the same time the approximation itself could be considerably distorted. Distortion, in this context, refers to a notable deviation of the approximations from the experimental data points in wide areas above the pole locations in the complex plane. This behaviour was connected to a strong weighting of the partial fractions in the evaluation process. The effects were found to be pronounced for the Tacoma Bridge girder section but not so significant for the East Great Belt Bridge girder section. In this respect, it should be noted that the The Hardanger Bridge girder section is a streamlined box-girder, similar to the case of the East Great Belt Bridge girder section.
Previous work on multimodal flutter analysis of the Hardanger Bridge [32] addressed the influence of the variable number of pole terms in the rational function approximations on the resulting flutter speed of the bridge. Deviations in flutter speed were found to be within 6% when using the above discussed function approximations, using four, six and nine pole terms.

Measurements and approximations of self-excited loads
Flutter derivatives representing lift, twisting moment and drag, a full set of 18 derivatives, were determined from tests in nominally smooth flow over a wide range of reduced frequencies (0. 3 -9). Details of the experimental set up for the forced vibration tests and the computational procedures for the identification of flutter derivatives are presented in Neuhaus et al. [20]. In the present investigation, the motion-dependent forces at zero mean wind angle of incidence are applied. Forcing amplitudes were up to 6% of the girder height in the vertical translatory motion and 1° in twisting motion. The tests were performed by forcing each section model degree of freedom separately and quantifying the associated load components in terms of the derivatives Hi * , Ai * and Pi * .
Indicial function approximations with the extended minimum-state method [23], have been used to define the aeroelastic loads in this study. The response analyses with the FE-model in time-domain are performed using the indicial function approximations with six lag terms. With respect to the response analyses, it is to be noted that the reduced torsional frequency of the bridge girder is about 1 at the design wind speed and nearly 0.4 at the estimated flutter speed. It is thus essential that the approximations are most accurate in the corresponding range of reduced frequencies.
The experimental data and the fitted approximations of the aeroelastic transfer functions are presented in figures 3 and 4 in terms of their real and imaginary parts.

Frequency-domain analysis
The results of the time-domain aeroelastic analysis will be compared to those obtained by performing a multi-modal analysis in the frequency.domain. By denoting the structural modal quantities Ms, Cs and Ks as well as the modal displacement vector Z and applying the load model in equation (7), the equation of motion is written: (10) Za in equation (10) defines the aerodynamic states. With the state vector set equal to Zst={Z Ż̇ Za} T , the bridge-wind system matrix is established with equation (11): where Cq and Kq incorporate both the structural and the aerodynamic modal properties of the system and coefficients Dij in D are equal to i(x) j(x)dx / i 2 (x)dx where i(x) j(x) are the combined modes.

Natural modes
Natural frequencies f and periods T for the still-air eigen-modes, estimated with the finite element model, are presented in Table 1. Vs-i, Hs-i and Ts-i, with i=1-3, denote the vertical, horizontal and the torsional direction. The frequency ratio between the first torsional and the first vertical mode is 2.6. The first three symmetrical modes in each direction are presented in figure 6.   A horizontal component of the first torsional mode, also identified in [18], is included in figure 7. The torsional response corresponding to the vertical motion of cable planes of ±1 m is accompanied by the horizontal displacement of 0.16 m at girder mid-span.    The motion-induced drag on the cables and hangers are implemented by quasi-steady assumptions in the form of equation (12): where CC is static shape coefficient, DC is representative cable diameter, LC the length of cable element and R C yd(t) is horizontal displacements relative to time-averages in the mean wind direction. The transient response, at design wind speed 38.6 m/s, at the girder mid-span, is shown in figure 8. The directions of structural displacements in the fixed global coordinate system are adopted such that the positive ryd is horizontal displacement in the mean wind direction, positive rd is the rotation of the bridge girder with the windward side moving upwards, and the positive rzd is the upwards vertical displacement.

Dynamic response with reduced set of motion-induced loads
In order to separate the effect of the motion-dependent wind load related to the horizontal bridge response, the characteristics of the aeroelastic response based on two different reduced sets of motioninduced loads is presented in the following. The first case concerns the time-domain analyses with 4 indicial functions Lz, L  , Mz and M  , and the second one involves an additional indicial function for the aeroelastic drag due to the horizontal motion, Dy. The two reduced sets of motion-dependent loads are denoted "4 i.f." and "5 i.f." The motion-induced loads on cables are included. Figure 12 presents 300 s long time-histories of responses at mid-span girder, for the two above specified load cases, both at the mean wind velocity of 76.3 m/s. The interaction between the vertical and the twisting motion is recognized in both load cases. In the case, based on the five indicial functions, the flutter oscillations develop within approximately 180 s. In the load case with four indicial functions, the lowest horizontal mode contributes significantly to the vertical response, even after 300 s. The impact of the motion-induced drag damping is demonstrated in figures 12(a) and 12(c). The horizontal motion in figure 12(c) is associated with the damping ratios of 2.5% and 10% for the load cases with 4 i.f. and 5 i.f. respectively.

Multimode flutter analysis
The stability of the aeroelastic system is also investigated in the frequency domain, by performing the eigenvalue analysis of the system matrix in equation (11). Eigenvalue solutions are studied by increasing the mean wind speed in steps of 0.1 m/s. Nine lowest symmetric modes are included in the analysis: Vsi (1-4), Ts-i (1-2) and Hs-i (1-3). The damping ratios of the in-wind bridge modes are presented in figure  13. Critical flutter speed is determined to be 71.9 m/s, the velocity at which damping for one of the modes becomes zero.
The flutter mode is illustrated in the complex plane in figure 14. For clarity, only the most significant contributions are given in the plot. Figure 14 shows that the three lowest vertical modes contribute in the flutter mode. The flutter frequency from the eigenvalue solutions is 0.245 Hz. Figure 15 shows the time-domain flutter response identified in figure 11 together with the eigenvector from the multimodal response analysis. The two flutter modes have a strong similarity.

Discussion
The wind-tunnel tests undertaken for verifications of the aerodynamic design of the Hardanger Bridge girder [16] indicated that the full-scale flutter velocity is 79.5 m/s in smooth flow and 70.3 m/s in turbulent flow, based on the testing of a 2 DOF section model in ambient vibrations. Flutter velocity estimated from a time-domain analysis using 2 DOF model with the set of 4 indicial functions and the methodology presented above was 76.9 m/s [33], or 3% lower than the experimental result reported in smooth flow. The present load model provides thus a critical velocity estimate in a very good agreement with the wind tunnel test results.
In the present study, the coupled flutter oscillations of the prototype bridge involve the three lowest vertical modes of vibration. In previous studies based on multimodal analysis, Jakobsen and Hjorth-Hansen [13] highlighted the contribution of the second vertical mode, but not the third one, to flutter oscillations. Their work, based on the motion-induced forces deduced from the ambient vibration data, also concerned a slightly different, preliminary design of the bridge.
A flutter velocity of 76.3 m/s, as estimated in the present analysis, is also consistent with the values computed by Øiseth et al. [14] for the same structure in the frequency domain, using aerodynamic derivatives (79 m/s) and modified quasi-steady load estimates (76 m/s). A similar value (78.3 m/s) was obtained by Øiseth et al [15] in a time-domain analysis with a finite element model, using the so-called aerodynamic degrees of freedom to model the motion-dependent forces on the bridge represented by a simply supported beam.

Conclusions
This paper presents the response analyses of a full-bridge FE-model of the Hardanger Bridge. The objective is to study bridge response in time-domain at high wind speeds considering all components of the motion-dependent loads.
The flutter velocity is found to be more or less unaffected by the motion-induced drag components as well as by the lift and twisting moment induced from the horizontal motion of the bridge girder. On the other hand, the motion-dependent drag loads due to the horizontal motion of the girder contributed to the horizontal and the vertical transient response of the bridge.
The critical flutter velocity was estimated as 75.6 m/s, using the FE-model with a complete set of aeroelastic load components. That is 5% higher than the estimate based on multimodal analysis. The flutter mode shapes in the time-domain and the frequency domain are consistent and involve contributions of the three lowest symmetric vertical modes, together with the lowest torsional mode.
The in-wind torsional frequencies in the time-domain and in the frequency-domain complied within 1% at the design wind-speed. The flutter frequency was 20% higher in time-domain analyses compared to the results from multi-modal analysis. Thus, at the flutter velocity, the apparent reduction in torsional stiffness is distinctly lower using the FE-model compared to the results from the multimodal analysis.
The methodology applied in this work, provides detailed response estimates and can be utilized in general studies of bridge response to increase understanding of the complex aeroelastic behaviour of long-span bridges.