Numerical investigation on the nonlinear flutter behaviors of a π-shaped sectional bridge deck

Abstract

Nonlinear flutter behaviors of a π-shaped sectional bridge deck are investigated numerically, focusing on the effects on the post-flutter limit cycle oscillation (LCO) amplitudes caused by various nonlinearities, including aeroelastic, structural damping, and geometric nonlinearity. The Multi-stage indicial functions (IFs) are adapted to describe the nonlinear aeroelastic effects in time-domain, which necessitate amplitude-dependent aeroelastic model parameters. The motion amplitude-dependent parameters are mainly extracted numerically with computational fluid dynamics (CFD) method and partially examined by experimentally extracted results. The CFD-based parameters turn out to be satisfactory compared with experimental results. Finally, the flutter thresholds and post-flutter limit cycle oscillations are simulated using the present time-domain multi-stage indicial functions model. The numerical results show that introduction of any one kind of nonlinearity (aeroelastic, structural damping, or geometric) in the numerical procedure can lead to LCO. These three kinds of nonlinearity are compared in detail in terms of the effects they have on the post-flutter LCO amplitude. It turns out that the damping nonlinearity, for the case concerned in this work, is the most important factor that determines the post-flutter LCO amplitude and the flutter onset wind speed as well. The time-domain multi-stage indicial functions procedure adopted in this work can serve as an efficient solving strategy to obtain nonlinear post-flutter performances of long-span bridges.

Keywords

Nonlinear aeroelasticity post-flutter limit cycle oscillation indicial function π-shaped bridge section

Introduction

Flutter thresholds have long been used to evaluate the aeroelastic performances of bridge structures. The aerodynamic response of a bridge structures duo to wind is assumed to diverge until structural failure for any oncoming wind speeds beyond a specific threshold. The critical wind speed dominated methodology has been extensively practiced, both numerically and experimentally, due to presumed catastrophic destiny of a flutter instability. Researches in the aerospace field indicate, however, that flutter in the form of a slender and thin airfoils usually end up with limit cycle oscillation (LCO) of small-to-moderate motion amplitudes as a result of structural nonlinearity, aeroelastic nonlinearity, or of a combination of both (Dowell 2015, Günter and Mai, 2019). It is a known fact that the original Tacoma Narrows Bridge experienced flutter instability at a wind speed of about 19 m/s, of which the deck rotation reached as large as about 70° and finally resulted in collapse of the whole span. However, the same mechanism experienced by the Golden Gate Bridge resulted in only slight structural motion without damage. (Morse-Fortier, 2005). Therefore, different deck configurations differ substantially in post-flutter properties.

Three approaches in general can be adopted to obtain linear or nonlinear flutter properties. The first one is sectional model testing (Amandolese et al., 2013; Tang and Hua, 2019), main features of which include the following: (1) With small motion amplitudes and nonlinearities of the model system being able to be neglected, the critical flutter wind speed can be accurately obtained. (2) When the motion of the model is large, the post-flutter characteristics of the sectional model are almost irrelevant to those of the full bridge. However, segment models can still be used to obtain motion-dependent nonlinear aeroelasticities. The second one is computational fluid dynamics (CFD) coupled with structural motion equations (Ying et al., 2017). This methodology does not use any semi-empirical or empirical model. Its reliability depends on the selected turbulence model and grid quality, and the numerical solving process may be very time-consuming. The third one involves a semi-empirical nonlinear model based on experimental results. Larsen et al. (2007) brought forward a semi-empirical model specifically used to describe the nonlinear aerodynamic force of wind turbine blades; Diana G et al. (2010) proposed a nonlinear aerodynamic model with transient wind angle of attack (AOA) and its first derivative as variables, which is able to describe the hysteresis effect of aerodynamic force at a given amplitude. Although the family of semi-empirical models deal with nonlinear hysteresis loops, they are not able to reflect the amplitude-dependent properties, which is the key to the nonlinear post-flutter issues.

In recent years, there are several studies have focused on the aerodynamic performance of bridge decks at large amplitudes. Wu T and Kareem A (2014a, 2015) proposed a method by using the Volterra theory to describe nonlinear bridge aerodynamics. Liu and Ge (2015) adopted the method of additional nonlinear differential equations and additional aerodynamic degrees-of-freedom to simulate the memory effect and nonlinear characteristics of aerodynamic self-excitation force. Time-domain self-excited force can be realized by indicial function (Wu T and Kareem A, 2014b). Therefore, Zhang (2018) put forward the multi-stage indicial function method, which describes amplitude-dependent aeroelasticity and hence is applicable to post-flutter issues. To reflect the nonlinear aeroelastic properties of bluff bridge decks. Gao and Zhu (2015) used third-order polynomials to describe the nonlinear characteristics of aerodynamic moments, and results obtained by this method are inconsistent with the experimental results. Zhang et al. (2019) tried a derivative-based nonlinear self-excited force model and verified with CFD results. Recently, machine learning technique is used to describe nonlinear aerodynamics (Li et al., 2020; Wu and Kareem, 2011), it is an important method that has been rapidly advanced.

The self-excited force model mentioned above is basically a time-frequency hybrid model, which may cause some false modes when applied to the full bridge calculation. At present, self-excited force time-domain model is rare, and it can consider various nonlinearity in the full bridge flutter analysis, such as geometric and material nonlinearity. Therefore, it is worthy to build the pure time-domain model to obtain the post-flutter characteristics of bridge deck.

In this paper, the LCO amplitude characteristics of a π-shaped deck are investigated by wind tunnel tests. A nonlinear time-domain aeroelastic model based on multi-stage IFs combined with amplitude-dependent structural damping is proposed. Its accuracy and reliability in calculating the flutter threshold and post-flutter motions are verified. The effects of a number of nonlinear factors on the post-flutter response are studied.

Time-domain aeroelastic model

Multi-stage IFs

