3.1 Flow governing equations

In this study, the flow analysis of the NACA0012 airfoil is performed by an in-house hybrid-unstructured flow solver which solves the URANS equations by a cell-centred finite volume approach. The integral form of the 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:

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}=[\unicode[STIX]{x1D70C}\;\;\unicode[STIX]{x1D70C}u\;\;\unicode[STIX]{x1D70C}v\;\;\unicode[STIX]{x1D70C}E\;\;\unicode[STIX]{x1D70C}\tilde{\unicode[STIX]{x1D708}}]^{\text{T}}$ . $\unicode[STIX]{x1D70C}$ is the density; $u,v$ are the $x$ and $y$ components of the velocity vector of the flow; $E$ denotes the specific total energy; and $\tilde{\unicode[STIX]{x1D708}}$ 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. $\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 advection upstream splitting method (AUSM) $+$ scheme is conducted to evaluate the inviscid flux, with a reconstruction technique based on the least-squares approach. High-order reconstruction schemes refer to Liu et al. (Reference Liu, Zhang, Jiang and Ye2016). The grid velocity $\boldsymbol{V}_{grid}$ is modified according to the geometric conservation law (Hyung & Yannis Reference Hyung and Yannis2006) 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 $\hat{\unicode[STIX]{x1D710}}$ is given by the standard form proposed by Spalart & Allmaras (Reference Spalart and Allmaras1992). Many reports have indicated that the S–A turbulence model with a second-order scheme is accurate enough to predict the flow instability (Crouch, Garbaruk & Magidov Reference Crouch, Garbaruk and Magidov2007; Sartor et al. Reference Sartor, Mettot and Sipp2015a ; Sartor & Timme Reference Sartor and Timme2017). The convective and source terms are discretized by the second-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. Details of the numerical method can be found in Zhang et al. (Reference Zhang, Gao, Liu, Ye and Jiang2015b ) and Gao, Zhang & Ye (Reference Gao, Zhang and Ye2016a ).

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 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.

A moving boundary is involved in the simulation due to the airfoil pitching motion when the frequency lock-in is simulated. Thus, a grid deformation method must be used to match the grid with the new wall boundary condition. The scheme of the computational grid deformation is based on radial basis function (RBF) interpolation (De Boer, Van der Schoot & Bijl Reference De Boer, Van der Schoot and Bijl2007; Wang, Mian & Ye Reference Wang, Mian and Ye2015). A compact Wendland’s C2 function (De Boer et al. Reference De Boer, Van der Schoot and Bijl2007) is chosen as the basis function.

3.2 Validation of the CFD method

Validation of CFD results is performed based on a transonic buffet wind tunnel test by Doerffer et al. (Reference Doerffer, Hirsch and Dussauge2010). They performed a wind tunnel test on a NACA0012 airfoil at Mach numbers from 0.7 to 0.8, at a Reynolds number of $3\times 10^{6}$ and within a range of angles of attack including the buffet onset. In this experimental test, the buffet frequency at $M=0.7$ and $\unicode[STIX]{x1D6FC}=5^{\circ }$ is 74 Hz, which is 0.176 in the non-dimensional reduced frequency scale. The reduced buffet frequency is defined as $k_{b}=\unicode[STIX]{x03C0}f_{b}c/U_{\infty }$ . $f_{b}$ is the buffet frequency; $c$ indicates the chord of the airfoil; and $U_{\infty }$ denotes the velocity of the free stream.

Figure 2.

Buffet onset boundaries calculated with different CFD codes at Mach numbers from 0.7 to 0.8, and $Re=3\times 10^{6}$ .

The computational hybrid-unstructured flow solver has 25, 361 nodes and 40 layers of structured viscous grid around the airfoil. The distance between the first layer and the wall in the perpendicular direction is $5\times 10^{-6}$ chords ( $y^{+}\sim 1$ ). The time step for the CFD simulation is $2.94\times 10^{-4}$  s. Figure 2 shows the comparison of the transonic buffet onset boundary for the NACA0012 airfoil. The circles are from the experiment of Doerffer et al. (Reference Doerffer, Hirsch and Dussauge2010) and the clean solid line is from the URANS method. It also shows the results from Soda & Voss (Reference Soda and Voss2005), in which the boundaries are larger than those of the experimental data. The experiment and URANS calculations are in very good agreement. At $M=0.7$ , the calculated buffet onset angle is $4.80^{\circ }$ , which is very close to the $4.74^{\circ }$ tested by Doerffer et al. (Reference Doerffer, Hirsch and Dussauge2010). Time histories of the lift coefficient at $4.6^{\circ }$ , $4.8^{\circ }$ and $5.5^{\circ }$ are shown in figure 3. It can be seen that the dramatic unsteady loads occur at $5.5^{\circ }$ . Therefore, further investigation of the lock-in phenomenon is performed at $M=0.7$ and $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ . Figure 4 shows reduced buffet frequencies at different angles of attack at $M=0.7$ . As the angle of attack increases, the buffet frequency increases smoothly and slightly. At $\unicode[STIX]{x1D6FC}=5^{\circ }$ , the reduced frequency is 0.180, which is very close to the experimental data. For the scale of Strouhal number, $St=f_{b}c/U_{\infty }$ , it is approximately 0.058. McDevitt & Okuno (Reference McDevitt and Okuno1985) obtained $St=0.045{-}0.070$ on the NACA0012 airfoil. These results indicate that the numerical method in this paper is accurate for the simulation of transonic buffet. The convergence of the grid and the time step are presented in appendix B. For more validation of the URANS method, one can refer to our earlier studies (Gao et al. Reference Gao, Zhang, Liu, Ye and Jiang2015; Zhang et al. Reference Zhang, Gao, Liu, Ye and Jiang2015b ).

The grid deformation method is important to the coupled CFD/CSD simulation. The validation of this method and the numerical algorithm are tested on a forced harmonic pitching NACA 0012 airfoil from the CT5 case by Landon (Reference Landon1982). Details of the case and validated results are both presented in our research in Gao et al. (Reference Gao, Zhang, Liu, Ye and Jiang2015) and Zhang et al. (Reference Zhang, Gao, Liu, Ye and Jiang2015b ). The method is appropriate for the present investigation on the frequency lock-in.

