Identification of Flutter Derivatives of Bridge Decks in Wind Tunnel Test by Stochastic Subspace Identification T. Janesupasaeree and V. Boonyapinyo

Problem statement: Flutter derivatives are the essential parameters i n the estimations of the flutter critical wind velocity and the response s of long-span cable supported bridges. These derivatives can be experimentally estimated from wi nd tunnel test results. Generally, wind tunnel test methods can be divided into free decay test and buf feting test. Compared with the free decay method, the buffeting test is simpler but its outputs appea r random-like. This makes the flutter derivatives extraction from its outputs more difficult and then a more advanced system identification is required. Most of previous studies have used deterministic sy stem identification techniques, in which buffeting forces and responses are considered as noises. Thes e previous techniques were applicable only to the free decay method. They also confronted some diffic ulties in extracting flutter derivatives at high wi nd speeds and under turbulence flow cases where the bu ff ting responses dominate. Approach: In this study, the covariance-driven stochastic subspace id entification technique (SSI-COV) was presented to extract the flutter derivatives of bridge decks fro m the buffeting test results. An advantage of this method is that it considers the buffeting forces an d responses as inputs rather than as noises. Numeri cal simulations and wind tunnel tests of a streamlined thin plate model conducted under smooth flow by the free decay and the buffeting tests were used to validate the applicability of the SSI-COV method. Then, wind tunnel tests of a two-edge girder blunt type of Industrial-Ring-Road Bridge deck (IRR) were conducted under smooth and turbulence flow. Results: The identified flutter derivatives of the thin plate model by the SSI-COV technique agree wel l ith those obtained theoretically. The results from the thin plate and the IRR Bridge deck validat ed he reliability and applicability of the SSI-COV technique to various experimental methods and condi ti s of wind flow. Conclusion/Recommendations: The SSI-COV was successfully employed to identify flutter derivatives of bridge decks with reliable results. I is a proven technique that can be readily applie d to identify flutter derivatives of other bridge decks either by the free decay or the buffeting tests.


