Numerical predictions of wind-induced buffeting vibration for structures by a developed pseudo-excitation method

Abstract

A numerical analysis method for wind-induced response of structures is presented which is based on the pseudo-excitation method to significantly reduce the computational complexity while preserving accuracy. Original pseudo-excitation method was developed suitable for adoption by combining an effective computational fluid dynamic method which can be used to replace wind tunnel tests when finding important aerodynamic parameters. Two problems investigated are gust responses of a composite wing and buffeting vibration responses of the Tsing Ma Bridge. Atmospheric turbulence effects are modeled by either k–ω shear stress transport or detached eddy simulation. The power spectral responses and variances of the wing are computed by employing the Dryden atmospheric turbulence spectrum and the computed values of the local stress standard deviation of the Tsing Ma Bridge are compared with experimental values. The simulation results demonstrate that the proposed method can provide highly efficient numerical analysis of two kinds of wind-induced responses of structures and hence has significant benefits for wind-induced vibration engineering.

Introduction

The increasing awareness of the harmful effect of wind-induced vibration on structures has led to extensive research. Such vibration of tall/slender structures may cause structural failure, discomfort to occupants, or malfunction of equipment. Hence determining wind-induced structural dynamic responses is important. It is very common to deal with wind loads as random functions of time and space, e.g. Sivakumar and Haran ¹ presented aircraft random vibration analysis with random runway mathematical profiles.

Traditional structural random vibration analysis methods are widely adopted for response prediction with the two most famous being the complete quadratic combination (CQC) and square root of the sum of squares (SRSS) methods. The former is theoretically accurate but involves heavy computation and the latter is only effective for structures for which all modes have small damping ratios. Also multiexcitation problem has been a challenging research difficulty. ² Lin ³ presented a pseudo-excitation method (PEM) (initially called a fast CQC algorithm) to deal with dynamic responses of structures subjected to random seismic excitation. PEM is an accurate high-speed computation method which can be employed to analyze single- or multi-excitation random seismic responses which can be stationary or nonstationary.^4,5 In this method, the determination of random response of a structure is converted into the determination of response of the structure to a series of harmonic loads, i.e. the so-called pseudo excitations.

Wind tunnel tests were always employed to predict the wind-induced structural responses. Useful aerodynamic parameters (e.g. aerodynamic force coefficients and flutter derivatives) were also important data which could be found through experiments. ⁶ Significant aeroelastic progress was achieved through wind tunnel tests, beginning with Scanlan and Tomko, ⁶ e.g. Ge et al. ⁷ investigated aerodynamic stabilization of central stabilizers for box girder suspension bridges. Over recent decades computational fluid dynamics (CFD) has been extensively used in solving fluid–structure interaction (FSI) problem due to advances in computer power and algorithms.^8,9 CFD aeroelastic analysis of slender structures involves the difficulties of how to account both for FSI and for turbulent flow. Zhao and Li ¹⁰ investigated coupling oscillations due to transient energy growth of flow disturbances with a simplified model based on Dirichlet boundary conditions, which developed advanced air flow analysis based on the usage of novel energy methods.^11,12 Bai et al.^13,14 proposed a numerical method for aerodynamic analysis of airfoil or bluff bodies. The FSI problem was solved by a block-iterative coupling method and different turbulence models were tested. Two-dimensional or three-dimensional (3D) CFD modeling was adopted to find aerodynamic parameters efficiently.

Evaluating gust buffeting response has been a main research topic in structure vibration for several decades. At the beginning, simplified algorithm was proposed. Solari ¹⁵ proposed a gust effect factor technique to determine the maximum displacements and internal forces of vertical structures subjected to gust-excited vibrations and investigated the limit value of this technique. Based on simplified algorithms and actual statistical excitations, gust response predictions continued to be implemented; Chaya et al. ¹⁶ investigated the coupled gust response in nonstationary wind.

An accurate and efficient structural dynamic analysis method has significant potential while investigating the gust-excited vibration of structures, because it can reduce the dependence on wind tunnel tests. The present work mainly focuses on proposing a valuable numerical method based on PEM and CFD to gust response analysis of a composite wing and to buffeting vibration response analysis of long-span bridges. The power spectral responses and variances of the wing are computed with a kind of atmospheric turbulence spectrum, whereas the computed values of the local stress standard deviation of the Tsing Ma Bridge are compared to the experimental values. During the computation process, turbulence effects are included through $k-\omega$ shear stress transport (SST) or detached eddy simulation (DES) turbulence models.

PEM for wind-induced vibration analysis

For complicated structures, wind-induced response computations may require the use of several modes of vibration, rather than just the first mode which suffices for most common structures. Use of the classical random-vibration-based mode superposition approach can be difficult in such cases, whereas PEM converts such problems to the determination of the response of the structure to a series of harmonic loads, i.e. the so-called pseudo excitations.

Nonstationary random processes are generally of short duration and their basic characteristic is that the statistical properties vary significantly with time. The basic principle of PEM for nonstationary random vibration analyses is illustrated by Figure 1⁵ as follows.

Figure 1.

Basic principle of pseudo-excitation method (nonstationary analysis).

Consider a structure subjected to an evolutionary random excitation

$$f\left(t\right)=g\left(t\right)x\left(t\right)$$

(1)

in which g(t) is a slowly varying modulation function, while x(t) is a zero-mean stationary random process with auto power spectrum density (auto-PSD) $S_{\mathit{xx}}(\omega )$ . The deterministic excitation functions g(t) and $S_{\mathit{xx}}(\omega )$ are both assumed to be given. The auto-PSD of f(t) is

$$S_{\mathit{ff}}\left(\omega ,t\right)=g^2\left(t\right)S_{\mathit{xx}}\left(\omega \right)$$

(2)

In order to compute the PSD functions of various linear responses caused by f(t), a form of the pseudo (denoted by ∼ below) excitation can be given as