Figure 3.

The lift and moment responses at $\unicode[STIX]{x1D6FC}=4.6^{\circ }$ , $\unicode[STIX]{x1D6FC}=4.8^{\circ }$ and $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ at $M=0.7$ and $Re=3\times 10^{6}$ .

Figure 4.

Buffet frequencies at different angles of attack, $M=0.7$ and $Re=3\times 10^{6}$ .

4.1 Reduced-order modelling by system identification

The application of the ROM to unsteady aerodynamics has been an active area of research in recent years (Lucia, Beran & Silva Reference Lucia, Beran and Silva2004; Ghoreyshi, Jirasek & Cummings Reference Ghoreyshi, Jirasek and Cummings2014). It is desirable to use a relatively simple but accurate approximation for unsteady aerodynamics by a ROM with a small number of orders (typically less than one hundred) rather than the full-order CFD model. Popular ROMs are mainly based on proper orthogonal decomposition (POD) and the system identification method. They have been successfully performed to study thrust modelling for a flapping airfoil (Li et al. Reference Li, Liu, Kou and Zhang2017), flutter in transonic flow (Hall, Thomas & Dowell Reference Hall, Thomas and Dowell2000; Thomas et al. Reference Thomas, Dowell and Hall2006; Zhang, Chen & Ye Reference Zhang, Chen and Ye2015a ), aeroservoelasticity (Zhang & Ye Reference Zhang and Ye2007; Zhang, Ye & Zhang Reference Zhang, Ye and Zhang2009; Huang et al. Reference Huang, Li, Hu and Zhao2015) and unstable flow control (Barbagallo, Sipp & Schmid Reference Barbagallo, Sipp and Schmid2009; Juillet, Schmid & Huerre Reference Juillet, Schmid and Huerre2013). To solve the above-mentioned stability problems, one needs to establish a dynamic linear model under a slight structural disturbance even if the flow contains some nonlinear features, such as the shock wave and separation. However, for the limit-cycle oscillation (LCO) in transonic flow, if larger shock motions are involved, the nonlinearity becomes the dominant characteristic. Zhang, Wang & Ye (Reference Zhang, Wang and Ye2012) established a nonlinear unsteady aerodynamic ROM using the RBF neural network, and accurately predicted LCOs with large shock motions in transonic flow. It became a hot topic to model nonlinear and unsteady aerodynamics (He, Yang & Gu Reference He, Yang and Gu2014; Mannarino & Mantegazza Reference Mannarino and Mantegazza2014; Huang et al. Reference Huang, Li, Hu and Zhao2015). Very recently, Gao, Zhang & Ye (Reference Gao, Zhang and Ye2016b ) studied the transonic buzz based on the ROM in pre-buffet flow. However, these ROMs were conducted by the forced perturbation based on stable flows. To the authors’ knowledge, they have not been documented in unstable transonic buffet flows.

We are interested in constructing a linear model to analyse the underlying mechanism of the lock-in phenomenon rather than a nonlinear model just to reproduce complex responses. This linear model is constructed based on the system identification method. The transonic buffet flow discussed in the present study is unstable, and this makes it a challenging construction for two reasons:

1. (i) The linear ROM must be set up based on the unstable steady-state base flow. However, it is a challenge to obtain the unstable steady solution of the transonic buffet flow.

2. (ii) In the training process, a competition exists between the unstable buffet flow itself and the unsteadiness caused by the forced pitching airfoil. Therefore, the other challenge is to design an appropriate forcing signal to avoid the training behaving nonlinearly.

The ROM on the basis of ARX (Auto Regressive with eXogenous input) is introduced in detail in our previous study (Zhang & Ye Reference Zhang and Ye2007), and it has been successfully performed in the investigation of VIV at supercritical Reynolds numbers (Zhang et al. Reference Zhang, Li, Ye and Jiang2015c ). For the unstable flow, the modelling process contains three steps – obtaining the unstable steady base flow, training and then performing identification. In this paper, the modelling of the transonic buffet flow is also based on the ARX model. Given that the unsteady loads are computed in a discrete domain, the single-input–single-output model can be expressed as follows:

where $y_{a}$ is the system output (pitching moment coefficient $C_{m}$ ) and $u$ is the system input (pitching angle $\unicode[STIX]{x1D703}$ ). $A_{i}$ and $B_{i}$ are constant coefficients to be estimated. Orders of the model chosen by the user are $na$ and $nb$ . The least-squares method is used to estimate unknown model parameters. To guarantee the mean data are zero, the constant levels need to be removed from the initial data before they are estimated.

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

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

where:

To couple with the structural motion equation, the discrete-time state-space form equation is turned into the continuous-time form by bilinear transformation, and the model in the state-space form is constructed as (4.5):

As we construct the state-space model of the unsteady aerodynamics, the instability problem of the flow is converted into the analysis of eigenvalues of the matrix $\unicode[STIX]{x1D63C}_{a}$ in (4.5).

4.2 Validation of the reduced-order model