INTRODUCTION
Long-span cable-supported bridges are highly susceptible to wind excitation because of their inherent flexibility and low structural damping. Wind loads play an important role in the design of these structures. A wind-induced aerodynamic force can be divided into two parts: a buffeting force that depends on the turbulence of incoming flow and an aeroelastic force that originates in the interaction between the airflow and the bridge motion. The motion-dependent forces feed back into the dynamics of the bridge as aerodynamic damping and stiffness; the effect is termed 'aeroelasticity' and is commonly described via 'flutter derivatives'. The problems of aerodynamic stability including vortex-induced vibrations, galloping, flutter and buffeting, may have serious effects on the safety and the serviceability of the bridges. Among these, flutter is the most serious wind-induced vibration of bridges and may destroy the bridges due to diverging motions either in single or torsion-bending coupled mode. Notorious examples by the flutter phenomenon are the failures of the Brighton Chain Pier Bridge in 1836 and the original Tacoma Narrow Bridge in 1940. The flutter derivatives depend primarily upon the conditions of wind, the cross-sectional shape and the dynamic characteristics of the bridges. Nevertheless, no theoretical values exist for these derivatives for various bridge shapes except only for a simple thin plate section. A major research tool in these studies is, therefore, a wind tunnel test, in which a geometrically and aerodynamically representative scale model of a length of a bridge deck is mounted in a wind tunnel. The flutter derivatives are non-dimensional functions of wind speed, geometry and frequency of vibrations; therefore they can be applied directly to full-scale bridge in a piecewise manner.
The experimental method used for a determination of flutter derivatives can be grouped under two types, i.e. forced [1] and free vibration methods [2][3][4][5] . Having less emphasis on elaborate equipment, time and work; the free vibration method seems to be more tractable than the forced method. In the determination of flutter derivatives by free vibration method, the system identification method is the most important part required to extract these parameters from the response output of the section model. The free vibration method depends on system identification techniques and can be classified into two types, i.e. free decay and buffeting tests. In the free decay method, the bridge deck is given an initial vertical and torsional displacement. The flutter derivatives are based on the transient (i.e., free decay) behavior that occurs when the bridge deck is released. The buffeting test, on the other hand, uses only the steady random responses (i.e., buffeting responses) of bridge deck under wind flow without any initial displacement given to the model. Compared with the free decay method, the buffeting test is simpler in the test methodology, more cost effective and more closely related to real bridge behaviors under wind flow, but with a disadvantage that the outputs appear randomlike. This makes the parameters extraction more difficult and a more advanced system identification is required.
In most of the previous studies, flutter derivatives were estimated by deterministic system identification techniques that can be applied to the free decay method only. Examples of previous deterministic system identification that were applied to the free decay method included Scanlan's method [2] , Poulsen's method [3] , Modified Ibrahim method (MITD) [4] and Unified Least Square method (ULS) [5] . In these system identification techniques, the buffeting forces and their responses are regarded as external noises, the identification process then requires many iterations [3][4][5] . It also confronted with difficulties at high wind speeds where the initial free decay is drowned by buffeting responses [3][4][5] . Moreover, at high reduced wind speed, the vertical bending motion of the structure will decay rapidly due to the effect of the positive vertical aerodynamic damping and thus the length of decaying time history available for system identifications will decrease. This causes more difficulties to the deterministic system identification techniques [4,5] . In case of turbulence flow, the presence of turbulence in the flow is equivalent to a more noisy-input signal to the deterministic system identification. This made the extraction process more complicated and most likely reduced the accuracy of the flutter derivatives identified [3,4] . In addition, due to test technique, the free decay method is impractical to determine flutter derivatives of real bridges in field.
On the other hand, the buffeting test uses random responses data of bridge motion from wind turbulence only. This mechanism is more closely related to a real bridge under wind flow and is applicable to real prototype bridges. The method costs less and simpler than the free decay since no operator interrupts in exciting the model. However, as wind is the only excited source, it results in low signal-to-noise ratio, especially at low velocity and therefore a very effective system identification is required. None of the aforementioned system identification techniques is applicable to the buffeting responses tests. System identification techniques can be divided into two groups, i.e., deterministic and stochastic. If the stochastic system identification technique [6][7][8][9] is employed to estimate the flutter derivatives of a bridge deck from their steady random responses under the action of turbulent wind, the above-mentioned shortcomings of the deterministic system identification technique can be overcome. The reason is that the random aerodynamic loads are regarded as inputs rather than noises, which are more coincident with the fact. Therefore, the signal-to-noise ratio is not affected by wind speed and the flutter derivatives at high reduced wind speeds are more readily available. These aspects give the stochastic system identification methods an advantage over the deterministic system identification.
Many stochastic system identification methods have been developed during the past decades, among which the Stochastic Subspace Identification (SSI in short) [7,8] has proven to be a method that is very appropriate for civil engineering. The merit points of SSI are: (1) The assumptions of inputs are congruent with practical wind-induced aerodynamic forces, i.e. stationary and independent on the outputs; (2) Identified modes are given in frequency stabilization diagram, from which the operator can easily distinguish structural modes from the computational ones; (3) Since the maximum order of the model is changeable for the operator, a relatively large model order will give an exit for noise, which in some cases can dramatically improve the quality of the identified modal parameters; (4) Mode shapes are simultaneously available with the poles, without requiring a second step to identify them. There are two kinds of SSI methods, one is data-driven and the other is covariance-driven.
In this study, the covariance-driven stochastic subspace identification method is used to estimate the flutter derivatives from random responses (buffeting) under the action of smooth and turbulent wind. Tests are also carried out with the free decay method (single and two-degree-of-freedom) in order to examine the robustness of the present technique that the results are not affected by test methods used. To validate the applicability of the present technique, first numerical simulations are performed then sectional-model tests of a quasi-streamlined thin plate model, which is the only section that theoretical flutter derivatives exist, are performed under smooth flow. Encouraged by the success in the evaluation process, the flutter derivatives of a real bridge are determined. The twoedge-girder type blunt section model of Industrial-Ring-Road Bridge (IRR in short), a cable-supported bridge with a main span of 398 m in Samutprakan province, Thailand, was tested both in smooth and turbulence flow. Tests were conducted in TU-AIT Boundary Layer Wind Tunnel in Thammasat University, the longest and the largest wind tunnel in Thailand.

