Wind flows in the atmospheric boundary layer are always turbulent. Any study of a wind-induced vibration problem must confront this issue, either by matching turbulence characteristics completely or by acknowledging uncertainties in the conclusions as a result of imperfect simulations. Not many studies in this field have focused clearly on the effects of turbulence on aeroelastic forces. Scanlan and Lin [1] are pioneers who used a truss bridge deck section model and concluded that there is an insignificant difference in flutter derivatives (FDs) between smooth and turbulent flows. However, Huston [2] conducted a test on a model of the Golden Gate Bridge deck section and obtained results that were different from those obtained by Scanlan and Lin [1] and Haan and Kareem [3]. In studies on the effect of turbulent flows on FDs, a free-vibration technique using a section model with the application of a system identification technique to extract FDs is widely used. Various system identification techniques have been developed by many authors: These include the extended Kalman filter algorithm [4], the modified Ibrahim time domain (MITD) [5], the unifying least-squares method [6], and the iterative least-squares method [7]. In these systems, buffeting force and its response are considered to be external noise, which causes increased difficulty at high wind velocity. Bartoli and Righi [8] used the combined system identification method (CSIM), which is based on Sarkar’s MITD, to simultaneously extract all FDs from a two degrees-of-freedom (2DOF) rectangular section model. Their conclusion was that the identification of FDs in a turbulent flow succeeded in spite of the difficulties caused by locally induced noise due to signature turbulence. The main reason for their success is that the CSIM is a deterministic system identification, and the effects of turbulence are regarded as a noisy-input signal to the system, which results in additional problems in the identification process. Nikitas et al. [9] extracted FDs from ambient vibration data from full-scale monitoring using a more elaborate stochastic identification technique (the covariance block Hankel matrix, CBHM) [10]; that study also illustrated the applicability of system identification techniques to extract valuable results from field measurement data. Kirkegaard and Andersen [11] compared three state-space systems: stochastic subspace system identification (or simply stochastic system identification (SSI)), the matrix block Hankel (MBH) stochastic realization estimator, and the prediction error method (PEM). The SSI was found to give a good result for estimated modal parameters and mode shapes, the MBH was found to give poor estimates of the damping ratios and the mode shapes compared with the other two techniques, and the SSI was found to be approximately ten times faster than the PEM.

In addition to the issues described above, there is a shortcoming involved in the free-vibration technique. The extraction of FDs cannot be executed accurately in a high wind speed range because the aerodynamic damping of the vertical mode is too high, and the vertical free-vibration data are rapidly damped out. Based on these considerations, the idea of applying SSI to estimate the FDs from the gust responses of a truss bridge deck section was developed. In this paper, the turbulence effects on the FDs of a truss bridge deck are investigated by employing a section model wind-tunnel test in a turbulent flow. Output-only SSI is applied to extract FDs from gust response. Tests were also carried out using the free-vibration method for comparison with the proposed method.

A wind-tunnel test was conducted in a closed-circuit wind tunnel at Yokohama National University. The working section is 1.8 m wide and 1.8 m high, and the investigated profile is a truss bridge deck section (Fig. 1). It was fabricated from wood with a scale of 1:80 to represent the cross-section of a long-span suspension bridge. The width and depth of the section model are 363 mm and 162.5 mm, respectively. The unit length mass is 8.095 kg·m⁻¹ and the unit length moment of inertia is 0.2281 kg·m²·m⁻¹. The first vertical frequency and damping ratio are 1.869 Hz and 0.00509, respectively, and the first torsional frequency and damping ratio are 3.296 Hz and 0.004186, respectively. The section model was attached to a rigid frame and supported at each corner by a linear spring with stiffness k. The mounting position of the spring was adjusted such that the elastic center and the gravity center of the cross-section coincided. The tests were carried out in both smooth and turbulent flows. The turbulent flows used in this study were generated with a biplane wooden grid, and the turbulent properties were controlled by changing the distance to the model.

For the purpose of matching the power spectrum of the turbulent flow of a model with that of the full-scale version, Katsuchi and Yamada [12] introduced the concept of reduced turbulence intensity (I_r), which can be written as follows:

where I_u is the along-wind turbulence intensity, defined by I_u = σ_u/U, for which σ_u is the standard deviation of the turbulent wind velocity and U is the mean wind velocity; D is the height of the section model; and $L_{\mathrm{u}}^{X}$ is the integral length scale for the turbulent component in the longitudinal direction. In this study, the integral length scale is defined by

where n_peak is the frequency at which the reduced power spectrum reaches the maximum.