As discussed in § 3.1, there are two difficulties in constructing a linear ROM with external input in the unstable transonic buffet flow. Firstly, the ROM must be trained by appropriate signals based on the steady-state base flow. Different from the time-averaged flow, the steady-state base flow strictly satisfies the governing equations and boundary conditions in the mathematical form. It plays an important role in the modelling of unsteady flow and stability analysis (Dowell & Hall Reference Dowell and Hall2001; Barkley Reference Barkley2006; Illingworth, Morgans & Rowley Reference Illingworth, Morgans and Rowley2012). For the symmetric cylinder wake flow, we can obtain the base flow by imposing a symmetric boundary condition or virtual damping (Zhang et al. Reference Zhang, Li, Ye and Jiang2015c ). While for transonic buffet flow, the approach is to stabilize the unsteadiness by some control methods, which should not cause extra changes on the initial flow conditions such as the angle of attack, airfoil shape and so on. In our recent work (Gao et al. Reference Gao, Zhang and Ye2016a ), the unstable steady-state solution of buffet flow is obtained by a closed-loop active control method. Under the optimal control law, the unsteadiness of the buffet flow is stabilized. The stabilized flow is the steady base flow for the given buffet state, which is used as the initial flow in the training process. Figure 5 shows pressure contours and streamlines of the time-averaged flow field and steady-state base flow field at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ . By comparing pressure contours, significant differences are observed between the two flows. For the time-averaged flow, the shock wave disappears; while for the steady-state base flow, it remains at 20 % after the leading edge of the airfoil. This directly leads to the difference in the pressure coefficient, as shown in figure 6, and consequently to the difference in the aerodynamic force, as shown in figure 7(a), the pitching moment coefficient of the unstable steady base flow being slightly larger than that of the time-averaged one. Figure 7(a) also presents the evolution of the buffet flow at $M=0.7$ and $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ based on the steady-state solution, and in figure 7(b), the response is plotted in the logarithm scale. When the non-dimensional time is less than 140, the amplitude of the moment coefficient is less than 0.015. At this stage, the evolution develops in an approximately linear way, which we define as the linear stage. When the non-dimensional time is larger than 140, the growing amplitude of the unstable aerodynamic force will ultimately give rise to a nonlinear, limit-cycling behaviour, which is certainly not linear. So the training process must be performed at the linear stage.

Figure 5.

Pressure contours and streamlines of flow fields at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ : (a) time-averaged flow; and (b) steady-state base flow.

Figure 6.

Comparison of pressure coefficients between the time-averaged flow and the steady-state base flow.

Figure 7.

Evolution of the pitching moment coefficients based on the steady-state base flow are plotted in (a) normal scale and (b) logarithm scale, where $C_{m0}$ indicates the mean pitching moment coefficient.

As to the second difficulty, the training signal must satisfy the following three requirements: (i) the frequency of the training signal should cover the buffet frequency; (ii) the amplitude of the training signal should be small, otherwise it will induce nonlinear responses; (iii) the training signal should not be too long and the training process must be finished before the amplitude of the moment reaches 0.015. As shown in figure 8(a), a chirp signal with an increasing frequency is designed as the input. Figure 8(b) shows the power spectrum density (PSD) analysis of the signal. Its dominant reduced frequencies are from 0.1 to 0.5, including the buffet reduced frequency. In addition, to reduce the impact of the airfoil sudden movement in the flow, a smoothing method with an average weight coefficient is used in the initial period of non-zero signals.

Figure 8.

(a) Time history and (b) PSD analysis of the training signal.

We set up the aerodynamic model at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ to verify the ROM. In order to accurately predict the unsteady aerodynamics, delay orders should be assigned with relatively large numbers, i.e. $na=nb=80$ . The identified results are compared with those of CFD simulations in figure 9; and a good agreement is observed between them. Therefore, the identification method based on the ARX model has a high numerical accuracy.

Figure 9.

Identified results compared with those of CFD simulations.

The ROM is then validated in the time domain by comparing with CFD simulations under the excitation of harmonic signals. The harmonic forcing signal is in the form of $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{am}\sin (2\unicode[STIX]{x03C0}nf_{b}t)$ . $\unicode[STIX]{x1D703}_{am}$ is the vibration amplitude of the airfoil; $f_{b}$ is the frequency of buffet; and $n$ is the non-dimensional coefficient. Figure 10 shows comparison of time histories of the pitching moment coefficient between ROM and CFD simulation at two typical forcing frequencies ( $n=0.7$ and $n=1.3$ ). ROM predictions exactly match those computed by the CFD simulation when the amplitude of the pitching moment coefficient is less than 0.015.

Figure 10.

Comparisons of time responses of a harmonic pitching airfoil at $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ for (a) $\unicode[STIX]{x1D703}=0.02\sin (1.4\unicode[STIX]{x03C0}f_{b}t)$ and (b) $\unicode[STIX]{x1D703}=0.02\sin (2.6\unicode[STIX]{x03C0}f_{b}t)$ .

By solving the eigenvalue of matrix $\unicode[STIX]{x1D63C}_{a}$ in (4.5), we can obtain the least stable fluid mode, which is associated with the instability of the transonic buffet flow. The real part of the eigenvalue indicates the damping of the fluid mode, and a positive one means that the fluid is unstable. The imaginary part of the eigenvalue, in the reduced frequency scale, indicates the reduced frequency of the fluid mode. Figure 11 displays the trajectory of the least stable eigenvalue at $M=0.7$ , which changes smoothly as the angle of attack increases. The flow is globally stable (eigenvalues lie in the left half-plane) at small angles of attack up to $\unicode[STIX]{x1D6FC}=4.70^{\circ }$ . Then, the eigenvalue crosses the imaginary axis between $\unicode[STIX]{x1D6FC}=4.70^{\circ }$ and $\unicode[STIX]{x1D6FC}=4.80^{\circ }$ . Note that the critical angle of the buffet onset obtained here corresponds to the angle computed by the CFD simulation, and it approximates to the $\unicode[STIX]{x1D6FC}=4.74^{\circ }$ observed in the experiment by Doerffer et al. (Reference Doerffer, Hirsch and Dussauge2010). At $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ , the buffet frequency from the present ROM is 0.196 and it is 0.198 from the CFD simulation (figure 4). This reveals that the linear buffet frequency by the ROM is almost the same as the nonlinear one by the CFD. Therefore, we will not differentiate them in the ensuing discussion, and they are both donated by $k_{b}$ . This is different from the cylinder wake flow, in which the linear frequency is lower than the nonlinear frequency (Zhang et al. Reference Zhang, Li, Ye and Jiang2015c ), and in which the larger the Reynolds number is, the bigger the difference is (Sipp et al. Reference Sipp, Marquet, Meliga and Barbagallo2010).

Figure 11.

Least stable eigenvalues loci obtained by the ROM at various angles of attack at $M=0.7$ .

