2.1. The CFD–CSD method

Through the CFD–CSD method, the data sources and the validation of the aeroelastic model were obtained. The flow analysis was performed by an in-house hybrid-unstructured flow solver using a cell-centred finite volume method to solve unsteady Reynolds-averaged Navier–Stokes (URANS) equations. The integral form of two-dimensional compressible URANS equations with the S-A turbulence model can be written for a cell of volume $V$ limited by a surface $S$. The equation can be expressed as

(2.1)\begin{equation} \frac{\partial }{{\partial t}}\int_V {{\boldsymbol{W}}\,{\rm d}V} + \oint_S {{{\boldsymbol{E}}^i}({\boldsymbol{W}},{{\boldsymbol{V}}_{grid}})} \boldsymbol{\cdot} {\boldsymbol{n}}\,{\rm d}S - \oint_S {{{\boldsymbol{E}}^v}({\boldsymbol{W}}) \boldsymbol{\cdot} {\boldsymbol{n}}\,{\rm d}S} = \int_V {{\boldsymbol{H}}\,{\rm d}V}, \end{equation}

where $\boldsymbol {n}$ represents the outer unit normal vector to the boundary $S$; $\boldsymbol {W}$ is a five-component vector of conservative variables, ${\boldsymbol {W}} = {[\rho {\rho u}\ {\rho v}\ {\rho E}\ {\rho \tilde {\nu }}]^{\rm T}}$, $\rho$ being the density; $u,v$ are $x$-wise and $y$-wise components of the velocity vector of the flow; $E$ denotes the specific total energy; $\tilde {\nu }$ denotes the working variable of the S-A turbulence model; ${{\boldsymbol {E}}^i}({\boldsymbol {W}},{{\boldsymbol {V}}_{grid}})$ and ${{\boldsymbol {E}}^v}({\boldsymbol {W}})$ are the inviscid flux and the viscous flux, respectively; and $\boldsymbol {H}$ is the source term.

The spatial discretization and time integration of the turbulence model equation and mean flow equations are carried out in a loosely coupled way. The second-order AUSM$+$ scheme is conducted to evaluate the inviscid flux, with a reconstruction technique based on the least-squares approach. For high-order reconstruction schemes, the reader is referred to Liu et al. (Reference Liu, Zhang, Jiang and Ye2016). The grid velocity ${{\boldsymbol {V}}_{grid}}$ is modified according to the geometric conservation law (Ahn & Kallinderis Reference Ahn and Kallinderis2006) when CFD–CSD simulations are performed. The viscous flux term is discretized by the standard central scheme. In the turbulence model, the transport equation for the working variable $\tilde {\nu }$ is given by the standard form proposed by Spalart & Allmaras (Reference Spalart and Allmaras1992). The convective and source terms are discretized by the first-order AUSM$+$ scheme, and the destruction and diffusion terms by the second-order central scheme in the present study. For unsteady computations, the dual time stepping method, proposed by Jameson (Reference Jameson1991), is used to solve the governing equations. At the sub-iteration, the fourth-stage Runge–Kutta scheme is used with local time stepping and residual smoothing to accelerate the convergence.

A no-slip wall boundary condition is applied to the airfoil surface. The location and velocity of the airfoil are updated at each real time step due to the pitching motion of the airfoil. Moreover, the far-field boundary is assigned with a non-reflective boundary condition based on Riemann invariants to ensure that the disturbance generated by the object does not return to the flow field. For more validation of the URANS method, one can refer to our earlier studies (Zhang et al. Reference Zhang, Gao, Liu, Ye and Jiang2015a; Gao, Zhang & Ye Reference Gao, Zhang and Ye2016b). The aforementioned partial results are presented in a condensed form in Appendix B.1.

The structural motion equation was used to characterize and solve the elastic deformation of the airfoil. The generalized structural motion equation in matrix form is as follows:

(2.2)\begin{equation} \boldsymbol{M}\ddot{\boldsymbol{\xi}} + \boldsymbol{G}\dot{\boldsymbol{\xi}} + \boldsymbol{K}\boldsymbol{\xi} = {\boldsymbol{Q}}, \end{equation}

where $\boldsymbol {\xi }$ is the structural generalized displacement; ${\boldsymbol{M}}$, ${\boldsymbol{G}}$ and ${\boldsymbol{K}}$ are the generalized mass, damping (${\boldsymbol{G}} = 0$ in this paper) and stiffness matrix, respectively; and ${\boldsymbol{Q}}$ is the generalized aerodynamic force provided by the CFD method.

Equation (2.2) can be solved by the fourth-order-accuracy hybrid linear multi-step scheme (Zhang, Jiang & Ye Reference Zhang, Jiang and Ye2007) with a predictor–corrector method. This approach helps to build an efficient, stable and loosely coupled algorithm to perform the CFD-based aeroelastic simulation in the time domain. A previous study (Zhang et al. Reference Zhang, Jiang and Ye2007) also investigated the influence of time steps in aeroelastic simulations, indicating that the numerical solution is stable/convergent as long as there are no fewer than 30 time steps in a period. A moving boundary was involved in the simulation of aeroelasticity. Thus, a scheme of computational grid deformation based on the radial basis function interpolation (Wang et al. Reference Wang, Mian, Ye and Lee2015) was employed to match the grid with the new wall boundary condition.

The CFD–CSD method has been applied to the study of transonic flutter (Zhang et al. Reference Zhang, Gao, Liu, Ye and Jiang2015a), transonic buzz (Gao et al. Reference Gao, Zhang and Ye2016a) and transonic buffeting (Gao et al. Reference Gao, Zhang, Li, Liu, Quan, Ye and Jiang2017b). Lacking a standard case for the 2DOF buffet, the coupled CFD–CSD simulation method is validated by classical flutter cases with two degrees of freedom. In our previous work, flutter boundaries of the BACT model (Zhang et al. Reference Zhang, Gao, Liu, Ye and Jiang2015a), the Isogai wing (Zhang et al. Reference Zhang, Jiang and Ye2007) and the LCO with NACA64010 airfoil (Zhang et al. Reference Zhang, Wang, Ye and Quan2012) agreed well with experimental data and reference results. The aforementioned partial results are presented in a condensed form in Appendix B.3.

2.2. Aeroelastic case and validation

The research aeroelastic case is a 2DOF NACA0012 airfoil in subcritical transonic buffeting (pre-buffet) flow, as depicted in figure 3. The free-stream state of the modelling case is $Ma = 0.7$, $Re = 3 \times {10^6}$, $\alpha = 4.5^\circ$. The computational hybrid-unstructured mesh has 25 361 nodes, and 40 layers of structured viscous grids around the airfoil, as depicted in figure 48 in Appendix B.2. The time step for the CFD is $3 \times {10^{-4}}$ s (physical), and the non-dimensional time step is 0.1. The convergence of the grid and the time step are presented in Appendix B.2. Figure 4 compares the transonic buffet onset boundary for the NACA0012 airfoil ($Ma = 0.7\unicode{x2013}0.8$, $Re = 3 \times {10^6}$). The circles are from the experiment of Doerffer et al. (Reference Doerffer, Hirsch, Dussauge, Babinsky and Barakos2010), and the thick solid line is from the URANS method. It also shows the results from Soda & Voss (Reference Soda and Voss2005), in which boundaries are larger than those of the experimental data. The experiment and URANS calculations are in satisfactory agreement.

Figure 3. Sketch diagram of a 2DOF airfoil system.

Figure 4. Comparison of CFD-calculated buffet boundaries with experimental results at $Ma = 0.7$–0.8 and ${Re} = 3 \times 10^{6}$.

Angle ${\alpha }$ is the equilibrium AOA, which is an important parameter affecting the stability of transonic flow, and the stability of the fluid mode is very important for the SDOF flutter. The buffet onset is $4.8^\circ$ for stationary NACA0012 airfoil at $Ma = 0.7$, $Re = 3 \times {10^6}$. Time histories of the lift coefficient at ${\alpha }=4.5^\circ$ and ${\alpha }=5.5^\circ$ are shown in figure 5. It can be seen that dramatic unsteady loads appear at ${\alpha }=5.5^\circ$. At ${\alpha }=4.5^\circ$, flow oscillates slowly and converges to a stable state after the initial excitation. This indicates that although the flow field is in a subcritical buffeting state, its flow stability margin is not high. Further, the dimensionless pressure ($p$) contour plots at typical moments are presented in figures 6 and 7, for the unsteady flow field at ${\alpha }=5.5^\circ$ and the convergence process at ${\alpha }=4.5^\circ$ flow. Vertical dashed lines indicate the position of shock feet. It is evident that in the supersonic buffeting flow, there are significant changes in the oscillation amplitude of shock waves, as well as in the size and position of the flow separation region. During the convergence process in subcritical buffeting flow, shock wave oscillation amplitudes are small, and the left-hand side of the flow separation region remains close to the shock foot, with minor variations on the right-hand side corresponding to shock oscillations. Figure 8 presents the pressure contour plot of the steady flow field at ${\alpha }=4.5^\circ$, while figure 9 provides the surface pressure coefficient ($C_p =2(p-p_{\infty })/(\rho {V^2})$) curve of the steady flow field at ${\alpha }=4.5^\circ$, compared with the dynamic pressure coefficient curve of the unsteady flow field at ${\alpha }=5.5^\circ$. It can be observed that the shock wave position in the subcritical flutter flow field is further aft, while the oscillation range of the shock wave in the supercritical flutter is entirely ahead of the shock wave position in the subcritical state. Because the impact of the fluid mode on the structural mode is too strong in the supercritical transonic buffeting flow, ${\alpha }$ is selected as $4.5^\circ$ in the present work.