MATERIALS AND METHODS
Theoretical formulation of covariance-driven SSI: The dynamic behavior of a bridge deck with two Degrees-Of-Freedom (DOF in short), i.e., h (bending) and α (torsion), in turbulent flow can be described by the following differential equations [9,10] : Where: m and I = The mass and mass moment of inertia of the deck per unit span, respectively ω i = The natural circular frequency ξ i = The modal damping ratio (i = h, α) L se and M se = The self-excited lift and moment, respectively L b and M b = The aerodynamic lift and moment The self-excited lift and moment are given as follows [11] : Where: ρ = Air mass density; B is the width of the bridge deck U = The mean wind speed at the bridge deck level k i = ω i B/U = The reduced frequency (i = h, α) H i * and A i * = The so-called flutter derivatives, (i = 1, 2, 3,4) which can be regarded as the implicit functions of the deck's modal parameters The alternate form of self-excited forces is as Eq. 2 but without the factor 1/2 [3] .
The aerodynamic lift and moment can be defined as [10] : Where: C L , C D and C M = The steady aerodynamic force coefficients C′ L and C′ M = The derivatives of C L and C M with respect to the attack angles, respectively u(t) and w(t) = The longitudinal and vertical fluctuations of wind speed, respectively χ L and χ M = The lift and moment aerodynamic admittances of the bridge deck By moving L se and M se to the left side and merging the congeners into column vectors or matrices, Eq. 1 can be rewritten as follows: Where: = The gross damping matrix, i.e., the sum of the mechanical and aerodynamic damping matrices = The gross stiffness matrix The fluctuations of wind speed u(t) and w(t) in Eq. 3 are random functions of time, so the identification of flutter derivatives of bridge decks can be simplified as a typical inverse problem in the theory of random vibration and thus can be solved by the stochastic system identification techniques. Let: then Eq. 4 is transformed into the following stochastic state equations: The discrete form of Eq. 7 can be written as: Defining: and combining Eq. 9 and 10 we obtain the following Lyapunov equations for the state and output covariance matrices: From (8) and (9), it can be deduced: Defining a block Toeplitz T 1|i as: one can infer from the definition of covariance matrix that T 1|i can be expressed as the product of two block Hankel matrices Y f and Y p : where, Y f and Y p are composed of the 'future' and 'past' measurements, respectively:  308 Next, applying the factorization property to T 1|i by the singular value decomposition yields: Where: U, S and V = Orthonormal matrices S = A diagonal matrix containing positive singular values in descending order The number of nonzero singular values indicates the rank of the Toeplitz matrix. The reduced diagonal matrix S 1 is obtained by omitting the zero singular values from the matrix S. Matrices U 1 and V 1 are obtained by omitting the corresponding columns from the matrices U and V respectively. Now realizations of the system matrices are almost achieved. Matrix A is realized by using factorization of a shifted Toeplitz matrix T 2|i+1 that has similar structure as of T 1|i but consists of covariance from lag 2-2i. In a manner similar to the classical Eigensystem Realization Algorithm (ERA in short), one can find: where, N is model order, i.e., the maximum number of modes to be computed. Thus, the modal parameters can be determined by solving the eigenvalue problem of the state matrix A. By now, the theoretical formulation of covariance-driven SSI has been achieved. According to Eq. 16-19, a different combination of i, j and N will give a different state matrix and thus a different pair of modal parameters. Therefore, modal parameters should be derived from a series of combinations, rather than a single combination. In the process of identification, N or i should be given in series for certain values of j in order to obtain a frequency stability chart. Solving the eigenvalue problem of the state matrix A by the pseudo-inverse method yields: Where: Ψ = The complex eigenvector matrix Φ = The mode shape matrix Λ = A diagonal matrix composed of the complex poles of the system Different combinations of i, j and N are employed to derive the modal parameters statistically [3,6] .
Once the modal parameters are identified, the gross damping matrix C e and the gross stiffness matrix K e in Eq. 4 can be readily determined by the pseudo-inverse method: † * e e 2 * * 2 * * where the superscript*denotes the complex conjugate of the corresponding term. Let: where, C 0 and K 0 are the 'inherent' damping and stiffness matrices, respectively. Thus, the flutter derivatives can be extracted from the following equations: Numerical simulation tests: In order to validate the applicability of the covariance-driven SSI technique in flutter derivatives estimation of bridge decks, numerical simulations of signals from different test methods are first carried out. The numerical tests included two syntheses but well controlled cases: two uncoupled degrees of freedom and two coupled degrees of freedom (simulated response including the motion induced aeroelastic terms). Both cases are first excited in the transient (i.e., free decay) motion and then by a white noise loading process. Measurement noises are also added by a white noise process with a standard deviation equal to 10% of the standard deviation of the original responses, in order to investigate the effect of measurement noise.
Two uncoupled degrees of freedom; free decay: Transient responses time-series were obtained by direct calculations of the displacement values for N = 4096 discrete time stations, with 'sampling' interval ∆t = 0.02 sec (fs = 50 Hz). Structural modal properties used in this simulation were chosen from the previously tested sectional model of the Great Belt Bridge [12] . The modal matrices are given per unit length as: i.e., f ho = 1.9472 Hz, f θ0 = 5.7573 Hz, ξ h0 = 0.0053, ξ θ0 = 0.0056, where damping ratios, ξ, are representatives for the range of small amplitudes. The damping ratios were then multiplied in turn with 5, 10, 20 and 40, in order to cover the values of total damping (structural + aerodynamic) which could be presented in vibration of model section under wind flow. Values as high as ξ = 0.2 could be expected for the vertical degree of freedom under wind flow. Frequency and damping ratio estimates are practically identical to the preset values (less than 0.5% for the highest damping case). The system matrices are also excellent even for the short useful signal case with only a few cycles of vibration motion. In the case where 10%-measurement noises were added, identified frequencies were changed at lesser than 0.8%. Damping ratios were changed at most by 2% except in the case of the lowest damping case which was 5.4%. The diagonal terms of the estimated system matrices (frequency and damping matrices) are also identical to the preset values. Estimates of diagonal terms are distorted within 1% except only for the case with lowest damping case in which values are within 2.82%.
Two coupled degrees of freedom; free decay and buffeting responses: The next step in the simulation was a simulation test with full effective stiffness and damping matrices (i.e., coupled degrees of freedom) and with lift and moment forces of the white noise type, as assumed in the SSI-method. For the mean-wind speed of 10.26 m sec −1 and the aerodynamic derivatives assumed according to the values reported for a similar bridge cross-section [12] , the effective structural matrices were pre-set at: The response time-series were simulated for both free decay and buffeting responses under turbulence wind with 10% turbulence intensity; then measurement white noises were superimposed on the simulated response. The free decay response time-series were computed by constant acceleration method and samples are as shown in Fig. 1. The SSI-COV method, applied to these responses data, returned the effective structural matrices with the deviation from the pre-set ones (C and K) in percentage as: % % 0.66 3.00 0.14 0.05 C , K 0.16 0.26 4.26 0.08 Superimposing 10% measurement white noise on the simulated response made the structural matrices differed from those of the noise-free cases within 3%. The response time-series were also simulated for the case of buffeting responses where wind turbulence is the only excited source. The effective stiffness and damping matrices were taken as in the case of transients; examples of response time-series are as shown in Fig. 2. Buffeting responses required longer data records (20,000 data points in the present study) as compared to that in the free decay case (4096 data points) to yield acceptable results. Estimates of the frequencies and damping ratios agree well with preset values where precisions are within 0.5 and 2%, respectively. The diagonal terms in stiffness and damping matrices also agree well with preset values where the differences are less than 1% except for the C 11 (related to vertical damping) where the difference is around 2.5%. The most differences in the off-diagonal terms are K 21 and C 21 which are related to A 4 * and H 2 * , respectively. In the case of 10%-measurement noise added, the deviation of the reconstructed matrices from the pre-set ones, in percentage, is:   The model was suspended by eight springs outside the wind tunnel (Fig. 4). To simulate a bridge section model with 2DOFs, i.e. vertical bending and torsion, piano wires were used to prevent the motion of the model in longitudinal direction; this can be shown from Fig. 5, the schematic diagram of the top view of the test setup. Two piezoelectric acceleration transducers were mounted at the mid length of the model to capture the acceleration signals. The responses of the models were captured by the acceleration transducers and then the vertical and torsional responses can be respectively obtained by: where, x 1 and x 2 are the measurements of transducers 1 and 2, respectively; l is the space between transducers.