5.1 Coupled CFD/CSD simulation method

The generalized structural motion equation in matrix form is as follows:

where $\unicode[STIX]{x1D743}$ is the structural generalized displacement; $\unicode[STIX]{x1D648}$ , $\unicode[STIX]{x1D642}$ and $\unicode[STIX]{x1D646}$ are the generalized mass, damping and stiffness matrices, respectively. $\boldsymbol{Q}$ is the generalized aerodynamic force provided by the CFD.

The structural motion equation can be non-dimensionalized by the time $\text{d}t=U_{\infty }/b\,\text{d}t_{physics}$ . For the dynamics of the SDOF pitching airfoil with zero damping, equation (5.1) can be written as:

where $\unicode[STIX]{x1D703}$ represents the airfoil pitching angle. $k_{s}=\unicode[STIX]{x1D714}_{s}b/U_{\infty }$ is the reduced frequency non-dimensionalized by the half-aerodynamic chord $b$ , natural frequency of the pitching spring $\unicode[STIX]{x1D714}_{s}$ and the free-stream velocity $U_{\infty }$ . Therefore, $k_{s}$ represents the natural frequency of the elastic airfoil in the ensuing discussion. $\unicode[STIX]{x1D707}=m/\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}b^{2}$ is the non-dimensional mass ratio, where $m$ and $\unicode[STIX]{x1D70C}$ represent the mass of the wing with unit span and the density of the free stream, respectively. $r_{\unicode[STIX]{x1D6FC}}$ is the gyration radius of the airfoil around the elastic axis at $a=0.224c$ . $C_{m}$ is the pitching moment coefficient. A sketch diagram of the SDOF pitching airfoil is shown in figure 12.

Figure 12.

Sketch diagram of the SDOF pitching airfoil system.

By defining the structural state vector $\boldsymbol{x}=[\unicode[STIX]{x1D703},\dot{\unicode[STIX]{x1D703}}]^{\text{T}}$ , the structural motion equation in the state-space form is as follows:

where $\unicode[STIX]{x1D63C}=\left[\begin{smallmatrix}\mathbf{0} & \unicode[STIX]{x1D644}\\ -\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D646} & -\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D642}\end{smallmatrix}\right]$ ; $\unicode[STIX]{x1D63D}=\left[\begin{smallmatrix}\mathbf{0}\\ \unicode[STIX]{x1D648}^{-1}\end{smallmatrix}\right]$ ; $\mathbf{0}$ is the zero matrix; $\unicode[STIX]{x1D644}$ is the identity matrix; and $Q$ is a function of $C_{m}$ , defined by $Q(C_{m})=2C_{m}/(\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}r_{\unicode[STIX]{x1D6FC}}^{2})$ . Equation (5.3) is a first-order differential equation system. A fourth-order accuracy hybrid linear multi-step method (Zhang, Jiang & Ye Reference Zhang, Jiang and Ye2007) is used to solve the aeroelastic equation (5.3) in the time domain.

Lack of a standard case for the SDOF flutter, the above-coupled CFD/CSD simulation method is validated by classical flutter cases with two degrees of freedom. In our previous work, flutter boundaries of the bench active control technology (BACT) model (Zhang et al. Reference Zhang, Gao, Liu, Ye and Jiang2015b ), the Isogai wing (Zhang et al. Reference Zhang, Jiang and Ye2007) and the LCO with NACA64010 airfoil (Zhang et al. Reference Zhang, Wang and Ye2012) agreed well with experimental data and reference results.

5.2 ROM-based aeroelastic model

Defining the structural state vector $\boldsymbol{x}_{s}=[\unicode[STIX]{x1D703},\dot{\unicode[STIX]{x1D703}}]^{\text{T}}$ for the present SDOF pitching airfoil, the structural motion equation (5.1) in the state-space form can be rewritten as:

where $\unicode[STIX]{x1D63C}_{\boldsymbol{s}}=\left[\begin{array}{@{}cc@{}}0 & 1\\ -k_{s}^{2} & 0\end{array}\right]$ , $\unicode[STIX]{x1D63D}_{\boldsymbol{s}}=\left[\begin{array}{@{}c@{}}0\\ 1\end{array}\right]$ , $\unicode[STIX]{x1D63E}_{\boldsymbol{s}}=\left[\begin{array}{@{}cc@{}}1 & 0\end{array}\right]$ , $\unicode[STIX]{x1D63F}_{\boldsymbol{s}}=[0]$ , $q=1/(\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}r_{\unicode[STIX]{x1D6FC}}^{2})$ .

By defining the state vector $\boldsymbol{x}_{\boldsymbol{a}\boldsymbol{e}}=[\boldsymbol{x}_{\boldsymbol{s}},\boldsymbol{x}_{\boldsymbol{a}}]^{\text{T}}$ and coupling the structural state equations (5.4) with the aerodynamic state equations (4.5), we obtain the state equation and output equation for the linear aeroelastic system as follows:

The ROM-based aeroelastic model of a SDOF pitching airfoil system is then obtained. The problem of aeroelastic stability is converted into the analysis of complex eigenvalues of $\unicode[STIX]{x1D63C}_{\boldsymbol{a}\boldsymbol{e}}$ . The process of the transonic aeroelastic analysis is:

1. (i) Use the URANS flow solver to compute the aerodynamic coefficient of the mode with the designed input signal in figure 8.

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

3. (iii) Couple the aerodynamic state-space equation (4.5) and the structural state-space equation (5.4) to obtain the ROM-based aeroelastic state-space equation (5.5).

4. (iv) Solve the complex eigenvalues of the state matrix $\unicode[STIX]{x1D63C}_{\boldsymbol{a}\boldsymbol{e}}$ in (5.5). We focus on the two most unstable eigenvalues, which correspond to the fluid and structural modes. With different structural parameters, natural frequencies of the elastic airfoil and mass ratios, we can obtain a set of eigenvalues. It is the root loci. The stability of the aeroelastic system changes when the root loci cross the imaginary axis.