Figure 5. The lift responses at ${\alpha }=4.5^\circ$ and ${\alpha }=5.5^\circ$ at $Ma= 0.7$, and ${Re} = 3 \times 10^{6}$.

Figure 6. Pressure contours and streamlines in one period at ${\alpha }=5.5^\circ$, $Ma= 0.7$ and ${Re} = 3 \times 10^{6}$. Here $T$ is the period of the shock oscillation; $X_{SW}$ is the shock position; and $p$ is the non-dimensional pressure: (a) $t_1= 0T$, (b) $t_2= 1/4T$, (c) $t_3= 1/2T$, (d) $t_4= 3/4T$.

Figure 7. Pressure contours and streamlines in one period at ${\alpha }=4.5^\circ$, $Ma= 0.7$ and ${Re} = 3 \times 10^{6}$. Here $L_{SF}$ is the length of the flow separation region on the wing surface: (a) $t_1= 0T$, (b) $t_2= 1/4T$, (c) $t_3= 1/2T$, (d) $t_4= 3/4T$.

Figure 8. Pressure contours and streamlines of the steady flow field at ${\alpha }=4.5^\circ$.

Figure 9. Pressure coefficient curve of supercritical/subcritical buffet flow field.

Below is the structural set-up for the aeroelastic case. For a 2DOF elastic system, with the dimensionless time $\tau = {\omega _\alpha } \times t$, the above matrices and vectors are evaluated as:

(2.3) \begin{equation} \left.\begin{gathered} {\boldsymbol{M}} = \left[{\begin{array}{@{}cc@{}} 1 & {{x_\alpha }}\\ {{x_\alpha }} & {r_\alpha ^2} \end{array}} \right] ,\quad {\boldsymbol{K}} = \left[ {\begin{array}{@{}cc@{}} {{{({\omega _h}/{\omega _\alpha })}^2}} & 0\\ 0 & {r_\alpha ^2} \end{array}} \right]\\ {\boldsymbol{Q}} = \frac{1}{{\rm \pi} }{V^*}^2\left\{{\begin{array}{@{}c@{}} { - ({C_l}-{C_{l0}})}\\ {2({C_m}-{C_{m0}})} \end{array}} \right\},\quad {\boldsymbol{\xi }} = \left\{ {\begin{array}{@{}c@{}} {h/b}\\ \theta \end{array}} \right\}. \end{gathered}\right\} \end{equation}

In the definition of dimensionless parameters, traditional conventions are followed: aerodynamic parameters such as the Reynolds number and force coefficients are non-dimensionalized using the chord length $c$ (or 2$b$), while structural parameters such as the reduced frequency, mass static margin, mass ratio and elastic axis position are non-dimensionalized using the half-chord length $b$. As shown in (2.3) and figure 3, $Ma$ and $Re={\rho }{U_\infty }c/{\mu }_f$ represent the Mach and Reynolds number of the free stream, where $\rho$, ${U_\infty }$ and ${\mu }_f$ respectively represent the density, the velocity and the viscosity coefficient of the free stream. Parameters $C_l$ and $C_m$ are the coefficients of the lift and pitching moment non-dimensionalized by ${C_f}=2F/(\rho {U_{\infty}^{2}})$, where $F$ represents aerodynamic concentrated force. Parameters $C_{l0}$ and $C_{m0}$ are the mean value of the lift and pitching moment coefficients.

The term $h / b$ represents the airfoil heaving displacement and $\theta$ represents the airfoil pitching angle. Frequency ${k} = \omega b/U_\infty$ is the reduced frequency non-dimensionalized by $b$, natural frequency ${\omega }$ and ${U_\infty }$. Therefore, $k_h$ represents the dimensionless natural heaving (plunging) frequency and $k_{\alpha }$ represents the dimensionless natural pitching frequency. Ratio $\mu = m /({\rm \pi} \rho {b^2})$ is the non-dimensional mass ratio, where $m$ represents the mass of the wing. Radius $r_a$ is the gyration radius of the airfoil around the elastic axis. Velocity ${V^*} = {U_\infty }/({\omega _\alpha } b {\mu ^{1/2}}) = 1/(k_{\alpha } {\mu ^{1/2}})$ is the non-dimensional velocity. Parameter $x_a$ is the airfoil static unbalance, which represents the dimensionless distance between the centre of gravity and the elastic axis. The position of the elastic axis and $x_a$ are non-dimensionalized by $b$. Distance $a$ is the dimensionless distance of the elastic axis after the midpoint of the airfoil, which will comprehensively affect the mutual coupling of each fluid/structure mode. The elastic axis is set at 0.224$c$ from the leading-edge point, which is the common position in the research of transonic buzz.

2.3. Unsteady aerodynamic models

The process of modelling is introduced in Step 1 of figure 2. An identification technique is used to construct the ROM. The time cost of the ROM-based flutter analysis mainly lies in training the model. Because the unsteady loads are computed in a discrete domain, the model is chosen for the two-input–two-output system, as shown in (2.4a,b):

(2.4a,b)\begin{equation} \boldsymbol{u} = \left[{\begin{array}{@{}c@{}} {h/b}\\ \theta \end{array}} \right],\quad \boldsymbol{y_f} = \left[ {\begin{array}{@{}c@{}} {{C_l-C_{l0}}}\\ {{C_m-C_{m0}}} \end{array}} \right]. \end{equation}

The ARX model (Zhang et al. Reference Zhang, Jiang and Ye2007) is used to build the aerodynamic ROM. The ARX model is based on difference equations, and it can be expressed as follows:

(2.5)\begin{equation} {\boldsymbol{y_f}}(k) = \sum_{i = 1}^{na} {{{\boldsymbol{A}}_i}{\boldsymbol{y_f}}(k - i)} + \sum_{i = 0}^{nb - 1} {{{\boldsymbol{B}}_i}{\boldsymbol{u}}(k - i)} , \end{equation}

where $\boldsymbol {y}$ is the system output and $\boldsymbol {u}$ is the input. Orders of the model chosen by the user are $na$ and $nb$. The output of the current step is always determined by the linear combination of the previous $na$ steps’ output and the previous $nb-1$ steps’ input. The delay parameters ($na$, $nb$) determine the prediction accuracy of the model in unsteady flow. Coefficients ${\boldsymbol {A}}_i$ and ${\boldsymbol {B}}_i$ are constant coefficients to be estimated. The least-squares method is used to estimate unknown model parameters. To achieve a zero mean, the constant levels need to be removed from the data before they are estimated. This is mainly for two reasons. One is to eliminate the influence of static aeroelastic deformation on the system as much as possible in numerical studies. The second is to improve the accuracy of the autoregressive architecture model, where the sample data need to meet the zero-input zero-output characteristic.

In order to obtain the state-space aeroelastic model, we define a state vector ${\boldsymbol {x_f}}(k)$ consisting of ($na+nb-1$) vector states as follows:

(2.6)\begin{equation} {\boldsymbol{x_f}}(k) = {[{y_f}(k - 1), \ldots {y_f}(k - na),u(k - 1), \ldots ,u(k - nb + 1)]^{\rm T}} . \end{equation}

The state-space form of the discrete-time aerodynamic model is as follows:

(2.7)\begin{equation} \left.\begin{array}{@{}l@{}} {{\boldsymbol{x_f}}}(k + 1) = {{{\boldsymbol{A}}}_a}{\boldsymbol{x_f}}(k) + {{{\boldsymbol{B}}}_a}{\boldsymbol{u}}(k)\\ {{\boldsymbol{y_f}}}(k) = {{{{\boldsymbol{C}}}}_a}{\boldsymbol{x_f}}(k) + {{{{\boldsymbol{D}}}}_a}{\boldsymbol{u}}(k), \end{array} \right\} \end{equation}

where