RESULTS AND DISCUSSION
Case 1: Thin plate model under smooth flow: A quasistreamlined thin plate (Fig. 6) was first selected for wind tunnel test. The width to height (thickness) ratio of the plate is about 22.5. Table 1 shows the main parameters of the model. The extraction of flutter derivatives of the thin plate, using the SSI-COV technique, were performed on results from three types of tests, namely, (a) Single-Degree-Of-Freedom (SDOF) motion tests [2] , (b) free decay coupledmotion test (2DOFs) and (c) buffeting coupled-motion test (2DOFs). Typical test results showing responses from the bridge model are in Fig. 7 and 8. The responses for the free decay and the buffeting tests are sampled at the rates of 1000 and 200 Hz, respectively.  The results are then removed trend and re-sampled at 250 and 50 Hz, respectively. The covariance-driven SSI technique is applied to identify modal parameters from these data and a pseudo-inverse method is applied to estimate the stiffness and damping matrices. The flutter derivatives are estimated by Eq. 23 and reported in the form of Eq. 2 but without the factor 1/2.
Comparisons between SDOF and 2DOF-coupledmotion tests: free decay method: Figure 9 and 10 compare the flutter derivatives of the thin plate that are estimated by the SSI-COV technique using the above mentioned three test methods together with the Theodorsen's theoretical values [13] . Unless otherwise noted, at any wind speed, H 1 * , H 4 * , A 1 * and A 4 * which are associated with the vertical motion were calculated using the frequency n 1 (lower). In addition, the derivatives H 2 * , H 3 * , A 2 * and A 3 * which are associated with the torsional motion were calculated using the frequency n 2 (higher).
The direct flutter derivatives H 1 * and H 4 * as found from the single-degree-of-freedom vertical-motion tests and A 2 * and A 3 * as found from the single-degree-offreedom torsional-motion tests are also plotted and compared with those from the coupled-motion tests. The results are shown in Fig. 9 and 10. The near perfect match shows that the direct-flutter derivatives are indeed not affected by the motion along the other degree of freedom, as predicted by theory i.e., those flutter derivatives associated with h motion are not affected by α motion and vice versa. It also demonstrates the reliability of both the coupled-motion tests and the system identification method (SSI-COV).