The model is then validated by the BACT case (Rivera et al. Reference Rivera, Dansberry, Bennett, Durham and Silva1992). Although this case is a two DOF system, procedures for the modelling are the same. Therefore, it is an alternative validation owing to lack of standard cases for the SDOF system. The results are shown in appendix C.

6.1 Characteristics of the frequency lock-in

Figures 13–15 show time history responses and PSD analyses at different natural frequencies ( $k_{s}$ ) of the elastic airfoil at $\unicode[STIX]{x1D707}=200$ . In the cases of $k_{s}=0.10$ (figure 13) and $k_{s}=0.60$ (figure 14), the amplitudes of the pitching angle responses are very small; and the amplitudes of the pitching moment coefficient are approximately equal to those of the rigid model (figure 7). PSD analyses show that coupling frequencies ( $k_{ae}$ ) are dominated by the buffet frequency ( $k_{b}=0.196$ ), which displays forced vibration under the unstable buffet loads. However, in the case of $k_{s}=0.35$ (figure 15), the basic coupling frequency is no longer dominated by the buffet frequency, but follows the natural frequency of the elastic airfoil. This is the frequency lock-in phenomenon. When it occurs, the amplitude of the oscillating pitching angle exceeds $8^{\circ }$ , and the amplitude of the moment coefficient is much larger than that of the rigid model. Such large-amplitude oscillations may cause serious flight accidents.

Furthermore, figure 16 shows the coupling frequency ( $k_{ae}$ ) of the aeroelastic system as a function of the natural frequency ( $k_{s}$ ). In most regions, the coupling frequency follows the buffet frequency. However, in the region of $0.21<k_{s}<0.46$ (the lower boundary is controversial, which will be discussed in next section), the coupling frequency locks onto the natural frequency. The dashed line $k_{ae}=k_{s}$ represents the lock-in phenomenon. It is found that lock-in can occur in a wider range at the ratio of natural frequency to buffet frequency ( $k_{s}/k_{b}$ ) from 1.07 to 2.34. Previous research attributed the frequency lock-in to resonance. However, the so-called ‘resonance region’ ( $0.21<k_{s}<0.46$ or $1.07<k_{s}/k_{b}<2.34$ ) does not display a symmetrical distribution against the buffet frequency (0.196), but it is entirely located in the higher frequency side. Is frequency lock-in in transonic buffeting flow really caused by resonance?

Figure 13.

Time histories at $k_{s}=0.10$ of (a) the pitching angle response; (b) the pitching moment coefficient response; and (c) the PSD analysis of the pitching moment coefficient response.

Figure 14.

Time histories at $k_{s}=0.60$ of (a) the pitching angle response; (b) the pitching moment coefficient response; and (c) the PSD analysis of the pitching moment coefficient response.

Figure 15.

Time histories at $k_{s}=0.35$ of (a) the pitching angle response; (b) the pitching moment coefficient response; and (c) the PSD analysis of the pitching moment coefficient response.

Figure 16.

Oscillating frequency of the aeroelastic system as a function of the natural frequency of the elastic airfoil.

6.2 Mechanism of the frequency lock-in

The mechanism underlying the lock-in will be discussed in this section using the ROM-based aeroelastic model. Figure 17 presents the eigenvalue loci of the coupled system varying with the natural frequency ( $k_{s}$ ) of the airfoil at $\unicode[STIX]{x1D707}=200$ . It can be found that the eigenvalue loci have two branches, which we define as branch A and branch B. For branch A, its mode is stable (eigenvalues lie in the left half-plane) at low natural frequencies; then as the natural frequency increases, it crosses into the right half-plane, and the coupling frequency of the system at the crossing point is $k_{ae}=0.196$ ; as the natural frequency increases furthermore, it approaches the uncoupled least stable fluid mode (refer to figure 11 at $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ ). However, for branch B, it is unstable at the initial stage; then as the increase of the natural frequency, it crosses into the left half-plane, and the coupling frequency at the crossing point is $k_{ae}=0.46$ .

Figure 17.

Eigenvalue loci of the aeroelastic system at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ and $\unicode[STIX]{x1D707}=200$ .

Figure 18.

Imaginary part of the eigenvalue loci as a function of the natural frequency.

Figure 18 shows the imaginary part of the eigenvalue loci in figure 17 as a function of the natural frequency $k_{s}$ . At low natural frequencies, the coupling frequency of branch A follows the natural frequency of the structural model (S mode); and the coupling frequency of branch B is dominated by the buffet frequency of the fluid mode (F mode). However, at high natural frequencies, the coupling frequency of branch A turns to be dominated by the buffet frequency; while that of branch B becomes to follow the natural frequency. The corresponding modes also exchange.

From the above analysis, it can be seen that there is mode veering between the two branches. Mode veering refers to the phenomenon associated with the imaginary part (frequency) of the eigenvalue loci of a system with a variable parameter: two branches approach each other and then veer away and diverge instead of intersecting (Perkins & Mote Reference Perkins and Mote1986; Chen & Kareem Reference Chen and Kareem2003; Vijayann & Woodhouse Reference Vijayann and Woodhouse2014). When mode veering occurs, eigenmodes (eigenvectors) associated with these two branches also exchange in a continuous way. Meliga & Chomaz (Reference Meliga and Chomaz2011) first found the veering phenomenon between F mode and S mode, which was defined as mode transition in their study, when they investigated the effect of mass ratio on the instability of a spring mounted circular cylinder. It is found that if the mass ratio is larger than the critical mass ratio, the instability is led by F mode; if it is smaller, the instability is led by S mode; and if it is close to the critical mass ratio, no clear distinction can be found between the two modes.