Table 1 depicts the results of three turbulent flow parameters, I_u, $L_{\mathrm{u}}^{X}$, and I_r, corresponding to different grid-to-model distances. I_r increases proportionally with I_u, but the inverse is true for $L_{\mathrm{u}}^{X}$.

The turbulence intensity and the integral length scale do not fully describe the properties of oncoming turbulent flows. Nakamura and Ozono [13] showed that small-scale turbulence affects flow fields and aerodynamic parameters more than large-scale turbulence. Therefore, the power spectral density (PSD) of the turbulence was also quantified for this research. Fig. 2(a) shows the non-dimensional PSD function of the along-wind turbulence, together with the von Karman and Eurocode 1 spectra. Compared with the von Karman spectrum, the measured data coincided well in the high-frequency range and was a little higher in the low-frequency range. Turbulent energy is generated in larger eddies (low frequency). For most structures, these low-frequency fluctuations give no significant response contribution. Fig. 2(b) shows three spectra obtained at different I_r. The values of the PSD function increase as I_r increases.

Tests were conducted in both smooth and turbulent flows. The objective of this test was to quantify the effect of oncoming turbulent flows on the dynamic responses of the section model.

Fig. 3 illustrates the vibration amplitudes of 2DOF (heaving and torsional modes) versus the reduced wind velocities (V_r = U/fB, where f is the frequency and B is the width of the section model) under smooth and various turbulent flows. In smooth flow, vertical vibration is limited when V_r ranges from 0 to 9, and then increases considerably; however, vertical divergent vibration does not occur in this test (Fig. 3(a)). On the other hand, the torsional displacement is very small until a sudden increase and flutter occur at a V_r of about 5.7 (Fig. 3(b)).

When the model is immersed in turbulent flows, vertical and torsional motions appear in spite of small reduced wind velocities. The vertical response increases proportionally with V_r, and when the I_r increases, the amplitude of vibration increases slightly. In the turbulent flows, vertical divergent vibration does not occur. In the cases of I_r = 6.97% and I_r = 11.09%, the torsional displacement gradually increases with V_r and flutter occurs at V_r = 7.2 and V_r = 7.7, respectively. The flutter speeds are higher than those in smooth flow (V_r = 5.7). Moreover, in the case of I_r = 21.02%, flutter occurs at V_r = 8.6. (The critical wind speed is defined at the amplitude of 0.5°.)

In general, turbulent flows on a section model induce larger vibration than smooth flow, and the vibration proportionally increases with the wind velocities; however, the vibration does not suddenly increase in turbulent flows as it does in smooth flow. The motion under turbulent flows is known as buffeting response; it affects the service state of bridges by causing issues such as fatigue problems.

In this study, the FDs of a bridge deck in turbulence are identified using an SSI method. The theories applied in this study are mainly based on the work of Peeters and Roeck [14].

When considering a 2DOF section model of a bridge deck in turbulent flow, fluctuating wind loads that act on the deck can be expressed by the linear superposition of a self-excited force and a buffeting force as follows:

where m and I are the mass and polar moment of the inertial per unit length, respectively; h and α are heaving and torsional displacement, respectively; (⋅) represents time differentiation; ω_h = 2πf_h and ω_α = 2πf_α are the natural circular frequencies of the heaving and torsional mode, respectively; ξ_h and ξ_α are the damping ratio of heaving and torsion, respectively; L_b and M_b are the buffeting forces in the vertical and torsional directions, respectively; and L_se and M_se are the self-excited lift and pitching moment, respectively, given by

where ρ is the air density, U is the mean wind velocity, B is the width of the bridge deck, K_k = ω_kB/U is the reduced frequency (where k = h, α), and Hi∗ and Ai∗ (where i = 1, 2, 3, 4) are the FDs.

By substituting Eq. (4) into Eq. (3) and moving the aerodynamic damping and stiffness terms to the left-hand side, Eq. (3) can be rewritten as follows:

where $ q(t)=\left[\begin{array}{ll}h(t) & \alpha(t)\end{array}\right]^{\mathrm{T}} $ is the generalized buffeting response, $ f(t)=\left[\begin{array}{ll}L_{\mathrm{b}} & M_{\mathrm{b}}\end{array}\right]^{\mathrm{T}} $ is the buffeting force, f(t) is factorized into matrix B₂ and input vector u(t), M is the mass matrix, C^e is the total damping matrix including the structural and aerodynamic damping, and K^e is the total stiffness matrix including the structural stiffness and aerodynamic stiffness.

The second-order differential equation, Eq. (5), can be transformed into a first-order state equation, Eq. (6):