$$\tilde{f}\left(\omega ,t\right)=g\left(t\right)\sqrt{\left[S_{\mathit{xx}}\left(\omega \right)\right]}\exp \left(i\omega t\right)$$

(3)

Now suppose that y(t) is an arbitrarily selected response vector (see Figure 1(b)) and that $\tilde{\mathbf{y}}(\omega ,t)$ is a corresponding transient response vector due to the pseudo excitation $\tilde{f}\hspace{0.17em}(\omega ,t)$ with the structure initially at rest (see Figure 1(b)). It has been proven that the desired PSD matrix of y(t) is

$${\mathbf{S}}_{\mathit{yy}}\left(\omega ,t\right)={\tilde{\mathbf{y}}}^{*}\left(\omega ,t\right){\tilde{\mathbf{y}}}^T\left(\omega ,t\right)$$

(4)

where the superscripts * and T denote matrix conjugate and transpose, respectively.

It has been found that when CQC, SRSS, and PEM are used to calculate power spectral responses, N-dimensional vectors are multiplied, respectively, $q^2$ times, q times, and once. ⁴ Moreover, especially for long-span structures, space modal coupling requires the number of participant modes to be of order 10², which hugely increases the computation time. ¹⁷ Of the three methods, PEM is theoretically accurate and computationally efficient and so can overcome the disadvantages of both the CQC and SRSS methods. Similar conclusions were drawn for nonstationary random excitation problems. ¹⁸ Therefore, PEM is chosen to be developed in this paper for numerical predictions of wind-induced responses for different structures.

Numerical algorithm for aerodynamic parameter computation

A 2D flat plate in incompressible flow is taken as an example (see Figure 2). Here U is the flow velocity, $\rho$ is the flow density, a is a dimensionless coefficient, and ab represents the distance between the rigid center and the chord midpoint of the section, so that if the rigid center coincides with the chord midpoint, a = 0 and b is the half length of the chord. Considering forced vibration of the form

$$h=h_0\exp \left(i\omega t\right),\alpha ={\alpha }_0\exp \left(i\omega t\right)$$

(5)

where h is the vertical structure motion and $h_0$ is the initial vertical forced vibration amplitude, $\alpha$ is the twist structure motion and ${\alpha }_0$ is the initial twist forced vibration amplitude, and $\omega =2\pi f$ with f being the forced vibration frequency.

Figure 2.

Schematic of the 2D thin plate.

For the same 2D flat plate of Figure 2, Scanlan and Tomko ⁶ proposed the expansions

$$L=\frac{\rho U^2}{2}\left(2b\right)\left[k{H}_1^{*}\frac{\dot{h}}{U}+k{H}_2^{*}\frac{b\dot{\alpha }}{U}+k^2{H}_3^{*}\alpha \right]$$

(6)

$$M=\frac{\rho U^2}{2}{\left(2b\right)}^2\left[k{A}_1^{*}\frac{\dot{h}}{U}+k{A}_2^{*}\frac{b\dot{\alpha }}{U}+k^2{A}_3^{*}\alpha \right]$$

(7)

where L is the aerodynamic lift force and M is the aerodynamic moment; $k=\frac{b\omega }{U}$ is the dimensionless reduced frequency of the motion; $H_1^{*},H_2^{*},H_3^{*},A_1^{*},A_2^{*}$ , and $A_3^{*}$ are the flutter derivatives which in general are functions of k; $\dot{h}$ is the time derivative of h; and $\dot{\alpha }$ is the time derivative of $\alpha$ .

To get the six flutter derivatives, a CFD method based on block-iterative coupling and turbulence modeling is used in this paper. ¹³ Wind–structure interaction effects can be included during the prediction process by using this method. Solvers based on this method were implemented by existing CFD and structural analysis codes, i.e. a CFD code CFX and a FEM code JIFEX (i.e. it was developed by Dalian University of Technology). ¹⁹

Gust analysis of a composite wing

The structure

Composite materials have many advantages in industrial aerospace applications, e.g. high strength, high stiffness, long fatigue life, and good damping characteristics. So the laminated NACA0012 (i.e. NACA represents National Advisory Committee for Aeronautics) wing structure consisting of composite materials is chosen in this paper (see Figure 3). The total span thickness (i.e. perpendicular to the section shown in Figure 3) of the wing is 5 m, the thickness of each layer is 0.02 m, and the wing is fixed at its root. The material constants of the laminate are longitudinal elastic modulus $E_1=138$  GPa, transverse elastic modulus $E_2=8.96$  GPa, positive axis shear modulus $G_{12}=7.1$  GPa, positive axis Poisson ratio ${\nu }_{12}=0.30$ , and $\rho =2000{\text{kg/m}}^3$ .

Figure 3.

The classical NACA0012 airfoil model.

Analysis method by developing PEM

When flying in a vertical simple harmonic gust wind, the wing vibrates harmonically. So equations (6) and (7) can be rewritten as

$$L\left(t\right)=\frac{1}{2}\rho U^{2}B\frac{w_B}{U}2\pi \left\{C\left({\omega }^{*}\right)\left[J_0\left({\omega }^{*}\right)-i{J}_1\left({\omega }^{*}\right)\right]\right\}+i{J}_1\left({\omega }^{*}\right)\exp \left(i\omega t\right)$$

(8)

$$M\left(t\right)=\left[\frac{x_a}{2}+\frac{1}{4}\right]L\left(t\right)B$$

(9)

where $C({\omega }^{\ast })$ (so that as $\mathit{Ma}\to \infty ,$ $C\left({\omega }^{\ast }\right)$  = 0) is an unsteady aerodynamic coefficient generated by gust wind, B is the chord length of the wing, $x_a$ is the distance between the rigid center and the midpoint of the chord, and $w_B$ is the vertical velocity component of the gust wind.