Higher-order aerodynamic components may exist for blunt body at large motion amplitude (Lin et al., 2019; Zhang et al., 2019). However, High-order force components have few effects on flutter response of bridge deck (Zhang et al.2019, 2021). Therefore, aerodynamic loads acting on the blunt body section of a bridge can be expressed by the Scanlan time-frequency hybrid model. The aerodynamic self-excited lift force and moment per unit length due to harmonic vibration are expressed as follows (Scanlan, 1993, 2000)

L_{se} = \frac{1}{2} ρ U^{2} (2 B) [{KH}_{1}^{*} \frac{h}{U} + {KH}_{2}^{*} \frac{B α}{U} + K^{2} H_{3}^{*} α + K^{2} H_{4}^{*} \frac{h}{B}]

(1)

M_{se} = \frac{1}{2} ρ U^{2} (2 B^{2}) [{KA}_{1}^{*} \frac{h}{U} + {KA}_{2}^{*} \frac{B α}{U} + K^{2} A_{3}^{*} α + K^{2} A_{4}^{*} \frac{h}{B}]

(2)

where ρ is fluid density; U is oncoming wind speed; B is the deck width; K = Bω/U is the reduced frequency; ω is circular frequency; $H_{i}^{*}$ (i = 1∼4), $A_{i}^{*}$ (i = 1∼4) are flutter derivatives, h and α are the plunges and torsional displacements, respectively; $h^{.}$ , $α^{.}$ are first derivatives of h and α with respect to time t, respectively.

The effects of along-wind motion on aeroelasticity are neglected. Equation (1) and (2) are time- and frequency-domain mixed expressions, which are not able to be used directly in time-domain. IF can be used to transform the mixed-domain expressions into time-domain. Taking the aerodynamic lift as an example, a general form as be expressed as

L (s) = \frac{1}{2} ρ U^{2} (2 B) C'_{L} α_{0} φ (s)

(3)

where $ϕ (s)$ is the indicial function; $C_{L}^{'} = d C_{L} / d C_{L} d α d α$ ; s =Ut/B is a dimensionless time.

The aerodynamic lift and torque due to an arbitrary combination of vertical and torsional motions can be obtained as follows

\begin{matrix} L_{se} (x, s) = L_{se α} (x, s) + L_{seh} (x, s) \\ = \frac{1}{2} ρ U^{2} (2 B) C'_{L} [\int_{- \infty}^{s} φ_{L α} (s - σ) α' (x, σ) d σ + \int_{- \infty}^{s} φ_{Lh} (s - σ) \frac{h " (x, σ)}{B} d σ] \end{matrix}

(4)

\begin{matrix} M_{se} (x, s) = M_{se α} (x, s) + M_{seh} (x, s) \\ = \frac{1}{2} ρ U^{2} (2 B^{2}) C'_{M} [\int_{- \infty}^{s} φ_{M α} (s - σ) α' (x, σ) d σ + \int_{- \infty}^{s} φ_{Mh} (s - σ) \frac{h " (x, σ)}{B} d σ] \end{matrix}

(5)

where $α'$ is the derivative of $α$ with respect to dimensionless time s, $h "$ is the second derivative of h concerning s, and x is the coordinate along the bridge deck. The motions themselves can vary with x coordinate along the bridge deck because of the deck’s rotations.

The IFs can be extracted directly from experiments or obtained indirectly from aerodynamic derivatives. According to the principle of spectral equivalence, the IFs are related to the aerodynamic derivative as follows

K (H_{1}^{*} - {iH}_{4}^{*}) = C'_{L} [φ_{Lh} (0) + \bar{φ}'_{Lh}]

(6)

K^{2} (H_{3}^{*} + {iH}_{2}^{*}) = C'_{L} [φ_{L α} (0) + \bar{φ}'_{L α}]

(7)

K (A_{1}^{*} - {iA}_{4}^{*}) = C'_{M} [φ_{Mh} (0) + \bar{φ}'_{Mh}]

(8)

K^{2} (A_{3}^{*} + {iA}_{2}^{*}) = C'_{M} [φ_{M α} (0) + \bar{φ}'_{M α}]

(9)

In this way, the interrelationships between the parameters of the IFs and the flutter derivatives are established.

The single-stage IF based method is derived on the basis of linear theory, where the aerodynamic forces are independent of the motion. If the motion amplitude-dependent properties are to be taken into account, however, multi-stage IFs need to be introduced as follows (Zhang, 2018)

ϕ_{FA} (s, A_{k}) = 1 - \sum_{i = 1}^{n} a_{ik} e^{- d_{ik} s}

(10)

The nonlinear aeroelasticity is determined by the groups of IFs φ_FA, describing the nonlinear aeroelasticity at a group of discrete motion states A_k (k =1, 2, ……, n).

When a bridge deck is experiencing flutter instability, its motion amplitude would vary with time until an LCO is formed. Corresponding to any intermediate state with torsional amplitude α* and vertical amplitude $h^{*}$ , suppose

α_{k} \leq α^{*} \leq α_{k + 1}, h_{n} \leq h^{*} \leq h_{n + 1}

(11)

then the group of IFs corresponding to this intermediate state is determined by interpolation as

\begin{matrix} ϕ_{L α} (s, α^{*}) = & \frac{α_{k + 1} - α^{*}}{α_{k + 1} - α_{k}} ϕ_{L α} (s, α_{k}) \\ + \frac{α^{*} - α_{k}}{α_{k + 1} - α_{k}} ϕ_{L α} (s, α_{k + 1}) \\ = η_{α, k} \cdot ϕ_{L α} (s, α_{k}) + η_{α, k + 1} \cdot ϕ_{L α} (s, α_{k + 1}) \end{matrix}

(12)

\begin{matrix} ϕ_{M α} (s, α^{*}) = \frac{α_{k + 1} - α^{*}}{α_{k + 1} - α_{k}} ϕ_{M α} (s, α_{k}) \\ + \frac{α^{*} - α_{k}}{α_{k + 1} - α_{k}} ϕ_{M α} (s, α_{k + 1}) \\ = η_{α, k} \cdot ϕ_{M α} (s, α_{k}) + η_{α, k + 1} \cdot ϕ_{M α} (s, α_{k + 1}) \end{matrix}