where

and where A_c is the designated state matrix with a size of 4 × 4, x(t) is the state vector, B_c is the input matrix, and I_u is a unit matrix.

The combination of the state equation and the observation equation fully describe the input and output behaviors of the structural system; as such, it is named the state-space system.

where y(t) is the output vector, C_c is the output matrix, and D_c is the direct transmission matrix in continuous time.

Eq. (7) is a deterministic state-space model in continuous time. Continuous time implies that the expression can be evaluated at each time instant. Deterministic implies that the input and output quantities can be measured exactly. This is not practical, because the measurements are mostly sampled at a discrete time. In addition, it is impossible to measure all DOFs, and measurements always cause disturbance effects. For all these reasons, the continuous deterministic system is converted into a suitable form—that of a discrete-time stochastic state-space model—as follows:

where $ x_{k}=x(k \Delta t)=\left(\begin{array}{ll}q_{k} & \dot{q}_{k}\end{array}\right)^{\mathrm{T}} $ is the discrete-time state vector containing the discrete sampled displacement q_k and velocity q̇k; w_k is the process noise due to disturbances and modeling error; v_k is the measurement noise due to sensor inaccuracy; and A and C are discrete state and output matrices, respectively. It is assumed that w_k and v_k are zero mean and that their covariance matrix is as follows:

where the indices p and q are time instants, E is the expectation operator, and δ_pq is the Kronecker delta. The correlations $ E\left[\begin{array}{ll}w_{p} & w_{q}^{\mathrm{T}}\end{array}\right] $ and $ E\left[\begin{array}{ll}v_{p} & v_{q}^{\mathrm{T}}\end{array}\right] $ are equal to zero in the case of different time instants. $ Q=E\left[\begin{array}{ll}w_{k} & w_{k}^{\mathrm{T}}\end{array}\right], R=E\left[\begin{array}{ll}v_{k} & v_{k}^{\mathrm{T}}\end{array}\right], \text { and } S=E\left[\begin{array}{ll}w_{k} & v_{k}^{\mathrm{T}}\end{array}\right]$. It is further assumed that the stochastic model x_k, w_k, and v_k are mutually independent and zero mean. It can be proven that the output covariance $ R=E\left[\begin{array}{ll}y_{k+1} & y_{k}^{\mathrm{T}}\end{array}\right]$ for any arbitrary time-lags iΔt can be regarded as the impulse response of the deterministic linear time-invariant system A, C, G, where $ G=E\left[\begin{array}{ll}x_{k+1} & y_{k}^{\mathrm{T}}\end{array}\right]$ is the next state-output covariance matrix, as shown in Eq. (10).

Eq. (10) is called the Lyapunov equation, and indicates that the output covariance can be regarded as impulse responses. Therefore, the theoretical application of the stochastic system can go back to an eigensystem realization algorithm (ERA) method [15].

The output measurement data obtained from l sensors (in this study, l = 2 for heaving and torsion) are as follows:

The output data are assembled in a block Hankel matrix (H) with 2i block rows and j columns. The Hankel matrix can be divided into two parts: The upper part is the past output and the lower is the future output, as follows:

Like the CBHM or ERA method, data-driven stochastic system identification (SSI_data) implements directly with the output of experimental data, without converting output data to correlation, covariance, or spectra [16]. The main step of SSI_data is a projection of the row space of the future outputs (Y_f) into the row of past outputs (Y_p). The orthogonal projection P_i is defined as follows:

The orthogonal projection is implemented by the QR factorization of the block Hankel matrix shown in Eq. (12). This is defined as follows:

where $ Q \in \mathbb{R}^{j \times j} $ is an orthogonal matrix $ Q^{\mathrm{T}} Q=Q Q^{\mathrm{T}}=I_{j} \text { and } R \in \mathbb{R}^{2 l i \times j}$ is a lower triangular matrix. Because 2li < j, it is possible to reject the zeros in R and the corresponding zeros in Q.

Substituting the QR factorization of the output Hankel matrix, Eq. (15), into Eq. (13) yields the simple expression of the projection:

The key of SSI is that the orthogonal projection P_i is factorized into the product of the observability matrix O_i and the Kalman filter state sequence $ \hat{X}_{i}$.

The observability matrix O_i and the Kalman filter sequence $ \hat{X}_{i}$ are obtained by applying single-value decomposition to the projection matrix:

Comparing Eq. (17) and Eq. (18) gives

where $ (\cdot)^{\dagger}$ represents the pseudo-inverse of a matrix.

If the past and future outputs of the Hankel matrix are shifted, a time-shift projection is achieved:

where

O_i−1 is obtained from O_i after deleting the last l rows. The shifted state sequence can be computed in Eq. (20) as follows:

From Eq. (19) and Eq. (22), the Kalman state sequences X̂i and X̂i+1 are obtained using only output data. The state and controllability matrices can be recovered from the over-determined set of linear equations, obtained by extending Eq. (8):

where Y_i|i is a Hankel matrix with only one block row. Since the Kalman state sequence and the outputs are known, and the residual $ \left[\begin{array}{ll}\rho_{w}^{\mathrm{T}} & \rho_{v}^{\mathrm{T}}\end{array}\right]^{\mathrm{T}}$ are uncorrelated with $ \hat{X}_{i}$, the set of equations can be solved for A and C using the least-squares method:

The modal parameters of the system can be obtained by solving the eigenvalue problem for the state matrix:

where Ψ is the complex eigenvector, Λ is the complex eigenvalue and the diagonal matrix, and Φ is the mode shape matrix. When the complex modal parameters are known, the total damping matrix C^e and total stiffness matrix K^e in Eq. (5) are determined by

Let $ \bar{C}^{\mathrm{e}}=M^{-1} C^{\mathrm{e}}, \bar{K}^{\mathrm{e}}=M^{-1} K^{\mathrm{e}}, \bar{C}=M^{-1} C^{0}, \bar{K}=M^{-1} K^{0}$, where C⁰ and K⁰ are the structural damping and stiffness matrix of the system, respectively, under still air conditions.

Thus, the FDs of a 2DOF system can be defined as follows:

The buffeting and decay responses were acquired at a sampling frequency of 100 Hz. Fig. 4 shows the time-history responses of the model at the mean wind speed of 6.7 m⋅s⁻¹.

In a high wind speed range, the aerodynamic damping of the heaving mode is too high and the vertical free response is too short, so the extraction of FDs cannot be accomplished with high accuracy. In addition, it is not practical to use a free-decay mechanism to describe real bridge behavior under field wind excitation. On the other hand, the extraction of FDs from the buffeting response more closely reflects full-scale bridge behavior in a turbulent wind field. The bridge deck section model will vibrate under the excitation of turbulent flow even at a low wind velocity.

This method is simpler than the free-vibration technique because no operator corrupts by exciting the section model. Fig. 5 and Fig. 6 show the FDs of the bridge deck extracted by the SSI_data method from both the free-decay and buffeting responses of 1DOF and 2DOF systems under turbulence flows (I_r = 11.09%). In general, most FDs are in good agreement with both the free-decay response and the buffeting response of 1DOF and 2DOF systems. The heaving-related damping FD ($H_{1}^{*}$) extracted from the buffeting response is slightly higher than that obtained from the free-decay response. The coupled terms ($H_{2}^{*}$ and $H_{3}^{*}$) extracted from buffeting responses are more scattered than those from free-decay responses, particularly at a highly reduced wind velocity. The trends for $A_{3}^{*}$ of the FDs are similar in both cases. In this study, the sectional profile was a truss bridge deck section where only torsional flutter occurred; therefore, $A_{2}^{*}$ is the most important derivative. The $A_{2}^{*}$ extracted from the buffeting response is scattered at low reduced wind velocities, but the scatter decreases more at a high reduced wind velocity for the buffeting response than for the free-decay response. The trends of the results from the buffeting response and free-decay response are closely coincident.

Fig. 7 and Fig. 8 illustrate the FDs of the heaving and torsional modes under smooth and turbulent flows with different I_r. It was found that under smooth flow, the FD $H_{1}^{*}$ decreases faster than that under turbulent flows. This is because the damping ratio of the heaving mode under smooth flow is higher than under turbulent flow. Turbulence has a very small effect on the vertical and torsional frequency terms $H_{4}^{*}$ and $A_{3}^{*}$, and these values extracted from the buffeting response are somewhat lower in turbulence than in smooth flow. The off-diagonal terms $H_{2}^{*}$, $H_{3}^{*}$, $A_{1}^{*}$, and $A_{4}^{*}$ fluctuate around zero, which means that in this experiment, coupled vibration did not appear. Of these derivatives, the torsional damping term $A_{2}^{*}$ plays an important role in torsional flutter stability, since its positive/negative value corresponds to the aerodynamic instability/stability of the torsional flutter. As shown in Fig. 8, under smooth flow, $A_{2}^{*}$ is positive at high reduced wind velocity V_r = 5.2 and coincides with the negative total torsional damping. The significant effects of turbulent flows on FDs are particularly illustrated for the aerodynamic torsional damping term $A_{2}^{*}$: The positive value corresponds to a V_r of around 6.5-7.8 under I_r = 6.97% and I_r = 11.09%, respectively, whereas in the case of I_r = 21.02%, flutter does not occur up to V_r = 8. On the other hand, the effects of different turbulent intensities on FDs are fairly modest. The influence of turbulence on FDs will depend on the section. Sarkar et al. [5] found a small effect for a streamlined section, whereas the test on a truss bridge deck showed an appreciable effect, which is clearly shown by the torsional damping term $A_{2}^{*}$.