(2.8) \begin{align} \left.\begin{gathered} {{{\boldsymbol{A}}}_a} = \left[ {\begin{array}{@{}cccccccccc@{}} {{{\boldsymbol{A}}_1}} & {{{\boldsymbol{A}}_2}} & \cdots & {{{\boldsymbol{A}}_{na - 1}}} & {{{\boldsymbol{A}}_{na}}} & {{{\boldsymbol{B}}_1}} & {{{\boldsymbol{B}}_2}} & \cdots & {{{\boldsymbol{B}}_{nb - 2}}} & {{{\boldsymbol{B}}_{nb - 1}}}\\ {\boldsymbol{I}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ \vdots & {\boldsymbol{I}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{I}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{I}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{I}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}}\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & {\boldsymbol{0}} & \cdots & {\boldsymbol{I}} & {\boldsymbol{0}} \end{array}} \right],\\ {{{\boldsymbol{B}}}_{a}} = \left[{{{{\boldsymbol{B}}}_{\boldsymbol{0}}}}\quad {\boldsymbol{0}} \quad{\boldsymbol{0}}\quad{\cdots}\quad{\boldsymbol{0}}\quad{\boldsymbol{I}}\quad{\boldsymbol{0}} \quad{\boldsymbol{0}}\quad{\cdots}\quad{\boldsymbol{0}}\right]^{\boldsymbol{{\rm T}}},\\ {{{\boldsymbol{C}}}_a} = \left[{{{\boldsymbol{A}}_1}}\quad{{{\boldsymbol{A}}_2}}\quad \cdots \quad{{{\boldsymbol{A}}_{na - 1}}}\quad{{{\boldsymbol{A}}_{na}}}\quad{{{\boldsymbol{B}}_1}}\quad{{{\boldsymbol{B}}_2}}\quad \cdots \quad{{{\boldsymbol{B}}_{nb - 2}}}\quad{{{\boldsymbol{B}}_{nb - 1}}} \right],\\ {{{\boldsymbol{D}}}_a} = \left[ {{{\boldsymbol{B_0}}}} \right]. \end{gathered}\right\} \end{align}

To couple the structural equations, the discrete-time form equation is turned into the continuous-time form, and the model in the state-space form is constructed as follows:

(2.9)\begin{equation} \left.\begin{array}{@{}l@{}} {{{\boldsymbol{{\dot x}_f}}}}(t) = {{\boldsymbol{A}}_a}{{\boldsymbol{x_f}}}(t) + {{\boldsymbol{B}}_a}{\boldsymbol{u}}(t)\\ {{\boldsymbol{y_f}}}(t) = {{\boldsymbol{C}}_a}{\boldsymbol{x_f}}(t) + {{\boldsymbol{D}}_a}{\boldsymbol{u}}(t). \end{array} \right\} \end{equation}

2.4. Aeroelastic models

The structural state vector is defined as ${{\boldsymbol {x_s}}} = {[\theta,{h}/{b},\dot {\theta },{{\dot {h}}}/{b}]^\textrm {T}}$ for the 2DOF case, and the structural motion equation (2.2) in state-space form can be rewritten as

(2.10)\begin{equation} \left.\begin{array}{@{}l@{}} {{{\boldsymbol{{\dot{x}}_s}}}}(t) = {{\boldsymbol{A}}_s}{{\boldsymbol{x_s}}}(t) + k_{\alpha} \boldsymbol{\cdot} {{\boldsymbol{B}}_s}{{\boldsymbol{y_f}}}(t)\\ {\boldsymbol{u}}(t) = {{\boldsymbol{C}}_s}{{\boldsymbol{x_s}}}(t) + k_{\alpha} \boldsymbol{\cdot} {{\boldsymbol{D}}_s}{{\boldsymbol{y_f}}}(t), \end{array} \right\} \end{equation}

where

(2.11a–d)\begin{align} {{\boldsymbol{A}}_{\boldsymbol{s}}} = \left[{\begin{array}{@{}cc@{}} \boldsymbol{0} & \boldsymbol{I}\\ {{\boldsymbol{M}^{ - 1}} \boldsymbol{\cdot} \boldsymbol{K}} & {{\boldsymbol{M}^{ - 1}} \boldsymbol{\cdot} \boldsymbol{G}} \end{array}}\right],\quad {{\boldsymbol{B}}_{\boldsymbol{s}}} = \left[{\begin{array}{@{}c@{}} \boldsymbol{0}\\ {{\boldsymbol{M}^{ - 1}}} \end{array}}\right],\quad {{\boldsymbol{C}}_{\boldsymbol{s}}} = \left[ \boldsymbol{I}\quad \boldsymbol{0} \right],\quad {{\boldsymbol{D}}_{\boldsymbol{s}}} = [ \boldsymbol{0} ]. \end{align}

Coupling the structural state equation (2.10) with the aerodynamic state equation (2.9), we can get the state equation and output equation of the aeroelastic system as follows:

(2.12)\begin{equation} \left\{ {\begin{array}{@{}c@{}} {{{{\boldsymbol{{\dot{x}}_s}}}}(t)}\\ {{{{\boldsymbol{{\dot{x}}_f}}}}(t)} \end{array}} \right\} = \left[{\begin{array}{@{}cc@{}} {{{\boldsymbol{A}}_s} + k_{\alpha} \boldsymbol{\cdot} {{\boldsymbol{B}}_s}{{\boldsymbol{D}}_a}} & {k_{\alpha} \boldsymbol{\cdot} {{\boldsymbol{B}}_s}{{\boldsymbol{C}}_a}}\\ {{{\boldsymbol{B}}_a}{{\boldsymbol{C}}_a}} & {{{\boldsymbol{A}}_a}} \end{array}} \right]\left\{{\begin{array}{@{}c@{}} {{{\boldsymbol{x_s}}}(t)}\\ {{{\boldsymbol{x_f}}}(t)} \end{array}}\right\} . \end{equation}

The analysis of the flutter characteristics is then used to solve the complex eigenvalues in (2.12). The process of the transonic 2DOF flutter analysis is as follows.

1. (1) Use URANS methodology flow solver to compute the mode aerodynamic coefficients based on the designed input signals.

2. (2) Use the identification technique to construct the input–output difference model (2.7) in the discrete-time domain, and then turn it into the continuous-time state-space form (2.9).

3. (3) Couple the aerodynamic state-space equations (2.9) with the structural state-space equations (2.10), to obtain the aeroelastic state-space equation (2.12).

4. (4) Solve the eigenvalues of the state matrix in (2.12) at different mass ratios and ratios of structural natural frequencies (${i_n}={k_h}/{k_{\alpha }}$); the stability of the aeroelastic system can be analysed by the root loci method. The real and imaginary parts of the eigenvalues are the growth rate and frequency of the aeroelastic system. The damping is usually defined as the opposite of the growth rate. When the growth rate is greater than zero, the system is unstable. On the contrary, the system is stable.

It is noteworthy that the aeroelastic model constructed in this paper is based on a linear model under small disturbances, possessing both time-domain and frequency-domain solving capabilities. The aeroelastic model can accurately and rapidly predict the instability boundary and type of the aeroelastic system, but it cannot precisely predict the nonlinear response of the aeroelastic system.

3.1. Construction and validation of the aerodynamic model

Training signal design is the key to dynamics modelling. The training signal should cover the frequencies of the structural modes that need to be excited. The training amplitude should not be too large; otherwise, the aerodynamic responses will be nonlinear, and the linear modelling method cannot be used. Based on similar considerations, random signals that more easily excite flow nonlinearity are not suitable for linear modelling; generally, sweep signals are used. At that pre-buffet state, the flow field still has some potential unsteady features, which may be excited by some very small disturbances. Therefore, we take the steady-state flow field as the initial solution in the training process. The time step of the training signal needs to be moderate: too small a time step can lead to excessively high model dimensions, while too large a step can cause the model to lose high-frequency flow features, potentially resulting in modelling failure. Based on past research experience, it is generally recommended that the entire cycle corresponding to the highest frequency in the training signal should ideally contain 25–30 time steps. In this paper, a chirp signal and its power spectrum density (PSD) are designed as depicted in figure 10. The reduced frequency of the signal is in a broad band of $[0.04, 0.80]$. Use of CFD yields the aerodynamic response as the output.

Figure 10. (a) The training signal with (b) its PSD.

The relative errors (3.1) of model training are shown in figure 11:

(3.1)\begin{equation} \text{Relative error} = \sum_{t = 0}^{{t_n}} {\frac{{\left| {{\boldsymbol{y}}_t - {\boldsymbol{y}}_{t,pre}} \right|}}{{\left| {{\boldsymbol{y}}_t - \overline {{\boldsymbol{y}}_t} } \right|}}}. \end{equation}

The error is up to $10\,\%$ when $na=nb=10$, while for $na=nb=30$ the error can be reduced to $2.05\,\%$, which satisfies the requirement of accuracy in the subsequent analysis. Therefore, in the present study, the orders of the ROM are $na=nb=30$, which are larger than those in subsonic or supersonic potential flow. The training results of the ROM are depicted in figure 12.

Figure 11. Identified errors with different $(na,nb)$ orders.

Figure 12. The training results of the ROM.

The model is then validated by comparing the time history of the lift and moment coefficient response between the ROM and CFD. The test signal is depicted in figure 13. The major non-dimensional frequency of signal A is 0.1 and that of signal B is 0.3. Figure 14 presents the time history of the lift and moment coefficient predicted by the ROM, which is in good agreement with that computed by CFD, further verifying the ROM in the present study.