(13)

\begin{matrix} ϕ_{Lh} (s, h^{*}) = \frac{h_{n + 1} - h^{*}}{h_{n + 1} - h_{n}} ϕ_{Lh} (s, h_{n}) + \frac{h^{*} - h_{n}}{h_{n + 1} - h_{n}} ϕ_{Lh} (s, h_{n + 1}) \\ = η_{h, n} \cdot ϕ_{Lh} (s, h_{n}) + η_{h, n + 1} \cdot ϕ_{Lh} (s, h_{n + 1}) \end{matrix}

(14)

\begin{matrix} ϕ_{Mh} (s, h^{*}) (s, h^{*}) = \frac{h_{n + 1} - h^{*}}{h_{n + 1} - h_{n}} ϕ_{Mh} (s, h_{n}) \\ + \frac{h^{*} - h_{n}}{h_{n + 1} - h_{n}} ϕ_{Mh} (s, h_{n + 1}) \\ = η_{h, n} \cdot ϕ_{Mh} (s, h_{n}) + η_{h, n + 1} \cdot ϕ_{Mh} (s, h_{n + 1}) \end{matrix}

(15)

and corresponding to this intermediate motion state, the self-excited aerodynamic lift and torque are obtained by equation (4) and (5) as

{\begin{matrix} L_{se} (s, α^{*}, h^{*}) \\ M_{se} (s, α^{*}, h^{*}) \end{matrix}} = [\begin{matrix} L_{se α} (s, α_{k}) & L_{se α} (s, α_{k + 1}) & L_{seh} (s, h_{h}) & L_{seh} (s, h_{n + 1}) \\ M_{se α} (s, α_{k}) & M_{se α} (s, α_{k + 1}) & M_{seh} (s, h_{h}) & M_{seh} (s, h_{n + 1}) \end{matrix}] {\begin{matrix} η_{α, k} \\ η_{α, k + 1} \\ η_{h, n} \\ η_{h, n + 1} \end{matrix}}

(16)

where

\begin{matrix} η_{a, k} = \frac{α_{k + 1} - α^{*}}{α_{k + 1} - α_{k}}, η_{a, k + 1} = \frac{α^{*} - α_{k}}{α_{k + 1} - α_{k}}, \\ η_{h, n} = \frac{h_{n + 1} - h^{*}}{h_{n + 1} - h_{n}}, η_{h, n + 1} = \frac{h^{*} - h_{n}}{h_{n + 1} - h_{n}} . \end{matrix}

(17)

Mean wind loads

Mean wind loads are functions of the wind angle of attack. When the deck rotation is small, they can be taken as

\bar{D} (x) \approx \frac{1}{2} ρ U^{2} (2 B) {C_{D} (α_{0}) + {\frac{\partial C_{D}}{\partial α} |}_{α_{0}} \cdot \bar{α} (x)}

(18)

\bar{L} (x) \approx \frac{1}{2} ρ U^{2} (2 B) {C_{L} (α_{0}) + {\frac{\partial C_{L}}{\partial α} |}_{α_{0}} \cdot \bar{α} (x)}

(19)

\bar{M} (x) \approx \frac{1}{2} ρ U^{2} (2 B^{2}) {C_{M} (α_{0}) + {\frac{\partial C_{M}}{\partial α} |}_{α_{0}} \cdot \bar{α} (x)}

(20)

where C_L, C_D, and C_M are the aerostatic coefficients; $\partial C_{D} / \partial α$ , $\partial C_{L} / \partial α$ , and $\partial C_{M} / \partial α$ are the derivatives of C_L, C_D, and C_M with respect to the wind attack angle α, respectively.

Time-domain process

In order to eliminate transient step response due to applied gravity and mean wind loads, an enlarged structural damping ratio is adopted within a number of initial time steps, as shown in Figure 1. The initial time duration with enlarged damping should be long enough for the structure reaching its equilibrium position, and after that, the intended damping ratio is restored and the self-excited loads are applied to the structure.

Figure 1.

Loading diagram.

Damping nonlinearity can be achieved in the time-domain solution process. Assuming that vertical and torsional damping is a function of motion amplitude; that is, $ξ_{α} (x_{i}, s) = f (α^{*} (x_{i}, s))$ , $ξ_{h} (x_{i}, s) = f (h^{*} (x_{i}, s))$ , which may lead to some errors in the process of flutter convergence because the amplitude in one motion period is used to calculate the damping of the next. However, it should have no effect on final results because the LCO amplitude would be stable in the last. The time-domain solution process is presented in Figure 2.

Figure 2.

Time-domain analysis process.

Sectional model tests

Experimental set-up

Sectional model tests are carried out in the HD-2 boundary layer wind tunnel at Hunan University. It is a closed circuit tunnel of 17.0 m (length)×3.0 m (width)×2.5 m (height) in dimension. The wind speed in the wind tunnel is adjustable between 0∼58 m/s. The turbulence intensity measured under uniform flow in the wind tunnel is less than 0.5% when the wind speed is higher than 4 m/s.

The standard section of a Π-shaped sectional model is shown in Figure 3. The length of the model is L=1.54 m and the width and depth are B=0.58 m and D=0.067 m, respectively. The model is elastically mounted in wind tunnel. Vertical stiffness (K_h) is provided by eight identical springs (with a stiffness of 3663.7 N/m). The torsional stiffness (K_α) is provided by a fixed distance x_SC between the model center and the spring plane. Elliptical plates with chamfers are installed at both ends of the model to eliminate end effects.

Figure 3.

The Π-shaped sectional bridge model (units: m).

Three laser displacement sensors (optoNCDT Model 2300) with an operating range of 100∼300 mm are used to measure vertical motions. The model rotation is calculated from vertical motions and the corresponding sensor-to-sensor distance. The coordinate system of the suspension system is shown in Figure 4, where the vertical displacement is positive upward and the torsional displacement is positive clockwise. The sampling frequency is 200 Hz. Figure 5(a) is the design drawing and Figure 5(b) shows the elastically mounted model in wind tunnel.