Comparisons of coupled-2DOF motion tests between the free decay and the buffeting tests:
The flutter derivatives found from both the free decay and the buffeting tests for the coupled-2DOF cases are compared in Fig. 9 and 10. The results show good agreement between the two methods. This validates the ability of the system identification method (SSI-COV) to apply with both the free decay and the buffeting tests although it was developed from a stochastic model (i.e., white noise loading assumption). However, when a relatively heavy model is excited at a very low reduced wind velocity, i.e., low wind energy, it becomes more difficult to extract the flutter derivatives from the buffeting responses.
The results also show that identified flutter derivatives agree well with the theoretical ones.  [2,3,5,10] .
Case 2: Section model of IRR Bridge: Encouraged by the success in the thin plate model, the flutter derivatives of IRR Bridge, a cable-supported bridge with 2-edge girder, as shown in Fig. 11, were estimated by the SSI-COV technique. The IRR Bridge has a main span of 398 m. The deck consists of a concrete deck slab and a web of steel girders. The deck is supported by two cable planes at outside edge girders. A 2-edgegirder bridge section with A-shape pylons has good cost performance, but at the same time the bridge crosssection is known to be aerodynamically unstable at high wind speed. Table 2    The difference of A 4 * identified from both tests, seems to be negligible as effect of this derivative is usually considered to be less significant. In smooth flow, the most important derivative A 2 * are steadily increased (more negative) up to the reduced wind velocity around 3 and then started to decrease. This sign reversal is the primary factor toward the SDOF-torsional instability ("stall flutter") for bluff type sections. The torsional flutter was found at the reduced wind speed around 4.7.