Figure 13. The validation signals. (a) Signal A at ${k = 0.1}$. (b) Signal B at ${k = 0.3}$.

Figure 14. Comparison of the time history of the aerodynamics responses. (a,c) The response caused by the excitation of signal A. (b,d) The response caused by the excitation of signal B.

3.2. Identification of the fluid mode

The fluid eigenmode represents the stability of the fluid itself. It significantly affects the coupled-mode flutter and plays a vital role in the SDOF flutter of a spring-mounted airfoil. Therefore, it is necessary to identify and characterize the fluid eigenmode.

The free-stream state is $Ma = 0.7$, $Re = 3 \times {10^6}$, $\alpha = 4.5^\circ$, $a=0.224c$, which is identical to the case presented in § 3.1. Figure 15 shows the distribution of the eigenvalues of the $\boldsymbol {A}_a$ matrix of the model built before. The real part of the complex eigenvalue indicates the growth rate ($g$) of the fluid eigenmode, while the imaginary part indicates the reduced frequency ($k$). Matrix $\boldsymbol {A}_a$ is a redundant matrix, and the larger the orders are, the greater is the redundancy of the matrix. Some of the eigenvalues may not represent the real fluid mode. Therefore, the higher damped eigenvalues have a larger scatter with changing orders. However, it can be seen that an eigenvalue almost stays at the point of ($-0.026, 0.155$) with order increasing. This is closely similar to the fluid mode observed in the previous transonic aerodynamic modelling of an SDOF airfoil, as indicated by the black filled circles in figure 15 (Gao et al. Reference Gao, Zhang and Ye2016a). It is also noted that this eigenvalue is the one closest to the imaginary axis, indicating the least stable eigenmode, which implies that its stability margin is the lowest. In the present study, we define this corresponding mode as ‘dominant fluid mode’ (Gao et al. Reference Gao, Zhang and Ye2016a), also called ‘least stable mode’ (Sartor et al. Reference Sartor, Mettot and Sipp2015) or ‘unstable global mode’ (Crouch et al. Reference Crouch, Garbaruk, Magidov and Travin2009), which is associated with the flow damping and frequency. The fluid mode mentioned later refers to the dominant fluid mode identified here.

Figure 15. Eigenvalues of the aerodynamic models for the 2DOF airfoil at various orders; ${Ma = 0.7, Re = 3 \times {10^6}, \alpha = 4.5^\circ }$ (Gao et al. Reference Gao, Zhang and Ye2016a).

This mode and its eigenvalue in the subcritical state can be explicitly displayed by the linear ROM. However, it is difficult to extract this mode in numerical simulations and wind tunnel experiments, because this is a steady state. Here, modal decomposition has been performed by the DMD method for a subcritical buffet flow in the converging process, corresponding to the red dashed line in figure 5. The first four dominant global modes characterized by the pressure contours are shown in figure 16, and the growth rates and reduced frequencies are summarized in table 1. All modes are conjugate modes except the first one; thus a total of seven modes are selected. It can be noticed that both the growth rate and the frequency are zero for the first mode. It is a static mode, close to the mean flow field. All the other modes reflect the oscillating features resulting from shock waves. From table 1, the growth rates of these modes are negative because the snapshots are recorded from the converging process. We also note that the reduced frequency of Mode 2 is 0.162, which is equal to the buffet frequency from the CFD simulation, as shown in figure 5. And both of these are close to the frequency of the dominant fluid mode obtained from the unsteady aerodynamic ROM. The growth rate of Mode 2 is $-0.045$, which differs from the $-0.025$ obtained using the unsteady aerodynamic ROM. The root of the discrepancy lies in the different flow conditions: the DMD modes are based on a stationary airfoil in the converging process, whereas the unsteady aerodynamic ROM uses aerodynamic responses under structural perturbations. The presence of structural perturbations reduces the stability margin of the pre-buffet flow, consistent with the current findings. The absolute values of the growth rates of Mode 3 and Mode 4 are higher than that of Mode 2. Therefore, Mode 2 is the most important coherent global mode for the present pre-buffet flow, which should be focused on to address the most relevant changes in the flow field and physical mechanisms.

Figure 16. First four dominant global modes from the DMD method of a subcritical buffet flow in the converging process: (a) Mode 1, (b) Mode 2, (c) Mode 3 and (d) Mode 4.

Table 1. Growth rates and reduced frequencies of first four dominant DMD modes.

Further validation of the linear model and fluid modes in this paper is presented in Appendix C.

4.1. Dynamic characteristics under typical frequency ratio

This section is dedicated to studying the stability of the aeroelastic system of 2DOF airfoil in the transonic pre-buffet flow and the influences of the related structural parameters, under typical frequency ratios. The aeroelastic model has been established in the state-space form. Solving the eigenvalues of the state matrix in (2.12) at different structural parameters, the stability of the aeroelastic system can be analysed by the root loci method.

The free-stream state and the aerodynamic model have been introduced in § 3.2. For the structural parameters $\mu =200$, ${r_a^2}=1.036$, $x_a=0.2$, ${k_h} / {k_{\alpha }} = 0.85$, the eigenvalue loci as the structural natural pitching frequency ${k_{\alpha }}$ decreases are depicted in figure 17(a). The horizontal (vertical) axes represent the real (imaginary) parts of the eigenvalue. The real parts mean the growth rate of the coupled aeroelastic system and the imaginary parts represent the frequency. In our aeroelastic model, negative growth rate (positive damping) represents that the system is stable, while positive growth rate (negative damping) indicates that the system is unstable. It can be found that there are a total of three branches, which we label herein as Branch 1 (the black circle line), Branch 2 (the blue triangle line) and Branch 3 (the red square line). The purple point represents the least stable fluid mode and the arrow points in the direction of decreasing natural frequency. The eigenvalues of Branch 1 march towards the imaginary axis and further cross into the right half-plane, resulting in instability (flutter) of the coupled system. The flutter boundary predicted by the aeroelastic model is at ${k_{\alpha }}=0.37$.

Figure 17. Eigenvalue loci of the aeroelastic system at ${{k_h} / {k_{\alpha }} = 0.85}$. (a) The root loci. (b) Real parts and (c) imaginary parts of the root loci change with the natural pitching frequency.

Figure 18 shows a comparison of the time history of the heaving displacement and the pitching angle between the coupled CFD and CSD methods and the aeroelastic model simulation in the vicinity of the flutter boundaries at ${k_h} / {k_{\alpha }} = 0.85$. The results computed by the aeroelastic model agree well with those obtained by the coupled CFD and CSD methods, further verifying the aeroelastic model in the present study. The flutter boundary at various frequency ratios has also been confirmed through the CFD–CSD method, as denoted in figure 27 of § 5 with red open squares.

Figure 18. Comparison of the time history of the heaving displacement and the pitching angle between the coupled CFD–CSD method and the aeroelastic model, $Ma=0.7, Re=3 \times {10^6}, \alpha =4.5^\circ, {k_h}/{k_{\alpha }}=0.85, \mu =200$, for (a,c) ${{k_{\alpha }}=0.36}$ and (b,d) ${{k_{\alpha }}=0.38}$.

The aeroelastic system consists of three main modes: fluid mode, heaving mode and pitching mode. However, the relationship between branches of the eigenvalue loci and coupled modes of the aeroelastic system cannot be directly reflected in figure 17(a). It is necessary to determine which modes these three branches represent at each value of the natural frequency. Indeed, the hidden coupled pitching mode can be represented by Branch 1 at one value of the natural frequency and switch to Branch 2 or Branch 3 at another value of the natural frequency (or vice versa). This process is described as mode veering by Gao et al. (Reference Gao, Zhang, Li, Liu, Quan, Ye and Jiang2017b).

For further analysis, the real and imaginary parts of the eigenvalue varying with the natural pitching frequency are displayed in figures 17(b) and 17(c), respectively. The grey and blue colour blocks respectively represent the region of heaving instability and that of pitching instability. Since the coupled modes switch between branches with the change of the natural frequency, the left and right sides of branches in figure 17(c) are marked with the coupled mode types: coupled fluid mode (FM), coupled heaving mode (HM), coupled pitching mode (PM). The grey, blue and orange dashed lines are used to mark the uncoupled heaving, pitching and fluid modes.

The type of the coupled mode represented by a branch can be determined by the imaginary parts of the branch. The frequency of the fluid mode is almost constant, and the heaving and pitching modes have a fixed frequency ratio. Based on this, the types of coupled modes are distinguished. Certainly, in fluid–structure coupled aeroelastic systems, the modes reflected by the system's eigenvalue are all coupled. The purpose of resolving these modes through the aforementioned method is to better understand the coupled system through a linear aeroelastic model.