Figure 4.

Two-dimensional flexibly mounted bridge model.

Figure 5.

Sectional model: (a) the design drawing; (b) testing in wind tunnel.

Free decay measurements of the model (α₀ being +3, 0, and −3) without wind are carried out for pure vertical and torsional vibrations. The vertical and torsional natural frequencies (f_h and f_α) of the model are obtained through spectral analysis. The structural damping ratio is calculated as

ξ = \frac{\ln \frac{v_{n}}{v_{n + 1}}}{2 π}

(21)

where v_n and v_n₊₁ are vertical or torsional amplitudes in two adjacent periods.

The damping characteristics of the system are of obvious nonlinear nature (Gao et al. 2018; Tang and Hua, 2019). The variation of damping ratio with respect to motion amplitude (α₀=+3°) is identified using local least-squares fitting (LLSF) as shown in Figure 6. For relevant LLSF detail of nonlinear damping, one can refer to the literature by Tang and Hua, 2019). Figure 6(b) shows that the torsional damping ratio versus the amplitude relationship is approximately linear from 0.5° to 0.95°; therefore, a linear fit is requisite, which can be used to take into account damping nonlinearity in time-domain solution process.

Figure 6.

(a) Free decay vibration for α₀=+3°; (b) Nonlinear damping ratio with respect to amplitude for α₀=+3°.

The final damping ratios (ξ_h and ξ_α) are obtained by the mean of 20 consecutive adjacent samples. ξ_h and ξ_α at different initial AOAs are about 0.7% for linear damping. The mass and moment of inertia (m and I_α) of the whole motion system are ca. d as follows

m = K_{h} / {({2 π f}_{h})}^{2}

(22)

I_{α} = \frac{K_{α}}{{(2 π f α)}^{2}} .

(23)

The main parameters of the suspension system are summarized in Table 1. The moment of inertia is calculated with regard to shear center (S.C) of the model section.

Table 1.

Main parameters of the suspension system (S.I. units).

m	I_α	K_h	K_α	f_h	f_α	x _SC
42.924	1.208	29,309.473	890.519	4.160	4.321	0.174

Test results

Based on the linear theory, the vibration of the sectional model would diverge at any oncoming wind speed over the flutter threshold (Ucr). However, the resulted vibration tends to be stable LCOs in the experiments. Figure 7 shows that the amplitudes of the model are LCOs rather than divergence and the LCO amplitude increases with the reduced wind velocity (U/f_αB).

Figure 7.

Response of flutter: (a) vertical; (b) torsional.

Numerical simulation of nonlinear flutter

FE model and CFD-based aerodynamic coefficients

The experimental model is simulated numerically to verify the feasibility of the post-flutter time-domain solution strategy. The ANSYS finite element software is used to simulate the model, where BEAM4 element is adopted to simulate the sectional model, LINK10 element is used to simulate the spring support, and MASS21 element is used to simulate the mass and mass inertia moment, as shown in Figure 8.

Figure 8.

Finite element configuration of sectional model test.

Rayleigh damping model is adopted. The damping matrix is expressed as a linear combination of the mass and stiffness matrices as

C = α M + β K

(24)

where C, M, and K are damping, mass, and stiffness matrices, and α and β are Rayleigh damping coefficients, respectively.

It has been shown that the elastically suspended model in the wind tunnel has obvious nonlinear stiffness characteristics (refer to Amandolese et al., 2013; Gao et al. 2018; Tang and Hua, 2019). Therefore, it is necessary to consider geometric nonlinearity in time-domain analysis. The main calculation parameters for linear damping are listed in Table 2

Table 2.

Dynamic calculation parameters for linear damping.

AOA (deg)	Time steps (s)	Integral method	Damping ratio	A	β
+3	0.01	Newmark-β	0.007	0.1864397	0.0002627
0
−3

The flutter derivatives are identified by using the single-degree-of-freedom forced vibration method. The forced vertical and torsional vibration modes are given as follows

h (t) = h_{0} \sin (2 π f_{h} t)

(25)

α (t) = α_{0} \sin (2 π f_{α} t)

(26)

where h₀ and α₀ are the vertical and torsional amplitudes, respectively. f_h and f_α are corresponding oscillating frequencies, respectively, and $t$ is time.

The f_h and f_α are set at 2 Hz. The vertical amplitude includes two schemes, h₀/B=1/116 and h₀/B=1/29, respectively. The torsional amplitude includes eight schemes, namely, 1°, 2°, 4°, 6°, 8°, 10°, 12°, and 14°. The wind velocity U varies from 2 to 12 m/s and U/fB varies from 1.7241 to 10.3448.

Figure 9 shows the size of the body section used in the numerical simulation, which is simplified to some extent compared with the experimental model. The cross-section is 0.58 m in width (B) and 0.067 m in height (D). The three-dimensional span-wise effect has not been considered in the CFD procedure. The model is subjected to single-degree-of-freedom vertical or torsional vibration, the same as the experimental set-up.

Figure 9.

Dimension of the CFD model.

The computational domain is shown in Figure 10. The left side is the velocity inlet boundary, and the right side is the pressure outlet boundary. Both upper and lower boundaries are defined as symmetric boundaries. The cross-section surface is set as a non-slip boundary. The upstream length of the computational domain is 12B and the downstream length is 30B. The height of the computational domain is 24B, which ensures a sufficiently small blocking rate ζ =D/24B≈0.00,481.

Figure 10.

Computational domain and boundary conditions for α₀=+3°.