There are two factors that may affect FDs in a turbulent flow: the spanwise turbulence coherence and the aerodynamic admittance. In Fig. 7 and Fig. 8, turbulence effects can be seen in some FD terms. Significantly different trends can be seen in the diagonal terms ($H_{1}^{*}$, $H_{4}^{*}$, $A_{2}^{*}$, and $A_{3}^{*}$) of the smooth and turbulent flows. On the other hand, this difference is not significant in Fig. 5 and Fig. 6, where different identification procedures (buffeting/free decay) were applied for the same turbulence. It may be concluded from these observations that aerodynamic admittance causing a buffeting response does not significantly affect FD identification, although the difference in spanwise turbulence coherence between smooth and turbulent flows does have a significant effect on FD identification.

In order to confirm the results of identified FDs under buffeting responses, the flutter critical wind velocity (V_cr) was obtained from an equation of motion of a 2DOF system:

where $ \begin{array}{l}M=\left[\begin{array}{cc}m & 0 \\0 & I\end{array}\right], C=\left[\begin{array}{cc}2 m \xi_{h} \omega_{h} & 0 \\0 & 2 I \xi_{\alpha} \omega_{\alpha}\end{array}\right], K=\left[\begin{array}{cc}m \omega_{h}^{2} & 0 \\0 & I \omega_{\alpha}^{2}\end{array}\right], F=\left[\begin{array}{cc}L_{h} & L_{\alpha} \\M_{h} & M_{\alpha}\end{array}\right], \text { and } q=\left[\begin{array}{l}h \\\alpha\end{array}\right].\end{array}$

For a stability check, only the self-excited force is considered, and L_h, L_α, M_h, and M_α are self-excited force components defined by

where L_hR, L_hI, L_αR, L_αI, M_hR, M_hI, M_αR, and M_αI are self-excited force coefficients (FDs) that can be compared with those using Scanlan’s format, as follows:

The FDs of the truss bridge deck section that are provided in this study were approximated polynomials of the results from Fig. 7 and Fig. 8.

Assuming sinusoidal motion q = q₀exp(iωt), and since structural damping of a long-span bridge can be negligibly small, the damping matrix in Eq. (28) can be dropped. The aerodynamically influenced equation of motion can then be written as follows:

Solving Eq. (31) as an eigenvalue problem gives the stability condition of the system. Fig. 9 shows the change in the aerodynamic damping of the torsional mode. The flutter critical wind velocity is defined as the cross point of the torsional aerodynamic logarithmic decrement and the equivalent torsional structural logarithmic decrement (δ = −0.0263). The flutter critical wind velocity found from the FDs coincided with that found using a wind-tunnel dynamic test (Fig. 3). In the case of I_r = 21.02%, the flutter critical wind velocity identified from the FDs (V_cr = 8.2) was slightly smaller than that obtained from the dynamic test (V_cr = 8.6).

This study investigated the effects of turbulence on the FDs of a truss bridge deck section using a wind-tunnel test and the SSI method. The following conclusions were obtained from this study:

• FDs can be successfully obtained from gust responses. The advantages of this method include the gust response being easy to obtain and this method being less time-consuming than the traditional method; in particular, the length of the time history makes it easy to meet the requirements for extracting FDs, even at a high wind speed—a situation that reflects full-scale bridge behaviors more closely than the traditional method.

• SSI_data shows good results even from gust response because this method holds the advantage of considering buffeting force and response as inputs instead of as noise.

• Turbulent flow significantly affects the dynamic response and FDs of the truss bridge deck section. $A_{2}^{*}$ became positive at V_r = 5.2 in the smooth flow and was delayed until V_r = 6.5-7.8 in the turbulent flows of I_r = 6.97% and I_r = 11.09%, respectively; in the case of I_r = 21.02%, a positive value for $A_{2}^{*}$ did not appear.

Hoang Trong Lam, Hiroshi Katsuchi, and Hitoshi Yamada declare that they have no conflict of interest or financial conflicts to disclose.