The coupled modes represented by the three branches undergo complex mode veering processes. As the natural pitching frequency $k_{\alpha }$ decreases (the non-dimensional velocity $V^*$ increases), Branch 1 veers from HM to PM; Branch 2 veers from PM to FM; and Branch 3 veers from FM to HM. The mode veering occurs within the instability region ($0.02, 0.37$) of Branch 1, so there are two kinds of flutter before and after the mode veering. As shown in figure 12, the mode veering phenomenon is actually a continuous process with $k_{\alpha } = (0.14,0.22)$, rather than a sharp transition. Within this transition range, it is difficult to determine the correspondence between the branches and the coupled modes. For simplification, we define the boundary at $k_{\alpha } = 0.22$, where the response frequency of Branch 1 is in the middle of the natural frequencies of the structural heaving and pitching modes. Consequently, the instability range of Branch 1 is qualitatively divided into the PM and HM ranges.

The grey block represents the SDOF heaving flutter region at $k_{\alpha } = (0.22, 0.37)$. In the initial stage of the region (right-hand side), the coupled frequencies of the heaving and pitching modes do not shift relative to their natural frequencies. So it should be defined as the SDOF heaving flutter. This region is far from the instability boundary of the SDOF heaving airfoil ($k_h<0.14$). On the contrary, it is within the instability boundary of the SDOF pitching airfoil ($0.13< k_{\alpha }<0.42$). This special phenomenon is caused by the coupling of the fluid mode, heaving mode and pitching mode. Therefore, if the pitching mode is decoupled from the 2DOF aeroelastic system, the SDOF heaving flutter phenomenon will no longer exist in this region, although the pitching mode does not obviously participate in the coupling.

Figure 18 also corroborates our judgment on SDOF heaving flutter. We present the wavelet transform diagrams of the time-domain responses calculated by CFD–CSD in figure 18(c), as shown in figure 19. It can be seen that, due to the initial disturbances applied to both the heaving and pitching modes, their responses initially vibrate in their natural modes. Over time, the natural modes transition to the coupled modes. It can be observed that the natural pitching mode gradually disappears and transitions to the coupled mode, which is almost identical to the natural heaving mode.

Figure 19. Wavelet transform diagrams for the responses of pitching mode calculated by CFD–CSD at ${k_h}/{k_{\alpha }}=0.85$ and $k_{\alpha }=0.36$.

The blue block at $k_{\alpha } = (0.02, 0.22)$ represents the coupled-mode flutter, in which the pitching mode is unstable. This is based on the fact that Branch 1 represents the coupled pitching mode after the mode veering. At the beginning of this region, the frequency of the coupled pitching mode is between the uncoupled natural pitching frequency and the natural heaving frequency. The structural natural frequency is low, which is different from the instability boundary of the SDOF pitching airfoil.

4.2. The effect of structural parameters

Regarding the previous paragraph on the determination of two distinct vibration types, further analysis and explanation are required. We analyse the influences of two structural parameters ($x_a, \mu$) on the aeroelastic system. Firstly, the effect of $x_a$ changes on the aeroelastic system is investigated. Parameter $x_a$ is the airfoil static unbalance (dimensionless distance between the centre of gravity and the elastic axis), which is an important parameter that affects the coupling of heaving and pitching modes. When $x_a$ is close to zero or negative, the critical dimensionless wind speed of the coupled-mode flutter is large. Conversely, when $x_a$ increases, the critical wind speed of the flutter will decrease. The value of $x_a$ is set as 0.0035 in the BACT model. In our simulations, the airfoil is fixed at $0.224c$ from the leading edge ($a=-0.554b$) in both the CFD–CSD simulation and the aeroelastic model. The change in $x_a$ actually reflects a change in the centre of gravity position, meaning that it only requires modifying the structural model rather than the unsteady aerodynamic model. Therefore, this study can conveniently investigate the impact of changes in $x_a$ on the aeroelastic system.

A typical frequency ratio of 0.85 is selected, with $x_a=(0.0035, 0.075, 0.1, 0.2)$. For convenience of understanding, only the unstable branches of the eigenvalue are presented in figure 20(a). It can be observed that as $x_a$ increases, the upper boundary of the instability region decreases, while the lower boundary rapidly increases. Around $x_a=0.1$, the coupling form of the aeroelastic system changes, and the instability region at low frequency increases from scratch. To further investigate this phenomenon, the root locus of the aeroelastic system at $x_a = 0.0035$ and $k_h / k_{\alpha } = 0.85$ has been added, as shown in figure 20(b,c). Branch 1 always represents the structural heaving mode. There is a transition between Branch 2 and Branch 3, where the fluid mode veers to the pitching mode. Only Branch 2 has an unstable region. At the right-hand boundary of the instability region, the response frequency is locked to the natural pitching frequency. At the left-hand boundary, the response frequency is equal to the flow frequency. Therefore, it should be defined as an SDOF pitching flutter induced by the fluid mode. This is a completely different coupling form from that of the aeroelastic system at $x_a=0.2$, as shown in figure 17(b,c). The coupled-mode flutter and SDOF heaving flutter that occurred at ${x_a=0.2}$ no longer manifest at $x_a=0.0035$. In summary, changes in $x_a$ can alter the region of instability, the mode of instability and the form of coupling. Increasing $x_a$ will enhance the coupling between structural modes and suppress the SDOF flutter.

Figure 20. (a) The unstable branch of the eigenvalue loci at ${{k_h} / {k_{\alpha }} = 0.85}$ with different ${x_a}$. (b) Real parts and (c) imaginary parts of the root loci change with the natural pitching frequency at ${{k_h} / {k_{\alpha }} = 0.85}$ and ${x_a=0.0035}$.

The following explores the influences of mass ratio on the aeroelastic system. Similarly, a frequency ratio of 0.85 and a kind of mass ratio ($50, 100, 200, 400, 800$) are selected, as depicted in figure 21(a). It can be observed that as the mass ratio increases, the upper boundary of the instability interval decreases and the lower boundary also decreases. Around a mass ratio of 400, the system undergoes some changes. When the response frequency is around 0.2, the system returns to a stable state. To further investigate this phenomenon, the root locus of the aeroelastic system at $\mu = 800$ and $k_h / k_{\alpha } = 0.85$ has been added, as shown in figure 21(b,c). Branch 2 only represents the pitching mode and the fluid mode veers to the heaving mode. It can be clearly seen here that the instabilities of coupled pitching and heaving modes are separate, supporting the aforementioned judgment. Thus, the conclusion is that increasing the mass ratio can eliminate the mode veering phenomenon as much as possible without changing the vibration properties of the system.

Figure 21. (a) The unstable branch of the eigenvalue loci at ${{k_h} / {k_{\alpha }} = 0.85}$ with different ${\mu }$. (b) Real parts and (c) imaginary parts of the root loci change with the natural pitching frequency at ${{k_h} / {k_{\alpha }} = 0.85}$ and ${\mu }=800$.

4.3. Various aeroelastic phenomena under different natural frequencies

In § 4.1, we analysed the flutter characteristics and types of the aeroelastic systems for frequency ratios of 0.85, as $k_{\alpha }$ varies. However, one frequency ratio is merely demonstrated in one line within the two-dimensional parameter space composed of $k_h$ and $k_{\alpha }$. In order to establish a more comprehensive understanding of the flutter mechanism, this section analyses the flutter characteristics within the two-dimensional parameter space ($k_h, k_{\alpha }$) at $\mu = 200$ and ${x_a=0.2}$. Based on the analysis method in § 4.2, multiple distinct frequency ratios were analysed within ${i_n}=k_h/k_{\alpha }=(0.2,5)$. Then, the flutter characteristics were obtained by connecting the boundaries of various aeroelastic phenomena under different frequency ratios, as depicted in figure 22.

Figure 22. Various aeroelastic phenomena within the two-dimensional parameter space (${k_h},k_{\alpha }$).

There are six distinct aeroelastic phenomena: SDOF heaving flutter, SDOF pitching flutter, heaving instability within coupled-mode flutter, pitching instability within coupled-mode flutter, forced vibration and stable state. The horizontal (vertical) axes represent the natural dimensionless frequencies of airfoil pitching (heaving) mode. The colour blocks and numbers in the two-dimensional graphs represent the type of vibration. Multiple vibration types are converted near the $45^\circ$ symmetry dashed line. The wide variety of distinct aeroelastic phenomena is the result of the joint coupling of the fluid mode, structural heaving mode and structural pitching mode.

The following provides a detailed introduction to each of the six types. It is noteworthy that the stability boundaries predicted by the aeroelastic model are quite accurate, as demonstrated in figure 18 and the open red squares in figure 27. However, the boundaries between different aeroelastic phenomena are determined empirically. In fact, the transition is continuous rather than discrete. Therefore, these internal boundaries are merely schematic representations.

I – SDOF pitching flutter. Two grey areas. For the upper area, due to the high frequency of the heaving mode, the instability boundary is highly similar to that of the SDOF pitching airfoil. For the lower area, there are many complex phenomena due to the participation of the heaving mode. Coupled-mode flutter on the left-hand side (region III/IV) and SDOF heaving flutter on the right-hand side (region II) have encroached upon the region of SDOF pitching flutter.