The SST k-ω model, which has good accuracy predicting test data and good robustness, is used for turbulence modeling. The turbulence of the inlet flow is set to 0.2% and the turbulent viscosity ratio is set to 5. The CFD calculation uses the commercial software ANSYS Fluent 15.0. The stationary regions are discrete using a structured mesh. The rigid and deformed regions are separated by hybrid mesh. The rigid mesh and main beam section together perform the rigid motion and the deformation of the mesh in the deformation region is realized by any Lagrangian dynamic mesh technique based on re-meshed. The second-order implicit strategy and the second-order upwind strategy, respectively, solve the time and space terms of the Reynolds average N–S equation. The SIMPLEC algorithm is used for decoupling of the velocity and pressure equations. The convergence accuracy of each physical quantity calculation is 10e-8, the time step at each wind speed satisfies Ut/B<0.01, and the maximum number of iterations in each time step is 40.

The final amount of meshes for numerical simulation is 915,540 after the grid and time-step independence are verified. The mesh details around the section are shown in Figure 11. The dimensionless height of the first layer grid, Y⁺, under the maximum and minimum Reynolds numbers is less than 1.

Figure 11.

Mesh near the deck section.

Figure 12 shows the flutter derivatives identified by CFD, which are fitted by the least-squares on the basis of equation (1) and (2). Coefficients $A_{2}^{*}$ , $A_{4}^{*}$ , $H_{1}^{*}$ and $H_{3}^{*}$ are in good agreement with those from experiments. Coefficients $A_{3}^{*}$ and $H_{2}^{*}$ under the reduced wind velocity of less than 5.16 (U/fB<5.16) are close to the tests while $A_{1}^{*}$ and $H_{4}^{*}$ deviate from the tested results. The IFs parameters, which are shown in Table 3, are identified based on equation (6) to (9).

Figure 12.

CFD- and EXP-based flutter derivatives.

Table 3.

CFD-based and EXP-based aeroelastic coefficients.

Type	AOA	IFs	Amp	a ₁	a ₂	a ₃	d ₁	d ₂	d ₃
—	+3°	φ_Lh φ_Mh	h₀/B = 1/29	−17.9620	15.5941	4.0319	0.1219	0.1000	0.1623
EXP-based aeroelastic coefficients		φ_Lh φ_Mh	h₀/B = 1/29	82.3208	−81.9532	−4.5081	0.1000	0.1065	2.1000
		φ_Lα φ_Mα	2°	−11.8813	4.4503	9.3374	0.1487	0.2184	0.1000
		φ_Lα φ_Mα	2°	−4.2306	12.9775	−12.6949	1.2219	0.1261	0.1744
	0°	φ_Lh	h₀/B = 1/29	18.4553	−21.7805	5.1317	0.1000	0.1295	0.2283
		φ_Mh	h₀/B = 1/29	4.1844	3.7025	−9.2533	2.1000	0.2014	0.6601
		φ_Lα	2°	−60.0390	1.4964	1.8791	75.7719	0.0500	1.6021
		φ_Mα	2°	9.6723	−12.3715	3.0760	2.1000	0.7625	0.1575
	−3°	φ_Lh	h₀/B = 1/29	−14.3485	13.1197	1.1300	1.7109	1.0659	0.0500
		φ_Mh	h₀/B = 1/29	−4.7826	2.2934	1.9553	0.4853	0.7141	0.1050
		φ_Lα	2°	15.4772	−16.6395	1.4174	1.2061	1.9791	0.0500
		φ_Mα	2°	−3.5636	3.2864	1.3485	0.5449	2.1000	0.0500
	+3°	φ_Lh	h₀/B = 1/116	15.8912	−14.6180	3.6320	0.1000	0.1719	13.0505
		φ_Lh	h₀/B = 1/29	1.9827	5.1346	−3.9260	8.5918	0.0500	0.3039
		φ_Mh	h₀/B = 1/116	−14.7396	12.3288	−94.2726	0.3833	0.5727	69.5526
		φ_Mh	h₀/B = 1/29	−4.2880	−206.2874	51.2280	0.2256	41.3135	11.9617
CFD-based aeroelastic coefficients		φ_Lα	1°	−9.5870	8.5268	3.9020	1.0973	1.8985	0.0803
			2°	−8.0970	7.5422	3.6095	1.1112	2.0446	0.0500
			4°	4.4420	3.6110	−5.4001	2.1000	0.0500	0.8870
			6°	4.2324	3.5026	−5.2486	1.8305	0.0500	0.8215
			8°	3.4810	−4.1933	3.2105	0.0500	0.6869	1.8940
		φ_Mα	1°	−6.6329	−25.8931	23.6811	2.1000	0.3734	0.4841
			2°	10.3719	−9.7545	−8.3850	0.6177	1.0757	0.2281
			4°	3.3478	−6.2785	−5.0457	0.6517	1.1992	0.1155
			6°	−10.6511	9.5060	−6.4126	0.1000	0.2164	0.7689
			8°	10.4441	−7.8570	−10.0518	0.2455	0.6860	0.1000

Zhang et al., 2017.

Flutter thresholds

The flutter thresholds for α₀=+3° are obtained and plotted in Figure 13. Figure 13(a) shows that generally, U_cr calculated by CFD-based flutter derivatives is higher than that by EXP-based flutter derivatives. Mean wind loads, geometric and aeroelastic nonlinearities have little effect on the flutter threshold for the Π-shaped section concerned in this work. However, damping nonlinearity, in comparison with the geometric one, influence flutter thresholds, lowering the onset wind speed from 9.3 m/s to 7.7 m/s for EXP-based aerodynamic coefficients and from 10.6 m/s to 9.9 m/s for CFD-based aerodynamic coefficients.

Figure 13.

Flutter results for α₀=+3°: (a) thresholds; (b) Relative errors.

The relative errors (Res) are defined to evaluate the accuracy for different cases

R e s = \frac{U_{c r, n u m} - U_{c r, \exp}}{U_{c r, \exp}} \times 100 %

(27)

Figure 13(b) shows that the results calculated by EXP-based aerodynamic coefficients are close to the experimental results and have a small Res. The flutter thresholds at different AOAs is shown in Figure 14, where nonlinear effects are not included and structural damping ratio is set to 0.007, the same as in the experiments. Using the EXP-based aerodynamic coefficients, the numerically obtained flutter thresholds for α₀=−3°, α₀=0°, and α₀=+3°are very close to that by experiments, and the relative errors are as low as −5.50%, −2.39%, and 6.49%, respectively.