Figure 19 shows the close-up view of the veering region, as well as the coupling frequency of the aeroelastic system from the coupled CFD/CSD simulation which corresponds to that of branch B. As the natural frequency ( $k_{s}$ ) of the elastic airfoil approaches to the buffet frequency ( $k_{b}$ ), the loci of branch A and branch B tend to veer and march along each other’s previous locus (dashed and solid lines in figure 19), thus avoiding an intersection; and the modes associated with these two frequencies exchange. We also plot the directions (the real part and the imaginary part in complex plane) of the complex eigenvectors. As shown in figure 19, only components of interest are displayed, black and blue arrows standing for branch A and branch B, respectively. The directions also exchange after the mode veering. Therefore, they can be used to track the mode veering region. A common approximate law is defined as: if the direction of an eigenvector displays a deflection more than $10^{\circ }$ compared with that of the initial (last) one, we define it as the onset (offset) of the veering region. Based on this law, the mode veering region (the grey area), $0.168<k_{s}<0.232$ or $0.83<k_{s}/k_{b}<1.18$ , can be marked in figure 19. Therefore, we can draw a panorama of the mode veering associated with branch B with the increasing natural frequency. Before the mode veering, $k_{\boldsymbol{s}}<0.168$ , branch B corresponds to F mode; while after that, $k_{\boldsymbol{s}}>0.232$ , it turns to correspond to S mode. In the veering region, $0.168>k_{s}>0.232$ , the mode property of branch B transforms from F to S. Due to the continuity of the transformation, it is difficult to distinguish S mode from F mode in the region.

Figure 19.

Mode veering of imaginary part (frequency) of the eigenvalue loci with the varying natural frequency of the elastic airfoil. The blue triangle symbol refers to branch B, and the black square symbol refers to branch A. The dashed line denotes the buffet frequency of the rigid airfoil, and the solid line denotes the natural frequency. Blue arrows denote directions of eigenvectors associated with branch B, and black arrows denote those of branch A. The grey area refers to the mode veering region based on the approximate law.

Figure 20 shows the real part (represents the damping characteristics) and imaginary part (represents the frequency characteristics) of the eigenvalue loci in figure 17, as well as the coupling frequency and the amplitude of the pitching angle from the coupled CFD/CSD simulation as functions of the natural frequency of the elastic airfoil. Since the real part and the imaginary part share the same veering range, the real part also displays the mode veering phenomenon (figure 20 b), which has been discussed in the imaginary part of the loci in figure 19. Before the veering ( $k_{s}<0.168$ ), branch B stands for F mode, but after that ( $k_{s}>0.232$ ), it turns to stand for S mode. From figure 20(b), it can also be seen that the damping (real part) of S mode is positive at $0.232<k_{s}<0.460$ . Based on the linear eigenvalue theory, a positive damping means the instability of S mode, namely, flutter in aeroelasticity. The upper flutter boundary (the crossing point) is at $k_{s}=0.460$ , which is exactly equal to the upper boundary of the lock-in region from the coupled CFD/CSD simulation.

Figure 20.

(a) Imaginary part of the eigenvalue loci, (b) real part of the eigenvalue loci from the ROM-based aeroelastic model; and (c) coupling frequency of the system, (d) oscillating amplitude of the pitching airfoil from the coupled CFD/CSD simulation varying with the natural frequency of the elastic airfoil.

Figure 21.

Comparison of the logarithmic increment between the coupled CFD/CSD simulation and the ROM-based aeroelastic model as a function of the natural frequency of the elastic airfoil.

Figure 22.

Time history of the pitching angle of the airfoil from the coupled CFD/CSD simulation and the sample for the calculation of the logarithmic increment at (a) $k_{s}=0.25$ , (b) $k_{s}=0.28$ and (c) $k_{s}=0.42$ .

Figure 21 shows the comparison of the logarithmic increment of the pitching angle between the coupled CFD/CSD simulation and the ROM-based aeroelastic model. The logarithmic increment from the ROM-based model is transformed from the real part (positive damping) of branch B in figure 20(b). In the coupled CFD/CSD simulation, it is calculated from the time history of the pitching angle at the initial stage with approximately exponential oscillation. It is defined as $\unicode[STIX]{x1D702}=Ln(Am^{j+1}(\unicode[STIX]{x1D703})-Am^{j}(\unicode[STIX]{x1D703}))/Am^{j}(\unicode[STIX]{x1D703})$ , where $Am^{j}(\unicode[STIX]{x1D703})$ represents the amplitude of the $j$ th cycle. Responses and samples for the calculation of the logarithmic increment at typical natural frequencies are shown in figure 22. It can be seen that the logarithmic increment from the coupled CFD/CSD simulation can reasonably match that predicted by the linear model. As shown in figure 22(a), at $k_{s}=0.25$ , the oscillation reaches a LCO rapidly because of the higher linear increment (positive damping). In this case, the nonlinearity dominates the system. The predicted increment from the sample of the simulation is also dominated by the nonlinearity. Therefore, errors at cases of $k_{s}<0.25$ are reasonable. In the case of $k_{s}=0.42$ , the time history of the pitching angle displays an exponential growth. Both logarithmic increments are in good agreement. Therefore, from the results in figures 19 and 21, both are convincing, the characteristics of the frequency and the damping in the lock-in region predicted by the linear ROM-based aeroelastic model. Besides, the frequency lock-in is dominated by linear dynamics.

In order to thoroughly understand the dynamics of the buffeting response, the natural frequency ( $k_{s}$ ) of the airfoil is divided into four sub-regions, as shown in figure 20(c). They are divided by three critical natural frequencies, namely, the lower and upper boundaries ( $k_{s}^{l}=0.168$ , $k_{s}^{m}=0.232$ ) of the frequency loci veering, and the upper flutter boundary ( $k_{s}^{u}=0.460$ ). In sub-regions 1 ( $k_{s}<0.168$ ) and 4 ( $k_{s}>0.460$ ), coupling frequencies follow the buffet frequency, namely, $k_{ae}=k_{b}$ . In fact, they are both forced vibration caused by the unsteady buffeting loads. Vibration amplitudes of the airfoil are comparatively small in these two sub-regions. In sub-region 2 ( $0.168<k_{s}<0.232$ ), namely, the mode veering region, the coupling frequency neither follows the buffet frequency nor the natural frequency. The oscillating amplitude of the airfoil in this sub-region is between that in sub-regions 1 and 3. It is a transitional stage. However, due to the continuity of the veering process, it is hard to identify which mode, S mode or F mode, dominates the system instability. In sub-region 3 ( $0.232<k_{s}<0.460$ ), the coupling frequency synchronizes with the natural frequency ( $k_{ae}=k_{s}$ ) – the phenomenon of the frequency lock-in, and the airfoil vibrates with an increasing amplitude. The peak oscillation amplitude is achieved at a distinctly high frequency ratio of $k_{s}/k_{b}=1.73$ rather than $k_{s}/k_{b}=1.0$ .