II – SDOF heaving flutter. The lower areas of region II fall within the instability boundary of the SDOF heaving airfoil and are separated by region I. This occurs because the positive growth rate of SDOF pitching flutter is much greater than that of SDOF heaving flutter. Additionally, the remaining area lies below the $45^\circ$ line, where the natural heaving frequency exceeds the instability boundary of the SDOF heaving airfoil.

III$\boldsymbol{/}$IV – coupled-mode flutter. The distribution is symmetrical around the $45^\circ$ line within the range of low natural frequencies, with the upper part representing the instability of the heaving mode and the lower part representing the instability of the pitching mode. Upon comparing the black solid line with the open/filled circles, minimal differences can be observed in the boundaries of coupled-mode flutter. Additionally, numerous other complex phenomena depicted in figure 22 primarily originate from the fluid mode.

V – forced vibration. The areas with grey stripes. Existing research has found that elastic structures can cause a reduction in the transonic buffeting boundary (${Ma=0.7}, {\alpha } =3.9^\circ$) (Gao, Zhang & Ye Reference Gao, Zhang and Ye2018). This phenomenon is also observed in this study. In the area where the structural modes are stable and the natural frequency is low, forced vibration caused by buffeting flow occurs. Detailed analysis of this phenomenon is presented in § 5.3.

VI – stable state. In the top right corner, all fluid/structure modes of the aeroelastic system remain stable due to the high structural frequency. In the lower left corner, all modes are stabilized as a result of the complex coupling and interference between these modes.

These complex aeroelastic phenomena can be analysed from two perspectives using auxiliary lines. The first perspective is the effect of the fluid mode on a 2DOF aeroelastic system, which is discussed in detail in § 4.4. The second perspective is the effect of the coupling with another structural mode on the SDOF flutter (transonic buzz). In figure 22, the vertical dash-dot lines represent the flutter boundaries ($k_{\alpha } = (0.13, 0.42)$) of the SDOF pitching airfoil case depicted in figure 23(a), while the horizontal dash-dot line represents the flutter boundaries ($k_h = (0, 0.138)$) of the SDOF heaving airfoil case depicted in figure 23(b). The black dashed lines represent the LCO amplitudes computed using the CFD–CSD method at different natural frequencies. By comparing the results of figure 22 with the two cases mentioned above, the influence of the coupling of other structural modes on SDOF flutter can be observed. There are two distinct phenomena. The first is that the SDOF heaving flutter occurs at a high natural frequency (${i_n}=(0.8,1), {k_h}=(0.2, 0.4)$), which is discussed in §§ 4.1 and 4.2. The second is that the SDOF flutter disappears by coupling with another structural mode, at ${i_n}=(0.3,0.8)$, ${k_{\alpha }}=(0.3, 0.42)$ and ${i_n}=(0.3,0.5)$, ${k_{\alpha }}=(0, 0.12)$. A detailed study and mechanistic analysis of this phenomenon are presented in § 5.2.

Figure 23. The stability boundaries and the LCO amplitudes of SDOF aeroelastic system of NACA0012 airfoil. (a) The SDOF pitching aeroelastic system and (b) the SDOF heaving aeroelastic system.

4.4. The effect of the fluid mode on the flutter boundary

In figure 22, the black line with open circles represents the coupled-mode flutter boundary of the airfoil under the same parameters at ${\alpha } = 0^\circ$. While the black line with filled circles represents the flutter boundary of the airfoil under Euler solution at ${\alpha } = 4.5^\circ$. Our goal is to eliminate the dominant fluid mode in the above two ways. The first approach involves staying away from the buffeting boundary as much as possible, and the second approach is to eliminate the fluid viscosity. The basis for these two approaches is the paper of Karnick & Venkatraman (Reference Karnick and Venkatraman2017), which found that the viscous effects manipulate the flow by influencing the strength and location of the shock such that the flutter boundary changes significantly. The eigenvalues of the $\boldsymbol {A}_a$ matrix for both cases are depicted in figure 24. The most dominant eigenvalue becomes the zero-frequency mode near the point of ($-0.12, 0$). Notably, the dominant fluid mode at the point of ($-0.026, 0.155$) in the case under NS solution at ${\alpha } = 4.5^\circ$ has been successfully eliminated using both approaches. A more detailed discussion on the variation of fluid modes with the AOA is presented in Appendix C.

Figure 24. Eigenvalues of the aerodynamic models under different conditions.

By comparing the three cases mentioned above, we can qualitatively analyse the influence of the fluid mode on the aeroelastic system. The flutter boundaries of these cases are shown in figure 25. Case A corresponds to ${\alpha } = 4.5^\circ$ under NS solution, Case B represents ${\alpha } = 0^\circ$ under NS solution and Case C is at ${\alpha } = 4.5^\circ$ under Euler solution. The lower horizontal axis displays the arctangent of $k_h / k_{\alpha }$, while the upper horizontal axis represents $k_h / k_{\alpha } =i_n$ directly. The vertical axis indicates the non-dimensional flutter velocity (${{V_f}^*} = 1/(k_{{\alpha }f} {\mu ^{1/2}})$). The flutter boundaries of Case B and Case C are similar and exhibit an approximately axisymmetric distribution. Their flutter velocities decrease as ${k_h} / {k_{\alpha }}$ approaches one and increase as it deviates from one. When ${k_h} / {k_{\alpha }}<1$, Case B has lower flutter velocities compared to Case C, while the opposite holds true when ${k_h} / {k_{\alpha }}>1$. Notably, the flutter velocity of Case A is significantly smaller than that of Case B and Case C. Figure 26 shows the root loci for these three cases at ${k_h} / {k_{\alpha }} = 0.85$, corresponding to the blue dashed line in figure 25. Case A has been thoroughly discussed in §§ 4.1 and 4.2. The flutter boundary is located at ${{V_f}^*} = 0.19$, and the flutter type is SDOF flutter. For Case B and Case C, there is no mode veering between branches after the dominant fluid mode disappears. Branch 1 represents the coupled heaving mode, while Branch 2 represents the coupled pitching mode, as shown in figures 26(c)–26( f). Unlike Case A, under the critical state of flutter, the frequencies of the coupled modes of Case B and Case C approach each other, representing a coupled-mode flutter. The non-dimensional flutter velocities of Case B and Case C are 0.35 and 0.38, respectively, which are higher than that of Case A. The presence of the fluid mode transforms the type of flutter from coupled-mode flutter to SDOF flutter at the boundary, significantly reducing the flutter boundary. This finding is consistent with the experimental observations conducted by NASA (Rivera et al. Reference Rivera, Dansberry, Bennett, Durham and Silva1992; Dansberry et al. Reference Dansberry, Durham, Bennett, Rivera, Silva, Wieseman and Turnock1993), where an increase in $Ma$/AOA leads to a decrease in transonic flutter velocity, with the response frequency locked to the uncoupled structural frequency (vacuo natural frequency). These experimental results align with the findings and perspectives presented in this study, mutually reinforcing their validity.

Figure 25. The flutter boundaries varying with ${{k_h / k_{\alpha }}}$ under different conditions.

Figure 26. Real part and imaginary part of the root loci change with the natural pitching frequency at ${{k_h} / {k_{\alpha }} = 0.85}$: (a,b) ${{Ma = 0.7, {\alpha } = 4.5^\circ }}$ under NS solution; (c,d) ${{Ma = 0.7, {\alpha } = 0^\circ }}$ under NS solution; (e, f) ${{Ma = 0.7, {\alpha } = 4.5^\circ }}$ under Euler solution.

5.1. Coupled-mode flutter versus SDOF flutter

In § 4.3, we analysed the flutter characteristics within the two-dimensional parameter space ($k_h, k_{\alpha }$), and discovered six different aeroelastic phenomena. Among them, there are four types of flutter, including pitching instability within coupled-mode flutter, heaving instability within coupled-mode flutter, SDOF pitching flutter and SDOF heaving flutter. They can be classified into two categories: coupled-mode flutter and SDOF flutter. The observed phenomena and their categorization are based on the frequency characteristics of unstable coupled modes predicted by the aeroelastic model. When the frequency of the unstable coupled mode is locked to its natural frequency, it is classified as an SDOF flutter. On the other hand, when the frequencies of the two coupled modes deviate from their natural frequencies and approach each other, it is considered a coupled-mode flutter. In order to verify the results of the aeroelastic model and to clarify the behaviour and characteristics of various flutters in the time-domain nonlinear response, we selected one state from each of the four types of flutter to demonstrate and analyse their time-domain responses obtained from CFD–CSD, as shown in figures 28–31, corresponding to points P1 to P4 in figure 27. The vertical axis range of figures 28–31 was set to the same value.

Figure 29. Heaving instability of coupled-mode flutter at P2: ${k_h} = 0.18$, ${k_{\alpha }} = 0.15$. (a) Force coefficients, (c) displacements, (b,d) PSD, ${k_{ae}=0.167=0.93k_h}$.