Effect of turbulence:
Most of the prototype bridges are submerged in turbulent wind, therefore, detailed investigations of the effects of turbulence on the flutter derivatives is significant. Almost all the wind tunnel tests for flutter derivatives have been generally carried out in smooth flow. Although few researchers have studied the problem using wind tunnel tests, results and the identification methods were individually proposed [4,14] and the results are still debatable and inconclusive. For streamlined section, tests showed little effect [4,14] , while tests on a rectangular box girder bridge showed galloping in smooth flow [15] .
From Fig. 12 and 13, it can be found that the influence of flow type on H 4 * and A 3 * , i.e., flutter derivatives related to direct aerodynamic stiffness, seems to be negligible. Though, the value of H 4 * from turbulence flow is somewhat lesser than that in the smooth flow case, it affected only the second decimal digit of the frequency value. The influence also has negligible effect on H 1 * and H 2 * i.e., direct and cross derivatives related to vertical and torsional aerodynamic dampings, respectively. On the other hand, the more important A 1 * A 2 * and H 3 * , show rather noticeable deviations from those in smooth flow, especially at high reduced wind speeds. The most important effect is that the reduced wind speed corresponding to the reversed sign of the torsional aerodynamic damping A 2 * increased in turbulent flow. It shows that turbulence tends to make bridges more aerodynamically stable by delaying torsional flutter. The deviations of flutter derivatives may reveal the fact that for those bridges with bluff type sections similar to IRR Bridge, the effects of turbulence can be significant. Hence, the wind tunnel tests of such bridges for flutter derivative estimation should be carried out in turbulent flow as well.

CONCLUSION
A theoretical model based on the covariance-driven SSI technique was proposed to extract the flutter derivatives of bridge deck sectional models from coupled two-degree-of-freedom system by free decay and buffeting responses. An advantage of the adopted SSI-COV technique is that it considers the buffeting forces and responses as inputs instead of as noises as typically assumed. The conclusions of this study are as follows: • derivative as obtained from the two test methods but they agree in trend. We also observed the sign reversal of the A 2 * derivative as the reduced wind speed reached the value of 4.7. This indicates that this bridge section is susceptible to flutter instability at high wind speed • The test result of bluff section model of the IRR Bridge under turbulence wind revealed that the most important and positive effect of the turbulence is that it tends to make the bridge more aerodynamically stable by delaying the sign reversal of the aerodynamic damping A 2 * . This may help explain that for those bridges with bluff type sections similar to the IRR Bridge, the effects of turbulence can be significant. Hence, the wind tunnel tests of such bridges for flutter derivatives estimation should be carried out in turbulent flow as well Applying the proposed SSI-COV technique to the buffeting test yields a straightforward, cost effective and reliable system identification process that can be used to identify flutter derivatives of various bridge decks. It, however, has some limitations. For example, it becomes more difficult to extract the flutter