Results from the above discussion indicate that the coupling between the S mode and F mode leads to the instability of the former. Besides, the instability range coincides with the lock-in region obtained from the coupled CFD/CSD simulation. It reveals that the frequency lock-in is caused by the linear coupled-mode flutter. However, the flutter in this study is not a classical mode coupling between structural modes but a coupling between one structural mode and one fluid mode. Therefore, the physical mechanism underlying the frequency lock-in phenomenon in transonic buffeting flow is SDOF flutter. But different from the classical SDOF flutter (e.g. transonic buzz) which occurs from an initial stable flow, the present one is in the unstable buffeting flow. Its particularity is the competition between two unstable modes. Therefore, the response of the airfoil system undergoes a conversion from the forced vibration to the self-sustained oscillation (flutter). The coupling frequency certainly should lock onto the natural frequency of the elastic airfoil. The occurrence of the frequency lock-in is dominated by the linear dynamics.

The conversion is then validated by the coupled CFD/CSD simulation at $\unicode[STIX]{x1D707}=200$ and $k_{s}=0.3$ with zero initial condition of the displacement and velocity. Based on the analysis of the linear dynamics, there are two unstable modes in this case, F mode and S mode. The complex response can be explained by the competition between the two unstable modes. Figure 23 shows the time history and the frequency content of the responses. The response history can be divided into three stages. In the first stage, $t<200$ , the pitching angle is very small and the pitching moment coefficient vibrates with the amplitude approximately equal to that of the rigid airfoil. The coupling frequency is followed by the buffet frequency ( $k_{b}=0.196$ ) (figure 23 c). Therefore, the system is dominated by the unstable F mode, displaying forced vibration under the unstable buffet loads in this stage. In the second stage, $200<t<440$ , there are two peak frequency components, the buffet frequency with less power, and the natural frequency with more power (figure 23 d). It means that, as the response time goes, the buffet frequency becomes weak, while the natural frequency becomes strong. That is, the unstable S mode gradually dominates the dynamic characteristics of the system. Finally, in the last stage, $t>440$ , the responses diverge. The oscillating frequency absolutely follows the natural frequency (figure 23 e). The frequency lock-in phenomenon occurs. The unstable S mode completely dominates the characteristics of the system. Therefore, the system in this stage displays self-sustained oscillations (flutter).

Figure 23.

Time history of (a) the pitching angle, (b) the pitching moment coefficient and (c–e) frequency contents of the aeroelastic system at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ and $\unicode[STIX]{x1D707}=200$ , $k_{s}=0.3$ with zero initial conditions.

With time going on ( $t>800$ ), responses of the pitching angle and the moment coefficient increase rapidly and eventually approach LCOs. Owing to a limited time range, LCO is not shown in figure 23, but it is exhibited in figure 15. It can be seen that the pitching moment does not display a simple harmonic response (figure 15 b). The PSD result shows two peaks (figure 15 c). The basic peak (0.35) is accurately located at the structural reduced frequency, and the second one equals twice that of the structural reduced frequency. This is a nonlinear response. We notice that the airfoil pitching angle exceeds $8^{\circ }$ (figure 15 a). Such a large-amplitude pitching motion certainly leads to a nonlinear flow, and consequently results in a double frequency response of the pitching moment. However, we want to specifically emphasize that the nonlinearity dominates the system only when the structural pitching angle is large enough. At the initial stage of the response with a small pitching angle, the system is dominated by the instability of the structural mode (flutter). In transonic flow, a divergent response (flutter) generally reaches a LCO finally, as a result of the shock wave. Therefore, the nonlinear response (LCO) is just the final performance of the coupling linear flutter.

Some researchers have investigated the frequency lock-in in unstable flows with harmonic forcing, such as the transonic buffet flow by Raveh & Dowell (Reference Raveh and Dowell2011) and the open-cavity flow by Sipp (Reference Sipp2012). These studies used uncoupled methods, and in particular, the lock-in range depends on the forcing amplitude. In essence, there exists a completion between the unstable flow and the unsteadiness caused by the forcing. As long as the forcing amplitude is large enough, the unsteadiness of the forcing flow will dominate the entire characteristics of the system.

Figure 24.

(a) Real part and (b) imaginary part of the eigenvalue loci at $\unicode[STIX]{x1D707}=60$ as functions of the natural frequency of the elastic airfoil.

6.3 Effect of the mass ratio

The investigation in §§ 6.1 and 6.2 is performed at a fixed mass ratio $\unicode[STIX]{x1D707}=200$ . This section will focus on the effect of the mass ratio on the dynamics of the buffeting response and frequency lock-in.