Figure 30. The SDOF pitching flutter at P3: ${k_h} = 0.1$, ${k_{\alpha }} = 0.2$. (a) Force coefficients, (c) displacements, (b,d) PSD, ${k_{ae}=0.2=k_{\alpha }}$.

Figure 31. The SDOF heaving flutter at P4: ${k_h} = 0.26$, ${k_{\alpha }} = 0.31$. (a) Force coefficients, (c) displacement, (b,d) PSD, ${k_{ae}=0.26=k_{h}}$.

Figure 28 shows the time-domain response of the pitching instability with coupled-mode flutter at P1. Figure 29 depicts the time-domain response of the heaving instability within coupled-mode flutter at P2. For coupled-mode flutter, the amplitudes of force coefficients and displacements rapidly increase. This linear growth persists until the limit of the CFD–CSD method employed in this study, with the final obtained amplitudes before computational failure being $A_{h/b} > 1$ and $A_{\theta } > 0.6$. Figures 30 and 31 depict the time-domain response of the SDOF pitching (heaving) flutter in LCOs at P3 (P4), respectively. For SDOF flutter, the response of the aerodynamic forces and displacement eventually reach a LCO.

To more intuitively display these four types of flutter, we present the pressure and streamline diagrams at four typical moments for P1 to P4, as shown in figures 32–35. The horizontal dashed lines indicate the $y=0$ position, assisting in illustrating the elastic deformation of the airfoil. The typical moments are divided into four equal parts of a full cycle of structural modal displacement, with P1 and P3 based on the pitching mode displacement, and P2 and P4 based on the heaving mode displacement.

Figure 32. Pressure contours and streamlines in one period of pitching mode at P1: (a) $t_1$, (b) $t_2$, (c) $t_3$, (d) $t_4$.

For P1 and P2 states, both the pitching and heaving mode amplitudes are large, but their phase angles differ. When the pitching mode displacement of P1 is at its maximum (figure 32a), the heaving mode displacement is close to zero; conversely, when the heaving mode displacement is at its maximum (figure 32b), the pitching mode displacement is close to zero. The pattern for the P2 state is different: when the heaving mode displacement is at its maximum (figure 33a), the pitching mode displacement is also significant. For the P3 and P4 states, there is a significant difference in the displacement amplitudes of the pitching and heaving modes. The P3 state is dominated by pitching mode vibrations (figure 34), while the P4 state is dominated by heaving mode vibrations (figure 35). This pattern is more clearly illustrated in Lissajous figure 36, where there are evident differences in the vibration forms among the four states.

Figure 33. Pressure contours and streamlines in one period of heaving mode at P2: (a) $t_1$, (b) $t_2$, (c) $t_3$, (d) $t_4$.

Figure 34. Pressure contours and streamlines in one period of pitching mode at P3: (a) $t_1$, (b) $t_2$, (c) $t_3$, (d) $t_4$.

Figure 35. Pressure contours and streamlines in one period of heaving mode at P4: (a) $t_1$, (b) $t_2$, (c) $t_3$, (d) $t_4$.

Figure 36. Lissajous figure of the pitching and heaving mode response and airfoil vibration range diagram: (a) P1, (b) P2, (c) P3, (d) P4.

The time-domain responses of coupled-mode flutter and SDOF flutter exhibit many different characteristics, with the displacement response curve being more concise and exhibiting better regularity. Below, we extract various features of the displacement response curve for comparison. These features primarily include ${\eta }_{ae}$, amplitude ratio $R_A$ and phase angle ${\varphi }_{\theta -h}$, as illustrated in (5.1):

(5.1)\begin{equation} \left.\begin{array}{@{}l@{}} {\eta _{ae}} = {{\left| {{k_{ae}} - {k_{\alpha}}} \right|} / {\left| {{k_h} - {k_\alpha }} \right|}},\\ {R_A} = {A_{h/b}}/{A_\theta },\\ {\varphi _{h/b - \theta }} = {{\left| {{t_{max (h/b)}} - {t_{max (\theta )}}} \right|} / {\left( {{T_{ae}}} \right)}}. \end{array} \right\} \end{equation}

The meanings of the parameters in (5.1) are shown in figure 37, where ${\eta }_{ae}$ donates the deviating degree of the response frequency compared with the natural frequency of the unstable mode. When ${\eta }_{ae}=1$, $k_{ae}={k_\alpha }$. When ${\eta }_{ae}=0$, $k_{ae}={k_h}$. The amplitude of the non-dimensional heaving displacement is denoted by $A_{h/b}$, and the amplitude of the pitching angle (in radians) is represented by $A_{\theta }$. The period of the response curve is denoted by ${T_{ae}}$, while $t_{max (\theta )}$ and $t_{max (h/b)}$ represent the non-dimensional time corresponding to the maximum pitching displacement and heaving displacement within the same period.

Figure 37. A schematic diagram illustrating the meanings of the parameters in (5.1). (a) Displacements and (b) PSD.

The features of the above four types of flutter are shown in table 2, where P1, P2, P3 and P4 represent pitching instability within coupled-mode flutter, heaving instability within coupled-mode flutter, SDOF pitching flutter and SDOF heaving flutter, respectively. For coupled-mode flutter, since the amplitude of the response curve keeps increasing, the features are extracted for each period within the range shown in figures 28 and 29, and then averaged. For SDOF flutter, these features are directly extracted from the LCOs shown in figures 30 and 31.

Table 2. Various features of the displacement response of six aeroelastic phenomena.

For ${\eta }_{ae}$ in table 2, it can be observed that at P1 and P2, their response frequencies lie between the natural frequencies of the structural modes, and they are closer to the unstable mode. On the other hand, the response frequencies at P3 and P4 are equal to the unstable structural mode natural frequencies. For $R_A$, both structural modes at P1 and P2 exhibit significant vibrations, with the ratio between them being close to one. Meanwhile, the values of $R_A$ at P3 and P4 are significantly far from 1. Regarding ${\varphi }_{\theta -h}$, there is a noticeable phase lag between the heaving displacement and pitching displacement under P1 and P2 states, which is approximately 15 % to 40 % of the vibration period. In contrast, there is almost no phase difference between the two modes under P3 and P4 states, with the difference being less than 5 %. These differences and characteristics fully demonstrate that the aeroelastic phenomenon under P1 and P2 states corresponds to coupled-mode flutter, while the aeroelastic phenomenon under P3 and P4 states belongs to SDOF flutter. It should be noted that P4 indicates a specific SDOF heaving flutter. Here, ${k_h}=0.26$, which is significantly higher than the flutter boundary for SDOF heaving airfoil ($0, 0.138$), as shown in figure 23. This special phenomenon has been analysed in § 4.1 by the aeroelastic model. The results of the CFD–CSD method confirm the accuracy of the aeroelastic model presented in this paper. The last two rows in table 2 also provide the eigenvalues at P7 and P6. State P7 represents a forced vibration state characterized by the dominant response frequency being equal to the frequency of the fluid mode. State P6 is a stable state where there is no periodic oscillation in the force coefficient and displacement response. The dynamical characteristics at P7 and P6 are discussed in detail in § 5.2.

The phenomenon of SDOF flutter that arises in the aeroelastic system of a 2DOF airfoil is worthy of further discussion. This phenomenon has been commonly observed in previous transonic aeroelastic problems. It can be categorized into two types. The first one is just like the SDOF flutter that occurs on transonic buzz and the frequency lock-in in transonic buffeting flow. The fluid mode is coupled with a structural mode, leading to the instability of that structural mode. The addition or removal of another structural mode has minimal impact on this SDOF flutter system. The second type occurs in a transonic flutter. As $Ma$/AOA increases, the critical velocity of transonic flutter rapidly decreases. The system's basic vibration frequency shifts from being between the natural frequencies of mode A and mode B to almost aligning with the natural frequency of mode A. Mode A and mode B mentioned here are just used as examples. Significantly different from the first type, when mode B is eliminated, the second type of SDOF flutter will no longer occur. Dowell conducted a study on the transonic flutter of 2DOF airfoils based on numerical simulation methods (Kholodar et al. Reference Kholodar, Dowell, Thomas and Hall2004). This study discovered and noted that: ‘The flutter frequency at that point is essentially the same as the coupled in vacuo natural frequency corresponding to a dominant plunge motion. Correspondingly, the flutter eigenvector is dominated by the plunge motion. This is an example of single-DOF flutter, but note that the critical aeroelastic mode is a mass coupled natural mode, albeit one that is plunge dominated’. In the present results shown in figure 22, region I and the lower right corner of region II belong to the first type of SDOF flutter. The upper section of region II is the second type. We found that the second type of SDOF flutter, which involves flow, heaving and pitching modes deeply coupled, is a typical result. This type of flutter will not occur in the absence of any mode.

5.2. Disappearance of SDOF flutter

