1. Introduction
Large-scale buffeting flow can result in excessive lift fluctuations on aircraft operating at higher lift coefficients, thus limiting the allowable flight envelope. In transonic flow conditions, this is associated with shock oscillations that are synchronized to oscillations of the downstream shear layer. In order to avoid this buffeting condition, a buffet-onset boundary must be determined. This is typically done by slowly increasing the lift coefficient (e.g. pulling a windup turn) until the airplane experiences unsteady normal accelerations. The buffet-onset boundary as a function of Mach number is then set by the airplane lift level where the peak-to-peak accelerations exceed some threshold value.
This sequence of increasing angle of attack leading to increasing lift coefficient, ultimately resulting in the onset of large-scale unsteadiness, is also observed on airfoils. As a result, the transonic buffet on airfoils has been used as a model problem to investigate basic mechanisms (McDevitt & Okuno Reference McDevitt and Okuno1985; Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009) – although swept-wing buffet is more complex and also involves additional instabilities. Transonic buffet on airfoils has been shown to result from a Hopf bifurcation, leading to an oscillatory flow instability (Crouch, Garbaruk & Magidov Reference Crouch, Garbaruk and Magidov2007; Crouch et al. Reference Crouch, Garbaruk, Magidov and Travin2009a,Reference Crouch, Garbaruk, Magidov, Jacquin and Brazab; Iorio, Gonzalez & Ferrer Reference Iorio, Gonzalez and Ferrer2014; Sartor, Mettot & Sipp Reference Sartor, Mettot and Sipp2015). The initial onset and flow structure of the instability is well predicted by the stability analysis, when compared with both experiments (Crouch et al. Reference Crouch, Garbaruk, Magidov and Travin2009a,Reference Crouch, Garbaruk, Magidov, Jacquin and Brazab; Crouch, Garbaruk & Strelets Reference Crouch, Garbaruk and Strelets2019; Fukumoto et al. Reference Fukumoto, Kouchi, Sugioka and Koike2023) and simulations (Crouch et al. Reference Crouch, Garbaruk, Magidov and Travin2009a; Moise et al. Reference Moise, Zauner, Sandham, Timme and He2023). Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009a) also suggested that the near-critical nonlinear amplitudes are consistent with a weakly nonlinear description of a supercritical bifurcation.
Global instabilities have also been shown to give rise to large-scale unsteady flow for both infinite swept wings (Crouch et al. Reference Crouch, Garbaruk and Strelets2019; Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019) and finite-span swept wings (Timme Reference Timme2020; Plante & Laurendeau Reference Plante and Laurendeau2023; Sansica & Hashimoto Reference Sansica and Hashimoto2023). For the swept wing, there is a travelling wave mode in addition to the oscillatory mode seen on airfoils. The travelling wave mode convects outboard along the wing, with a dominant spanwise wavelength approximately equal to the wing chord (Iovnovich & Raveh Reference Iovnovich and Raveh2014; Crouch et al. Reference Crouch, Garbaruk and Strelets2019; Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019; He & Timme Reference He and Timme2021). The relative role of these modes is dependent on the configuration and the flow conditions, but in any case, some insight into the finite-amplitude behaviour is needed to assess the buffet boundary.
Earlier studies using unsteady simulations combined with proper orthogonal decomposition or dynamic mode decomposition have shown that the post-critical buffeting flow can be well represented using low-order modelling (Ohmichi, Ishida & Hashimoto Reference Ohmichi, Ishida and Hashimoto2017; Poplingher & Raveh Reference Poplingher and Raveh2018; Fukumoto et al. Reference Fukumoto, Kouchi, Sugioka and Koike2023). These studies also show that the dominant flow structure is qualitatively the same as predicted by the linear global stability analysis (GSA). Sansica et al. (Reference Sansica, Loiseau, Kanamori, Hashimoto and Robinet2022) used the sparse identification of nonlinear dynamics technique in conjunction with dynamic mode decomposition to extract a low-order representation of the entire flow field from unsteady simulations. The underlying dynamics are linked to the Landau equation, consistent with a supercritical bifurcation.
Here, we apply weakly nonlinear analysis to predict the post-critical buffeting flow as an extension of the global stability theory, focused on the oscillatory mode of a two-dimensional airfoil. The Landau equation is used to model the nonlinear amplitude variation, with the constants determined by various combinations of numerical simulation and stability analysis. The problem formulation provides a description of the numerical formulations used for the unsteady simulations and the GSA, followed by the weakly nonlinear analysis. Results are then presented for two different airfoils, followed by discussion and conclusions.
2. Problem formulation
The study is based on an airfoil at fixed chord Reynolds number 
$Re_C= U_{\infty } c / \nu$ and Mach number 
$M$. The airfoil angle of attack 
$\alpha$ is generally fixed, but may also vary slowly in time. Results are presented for the OAT15A airfoil at 
$Re_C=3 \times 10^6$ and the NACA0012 airfoil at 
$Re_C= 10^7$. The buffet instability growth rate 
$\gamma$ and angular frequency 
$\omega$ are non-dimensionalized by 
$c/U_{\infty }$, giving a Strouhal number of 
$St= f c/U_{\infty } = \omega /2 {\rm \pi}$.
The flow is governed by the unsteady Reynolds-averaged Navier–Stokes (URANS) equations, using the Spalart–Allmaras turbulence model (Spalart & Allmaras Reference Spalart and Allmaras1994) for closure of the Reynolds stresses. The turbulence model includes the compressibility correction (Spalart Reference Spalart2000), with the constant 
$C_5 = 3.5$.
The URANS equations are solved numerically using two independent in-house finite-element codes. The unsteady simulations are based on the code FELight, which is run on a fixed unstructured grid. This code follows the discretization and solution methods of the HOMA solver (Ahrabi & Mavriplis Reference Ahrabi and Mavriplis2020). The GSA uses a solution-adaptive scheme (loosely) coupling the GGNS-T1 solver with the EPIC grid adaption module (Michal et al. Reference Michal, Kamenetskiy, Krakos, Mani, Glasby, Erwin and Stefanski2018). EPIC is an anisotropic grid generation and adaptation code which takes a Riemannian metric and target grid complexity (usually, number of grid nodes or cells) and produces an unstructured grid approximately conforming to the input metric, following the paradigm of grid/metric duality (Loseille & Alauzet Reference Loseille and Alauzet2011). To be meaningful, the Riemannian metric is tied to some sort of error estimation in the localized form. For this study, we use a metric arising from the goal-oriented error-estimation technique described in detail in Kamenetskiy et al. (Reference Kamenetskiy, Krakos, Michal, Clerici, Alauzet, Loseille, Park, Wood, Balan and Galbraith2022). The adaptive grid approach is very efficient and robust for the GSA, but not suitable for the unsteady simulations.
In both the FELight and GGNS-T1 codes, the URANS equations are spatially discretized using streamline upwind/Petrov–Galerkin (SUPG), which is a stabilized finite-element discretization utilizing globally continuous finite elements (Brooks & Hughes Reference Brooks and Hughes1982). For the current study, piecewise linear elements are used. For the mixed types of cells (such as prisms, pyramids and hexes) used in the fixed-grid unsteady simulations, we use standard globally continuous finite-element spaces, which have locally poly-linear elements as well. The resulting discretization is second order for all equations, including the turbulence model equation governing 
$\tilde \nu$. Written in terms of the primitive variables 
$\boldsymbol {q} =[\rho, u, v, T, \tilde \nu ]^{\rm T}$, the semi-discretized system of ordinary differential equations has the form
\begin{equation} \boldsymbol{M}\frac{\partial \boldsymbol{q}}{\partial t}=R(\boldsymbol{q}), \end{equation}
where 
$R(\cdot )$ is the discrete residual and 
$\boldsymbol {M}$ is the finite-element mass matrix independent of 
$\boldsymbol {q}$. Other variable choices can also be used for the GSA to streamline the extraction of the eigenfunctions.
The boundary conditions are implemented either as essential (strong) conditions (e.g. non-slip viscous walls, where the velocity components and 
$\tilde \nu$ are set to zero) or as flux boundary conditions (e.g. in the far field and on free-slip inviscid walls). For this transonic flow, additional shock-capturing artificial dissipation is used, as in Holst et al. (Reference Holst, Glasby, Erwin, Stefanski, Prosser, Anderson and Wood2021).
To solve this system of equations, we employ a pseudo-transient continuation method (Kelley & Keyes Reference Kelley and Keyes1998), in which the temporal term is augmented as follows:
\begin{equation} \boldsymbol{M}\left( \frac{\partial \boldsymbol{q}}{\partial t} + \frac{\partial \boldsymbol{q}}{\partial \tau} \right) =R(\boldsymbol{q}). \end{equation}
In this equation, 
$\tau$ serves as the pseudo-time, distinct from 
$t$, which denotes physical time.
For steady-state computations, the physical time term is set to zero. The pseudo-transient term is discretized using the first-order backward difference formula. The resulting implicit system is linearized exactly and solved in a marching manner. To expedite convergence, the pseudo-time is computed for each grid node, resulting in a local time-stepping method. The local time steps are amplified by a global Courant–Friedrichs–Lewy (CFL) number to adjust the pace of nonlinear advancements. The linearized system is solved using an inexact Newton method, wherein the linear system is solved via an ILU-preconditioned GMRES method (Saad & Schultz Reference Saad and Schultz1986). Additionally, a line search process is employed to determine an optimal relaxation factor for the Newton updates, derived from the solution of the linearized system. A CFL number control mechanism is also integrated to regulate the pace of nonlinear advancement. The algorithm employed closely resembles those outlined in Kelley & Keyes (Reference Kelley and Keyes1998).
2.1. Unsteady simulation approach
The unsteady simulations are initialized with steady equilibrium states, which are solutions to the steady RANS equations. For the unsteady problem, the temporal term is discretized using the second-order backward difference formula. Treating the pseudo-time term similar to that of the steady-state scenario leads to a dual time-stepping technique (Jameson Reference Jameson1991). In this technique, the equations for each implicit time step are treated as modified steady-state problems, which are solved by iteratively progressing through the pseudo-time variable. The application of the Newton method with ILU-preconditioning and the inclusion of a line search process facilitates achieving deep levels of convergence (down to machine precision) for both steady-state and unsteady problems within each time step.
For the unsteady simulations, the initial steady state and the unsteady evolution are calculated on a fixed grid, as shown in figure 1(a). The grid has been refined in the areas of the shock and downstream shear-layer motions. The grid resolution and the physical time step are evaluated based on their impact on the buffet-instability growth rate and frequency, as discussed in the results.
Figure 1. Representative grids used for the (a) unsteady simulations and the (b) global stability calculations. Background colour for (b) shows the Mach number distribution.
2.2. Global stability analysis
For the GSA, we consider a small (linearized) perturbation around a steady-state ‘fixed-point’ solution (
$\boldsymbol {q}=\boldsymbol {q}^{(B)}+\boldsymbol {q}^\prime$), yielding
\begin{equation} \boldsymbol{M}\frac{\partial \boldsymbol{q}^\prime}{\partial t}= \boldsymbol{J} \boldsymbol{q}^\prime, \end{equation}
with 
$\boldsymbol {J}= {\partial R}/{\partial \boldsymbol {q}}|_{\boldsymbol {q}=\boldsymbol {q}^{(B)}}$. Expanding the linear perturbation in the normal-mode form 
$\boldsymbol {q}^\prime =\hat {\boldsymbol {q}} {\rm e}^{\lambda t}$, with 
$\lambda = \gamma + {\rm i} \omega$ gives
\begin{equation} \boldsymbol{M}^{{-}1} \boldsymbol{J} \hat{\boldsymbol{q}}= \lambda \hat{\boldsymbol{q}}. \end{equation}
Exact Jacobians 
$\boldsymbol {J}$ are obtained via the automatic differentiation approach, which provides an exact linearization at the converged steady state. This is essential for the current ‘discretize-then-linearize’ approach for GSA – in contrast to the ‘linearize-then-discretize’ approach of Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009a).
For the stability analysis, both the initial steady state and the global stability equations are solved using an adaptive grid with the drag coefficient chosen as the target functional output. In the process of adaptation, the curved geometry is resolved – that is, new grid nodes are projected to the exact geometry as driven by the error estimate. An example of the solution adaptive grid is shown in figure 1(b). Note that as the goal-oriented adaption approach utilizes information from both the primal and adjoint solutions, both primal and adjoint ‘flow features’ can be clearly seen as being resolved by the grid.
For each grid in the solution-adaptive sequence, as the solver reaches machine-zero convergence for the discrete residuals, the exact Jacobian 
$\boldsymbol {J}$ of the steady-state residual is computed and stored. The same is done for the mass matrix, which is not a trivially diagonal matrix for the considered case of the continuous finite-element spaces. Then, a separate code is used to solve the eigenvalue problem, based on the SLEPc package (Hernandes, Roman & Vidal Reference Hernandes, Roman and Vidal2005). The shift-and-invert approach with prescribed (
$\gamma ^*$, 
$\omega ^*$), and 
$\gamma ^* >0$, is used to capture the least-stable eigenvalues over a range of frequencies.
2.3. Weakly nonlinear analysis
Transonic buffet on high-Reynolds-number airfoils results from a Hopf bifurcation, similar to the onset of vortex shedding behind a circular cylinder at moderately low Reynolds numbers (Crouch et al. Reference Crouch, Garbaruk and Magidov2007). As the angle of attack is increased (for a fixed Mach number), the steady base flow loses stability and an oscillatory mode begins to grow. In general, a Hopf bifurcation could be supercritical or subcritical, resulting in different finite-amplitude behaviour. In the neighbourhood of the bifurcation, the amplitude of the variation can be described by the Landau equation (Drazin & Reid Reference Drazin and Reid1981):
\begin{equation} \frac{\partial A }{\partial t} = \gamma _0 A + \gamma _1 A^3, \end{equation}
where 
$\gamma _0$ is the linear growth rate and 
$\gamma _1$ is the Landau constant, also referred to as the first Lyapunov coefficient. The amplitude 
$A$ is a characteristic measure of the strength of the oscillation; here we define 
$A$ as one-half of the peak-to-peak variation of the normalized lift fluctuation 
$C_L'(t)= (C_L(t)-C_L^{(B)}) / C_L^{(B)}$, where 
$C_L(t)$ is the airfoil lift coefficient and 
$C_L^{(B)}$ is the lift coefficient of the steady base state. A lift-based amplitude definition provides a direct link to a buffet boundary, but an alternative amplitude definition such as the normalized shock displacement could also be used.
Based on experimental observations showing a smooth increase in 
$C_L(t)-C_L^{(B)}$ with angle of attack (McDevitt & Okuno Reference McDevitt and Okuno1985), the bifurcation is expected to be supercritical. For a supercritical bifurcation, the Landau constant 
$\gamma _1$ is negative and the supercritical oscillation is characterized by a finite-amplitude limit cycle. For these equilibrium states, the Landau constant can be determined from the equilibrium amplitude 
$A_e$ and the linear growth rate 
$\gamma _0$:
\begin{equation} \frac{\partial A_e }{\partial t} = \gamma _0 A_e + \gamma _1 {A_e}^3 = 0 \quad \rightarrow \quad \gamma_1 ={-}\gamma_0 / {A_e}^2. \end{equation}
Alternatively, the equilibrium amplitude can be determined from known values of the growth rate 
$\gamma _0$ and Landau constant 
$\gamma _1$.
In general, the Landau constant describing the finite-amplitude correction is complex, and includes a correction to the frequency 
$\omega _0$. In the current study, this frequency correction is less than 2 % (i.e. the equilibrium frequency is essentially the same as the linear frequency) and is not considered in the detailed results. In a perturbation analysis, two additional finite-amplitude corrections would be obtained at 
$O(A^2)$: a mean-flow correction and a second harmonic. These are briefly discussed in the results.
3. Results
The numerical simulations are initiated by perturbing a steady RANS solution with a small increase in angle of attack 
$\Delta \alpha$. Figure 2 shows the variation of the lift coefficient 
$C_L$ and the normalized lift fluctuation 
$C_L'$ for different values of 
$\Delta \alpha$. The resulting angle of attack 
$\alpha + \Delta \alpha$ is slightly different in each case, but this has negligible effects on the results of interest. Figure 2(b) shows that the initial amplitude is proportional to 
$\Delta \alpha$, but the linear growth rate and the equilibrium amplitude are independent of the initial amplitude.
Figure 2. Variation of (a) lift coefficient 
$C_L$ and (b) normalized lift fluctuation 
$C_L'$ with different levels of initial perturbation 
$\Delta \alpha$. OAT15A airfoil with 
$M=0.73$, 
$\alpha =3.3^{\circ }$.
The linear growth rate 
$\gamma _0$ and frequency 
$\omega _0$ are used to assess the grid resolution and time-step sensitivity for the unsteady simulations. These linear growth parameters show greater sensitivities than the equilibrium values. Table 1 shows the calculated instability characteristics for three different grids and time steps. The values for the medium grid with a 
$0.02 c/U_{\infty }$ time step are within 
$0.3\,\%$ of the finest-grid and minimum time-step results – thus, all subsequent results are based on these intermediate values.
Table 1. Variation of linear growth rate and frequency (
$\gamma _0$, 
$\omega _0$) from simulations as a function of time step and grid size. OAT15A airfoil with 
$M=0.73$, 
$\alpha =3.3^{\circ }$.
Figure 3 shows the growth rates and frequencies for different values of Mach number and angle of attack for the OAT15A airfoil. The extracted values from the numerical simulations (all based on 
$\Delta \alpha = 10^{-5}$) are in very good agreement with the GSA, even though these are from different codes using different grids and different numerical solvers. The stability results show slightly higher growth rates compared with the numerical simulations, consistent with the observations of Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009a). The symbols (
$\times$) show the buffet-onset angles of attack (for 
$M=0.72$, 
$0.73$) and the buffet frequencies (at 
$\alpha = 3.5^{\circ }$) from the experiments of Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009), which are in good agreement with the predictions.
Figure 3. (a) Growth rate 
$\gamma _0$ and (b) frequency 
$\omega _0$ as determined by the stability analysis (blue diamonds), the unsteady simulations (red squares), the earlier approach of Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009a) (grey dash) and the experiments of Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009) (
$\times$ symbols) for the OAT15A airfoil with 
$M=0.72$, 
$0.73$ and 
$0.74$.
For the 
$M=0.73$ case, results are also provided from the approach of Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009a) for comparison. Both the critical angle of attack and the frequency are in generally good agreement. Analyses of Sartor et al. (Reference Sartor, Mettot and Sipp2015) and Sansica et al. (Reference Sansica, Loiseau, Kanamori, Hashimoto and Robinet2022) both predict a critical angle of attack closer to 
$3.5^{\circ }$ for this Mach number, notably higher than the current results. This is attributed primarily to differences in the turbulence model used in those studies. The critical angle of attack from Paladini et al. (Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019) is 
$3.3^{\circ }$, only slightly higher than the 
$3.1^{\circ }$ reported here – the version of the turbulence model used in that study being similar to that of the current work. Each of these earlier works predicts a frequency near 
$\omega \approx 0.45$, in very good agreement with the current predictions.
The numerical simulation cases are each run sufficiently long to enable the extraction of an equilibrium amplitude 
$A_e$, similar to the example of figure 2. Equilibrium amplitudes for different Mach numbers are plotted as a function of 
$\alpha$ in figure 4. The results confirm that the Hopf bifurcation is supercritical, consistent with Sansica et al. (Reference Sansica, Loiseau, Kanamori, Hashimoto and Robinet2022). Although the results of Sartor et al. (Reference Sartor, Mettot and Sipp2015) and Sansica et al. (Reference Sansica, Loiseau, Kanamori, Hashimoto and Robinet2022) identify a different critical angle of attack, for an 
$\alpha$ increment of 0.5 above the critical value, their results give an equilibrium amplitude of 
$A_e \approx 0.1$ for 
$M=0.73$ – in good agreement with figure 4.
Figure 4. Equilibrium amplitudes as determined from the unsteady simulations for the OAT15A airfoil with (a) 
$M=0.72$, (b) 
$M=0.73$ and (c) 
$M=0.74$.
Using (2.6), the associated Landau constant is determined for each case based on the simulation amplitude 
$A_e$ and growth rate 
$\gamma _0$. The Landau constants are shown in figure 5, plotted against the 
$\alpha$ increment from the critical value 
$\alpha _{crit}$. Neglecting the value of 
$\gamma _1$ nearest to the onset of instability (which is highly sensitive to the estimated 
$A_e^2$ in (2.6)), the Landau constant shows very little change with 
$\alpha$ and Mach number, ranging between 
$-7 < \gamma _1 < -4.5$. This suggests that the buffet dynamics for a given airfoil can be represented by a simple function 
$\gamma _1 = f(\alpha - \alpha _{crit})$, or by a constant value – a linear fit is shown in figure 5 for comparison. Note that a meaningful quantitative comparison with the Landau constant calculated in Sansica et al. (Reference Sansica, Loiseau, Kanamori, Hashimoto and Robinet2022) is not possible, due to differences with the turbulence model used in that study.
Figure 5. Variation of the Landau constant 
$\gamma _1$ as determined from the unsteady simulations for the OAT15A airfoil with 
$M=0.72$, 
$0.73$ and 
$0.74$.
As noted in § 2.3, there is also an 
$O(A^2)$ correction to the average lift level resulting from the nonlinearity. Using results similar to figure 2, the adjustment to the average 
$C_L$ is estimated to be 
$C_L^{ave} \approx C_L^{(B)}(1-A^2)$. This correction estimate does not affect the estimated buffet amplitude, but it does provide a more complete description of the nonlinear dynamics.
Having determined the airfoil Landau constant from a few isolated simulations, the bifurcation diagrams can be efficiently constructed using the growth rates derived from the GSA. Figure 6 shows the global-stability-based predictions (blue curves using 
$\gamma _1 =-7.5 + 7(\alpha - \alpha _{crit})$) in comparison with the simulation results (red symbols). The comparisons show very good agreement over the range of amplitudes, extending to rather large peak-to-peak values greater than 
$20\,\%$.
Figure 6. Equilibrium amplitudes from stability analysis (blue lines) and simulation (red symbols) for the OAT15A airfoil with (a) 
$M=0.72$, (b) 
$M=0.73$ and (c) 
$M=0.74$.
To assess the generalizability of the results, we also consider the NACA0012 airfoil. Following a similar process to determine the Landau constant for two different Mach numbers, we obtain 
$\gamma _1 = -1.4 + 1.2 (\alpha - \alpha _{crit})$. The nominal value for the Landau constant 
$-1.3 < \gamma _1 < -0.7$ is smaller than for the OAT15A airfoil. Using the linear fit for 
$\gamma _1$, in combination with the growth rates from the GSA, the predicted bifurcation curve (blue line) is compared to the simulation (red symbols) in figure 7. The overall agreement is good, and is very good up to peak-to-peak amplitudes of 
$30\,\%$.
Figure 7. Equilibrium amplitudes from stability analysis (blue lines) and simulation (red symbols) for the NACA0012 airfoil with (a) 
$M=0.74$ and (b) 
$M=0.76$.
We now consider the case of an airfoil undergoing a continuous increase in angle of attack, analogous to the windup turn used for aircraft buffet determination. Figure 8 shows the variation of the lift coefficient as a function of the angle of attack 
$\alpha$ for two different pitch rates, 
$0.0001^{\circ }/(c/U_{\infty })$ and 
$0.001^{\circ }/(c/U_{\infty })$. In both cases, the 
$C_L$ values for 
$\alpha < 3.1$ are consistent with the steady-state calculations. Beyond 
$\alpha \approx 3.2$, the 
$C_L$ value for the lower pitch rate oscillates rapidly, consistent with the stability analysis. On the scale of this plot, the oscillations (in terms of cycles per degree of 
$\alpha$) are much more frequent than can be seen in the figure. The grey zone in the figure shows the trace of these oscillations. For the higher pitch rate, the 
$C_L$ oscillation is easily observed beyond 
$\alpha \approx 3.4$. While the discernable oscillations initially occur at higher 
$\alpha$ for the higher pitch rate, the envelope of the oscillations is independent of the pitch rate. The delay in the initial oscillation for higher pitch rates is a result of less time for the instability to grow over a given interval of 
$\alpha$ variation. The results show that empirical estimates of the buffet onset from a continuous-pitch study can be influenced by the pitch rate.
Figure 8. Variation of lift coefficient 
$C_L$ for a continuously pitched OAT15A airfoil at 
$M=0.73$. Time-varying simulations with pitch rates 
$0.0001^{\circ }/(c/U_{\infty })$ (grey line) and 
$0.001^{\circ }/(c/U_{\infty })$ (red line). Amplitude envelope from stability analysis (blue line).
The blue curves in figure 8 show the amplitude envelope as determined from the GSA using 
$C_L = C_L^{(B)} (1 - A_e^2 \pm A_e)$. Instability occurs at 
$\alpha \approx 3.1$, and for higher 
$\alpha$ values the limits on 
$C_L$ are based on the equilibrium amplitudes. The GSA provides a good estimate for the limiting behaviour, even in the continuous-pitch case. While the reduction in average 
$C_L$ is somewhat greater than the simple 
$A_e^2$ estimate, the fluctuation content of primary interest is well captured.
4. Conclusions
Weakly nonlinear analysis provides an effective framework for efficiently modelling the finite-amplitude oscillations linked to airfoil buffeting. Building on the global stability approach, a Landau equation is used to capture the post-critical behaviour of the dominant unsteadiness. Results show that the value of the Landau constant changes only weakly for a given airfoil, with varying Mach number and angle of attack. In practice, the Landau constant can be determined from a perturbation analysis, or from a small number of unsteady simulations. This can then be combined with linear growth rates from stability analysis to calculate equilibrium amplitudes. This provides a very efficient approach to develop finite-amplitude buffet boundaries.
The weakly nonlinear approach has been used to evaluate the amplitude behaviour resulting from a continuously pitched airfoil, analogous to an airplane wing during a windup manoeuvre. The stability analysis provides an envelope for the lift fluctuations, which is compared with URANS simulations using different pitch rates. Higher pitch rates result in delayed onset of finite-amplitude buffet oscillations, but the limiting values for each angle of attack are unchanged, and are well predicted by the GSA estimations.
The relative simplicity of the current approach – using only GSA and a limited number of simulations to capture the general nonlinear dynamics – makes it a good candidate for other problems where global stability results already exist. The effectiveness of the approach for the more general swept-wing buffet problem will be dependent on the relative role, or dominance, of different instabilities over a relevant range of angles of attack.
Declaration of interests
The authors report no conflict of interest.