Considering the case in which the rigid center and the midpoint of the chord coincide, the wing has two degrees of freedom, i.e. vertical translation and rotation about its center. The equation of motion of the wing after discretization can be written as

$$\mathbf{M}\ddot{\mathbf{y}}+\mathbf{C}\dot{\mathbf{y}}+\mathbf{Ky}={\mathbf{f}}_{\mathrm{ae}}+{\mathbf{f}}_g$$

(10)

where ${\mathbf{f}}_{\mathrm{ae}}$ and ${\mathbf{f}}_g$ are the aerodynamic forces due to, respectively, motion of the wing and the harmonic gust wind. Reducing the orders of the equation by using the first q modes of ${\Phi }_q$ (i.e. $(q<<n)$ ) in the mode superposition method ²⁰ gives

$${\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{M}{\Phi }_q\ddot{\mathbf{y}}+{\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{C}{\Phi }_q\dot{\mathbf{y}}+{\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{K}{\Phi }_q\mathbf{y}={\Phi }_{\mathit{q}}^{\mathit{T}}\left({\mathbf{f}}_{\mathit{ae}}+{\mathbf{f}}_g\right)\Phi$$

(11)

$$\begin{array}{l}{\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{M}{\Phi }_q={\mathbf{I}}_q\\ {\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{C}{\Phi }_q=2\Theta \\ {\Phi }_{\mathit{q}}^{\mathit{T}}\mathbf{K}{\Phi }_q={\Omega }^2\end{array}$$

(12)

where ${\mathbf{I}}_q$ is the unit matrix of order q; $\Theta$ is a diagonal matrix with ${\xi }_j{\omega }_j$ as its jth diagonal element; ${\xi }_j$ and ${\omega }_j$ are the damping ratio and circular frequency corresponding to the jth mode, respectively; ${\Omega }^2$ is a diagonal matrix with ${\omega }_j^2$ as its jth diagonal element; and the superscript T represents transposition. Then equation (10) can be rewritten as

$${\mathbf{M}}^{*}\ddot{\mathbf{u}}+\left({\mathbf{C}}^{*}+{\mathbf{C}}_{\mathrm{ae}}^{*}\right)\dot{\mathbf{u}}+\left({\mathbf{K}}^{*}+{\mathbf{K}}_{\mathrm{ae}}^{*}\right)\mathbf{u}={\mathbf{f}}_g^{*}$$

(13)

where ${\mathbf{C}}_{\mathrm{ae}}^{*}$ and ${\mathbf{K}}_{\mathrm{ae}}^{*}$ are, respectively, the aerodynamic damping and stiffness matrices after reducing the order of the aerodynamic force due to the simple harmonic vibration of the structure and are both functions of ω. ${\mathbf{f}}_g^{*}=\mathbf{p}(\omega )\}{w}_B$ can be obtained by reducing the order of equations (8) and (9). $w_B$ also represents the vertical velocity component of atmosphere turbulence, so PEM can be adopted to predict the PSD for atmosphere turbulence response to give

$$\tilde{\mathbf{f}}\left(\omega ,t\right)=\mathbf{p}\left(\omega \right)\sqrt{S_{\mathit{xx}}\left(\omega \right)}\exp \left(i\omega t\right)$$

(14)

where $\sqrt{S_{\mathit{xx}}(\omega )}$ is the velocity spectrum of the vertical atmospheric turbulence.

Rewriting the pseudo excitation $\tilde{\mathbf{f}}$ as $\tilde{\mathbf{f}}={\mathbf{p}}_1\exp (i\omega t)+i{\mathbf{p}}_2\exp (i\omega t)$ and substituting $\mathbf{u}=\mathbf{a}+\mathit{i}\mathbf{b}$ into equation (9) yields

$$\{\begin{array}{l}\mathbf{E}\mathbf{a}+\mathbf{Db}={\mathbf{p}}_1\\ -\mathbf{Da}+\mathbf{Eb}={\mathbf{p}}_2\end{array}$$

(15)

where

$$\begin{array}{l}\mathbf{E}=\left({\mathbf{K}}^{*}+{\mathbf{K}}_{\mathit{ae}}^{*}\right)-{\omega }^2\mathbf{M}\\ \mathbf{D}=-\omega \left({\mathbf{C}}^{*}+{\mathbf{C}}_{\mathit{ae}}^{*}\right)\end{array}$$

Then, equation (12) can be solved to give

$$\mathbf{u}={\left(\mathbf{E}{\mathbf{D}}^{-1}\mathbf{E}+\mathbf{D}\right)}^{-1}\left\{\left(\mathbf{E}{\mathbf{D}}^{-1}{\mathbf{p}}_1-{\mathbf{p}}_2\right)+i\left({\mathbf{p}}_1+\mathbf{E}{\mathbf{D}}^{-1}{\mathbf{p}}_2\right)\right\}$$

(16)

or

$$\mathbf{u}={\left(\mathbf{D}{\mathbf{E}}^{-1}\mathbf{D}+\mathbf{E}\right)}^{-1}\left\{\left({\mathbf{p}}_1-\mathbf{D}{\mathbf{E}}^{-1}{\mathbf{p}}_2\right)+i\left({\mathbf{p}}_2+\mathbf{D}{\mathbf{E}}^{-1}{\mathbf{p}}_1\right)\right\}$$

(17)

With the initial conditions $\mathbf{u}(0)=\dot{\mathbf{u}}(0)=\mathbf{0}$ , all the discrete numerical solutions of each dynamic response can be obtained by solving equation (13). Then the discrete numerical solutions of PSD follow from equation (4).

With the presented numerical method, aerodynamic forces due to the harmonic vertical and twist vibration can be converted into the aerodynamic damping and stiffness matrices of the dynamic equations, and the unsteady aerodynamic forces due to the harmonic gust wind can be taken as the external excitation to be applied. Then the random atmosphere turbulence excitation can be replaced by deterministic harmonic vertical gust wind excitation based on PEM, which enables the response PSD and variance of the wing to be found.

Atmosphere turbulence spectrum and corresponding numerical modeling

The Dryden and Von Karman models are the two main representations of the atmosphere turbulence spectrum ²⁰ and the former is adopted in this paper because of its simplicity of expression and its convenient computation. Its vertical and transverse spectrum functions are, respectively

$${\phi }_{\mathit{uu}}\left(\Omega \right)={\sigma }_u^2\frac{R}{\pi }\frac{1}{1+R^2{\Omega }^2}$$

(18)

$${\phi }_{\mathit{vv}}\left(\Omega \right)={\phi }_{\mathit{ww}}\left(\Omega \right)={\sigma }_w^2\frac{R}{\pi }\frac{1+3R^2{\Omega }^2}{{\left(1+R^2{\Omega }^2\right)}^2}$$

(19)

where $R={\displaystyle {\int }_0^{\infty }}{e}^{-\xi /L}\text{d}\xi$ is the turbulence scale where $\xi$ is a relevant function based on measurement and statistical data for the atmosphere, and ${\sigma }_u$ and ${\sigma }_w$ are the mean square deviations of the turbulence velocity. The asymptotic property of the Dryden model is

$$\underset{\Omega \to 0}{\lim }\Phi (\Omega )=\mathit{const},\underset{\Omega \to \infty }{\lim }\Phi (\Omega )\propto {\Omega }^{-2}$$

(20)

Unsteady Reynolds-averaged Navier–Stokes models are popular for CFD simulations. In this paper, the k–ω SST model is employed for CFD simulations, because it allows direct integration through the boundary layer and also benchmark testing shows that it is particularly superior for wall layer simulation. ²¹

Numerical wing results comparison between CFD simulation and Theodorsen theory

With a mesh generation method proposed by Bai et al., ¹³ 2D CFD modeling of the airfoil of the wing is given in Figure 4. The mesh region was divided into three parts: rigid, wake, and buffer region. It was cylindrical and centered on the airfoil and has radius $R_2=16B$ . The meshes used for the rigid, wake, and remaining buffer region were, respectively, (392 × 113), (64 × 160), and (104 × 39) (i.e. total 58,592 cells), where the first and second numbers are the number of cells in the tangential and radial directions, respectively. The viscous boundary layer over the structure surface is well resolved by the fine mesh with the overall y+ less than 2. By this way, uniform mean wind speeds have been imposed with turbulence model k–ω SST used and a pressure condition is applied at the outlet. Mesh independence test is shown in Table 6 in Appendix 1. Then flutter derivatives which were computed first for gust response analysis can be obtained.

Figure 4.

Mesh generation of the airfoil of the proposed wing.

Table 1 gives the first 10 eigenvalues of the wing structure introduced in the beginning of this section. For the computation, the first-order natural vibration frequency was included, the wind velocity U = 10 m/s, the Dryden atmospheric turbulence spectrum was used, the frequency integration range was $\omega \in [4,32]$ (rad/s), and the interval of the adjacent frequency points was $\Delta \omega =0.1$ (rad/s). After finding the aerodynamic forces of each frequency point by both CFD simulation and Theodorsen theory, the PSD of vertical displacement for the wing was obtained. Figure 5 presents the PSD of vertical displacement response of the leading edge of the wingtip, from which it can be seen that the response increases significantly near the natural vibration frequency of the wing. Additionally, the computational values given by CFD simulations are less than those from Theodorsen theory, which happens because the CFD simulations include actual physical factors which can induce such a decrease of the structure response, e.g. aerodynamic viscosity and the thickness of the wing. Hence the presented CFD simulation is a valuable method for predicting realistic wind-induced vibration behavior of structures.

Figure 5.

PSD of the vertical displacements of the leading edge of the wingtip. The solid and dashed lines represent the computed results based on CFD and Theodorsen theory, respectively. PSD: power spectrum density.

Table 1.

Natural vibration characteristic of the wing.

Mode	Frequency (rad/s)	Period (s)
1	25.66	0.24472
2	158.53	0.03961
3	220.83	0.02844
4	241.74	0.02598
5	438.80	0.01431
6	733.78	0.00856
7	846.73	0.00742
8	1188.60	0.00528
9	1249.11	0.00503
10	1370.25	0.00458

Results and discussion for the wing

This section presents gust wind response predictions for the wing based on the proposed method. The first 10 modes were still used for computation and the modal damping ratio of the wing was 0.02, the flow velocity $U_{\infty }<<0.3\text{Ma}$ , and incompressible flow was assumed. Unsteady aerodynamic force was applied at the quarter point of the chord. The Dryden atmospheric turbulence spectrum was still used, with R=750 m and ${\sigma }_w=1\text{m/s}$ . The alternative horizontal wind velocities of 10, 20, 30, and 40m/s were used. The frequency integration range used was $\omega \in [0,300]$ (rad/s) with the interval of the adjacent frequency points being $\Delta \omega =0.1$ (rad/s). The aerodynamic parameters were still obtained by the CFD method.

The computed vertical displacement variances at four locations on the wingtip are shown in Table 2, and Figure 6 shows the related PSD variation across the wingtip for wind velocity 40 m/s. Figure 7 shows how the PSD of the vertical displacement responses varies with velocity for the leading edge (left-hand picture) and midpoint of the wingtip. It is readily seen that the values of PSD fluctuate near the natural vibration frequencies of the structure. The PSD of the response is very small in the high frequency range, due to the velocity spectrum of the gust wind attenuating quickly at high frequencies.

Figure 6.

Response of the vertical displacement PSD at different positions on the wingtip when U = 40 m/s. PSD: power spectrum density.

Figure 7.

Variation with velocity of the vertical displacement PSD at the leading edge and at the center of the wingtip. PSD: power spectrum density.

Table 2.

Vertical displacement variance at the positions on the wingtip shown.

Velocity (m/s)	Leading edge	Quarter length from the leading edge	Centre	Trailing edge
10	1.025E-04	1.044E-04	1.036E-04	1.020E-04
20	2.137E-03	2.122E-03	2.106E-03	2.072E-03
25	5.849E-03	5.806E-03	5.762E-03	5.673E-03
30	1.347E-02	1.337E-02	1.327E-02	1.307E-02
35	2.737E-02	2.717E-02	2.697E-02	2.656E-02
40	5.050E-02	5.015E-02	4.978E-02	4.902E-02

Table 3 compares the computation times needed when using PEM, CQC, and SRSS. Clearly PEM dramatically improves the computational efficiency, even though it retains accuracy.

Table 3.

Efficiency comparison.

Method	CQC	SRSS	PEM
Time	86 min	55 min	40 s

CQC: complete quadratic combination; PEM: pseudo-excitation method; SRSS: square root of the sum of squares.

Buffeting analysis of Tsing Ma Bridge

The structure with FEM information

The Tsing Ma Bridge in Hong Kong is a suspension bridge with a main span of 1377 m between the Tsing Yi tower in the east and the Ma Wan tower in the west (see Figure 8).

Figure 8.

Tsing Ma Bridge in Hong Kong.

The height of the two reinforced concrete towers is 206 m and the two main cables are of 1.1 m diameter, 36 m apart, and supported by saddles located at the top of the four tower legs. The bridge deck is a hybrid steel structure consisting of a series of cross-frames supported on two outer-longitudinal trusses and two inter-longitudinal trusses acting compositely with stiffened steel plates, as shown in Figure 9. The bridge deck is supported by suspenders on the main span and on the Ma Wan side span and by three piers on the Tsing Yi side span. It carries a dual three-lane highway on its upper level and two railway tracks and two carriageways within its lower level. For structural health monitoring, a WASHMS was installed in the bridge by the Highways Department of Hong Kong in 1997. There are nine types of sensor in the WASHMS, including 6 anemometers and 110 strain gauges.

Figure 9.

A typical 18 m deck section of the main span of the Tsing Ma Bridge.

The finite element model of the Tsing Ma Bridge was generated based on the following principles: the geometry of the finite element model was consistent with the actual structure, there was a one-to-one correspondence between the structural elements of the model and those of the actual structure, the equivalent treatment was adopted for the simulations of the structural stiffness and density, and the boundary conditions were realistic. The global finite element model indicated in Figure 10 includes 12,898 nodes, 21,946 elements (2906 plate elements and 19,040 beam elements) with 4788 multipoint connection elements.

Figure 10.

Three-dimensional refined finite element model of the Tsing Ma Bridge.

Buffeting analysis method

The 3D integral equations of motion of the Tsing Ma Bridge can be represented as

$${\mathbf{M}}^{\text{s}}\ddot{\mathbf{x}}\left(t\right)+{\mathbf{C}}^{\text{s}}\dot{x}\left(t\right)+{\mathbf{K}}^{\text{s}}x\left(t\right)={\mathbf{f}}^{\text{bf}}+{\mathbf{f}}^{\text{se}}$$

(21)

where $\mathbf{x}(t)={\left\{{\mathrm{x}}_1^{\text{T}}(t),{\mathrm{x}}_2^{\text{T}}(t),\mathrm{\dots },{\mathrm{x}}_N^{\text{T}}(t)\right\}}^{\text{T}}$ is the nodal displacement vector of the structure finite element model; ${\mathbf{x}}_j(t)$ (j = 1,2,…,N) is the jth displacement column vector comprised of three translation displacements and three rotational displacements; N is the total number of nodes; ${\mathbf{M}}^{\text{s}}$ , ${\mathbf{C}}^{\text{s}}$ , and ${\mathbf{K}}^{\text{s}}$ are the 6N × 6N structural mass, damping, and stiffness matrices, respectively; and ${\mathbf{f}}^{\text{bf}}$ and ${\mathbf{f}}^{\text{se}}$ are the 6N order nodal buffeting and self-excited force vectors, respectively.

In this paper, the power spectral method is adopted to solve the motion equations in the frequency domain and PEM is used with a highly detailed model in order to get solutions accurately and efficiently.

The nodal displacement of the bridge was represented by the generalized displacement vector as

$$\mathbf{x}\left(t\right)=\mathbf{\Phi }\left(t\right)\mathbf{q}\left(t\right)$$

(22)

where the generalized displacement vector $\mathbf{q}(t)={\left\{{\mathrm{q}}_1(t),{\mathrm{q}}_2(t),\mathrm{\dots },{\mathrm{q}}_{N_m}(t)\right\}}^{\text{T}}$ , N_m is the total number of modes used, and the modal matrix $\mathbf{\Phi }\left(t\right)=\left[{\mathbf{\Phi }}_1\left(t\right),{\mathbf{\Phi }}_2\left(t\right),\mathrm{\dots },{\mathbf{\Phi }}_{N_m}\left(t\right)\right]$ has 6N×N_m dimension. Then equation (21) can be written as

$$\overline{\mathbf{M}}\ddot{\mathbf{q}}+\overline{\mathbf{C}}\dot{\mathbf{q}}+\overline{\mathbf{K}}\mathbf{q}={\overline{\mathbf{f}}}^{\text{bf}}+{\overline{\mathbf{f}}}^{\text{se}}$$

(23)

where $\overline{\mathbf{M}}$ , $\overline{\mathbf{C}}$ , and $\overline{\mathbf{K}}$ are N_m×N_m generalized mass, damping, and stiffness matrices, respectively; $\overline{\mathbf{M}}=\hspace{0.17em}{\mathbf{\Phi }}^{\text{T}}{\mathbf{M}}^{\text{s}}\mathbf{\Phi }$ , $\overline{\mathbf{C}}=\hspace{0.17em}{\mathbf{\Phi }}^{\text{T}}{\mathbf{C}}^{\text{s}}\mathbf{\Phi }$ , and $\overline{\mathbf{K}}=\hspace{0.17em}{\mathbf{\Phi }}^{\text{T}}{\mathbf{K}}^{\text{s}}\mathbf{\Phi }$ ; and ${\overline{\mathbf{f}}}^{\text{bf}}={\mathbf{\Phi }}^{\text{T}}{\mathbf{f}}^{\text{bf}}$ and ${\overline{\mathbf{f}}}^{\text{se}}={\mathbf{\Phi }}^{\text{T}}{\mathbf{f}}^{\text{se}}$ are the generalized buffeting and self-excited forces, respectively.

The power spectrum matrix of the motion equations due to random excitation is a nonnegative definite Hermitian matrix, so it can be decomposed with complex ${\mathbf{L}}^{\mathbf{*}}\mathbf{D}{\mathbf{L}}^{\text{T}}$ , which is similar to real $\mathbf{LD}{\mathbf{L}}^{\text{T}}$ decomposition. The elements of the diagonal matrix $\mathbf{D}$ are nonnegative.

Now using the corresponding column vectors of the diagonal matrix $\mathbf{D}$ and the lower tridiagonal matrix $\mathbf{L}$ gives a series of simple harmonic pseudo excitations. Then pseudo responses were determined by applying the pseudo excitations to the motion equation. Assuming $r_j$ is the pseudo displacement or internal force response induced by the $j$ th pseudo excitation, the final response spectrum matrix of the structure is

$${\mathbf{S}}_{\mathbf{rr}}\left(\omega \right)={\displaystyle \sum _{j=1}^{m_p}{\mathbf{r}}_j^{*}\left(\omega \right){\mathbf{r}}_j^{\text{T}}\left(\omega \right)}$$

(24)

Self-excited forces arise from wind–structure interaction, so the computation of the nodal self-excited force must consider the displacement relationship of each main span section node. Finally, equation (23) becomes

$$\overline{\mathbf{M}}\ddot{\mathbf{q}}+\left(\overline{\mathbf{C}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{D}\mathbf{\Phi }\right)\dot{\mathbf{q}}+\left(\overline{\mathbf{K}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{S}\mathbf{\Phi }\right)\mathbf{q}={\overline{\mathbf{f}}}^{\text{bf}}$$

(25)

where ${\overline{\mathbf{f}}}^{\text{bf}}={\mathbf{\Phi }}^{\text{T}}{\mathbf{f}}^{\text{bf}}$ is the generalized buffeting force; the generalized mass matrix $\overline{\mathbf{M}}=\hspace{0.17em}\mathbf{I}$ ; the generalized stiffness matrix $\overline{\mathbf{K}}=\hspace{0.17em}\text{diag[}{\omega }_{\text{1}}^{\text{2}},{\omega }_2^2,\mathrm{\dots },{\omega }_{N_m}^2]$ , where ${\omega }_i$ ( $i=1,2,\mathrm{\dots },N_m$ ) is the $i$ th natural frequency of the structure; and the generalized damping matrix $\overline{\mathbf{C}}=\hspace{0.17em}\text{diag[}2{\varsigma }_1{\omega }_1,2{\varsigma }_2{\omega }_2,\mathrm{\dots },2{\varsigma }_{N_m}{\omega }_{N_m}]$ where ${\varsigma }_i$ is the $i$ th modal damping ratio.

Assuming that the time-history of the fluctuating wind velocity is a stationary random process and based on random vibration theory, the spectrum matrix of the generalized nodal displacement $q(t)$ of the buffeting response was represented as

$${\mathbf{S}}_{\mathit{qq}}\left(\omega \right)={\overline{\mathbf{H}}}^{*}\left(\omega \right){\mathbf{S}}_{\overline{F}\overline{F}}\left(\omega \right){\overline{\mathbf{H}}}^{\text{T}}\left(\omega \right)$$

(26)

where ${\mathbf{S}}_{\overline{F}\overline{F}}(\omega )$ is the N_m×N_m spectrum matrix of the generalized buffeting force and $\overline{\mathbf{H}}(\omega )$ is the generalized frequency response matrix. Expressions for them are, respectively

$${\mathbf{S}}_{\overline{F}\overline{F}}\left(\omega \right)={\mathbf{\Phi }}^{\text{T}}{\mathbf{Q}}^{\text{bf}}{\hspace{0.17em}}^{*}{\mathbf{S}}_{\mathit{aa}}\left(\omega \right){\mathbf{Q}}^{\text{bf}}{\hspace{0.17em}}^{\text{T}}\mathbf{\Phi }$$

(27)

$$\overline{\mathbf{H}}\left(\omega \right)=\frac{1}{\left(\overline{\mathbf{K}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{S}\mathbf{\Phi }\right)-{\omega }^2\overline{\mathbf{M}}+i\omega \left(\overline{\mathbf{C}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{D}\mathbf{\Phi }\right)}$$

(28)

where the power spectrum ${\mathbf{S}}_{\mathit{aa}}(\omega )$ of the turbulence fluctuating wind is a nonnegative definite Hermitian matrix to which ${\mathbf{L}}^{\mathbf{*}}\mathbf{D}{\mathbf{L}}^{\text{T}}$ decomposition was applied. Hence

$${\mathbf{S}}_{\mathit{aa}}\left(\omega \right)={\mathbf{l}}^{*}\left(\omega \right)\mathbf{d}\left(\omega \right){\mathbf{l}}^{\text{T}}\left(\omega \right)={\displaystyle \sum _{j=1}^{m_p}{\mathbf{S}}_{\mathit{aa},j}\left(\omega \right)}$$

(29)

$${\mathbf{S}}_{\mathit{aa},j}\left(\omega \right)=d_j{\mathbf{l}}_j^{*}\left(\omega \right){\mathbf{l}}_j^{\text{T}}\left(\omega \right)$$

(30)

where $m_p$ is the rank of the matrix ${\mathbf{S}}_{\mathit{aa}}(\omega )$ ; $\mathbf{d}(\omega )$ is a 3M × 3M real diagonal matrix with $d_j$ as its $j$ th nonzero diagonal element (i.e. $j=1,2,\mathrm{\dots },m_p$ ); and $\mathbf{l}(\omega )$ is a 3M × 3M lower triangular matrix with unity as its diagonal elements, so that ${\mathbf{l}}_j(\omega )$ is the column vector consisting of the lower triangular elements belonging to the $j$ th column.

For PEM, the deterministic harmonic pseudo excitation applied to the global bridge which corresponds to the sub-spectrum ${\mathbf{S}}_{\mathit{aa},j}$ was represented as

$${\mathbf{a}}_{p,j}\left(\omega ,t\right)={\mathbf{l}}_j\left(\omega \right)\sqrt{d_j\left(\omega \right)}\exp \left(i\omega t\right)$$

(31)

Under the action of the $j$ th pseudo excitation ${\mathbf{a}}_{p,j}(\omega ,t)$ , the generalized bridge motion equation (25) is rewritten as

$$\overline{\mathbf{M}}{\ddot{\mathbf{\eta }}}_{p,j}\left(\omega ,t\right)+\left(\overline{\mathbf{C}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{D}\mathbf{\Phi }\right){\dot{\mathbf{\eta }}}_{p,j}\left(\omega ,t\right)+\left(\overline{\mathbf{K}}-{\mathbf{\Phi }}^{\text{T}}\mathbf{S}\mathbf{\Phi }\right){\mathbf{\eta }}_{p,j}\left(\omega ,t\right)={\mathbf{\Phi }}^{\text{T}}{\mathbf{Q}}^{\text{bf}}{\mathbf{a}}_{p,j}\left(\omega ,t\right)$$

(32)

where the generalized pseudo displacement response is

$${\mathbf{\eta }}_{p,j}\left(\omega ,t\right)=\overline{\mathbf{H}}\left(\omega \right){\overline{\mathbf{l}}}_j\left(\omega \right)\sqrt{d_j\left(\omega \right)}\exp \left(i\omega t\right)$$

(33)

$${\overline{\mathbf{l}}}_j\left(\omega \right)=\hspace{0.17em}{\mathbf{\Phi }}^{\text{T}}{\mathbf{Q}}^{\text{bf}}{\mathbf{l}}_j\left(\omega \right)$$

(34)

and the pseudo bridge displacement response to the jth pseudo excitation ${\mathbf{a}}_{p,j}(\omega ,t)$ is

$${\mathbf{x}}_{p,j}\left(\omega ,t\right)=\mathbf{\Phi }\overline{\mathbf{H}}\left(\omega \right){\overline{\mathbf{l}}}_j\left(\omega \right)\sqrt{d_j\left(\omega \right)}\exp \left(i\omega t\right)=\sqrt{d_j\left(\omega \right)}{\overline{\mathbf{x}}}_j\left(\omega \right)\exp \left(i\omega t\right)$$

(35)

The pseudo velocity and acceleration response are, respectively

$${\dot{\mathbf{x}}}_{p,j}\left(\omega ,t\right)=i\omega \sqrt{d_j\left(\omega \right)}{\overline{\mathbf{x}}}_j\exp \left(i\omega t\right)$$

(36)

$${\ddot{\mathbf{x}}}_{p,j}\left(\omega ,t\right)=-{\omega }^2\sqrt{d_j\left(\omega \right)}{\overline{\mathbf{x}}}_j\exp \left(i\omega t\right)$$

(37)

The pseudo element stress determined from the element modal stress was represented as

$${\mathbf{\sigma }}_{p,j}=\mathbf{\Gamma }\overline{\mathbf{H}}\left(\omega \right){\overline{\mathbf{l}}}_j\left(\omega \right)\sqrt{d_j\left(\omega \right)}\exp \left(i\omega t\right)$$

(38)

where $\Gamma$ represents the modal stress matrix. $\Phi$ and $\Gamma$ were obtained from the modal analysis of the structure. Now the structural displacement and stress response spectrum matrix were, respectively, determined by equation (24)

$${\mathbf{S}}_{\mathbf{xx}}\left(\omega \right)={\displaystyle \sum _{j=1}^{m_p}{\mathbf{x}}_{p,j}^{*}{\mathbf{x}}_{p,j}^{\text{T}}}={\displaystyle \sum _{j=1}^{m_p}{d}_j{\overline{\mathbf{x}}}_j^{*}{\overline{\mathbf{x}}}_j^{\text{T}}}$$

(39)

$${\mathbf{S}}_{\mathrm{\sigma \sigma }}\left(\omega \right)={\displaystyle \sum _{j=1}^{m_p}{\mathbf{\sigma }}_{p,j}^{*}{\mathbf{\sigma }}_{p,j}^{\text{T}}}$$

(40)

Finally, the standard deviation of the bridge response was found by integrating the response spectrum matrix in the frequency domain

$${\left(\text{RMS}\right)}_{\mathbf{a}}=\sqrt{2{\displaystyle {\int }_0^{+\infty }{S}_{{\mathbf{a}}_i}\left(\omega \right)\text{d}\omega }}$$

(41)

where a represents either nodal displacement x or stress response $\mathbf{\sigma }$ .

Results and discussions for the Tsing Ma Bridge

First, aerodynamic parameters (i.e. flutter derivative) of the Tsing Ma Bridge were obtained by both the 3D CFD method and wind tunnel experiments. ²² The turbulence of the wind was modeled with DES, which has been proved to be effective for simulating 3D turbulent flow. ¹⁴ Mesh generation of 3D modeling is shown in Figure 11. With the similar algorithm of airfoil, total cell number of the deck section of Tsing Ma Bridge is 3,283,496 (i.e. the thicknesses of the 3D model was 1 m). Mesh independence test is shown in Table 7 in Appendix 1.

Figure 11.

CFD modeling of the deck of Tsing Ma Bridge.

The comparison between the two methods is shown in Figure 12. It can be seen that the computed values of the 3D CFD method are close to those of the wind tunnel experiments.

Figure 12.

Flutter derivatives of the main span of the Tsing Ma Bridge. CFD: computational fluid dynamics.

Figure 13.

Global positions of 10 strain gauges used by WASHMS.

Then a numerical prediction of the buffeting response of the Tsing Ma Bridge was carried out. The first 80 modes of the structure were used with the highest natural frequency being 1.1 Hz and the modal damping ratio being 0.005, so that the computational range was 0–1.1 Hz. The records of the structural health monitoring system WASHMS when the Tsing Ma Bridge went through the typhoon named York in 1999 were taken as the actual stress data. The record values of 10 strain gages were chosen for comparison with the computed values. Their positions are shown in Figure 13. The notation scheme for the strain gages is as follows: SS-TLN-xx and SS-TLS-xx represent a uniaxial strain gage in the north and south longitudinal trusses of the section L, respectively; and SP-TLN-xx and SP-TLS-xx represent a pair of uniaxial strain gages in the north and south longitudinal trusses of the section L, respectively.

Tables 4 and 5 show the stress response values of the computations based on the aerodynamic parameters from the CFD method and from the wind tunnel experiment. The measured values from WASHMS are also shown for comparison in Table 4. It can be seen that for seven stress responses, the difference between the computed values achieved by the presented method and the measured values is less than 30%; the differences of the other three are in the range of 30–40%. All the computational results were less than the measured values, and the computational results can reveal the buffeting response characteristic of the Tsing Ma Bridge. According to the engineering practice, this difference can be brought due to manufacturing difference between the design model and the actual bridge, local model of the actual bridge, etc. This difference can be acceptable for engineering applications if the relative error can be estimated.

Table 4.

Standard deviation comparison of measured and computed local stresses.

Strain gage number	Stress standard deviation (MPa)		Relative error (%)
	Measured stress (WASHMS)	Computed stress (CFD)
SP-TLN-01	1.487	1.321	11.2
SS-TLN-01	1.426	1.339	6.1
SP-TLN-05	1.697	1.383	18.5
SS-TLN-03	1.907	1.367	28.3
SP-TLS-12	1.659	1.432	13.7
SS-TLS-09	1.832	1.212	33.8
SP-TLN-02	1.441	1.155	19.8
SP-TLN-03	1.409	1.089	22.7
SP-TLS-02	1.351	0.887	34.3
SP-TLS-03	1.385	0.837	39.6

CFD: computational fluid dynamics; WASHMS: Wind and Structural Health Monitoring System.

Table 5.

Standard deviation comparison of local stresses computed using flutter derivatives from the wind tunnel test and CFD numerical simulation.

Strain gage number	Computed stress standard deviation (MPa)
	Based on experiment	Based on CFD
SP-TLN-01	1.61	1.321
SS-TLN-01	1.59	1.339
SP-TLN-05	1.875	1.383
SS-TLN-03	1.904	1.367
SP-TLS-12	1.518	1.432
SS-TLS-09	1.777	1.212
SP-TLN-02	1.431	1.155
SP-TLN-03	1.249	1.089
SP-TLS-02	0.904	0.887
SP-TLS-03	1.062	0.837

CFD: computational fluid dynamics.

From Table 5, it can also be seen that the difference between the computed stress response values based on the aerodynamic parameters achieved by wind tunnel experiment and the measured values from WASHMS is similar to that between the computed values achieved based on the presented method and the measured values. Hence the usage of the presented method which was developed based on PEM and CFD can certify the efficiency of the complex structural buffeting analysis of the Tsing Ma Bridge.

Conclusions

This paper presents numerical predictions of wind-induced responses for different structures by proposing a developed numerical method based on PEM and CFD. PEM converts the nonuniform random response to dynamic response analysis under deterministic external loads, and CFD computes the aerodynamic parameters accurately. A gust response analysis method for the composite wing and a buffeting response analysis method for the Tsing Ma Bridge are proposed.

For the initial investigation of a composite wing, the computational values based on CFD simulations were less than those based on Theodorsen theory. This occurs because the presented CFD simulations include actual physical phenomena, e.g. aerodynamic viscosity and the thickness of the wing, which can reduce the structure response. Then random response analysis of the composite wing has shown that the changes of the PSD of the wing due to different velocities are similar and exhibit increases near the natural frequencies of the wing. The efficiency was hugely improved without losing accuracy, because the use of PEM and CFD does not cause any approximations or simplifications. Hence the method has significant application value in the field of numerical analysis of wing structures.

Numerical method of buffeting analysis of Tsing Ma Bridge was also developed based on PEM and 3D CFD. The stress response of the Tsing Ma Bridge in Hong Kong was also computed. It was found that the difference between the computed values achieved by the presented method and the measured values from WASHMS can be sufficiently acceptable for engineering applications. Hence the 3D CFD numerical method can be used to replace wind tunnel experiments in order to find aerodynamic parameters effectively. Furthermore, using PEM with the 3D CFD method provides accurate and efficient dynamic response analysis of complex bridge structures.

The simulation results presented demonstrate that the proposed method can provide highly efficient numerical analysis of wind-induced response for structures and so has significant benefits in wind engineering applications.