In figure 22, we observed that the SDOF flutter disappears due to coupling with another structural mode. This phenomenon can be categorized into two regions. The first one is the reduction in SDOF-pitching-flutter boundaries due to coupling with the heaving mode, occurring at (${i_n}=k_h/k_{\alpha }=(0.3,0.8), {k_{\alpha }}=(0.3, 0.42)$). The second is that the SDOF heaving flutter disappears by the coupling with the pitching mode, happening at (${i_n}=(0.3,0.5), {k_{\alpha }}=(0, 0.12)$). It is evident that the line with $k_h/{k_{\alpha }} = 0.5$ intersects these two regions, as depicted in figure 27. The subsequent analysis focuses on the root locus at ${k_h} / {k_{\alpha }} = 0.5$ to elucidate the underlying mechanism of this phenomenon, as illustrated in figure 38. The mode veering occurs between Branch 1 and Branch 3. Branch 2 consistently represents the pitching mode, and only Branch 2 exhibits an unstable region. At the instability region, the frequency of the coupled pitching mode locks onto the natural pitching frequency. Thus, it is defined as SDOF pitching flutter. The right flutter boundary of the 2DOF aeroelastic system is located at $k_{\alpha }=0.35$. This is significantly lower than that ($k_{\alpha }=0.42$) for an SDOF pitching airfoil, as shown in figure 39(a,b). In other words, at ${k_h} / {k_{\alpha }} = 0.5$, there is a reduction in the SDOF-pitching-flutter boundaries in the 2DOF aeroelastic system compared with an SDOF pitching airfoil. The left boundary of the 2DOF aeroelastic system is located at $k_{\alpha }=0.12$, $k_h=0.06$. For an SDOF heaving airfoil, the region of SDOF flutter is $k_h=(0, 0.138)$, as shown in figure 39(c,d). In other words, at ${k_h} / {k_{\alpha }} = 0.5$, compared with an SDOF heaving airfoil, the SDOF heaving flutter disappears in the 2DOF aeroelastic system. Since the frequency ratio of the 2DOF aeroelastic system is 0.5, the range of the horizontal axis in figure 39(a,b) has been set to $k_{\alpha }=(0, 0.5)$. The range of the horizontal axis in figure 39(c,d) has been set to $k_{\alpha }=(0, 0.25)$. This allows for a direct comparison between figures 39 and 38(b,c).

Figure 38. (a) Eigenvalue loci of the 2DOF aeroelastic system at ${{k_h} / {k_{\alpha }} = 0.5}$. (b) Real parts and (c) imaginary parts of the root loci changes with the natural pitching frequency.

Figure 39. Real part and imaginary part of the root loci of the SDOF airfoil. (a,b) The SDOF pitching airfoil. (c,d) The SDOF heaving airfoil.

For the SDOF airfoils shown in figure 39, (2.2) and (2.3) degenerate into the following forms:

(5.2) \begin{equation} \left.\begin{array}{@{}ll@{}} \text{SDOF heaving case:} & {\dfrac{{\ddot{h}}}{b} + {k_h^2}\dfrac{h}{b} = \dfrac{1}{{{\rm \pi} \mu }}( - {C_l})} \\ \text{SDOF pitching case:} & {\ddot{\theta} + {k_\theta^2}\theta = \dfrac{1}{{{\rm \pi} \mu r_\theta ^2}}(2{C_m})}, \end{array} \right\} \end{equation}

where $\mu =200$, ${r_a^2}=1.036$. The free-stream conditions are $Ma=0.7$, $Re=3 \times {10^6}$, $\alpha =4.5^\circ$.

The instability of aeroelastic systems can be attributed to the coupling between different modes. Specifically, in SDOF airfoils, the fluid mode is only coupled with a structural mode, leading to instability. On the contrary, 2DOF aeroelastic systems exhibit coupling among multiple modes including fluid, pitching and heaving modes. Due to the presence of an additional structural mode, the energy of the fluid mode is distributed across two structural modes, delaying the onset of instability and reducing its unstable region. For instance, in figure 39(a,b), the coupling between pitching and fluid modes results in instability of the coupled pitching mode at $k_{\alpha }=0.42$. Meanwhile, in figure 39(c,d), near $k_h=0.21$, corresponding to $k_{\alpha }=0.42$, the natural frequencies of heaving mode and fluid mode become close, and mode veering occurs. In the 2DOF system, the mode veering between fluid mode and heaving mode weakens the ability of the fluid mode to induce instability of the pitching mode, thereby delaying the occurrence of SDOF pitching flutter. In figure 38(b,c), an ongoing SDOF pitching flutter inhibits the occurrence of SDOF heaving flutter in the range of $k_{\alpha }=(0.12, 0.35)$ and $k_h=(0.06, 0.175)$. In the range of $k_{\alpha }=(0, 0.12)$ and $k_h=(0, 0.06)$, the coupled pitching mode regains a stable state. However, the natural frequency of the pitching mode is closer to that of the fluid mode, which still suppresses the instability of the heaving mode, preventing SDOF heaving flutter. Therefore, under specific conditions, the release of other degrees of freedom can enhance aeroelastic stability. In conclusion, we have demonstrated that the coupling of different modes plays a crucial role in determining the stability of aeroelastic systems.

To further illustrate this point, we have performed CFD–CSD simulations for typical states denoted by the red vertical line (P5, P6) in figures 38 and 39, corresponding to figures 40 and 41. Compared with the flutter, the aerodynamic and displacement are stable at P5, as shown in figure 40. For the SDOF pitching airfoil, the pitching mode is unstable at $k_{\alpha } = 0.4$, and the vibration amplitude is large due to the SDOF pitching flutter, as depicted by the grey dotted line in figure 40. However, for the 2DOF aeroelastic system, the structural pitching mode remains stable due to the release of the heaving mode. Figure 41 shows the time-domain response at P6. The response amplitude is mainly due to the residual energy of the initial disturbance. For the SDOF heaving airfoil, the heaving mode is unstable at $k_h = 0.05$, and the vibration amplitude is large due to the SDOF heaving flutter, as depicted by the grey dotted line in figure 41. However, for the 2DOF aeroelastic system, the structural heaving mode remains stable due to the release of the pitching mode.

Figure 40. Stable state at P5: ${k_h} = 0.2$, ${k_{\alpha }} = 0.4$. (a) Force coefficients, (c) displacements, (b,d) PSD.

Figure 41. Stable state at P6: ${k_h} = 0.05$, ${k_{\alpha }} = 0.1$. (a) Force coefficients, (c) displacements, (b,d) PSD.

5.3. Fluid mode instability induced by the elastic structure

This section analyses the mechanism of phenomenon V – forced vibration – shown in figure 22. Existing research has found that elastic structures can cause a reduction in the transonic buffeting boundary ($Ma=0.7, {\alpha } =3.9^\circ$) (Gao et al. Reference Gao, Zhang and Ye2018). A stability analysis based on a linear ROM was conducted on an SDOF pitching airfoil at $\mu = 60$, $Ma=0.7$, ${\alpha } =4.5^\circ$, resulting in the root loci of the aeroelastic system as a function of the natural pitching frequency, as shown in figure 42(a,b). Only Branch 1 becomes unstable. Based on the relationship between the coupled frequency of Branch 1 and the natural frequencies of the structural/fluid modes, the dynamic modes are distinguished. At low pitching natural frequencies, it is forced vibration due to the instability of the dominant fluid mode; at high pitching natural frequencies, it is flutter due to the instability of the structural mode. This unique phenomenon is summarized as a reduction in the buffet boundary caused by the elastic structure.

Figure 42. Real part and imaginary part of the root loci change with the natural pitching frequency. (a,b) An SDOF pitching airfoil at $\mu = 60$, $Ma=0.7$, ${\alpha } =4.5^\circ$ (Gao et al. Reference Gao, Zhang and Ye2018). (c,d) A 2DOF airfoil at $\mu = 200$, ${k_h} / {k_{\alpha }} = 3$.

A similar phenomenon was observed in the study of a 2DOF airfoil, as shown in figure 42(c,d). In the root locus at $\mu = 200$, ${k_h} / {k_{\alpha }} = 3$, only Branch 2 becomes unstable. Based on the relationship between the coupled frequency of Branch 2 and the natural frequencies of the structural/fluid modes, three different types of dynamic modes were identified for this frequency ratio. The blue region at $k_{\alpha } =(0.15,0.42)$ represents pitching mode instability, the red region at $k_{\alpha } =(0.05,0.15)$ represents dominant fluid mode instability and the grey region at $k_{\alpha } =(0,0.05)$ represents heaving mode instability. Consequently, it can be concluded that in the red region, the instability of the dominant fluid mode occurs. To further illustrate this point, we have performed CFD–CSD simulations for typical states denoted by the red vertical line (P7) in figure 42(c,d), as shown in figure 43. It can be observed that under this condition, the coupled frequency $k_{ae}$ of the system is dominated by the buffeting frequency $k_b$, with the natural frequency of the structural mode completely dissipating, which verifies the aforementioned conclusion.

Figure 43. Forced vibration at P7: ${k_h} = 0.24$, ${k_{\alpha }} = 0.08$. (a) Force coefficients, (c) displacements, (b,d) PSD, ${k_{ae}=0.16=k_{b}}$.