Figure 14.

Flutter thresholds: (a) Comparison between tested and simulated; (b) Influence of structural damping.

The structural damping is an important factor affecting flutter threshold (Chen and Kareem 2003; Ying et al., 2017), and its effects on flutter threshold of bridge decks depend on the properties of aerodynamic damping. As the wind speed increases, if the growth rate of aerodynamic damping is large, the flutter threshold would not change significantly as the structural damping increases. If, however, the growth rate of the aerodynamic damping is small, meanwhile, aerodynamic and structural damping are of the same order of magnitude, the flutter threshold would be affected significantly by the structural damping. For the Π-shaped section concerned in this work, the evolution of flutter thresholds with respect to the damping ratio are shown in Figure 14(b). It can be seen the flutter threshold varies almost linearly with the damping ratio for α₀=+3°, and increases more drastically at α₀=0° and α₀=−3° than at α₀=+3°. It is noted that this important property is not considered in the currently Chinese Wind-resistance Design Specification for Highway Bridges, which is based on the ideal plate related aerodynamic theories.

LCO phenomena

With the time-domain nonlinear aeroelastic model presented in this work, the aerodynamic responses for α₀=+3° are obtained taking into account of damping or aeroelastic nonlinearities only (see Figure 15).

Figure 15.

Time histories of post-flutter response for α₀=+3° (CFD-based aerodynamic coefficients; mean wind loads included): (a) linear; (b) aeroelastic nonlinearity only; (c) damping nonlinearity only.

The linear model shows an indefinite increase in the motion amplitude when the oncoming wind speed exceeds U_cr, as shown in Figure 15(a). When aeroelastic nonlinearity is introduced, the flutter instability evolves into a stable LCO instead of a divergent motion, as illustrated in Figure 15(b). It can be found that damping nonlinearity also result in a LCO (see Figure 15(c)). While the damping nonlinearity leads to a much smaller LCO amplitude than that resulted from the aeroelastic nonlinearity.

The governing equations on the basis of equation (2) for self-excited vibration of the torsional motion in smooth flow can be expressed as

I_{α} \overset{\cdot\cdot}{α} + 2 I_{α} ω_{α} (ξ_{struct_α} + ξ_{aero_α}) \overset{\cdot}{α} + (K_{α} - ρ ω^{2} B^{4} A_{3}^{*}) α = ρ U^{2} B^{2} [{KA}_{1}^{*} \frac{\overset{\cdot}{h}}{U} + K^{2} A_{4}^{*} \frac{h}{B}]

(28)

where ξ_{struct_α} is structural damping ratio caused by torsional motion; ξ_{aero_α} is the torsional aerodynamic damping ratio, and it can be calculated by $- ρ ω B^{4} A_{2}^{*} / - ρ ω B^{4} A_{2}^{*} 2 I_{α} ω_{α} 2 I_{α} ω_{α}$ .

It can be seen from equation (28) that vertical motion will give rise to coupled aerodynamic damping ratio (ξ_couple). In the process of flutter amplitude increasing, ξ_{struct_α}, ξ_{aero_α} and ξ_couple change with the amplitude. If the sum of ξ_{struct_α}, ξ_{aero_α} and ξ_couple is close to zero, the stable LCO will form. Figure 16 illustrates the flow diagram of identifying total torsion damping ratio (ξ_{total_α}). At first, extracting maximum value and corresponding time and using smooth function in MATLAB Toolbox smooth amplitude curve with time. Then, using two adjacent amplitude points to calculate the damping ratio corresponding to the midpoint of the two amplitudes. Finally, the ξ_couple can be calculated by the following formula

ξ_{couple} = ξ_{total_α} - ξ_{struct_α} - ξ_{aero_α}

(29)

Figure 16.

Flow diagram for identification of total torsional damping ratio.

The LCO mechanism can be explained from Figure 17. For the case without any nonlinearity illustrated in Figure 17(a), all damping ratios are basically constant and independent of the amplitude during the flutter development, which shows that the flutter instability is linear behavior. When only aeroelastic nonlinearity is considered, ξ_{aero_α}, ξ_{total_α} and ξ_couple vary with the amplitude. It can be seen from Figure 17(b) that ξ_{aero_α} and ξ_couple decrease with amplitude for the LCO amplitude greater than 4°, which can be seen from Figure 18 that A₂^*coefficient decreases as LCO amplitude increases at U/fB=4.619 when LCO amplitude is between 4° and 6°, which will cause the aerodynamic damping to decrease with the increase of amplitude. A new LCO will be formed as long as the ξ_{total_α} is zero. In addition, ξ_couple is reduced more than ξ_{aero_α} during the formation of LCO, which shows that ξ_couple plays a greaVter role in promoting the formation of LCO. For the case that only includes damping nonlinearities, ξ_{aero_α} and ξ_couple are basically independent of the torsional amplitude. As shown in Figure 17(c), the increasing trend of ξ_{total_α} with the amplitude is basically the same as that of ξ_{struct_α}. the ξ_{aero_α} and ξ_{struct_α} are of the same order of magnitude. Therefore, it is reasonable that the flutter threshold is affected by structure damping in the flutter thresholds. As long as the ξ_{struct_α} is increased to make the ξ_{total_α} zero, a new LCO can also be formed.

Figure 17.

The evolution of damping ratio versus amplitude for α₀=+3° (CFD-based aerodynamic coefficients, U=11.1 m/s): (a) linear; (b) aeroelastic nonlinearity only; (c) damping nonlinearity only.

Figure 18.

Torsional motion-related flutter derivatives U/fα_B=4.619, α₀=3°: (a) H₂*, (b) A₂*.