Figures 24 and 25 show the complex eigenvalue loci as functions of the natural frequency ( $k_{s}$ ) of the elastic airfoil at $\unicode[STIX]{x1D707}=60$ and $\unicode[STIX]{x1D707}=1000$ , respectively. Generally, the analysis procedure is similar to that of the eigenvalue loci at $\unicode[STIX]{x1D707}=200$ . We first recognize properties of the two branches and boundaries of the veering region by directions of eigenvectors based on the imaginary part of the eigenvalue loci (figures 24 b and 25 b). Then, we determine the upper flutter boundary (the crossing point of branch B) in the real part of the eigenvalue loci (figures 24 a and 25 a). Therefore, we can get three critical frequencies, the lower and upper boundaries ( $k_{s}^{l}$ and $k_{s}^{m}$ ) of the veering region and the upper flutter boundary ( $k_{s}^{u}$ ). In the case of $\unicode[STIX]{x1D707}=60$ , $k_{s}^{l}=0.126$ , $k_{s}^{m}=0.253$ and $k_{s}^{u}=0.460$ ; while in the case of $\unicode[STIX]{x1D707}=1000$ , $k_{s}^{l}=0.190$ , $k_{s}^{m}=0.210$ and $k_{s}^{u}=0.460$ . According to the three critical points, the natural frequency can be also divided into four sub-regions, namely, the forced vibration regions as sub-region 1 and 4, the mode veering region as sub-region 2 and the frequency lock-in region as sub-region 3. As discussed in § 5.2, the responses are forced vibrations in sub-regions 1 and 4 with smaller oscillating amplitudes, and the coupling frequencies always follow the buffet frequency. In sub-region 3, S mode is unstable, and the unstable region coincides with the lock-in region. It provides further evidence for the fact that the frequency lock-in phenomenon in transonic buffeting flow is caused by the linear coupled flutter. Furthermore, we notice that the upper flutter boundary (the crossing point, $k_{s}^{u}$ ) is substantially fixed at $k_{s}^{u}=0.460$ at different mass ratios. This reveals that the mass ratio has a slight influence on this boundary. This is different from the conclusion in VIV (Meliga & Chomaz Reference Meliga and Chomaz2011; Zhang et al. Reference Zhang, Li, Ye and Jiang2015c ), in which the upper boundary of the lock-in is significantly affected by the mass ratio, especially when the mass ratio is rather small. In sub-region 2, veering boundaries greatly rely on the mass ratio. The veering range changes from (0.126, 0.253) to (0.190, 0.210) as the mass ratio increases from $\unicode[STIX]{x1D707}=60$ to $\unicode[STIX]{x1D707}=1000$ . Thus, the larger the mass ratio is, the narrower the veering region is and the closer veering boundaries are to the buffet frequency.

Figure 25.

(a) Real part and (b) imaginary part of the eigenvalue loci at $\unicode[STIX]{x1D707}=1000$ as functions of the natural frequency of the elastic airfoil.

Figure 26 shows the comparison of the coupling frequency ( $k_{ae}$ ) between the coupled CFD/CSD simulation and the ROM-based aeroelastic model at both mass ratios. The coupling frequency obtained by the simulation reasonably matches that of the model. Especially, the linear ROM-based aeroelastic model can accurately predict the upper flutter boundary. Therefore, it further proves that the frequency lock-in in transonic buffeting flow is dominated by linear dynamics. However, the characteristics of veering, to a certain extent, are associated with nonlinear dynamics when the mass ratio is small. In this region, the interaction between fluid and structure is strong, and the linear stage is very short, which is believed to be responsible for the gap between the linear model and the nonlinear coupled CFD/CSD simulation at $\unicode[STIX]{x1D707}=60$ . A full weakly nonlinear model, such as that reported by Meliga & Chomaz (Reference Meliga and Chomaz2011), may help to analyse cases with very small mass ratios.

Figure 26.

Comparison of the coupling frequency between the ROM-based aeroelastic model and CFD/CSD simulation at (a) $\unicode[STIX]{x1D707}=60$ and (b) $\unicode[STIX]{x1D707}=1000$ .

Figure 27.

Eigenvalue loci of the coupled system at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ and $\unicode[STIX]{x1D707}=5000$ .

We further enlarge the mass ratio to $\unicode[STIX]{x1D707}=5000$ . The eigenvalue loci are shown in figure 27. It can be seen that the loci display two separated and clean clusters. One fluctuates in the vicinity of the least stable F mode (figure 11), and forms a circle with the increasing natural frequency ( $k_{s}$ ) of the elastic airfoil. The other fluctuates in the vicinity of the imaginary axis, first crossing the axis from the left half-plane to the right half-plane at $k_{s}=0.196$ and then crossing from the right to the left at $k_{s}=0.460$ . This is different from the eigenvalue loci in figure 17. In this case, $\unicode[STIX]{x1D707}=5000$ , the branch in the vicinity of the least stable F mode is a coupled F mode. And the other branch near the imaginary axis is a coupled S mode. Figure 28 shows the real and imaginary parts of complex eigenvalue as functions of the natural frequency of the elastic airfoil. From the imaginary part in figure 28(b), it can be seen that the coupled-mode frequencies always follow their uncoupled frequencies. When $k_{s}=0.196$ , the two loci intersect but do not repeal. This is the phenomenon of loci crossing. Different from the loci veering, eigenmodes associated with the two branches do not exchange in loci crossing. As shown in figure 28(a), the damping of S mode is positive when the natural frequency is within the range of $0.196<k_{s}<0.460$ , indicating the instability of S mode (flutter). This range is also the lock-in region from the coupled CFD/CSD simulation. This further verifies that the frequency lock-in in transonic buffeting flow is caused by the linear coupled-mode flutter.

Figure 28.

(a) Real part of the eigenvalue loci and (b) imaginary part of the eigenvalue loci as functions of the natural frequency of the elastic airfoil at $M=0.7$ , $\unicode[STIX]{x1D6FC}=5.5^{\circ }$ and $\unicode[STIX]{x1D707}=5000$ .

It can be seen that with the increase of the mass ratio, the imaginary part (frequency) of the eigenvalue loci changes from veering to crossing, which suggests the existence of a critical mass ratio. It is closely related to characteristics of the least stable F mode. That is, free-stream parameters, namely, Mach number, angle of attack and Reynolds number determine the critical mass ratio. In the present study, we calculate that the critical mass ratio is $\unicode[STIX]{x1D707}=3350$ . Figure 29 shows the mode veering region and the frequency lock-in region at different mass ratios. Good agreement is observed between boundaries calculated by the coupled CFD/CSD simulation and those obtained by the ROM-based aeroelastic model. As the mass ratio increases, the mode veering region becomes smaller. At $\unicode[STIX]{x1D707}=3350$ , the mode veering region vanishes and the lower boundary of the lock-in region equals the buffet frequency.

Figure 29.

Mode veering and frequency lock-in boundaries as functions of the mass ratio.