At present, the pure time-domain post-flutter algorithm is still rare, so the influence of geometric nonlinearity on post-flutter has not been paid attention to. The geometric nonlinear effect reduces the stiffness of structure under large amplitude, which further affects the dynamic response characteristics of the structure. On the basis of Figure 15(b), once geometric nonlinearity is introduced, the LCO amplitude decreases from 5.63° to 1.01° (Figure 19). In theory, the geometric nonlinearity effect is weak at small LCO amplitude, and it should not have such a large impact on the LCO amplitude. It may be caused by the following two reasons: (i) The Rayleigh damping model, determined by the stiffness and mass matrix, may get rise to additional nonlinear damping if structural stiffness is changed. (ii) The geometric nonlinearity causes the stiffness to decrease, which make the flutter frequency (f) decrease and U/fB increase. That is, geometric nonlinearity has influence on aerodynamic damping. As illustrated in Figure 19(b), it is stranger that the LCO will be formed when only geometric nonlinearity is considered. There is no doubt that the LCO is caused by the Rayleigh damping model.

Figure 19.

Influence of geometric nonlinearity for post-flutter LCO in time-domain analysis (CFD-based aerodynamic coefficients, U = 11.1 m/s): (a) Aeroelastic and geometric nonlinearity included, ξ_α = ξ_h = 0.007; (b) geometric nonlinearity only, ξ_α = ξ_h = 0.007; (c) geometric nonlinearity only, ξ_α = ξ_h = 0.

In order to eliminate the influence of additional nonlinear damping caused by the Rayleigh damping model, the damping ratio of the structure is artificially set to zero and the stiffness of the structure will not affect numerical damping ratio. It can be seen that no LCO appears at any wind speed at ξ_α=ξ_h=0, which is sufficient to show that Rayleigh Damping model will overestimate structural damping when geometric nonlinearity is considered. Therefore, the post-flutter time-domain analysis should carefully use the combination of geometric nonlinearity and Rayleigh damping model.

It has also been discussed above that geometric nonlinearity has effect on aerodynamic parameters. Since CFD-based coefficients combined with geometric nonlinearity can no produce LCO under zero damping. The EXP-based coefficient can be used to verify whether the geometric nonlinear can form a LCO under zero damping. It is interesting that the LCO is formed when EXP-based coefficients (+3°, 0° and −3°) are adopted. Figure 20 shows that the LCO amplitude caused by geometric nonlinearity is larger than that caused by other nonlinear factors, which is also reasonable because the geometric nonlinearity is obvious only at large LCO amplitude. Figure 21 illustrates the flow diagram for obtaining the transient frequency. Firstly, all maximum values will be identified, then the data of 10 cycles before and after A_n maximum value will be selected. In order to obtain the frequency value with higher resolution, original data will be filled in 0 to make the data volume of 20 cycles be 300,000. Finally Fast Fourier transform will be used to obtain the transient frequency (f_n) corresponding to A_n. Assuming that there are k maximum values, the same operation is performed on the 10th to (k-10)th maximum values in turn, and the transient frequency evolution curve with amplitude can be obtained.

Figure 20.

Influence of geometric nonlinearity for post-flutter LCO in time-domain analysis (EXP-based aerodynamic coefficients, geometric nonlinearity only, ξ_α = ξ_h = 0): (a) α₀ = +3°, U = 7.0 m/s; (b) α₀ = 0°, U = 8.6 m/s; (c) α₀ = −3°, U = 9.0 m/s.

Figure 21.

Flow diagram for identification of transient frequency with respect to amplitude.

It is obvious that flutter transient frequency (f) decreases with respect to torsional amplitude, and the frequencies of the +3°, −3°, and 0° are reduced by 0.46%, 0.47%, and 0.52%, respectively, during the formation of LCO (see Figure 22), which will give rise to the increase of U/fB and change the flutter derivatives shown in Figure 23. Therefore, the LCO mechanism caused by geometric nonlinearity is essentially aeroelastic nonlinearity, which different from the aeroelastic nonlinearity caused by LCO amplitude shown in Figure 18.

Figure 22.

Transient frequency and reduced wind speed with respect to amplitude.

Figure 23.

Flutter derivatives with respect to reduced wind speed.

The post-flutter LCO amplitudes (α₀=+3°), under various combinations of nonlinearities, are obtained and pictured in Figure 24, where all cases have taken into account the mean wind loads. It is noted when the experimentally obtained aeroelasticity is used, the LCO amplitudes resulted from the damping nonlinearity only is in good agreement with experimental results, especially the vertical amplitude (see Figure 24(a)). The computed vertical and torsional LCO amplitudes are lower than the experimental results when the geometric and damping nonlinearities are included (see Figure 24(a) and (b)), which is reasonable because of that Rayleigh Damping model will overestimate structural damping when geometric nonlinearity is considered (see Figure 19).

Figure 24.

Post-flutter LCOs (α₀=+3°; mean wind loads included): (a) vertical amplitude A_{v_amp}/D, EXP-based aeroelasticity; (b) torsional amplitude A_{α_amp}, EXP-based aeroelasticity; (c) vertical amplitude A_{v_amp}/D, CFD-based aeroelasticity; (d) torsional amplitude A_{α_amp}, CFD-based aeroelasticity.

When CFD-obtained aeroelastic properties are used, the LCOs obtained from considering aeroelastic nonlinearity only are obviously larger than that from the tested results (see Figure 24(c) and (d)). Moreover, the torsional LCO amplitude obtained from aeroelastic nonlinearity only increases much more rapidly than other cases (see Figure 24(d)). The vertical LCOs, when both the aeroelastic and damping nonlinearities being taken into account, are close to the tested results while the torsional LCOs deviates to some extent from the experimental result when U/f_αB > 5.20 (see Figure 24(c) and (d)). The responses, with all nonlinearities being considered, are in good tendency with experimental results (see Figure 24(c) and (d)), but its value is lower than the experimental results, which may be caused by the Rayleigh damping model overestimating structural damping in post-flutter time-domain analysis. In general, the numerical results based on EXP- and CFD-based aerodynamic coefficients are quite precise in terms of predicting the experimental results.

The discrepancy between the numerical and the experimental results can be attributed to the following aspects: (i) Identification of flutter derivatives and ensuing indicial functions based on the CFD simulations would give rise to certain amount of error. (ii) The identified structural damping of the sectional model using LLSF must result in unavoidable error, which may have a non-negligible impact on the post-flutter response.

The post-flutter (α₀=+3°) LCOs with respect to the reduced wind velocity (U/f_αB) under different damping ratios are also obtained and presented in Figure 25, by using the present time-domain model. It is obvious that the structural damping ratio affects remarkably the LCO amplitudes.

Figure 25.

LCO calculated from for CFD-based aeroelasticity and different damping ratios (α₀=+3°, geometric and aeroelastic nonlinearities included).

Concluding remarks

The flutter thresholds and post-flutter response of a π-shaped sectional bridge model are investigated using a self-excited time-domain model. The time-domain model can serve as a theoretical basis and solving strategy for the nonlinear flutter analysis. The main conclusions are as follows:

Damping, aeroelastic and geometric nonlinearities are realized in post-flutter time-domain analysis. The numerically obtained flutter thresholds and LCO amplitude of a π-shaped bridge deck are in good agreement with the experimental results, indicating applicability of the presented time-domain model.

Damping, aeroelastic, and geometric nonlinearities can give rise to LCO; the LCO mechanism caused by geometric nonlinearity is essentially aeroelastic nonlinearity. The Rayleigh damping model should be used cautiously in time-domain of post-flutter considering geometric nonlinearity because of it will change with structural stiffness and may overestimate the structural damping.

While the flutter threshold of the model is independent on the aeroelastic and geometric nonlinearities, it is sensitive to the structural damping. Structural damping can also influence significantly the amplitude of LCO.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (Grant Number 51578233 and 51938012).

ORCID iD

Zhixiong Wang

References

Amandolese

Michelin

Choquel

(2013) Low speed flutter and limit cycle oscillations of a two-degree-of-freedom flat plate in a wind tunnel. Journal of Fluids and Structures 43: 244–255.

Chen

Kareem

(2003) Efficacy of tuned mass dampers for bridge flutter control. Journal of Structural Engineering 129(10): 1291–1300.

Diana

Rocchi

Argentini

, et al. (2010) Aerodynamic instability of a bridge deck section model: linear and nonlinear approach to force modeling. Journal of Wind Engineering and Industrial Aerodynamics 98: 363–374.

Dowell

(2015) A Modern Course in Aero-Elasticity. Fifth Revised and Enlarged. Edition. Dordrecht: Springer International Publishers.

Gao

Zhu

(2015) Nonlinearity of mechanical damping and stiffness of a spring-suspended sectional model system for wind tunnel tests. Journal of Sound and Vibration 355: 369–391.

Gao

Zhu

Han

, et al. (2018) Nonlinear post-flutter behavior and self-excited force model of a twin-side-girder bridge deck. Journal of Wind Engineering and Industrial Aerodynamics 177: 227–241.

Günter

Mai

(2019) Experiments on transonic limit-cycle-flutter of a flexible swept wing. Journal of Fluids and Structures 84: 153–170.

Liu

(2015) Nonlinear dynamic subsystem model for large-amplitude motion-induced aerodynamic forces of bridge decks. Journal of Harbin Institute of Technology 47(9): 73–78. (in Chinese).

Lin

Wang

Nikitas

, et al. (2019) Effects of oscillation amplitude on motion-induced forces for 5:1 rectangular cylinders. Journal of Wind Engineering and Industrial Aerodynamics 186: 68–83.

10.

Larsen

Nielsen

SRK

Krenk

(2007) Dynamic stall model for wind turbine airfoils. Journal of Fluids and Structures 23: 959–982.

11.

Morse-Fortier

(2005) Professor Robert H. Scanlan and the Tacoma Narrows Bridge. Structures congress.

12.

Scanlan

(1993) Problematics in formulation of wind force models for bridge decks. Journal of Engineering Mechanics 119: 1353–1375.

13.

Scanlan

R.H

(2000) Motion-related body-force functions in two-dimensional low-speed flow. Journal of Fluids and Structures 14: 49–63.

14.

Tang

Hua

Chen

, et al. (2019) Experimental investigation of flutter characteristics of shallow Π section at post-critical regime. Journal of Fluids and Structures 88: 275–291.

15.

Liu

(2020) Nonlinear unsteady bridge aerodynamics: reduced-order modeling based on deep LSTM networks. Journal of Wind Engineering and Industrial Aerodynamics 198: 104116.

16.

Kareem

(2011) Modeling hysteretic nonlinear behavior of bridge aerodynamics via cellular automata nested neural network. Journal of Wind Engineering and Industrial Aerodynamics 99: 378–388.

17.

Kareem

(2014a) Simulation of nonlinear bridge aerodynamics: a sparse third-order volterra model. Journal of Sound and Vibration 333(1): 178–188.

18.

Kareem

(2014b) Revisiting convolution scheme in bridge aerodynamics: a comparison of step and impulse response functions. Journal of Engineering Mechanics 140(5): 70–75.

19.

Kareem

(2015) A nonlinear analysis framework for bluff-body aerodynamics: a volterra representation of the solution of navier-stokes equations. Journal of Fluids and Structures 54: 479–502.

20.

Ying

Zhang

, et al. (2017) Numerical explorations of the limit cycle flutter characteristics of a bridge deck. Journal of Wind Engineering & Industrial Aerodynamics 170: 226–237.

21.

Zhang

(2018) Multistage indicial functions and postflutter simulation of long-span bridges. Journal of Bridge Engineering 23(4): 04018010.

22.

Zhang

, et al. (2019) Energy budget analysis and engineering modeling of post-flutter limit cycle oscillation of a bridge deck. Journal of Wind Engineering and Industrial Aerodynamics 188: 410–420.

23.

Zhang

Liu

Wang

, et al. (2021) A practical method for predicting post-flutter behavior of a rectangular section. Journal of Wind Engineering and Industrial Aerodynamics 216: 104707.