2.1.1. Flow solver

The numerical simulations were carried out using the in-house code, SBLI (Yao et al. Reference Yao, Shang, Castagna, Johnstone, Jones, Redford, Sandberg, Sandham, Suponitsky and De Tullio2009), which is a scalable compressible flow solver with multi-block and shock-capturing capabilities and has been used previously to study buffet (Zauner et al. Reference Zauner, De Tullio and Sandham2019; Zauner & Sandham Reference Zauner and Sandham2020a,Reference Zauner and Sandhamb). SBLI solves for the compressible Navier–Stokes equations which govern the flow evolution in a dimensionless form (see Zauner & Sandham Reference Zauner and Sandham2020b, pp. 3–4). The aerofoil chord, the freestream density, streamwise velocity (non-swept) and temperature are used as reference scales implying that their corresponding dimensionless equivalents are given by $c = \rho _\infty = U_\infty = T_\infty = 1$, respectively. Fourth-order finite difference schemes (central at interior and the Carpenter scheme Carpenter, Nordström & Gottlieb (Reference Carpenter, Nordström and Gottlieb1999) at boundaries) are employed for spatial discretisation, and a low-storage third-order Runge–Kutta scheme is used for time discretisation. To capture features of shock waves, a total variation diminishing scheme is employed, details of which can be found in Sansica (Reference Sansica2015).

2.1.2. LES

To perform the LES, we adopt the spectral-error-based implicit LES approach which has been validated against DNS for flows where buffet is observed (Zauner & Sandham Reference Zauner and Sandham2020b). This approach utilises the error estimator proposed in Jacobs et al. (Reference Jacobs, Zauner, De Tullio, Jammy, Lusher and Sandham2018) for identifying regions of insufficient grid resolution in DNS. Based on this estimator, a low-pass filter is used on all conserved variables to locally correct spectral deviations whenever and wherever they occur. Sixth-order compact finite difference schemes for filtering applications (Lele Reference Lele1992) are used for this purpose. A blending function is used to reduce the impact of filtering at higher wavenumbers. This is given by $q_{{updated}} = q_{{unfilt.}} - a_{{lim}} (q_{{unfilt.}}- q_{{filt.}})$, where the updated flow field is computed as an affine combination of the unfiltered and filtered flow fields (Bogey & Bailly Reference Bogey and Bailly2004) with the constant, $a_{{lim}} = 0.4$. The parametric values used here are the same as those used in Zauner & Sandham (Reference Zauner and Sandham2020b) (see their table 1).

2.1.3. Geometry, grid and boundary conditions

Dassault Aviation's supercritical, laminar V2C profile with a blunt TE (thickness 0.5 % chord), as used in the TFAST project (Billard et al. Reference Billard2021), is employed in the simulations. To accurately model the blunt TE, maintain grid smoothness and ensure flexibility in the distribution of grid points, C-H multiblock structured grids were generated using an in-house, open-source code (Zauner & Sandham Reference Zauner and Sandham2018). A grid was generated for each incidence angle so as to ensure that the aerofoil wake remains consistently well resolved as $\alpha$ changes. Features of a typical case of $\alpha = 4^{\circ }$ (‘reference grid’) are highlighted in figure 1(a). The aerofoil (cyan) is treated as an isothermal wall with the temperature set equal to that of the freestream. Integral and zonal characteristic boundary conditions (CBCs) (Sandhu & Sandham Reference Sandhu and Sandham1994; Sandberg & Sandham Reference Sandberg and Sandham2006) are used on the freestream (green) and outflow boundaries (brown), respectively. A localised filter is adopted to handle the singular points at the corners of the blunt TE by employing the strategy proposed in Jones, Sandberg & Sandham (Reference Jones, Sandberg and Sandham2006). For the unswept cases, $x$ and $z$ represent the dimensionless streamwise and spanwise Cartesian coordinates based on the freestream, respectively. The coordinate orthogonal to the two is represented by $y$. The chordwise coordinate is represented by $x'$. The curvilinear circumferential and radial coordinates are $\xi$ and $\eta$, respectively. The 3D grid was obtained by extruding in the $z$-direction with a uniform grid spacing of $10^{-3}$ for two different spanwise lengths, $L_z = 0.05$ (‘narrow domain’) and $L_z = 1$ (‘wide domain’ (WD), see Appendix A.1), with a total of approximately 75 million and 1.5 billion grid points present, respectively. For all cases, the spanwise boundaries are assumed to be periodic (to model an infinite unswept wing). The narrow domain is employed for all cases other than those used to examine the effect of sweep. Note that when $\varLambda = 0^{\circ }$, the flow configuration is of an unswept wing and the $x$-coordinate is based on the streamwise direction. Thus, the freestream $x$-velocity component is unity whereas other components are zero.

Figure 1. Multi-block grid shown by plotting every $15\mathrm {th}$ grid point in $\xi$ and $\eta$ directions for the case of $\alpha = 4^{\circ }$: (a) entire domain and (b) vicinity of aerofoil. Characteristic boundary conditions (CBCs) are applied on the inflow and outflow boundaries, while isothermal no-slip conditions are applied on the wall. The pink dashed curve is positioned at a normal distance of 0.05 from the aerofoil surface and is used to monitor shock wave features.

The reference grid's features in the aerofoil's vicinity are highlighted in figure 1(b). The blunt TE contains 30 grid points. A curve at a constant wall-normal distance of $0.05c$ from the aerofoil's surface, is also shown (dashed curve). The latter is used for monitoring shock wave motion and will be referred to as C5 henceforth. The grid clustering is relatively denser close to the aerofoil surface, in its wake, and in the region where the shock wave is expected. The wall-normal and wall-parallel spacings at the wall vary between $1\times 10^{-4}$ to $1.7\times 10^{-4}$ and $4\times 10^{-4}$ to $2\times 10^{-3}$, respectively.

2.1.4. Flow parameters

The fluid is assumed to be a perfect gas with a specific heat ratio, $\gamma = 1.4$ which satisfies Fourier's law of heat conduction (Prandtl number, $\textit {Pr} = 0.72$). It is also assumed to be Newtonian, with its viscosity variation with temperature governed by the Sutherland's law (Sutherland coefficient, $C_{Suth} = {110.4}/{268.67} \approx 0.41$). We examine the effect of the flow parameters, $M$, $\alpha$, $\textit {Re}$ and $\varLambda$ on buffet by starting from a baseline ‘reference’ case with $M = 0.7$, $\alpha = 4^{\circ }$, $\textit {Re} = 5\times 10^{5}$ and $\varLambda = 0^{\circ }$ and varying the value of only one of these parameters while keeping the other parameters the same as that of the reference. In addition to these, a case at $M = 0.8$ and $\alpha = 0^{\circ }$, with other parameters at reference value, was also simulated to examine buffet features at zero incidence. This is referred to as the A0M8 case (see § 3.2.3). A complete list of all cases studied (except those on grid convergence and domain extent) and the parametric values for each are provided in table 1. Note that the grid resolution close to the leading edge is relatively coarser than at other regions of the aerofoil (figure 1b). This was found to introduce minor grid-level oscillations in the flow field only for the highest incidence angle considered ($\alpha = 6^{\circ }$, see table 1). This occurred in parts of the buffet cycle when the shock wave moved in close proximity to the leading edge due to the high amplitude of buffet at this $\alpha$. However, this was only during a small part of the buffet cycle and no significant effect on the buffet flow features could be discerned. Nevertheless, higher angles of attack ($\alpha > 6^{\circ }$) were not examined for this reason. The effect of varying sweep angle, $\varLambda$ is also reported here, but due its limited scope, we present it in Appendix A.

Table 1. Parameter values for various cases simulated (cases with buffet in boldface).

2.1.5. Temporal features

For all cases except those where $\alpha$ is varied, the simulations were initialised with freestream conditions used at all grid points. For the study on varying $\alpha$, because it was originally initiated as part of a different study, the initial conditions were chosen as the fully evolved solutions of the reference case. Note that the choice of initial conditions have been shown to not affect flow features that occur past the transients for supercritical aerofoils (Xiao et al. Reference Xiao, Tsai and Liu2006). A constant dimensionless time step of $3.2\times 10^{-5}$ is used (approximately $3\times 10^{5}$ time steps per buffet cycle).

We denote the time taken for the flow to evolve past initial transients as $t_0$. When oscillations resembling buffet are present in $C_L(t)$ (with $St \sim 0.1$), this approximate transient time, $t_0$, is chosen such that it coincides with the high-lift phase and such that $C_L(t)$ is periodic and the oscillations persist (i.e. do not monotonically dampen) for $t \geqslant t_0$ when visually inspected (e.g. compare figures 2 and 5(a) for the reference case, where $t_0 = 18$). Beyond the transients, the simulations were run until at least 10 low-frequency cycles were completed (or for 20 time units when such oscillations are absent) for the non-swept cases so as to improve statistical convergence of mean quantities and the accuracy of modal decomposition. To confirm if buffet is present, the power spectral density (PSD) of the fluctuating component of the span-averaged lift coefficient, PSD($C_L'$), was computed based only the data past the initial transient time (i.e. $C_L'(t)$ for $t > t_0$). As the signals are strongly coherent, a periodogram of $C_L'$ was computed with frequency resolution $\delta f \approx 0.01$, sampling frequency, $f_S = 128$ and a Hamming window (to reduce spectral leakage). Through visual inspection, buffet is confirmed to have occurred when there is a time past which $C_L(t)$ exhibits sustained (non-decaying) low-frequency oscillations and when a discernible local maximum is present in the power spectrum in the low-frequency range $0 < St < 0.5$. In addition to the PSD, a continuous wavelet transformation of the same signal ($C_L'$) was also computed using MATLAB. For this, the commonly used Morse wavelet, with symmetry parameter as 3 and time-bandwidth product as 60, was employed.

Figure 2. Temporal variation of span-averaged aerofoil coefficients, lift, pressure drag, friction drag and TE pressure (corner on suction surface) past initial transients for the reference case. Dashed lines correspond to high-lift (red), low-lift (blue), low-friction-drag (green) and high-friction-drag (brown) phases.

2.2.1. Considerations on domain extent

The validation of the LES with DNS in Zauner & Sandham (Reference Zauner and Sandham2020b) was performed in a narrow domain of $L_z = 0.05$. To examine the effect of spanwise width, an LES for a WD of $L_z = 1$ was also performed. It was observed that some differences exist between the two, including increases in the buffet amplitude (${\approx }40\,\%$ difference in PSD of lift fluctuations, compare reference and $\varLambda = 0^{\circ }$ (WD) cases in table 2) and regularity of the oscillations for the WD, but the major flow features, including the development of multiple shock wave structures, their propagation features and buffet frequency (variation of ${\approx }6\,\%$) remain similar.

Table 2. Mean aerofoil coefficients, separation length and buffet measures for various cases simulated. Here $St_b$ is the buffet frequency and PSD($C_L',St_b$) is the PSD of the fluctuating component of lift at $St_b$.

As the focus of this study is on examining qualitative aspects of buffet, we have chosen $L_z = 0.05$. Results from other previous studies on turbulent buffet also support this choice. Experiments reported in Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009) have shown that although weak 3D features associated with flow separation are present in the flow, buffet on unswept wings is characterised by a 2D mode. Wall-modelled LES in a domain of span 0.065$c$ was shown in Fukushima & Kawai (Reference Fukushima and Kawai2018) to accurately capture buffet features observed in experiments. Similarly, in Grossi et al. (Reference Grossi, Braza and Hoarau2014), 3D structures with the spanwise wavelength between $0.029c$ and $0.04c$ were observed in DES.

The flow field in the vicinity of the inflow and outflow boundaries was checked a posteriori for gradients. These were found to be insignificant for all cases simulated. To further investigate the effect of domain extent, a new grid was generated with the domain boundaries extended further by $40\,\%$. Simulations were carried out for the same settings as $M= 0.85$ and $M = 0.9$. The latter is the case for which the supersonic region's extent was found to be largest, whereas the former is the case for which it is largest and for which buffet occurs (see, e.g., figures 11 and 13). We observed that for $M = 0.85$, the lift coefficient variation and buffet frequency are not significantly different ($St_b\approx 0.33$, $\bar {C}_L = 0.16$ and the root-mean-square value of $\bar {C}_L' = 0.025$ for both cases) suggesting that for this $M$ and, by extension, for lower $M$ (where supersonic regions are smaller), the effect of extending the domain further on buffet is negligible. For $M = 0.9$, $\bar {C}_L$ was found to increase from 0.35 to 0.39 and the extent of the supersonic region also expanded. Thus, the results of this case alone must be interpreted with caution. However, because the flow remained approximately steady with no buffet present for either domain extent, this case of $M = 0.9$ was not explored further.

2.2.2. Grid resolution

In the present simulations, the main feature of the flow is a laminar BL that develops from the leading edge up to the shock foot. To examine if this is accurately captured, a new grid with 25 % more points was generated (${\approx }100$ million grid points), with the grid spacing in $\xi$ and $\eta$ directions in the laminar BL region ($0 \leqslant x \leqslant 0.5$) approximately halved. The simulation was carried out at the highest Reynolds number studied here ($\textit {Re} = 1.5\times 10^{6}$). The results, including instantaneous lift coefficient evolution and mean surface coefficients, match well for the two grids (not shown, provided as supplementary material available at https://doi.org/10.1017/jfm.2022.471), confirming the adequacy of the original grid in capturing the laminar BL.

In the regions where the mean flow is attached and turbulent, for most cases studied, the grid spacing in wall units at the aerofoil surface computed based on the mean wall shear stress, $\Delta \xi ^{+}$, $\Delta \eta ^{+}$ and $\Delta z^{+}$, were found to be approximately 10, 1 and 10, respectively. This indicates that the grid resolution is more than adequate for LES (Garnier & Deck Reference Garnier and Deck2010; Dandois et al. Reference Dandois, Mary and Brion2018). For a few cases close to buffet onset ($\alpha = 3^{\circ }$ and $M \leqslant 0.69$), the $\Delta \eta ^{+}$ value is higher (${\approx }2$), but only in a small region downstream to where shock waves occur (approximately $0.85 \leqslant x \leqslant 0.9$, see figure 6b). Additionally, for $\textit {Re} > 5\times 10^{5}$, we found $\Delta \xi ^{+}$, $\Delta z^{+} \approx 30$ and $\Delta \eta ^{+} \approx 5$ in some regions of the flow. Although the $\Delta \xi ^{+}$ and $\Delta z^{+}$ are still adequate (for comparison, Garnier & Deck (Reference Garnier and Deck2010) used $\Delta \xi ^{+} \approx 50$ and $\Delta z^{+} \approx 20$ for LES), the $\Delta \eta ^{+}$ values are higher than recommended. We do not expect this to play a significant role in affecting buffet dynamics because we have a laminar BL upstream to the shock wave and, downstream, a separation bubble which reattaches only in parts of the buffet cycle to form a turbulent BL which persists only a short distance before reaching the TE. Care has been taken to ensure that the regions associated with separated flow are well-resolved, based on visual checks of instantaneous flow fields and using the spectral-error indicator to check filter activity. In addition, as noted previously, a grid study at $\textit {Re} = 1.5 \times 10^{6}$, did not show any significant changes in buffet features. Thus, the adequacy of the grid is supported by a combination of wall-unit checks, grid sensitivity studies and careful monitoring of simulations using the spectral error detector.

2.3.1. Formulation

Given an ensemble of spatiotemporal realisations that represent a stochastic process, $\{\boldsymbol {q}_\zeta (\boldsymbol {x},t)\}$, with $\zeta$ denoting the realisation's index, a set of orthogonal basis functions, $\{\boldsymbol {\phi }_i(\boldsymbol {x},t)\}$, that ideally represents the coherent aspects of the process is determined by maximising a utility functional (Lumley Reference Lumley1970; Towne et al. Reference Towne, Schmidt and Colonius2018),

(2.1)\begin{equation} J(\boldsymbol{\phi}) = E_\zeta\left[\frac{|\langle\boldsymbol{q}_\zeta,\boldsymbol{\phi}\rangle_{(\boldsymbol{x},t)}|^{2}}{ \langle\boldsymbol{\phi},\boldsymbol{\phi}\rangle_{(\boldsymbol{x},t)}}\right]. \end{equation}

Here $E_\zeta []$ is the expectation operator, and $\langle,\rangle _{(\boldsymbol {x},t)}$ denotes the appropriate inner product over space, $\boldsymbol {x}$, and time, $t$. This allows for choosing each basis function such that the projection of $\boldsymbol {q_\zeta }$ on $\boldsymbol {\phi }_i$, (i.e. $\langle \boldsymbol {q}_\zeta,\boldsymbol {\phi }_i\rangle _{(\boldsymbol {x},t)}$), is maximised over $\zeta$ in a least-square sense. It is possible to show (e.g. by assuming that $\boldsymbol {\phi }_i$ are normal and reformulating using Lagrange multipliers) that the extrema of this functional satisfy the eigenvalue problem

(2.2)\begin{equation} \langle {\boldsymbol{\mathsf{C}}}(\boldsymbol{x},\boldsymbol{x}',t,t'),\boldsymbol{\phi}_i(\boldsymbol{x}',t')\rangle_{(\boldsymbol{x}',t')} = \lambda_i\boldsymbol{\phi}_i(\boldsymbol{x},t), \end{equation}

where ${\boldsymbol{\mathsf{C}}}$ is the two-point spatiotemporal correlation tensor based on any two points $\boldsymbol {x}$ and $\boldsymbol {x}'$ in the spatial domain and $t$ and $t'$ in time, and $\boldsymbol {\phi }_i$ and $\lambda _i$ are the $i$th eigenfunction and eigenvalue, respectively. Note that, as is common for eigenvalue problems, the index $i$ on the right-hand side does not denote a summation. Here, $\lambda _i = J(\boldsymbol \phi _i)$, implying that the maximum of $J$ is given by the maximum eigenvalue. The expansion of any realisation, $\boldsymbol {q}_\zeta$, using the basis $\boldsymbol {\phi }$ is given by

(2.3)\begin{equation} \boldsymbol{q}_\zeta(\boldsymbol{x},t) = \sum_i \langle\boldsymbol{q}_\zeta,\boldsymbol{\phi}_i\rangle_{(\boldsymbol{x},t)} \boldsymbol{\phi}_i(\boldsymbol{x},t) = \sum_i a_{i,\zeta} \boldsymbol{\phi}_i(\boldsymbol{x},t), \end{equation}

where $a_{i,\zeta }$ are the expansion coefficients. It follows from the above equations and the orthonormality of the basis functions ($\langle \boldsymbol {\phi }_j, \boldsymbol {\phi }_k\rangle _{(\boldsymbol {x},t)} = \delta _{jk}$) that

(2.4)\begin{align} E_\zeta[a_{j,\zeta}a_{k,\zeta}^{*}] &= \langle\langle {\boldsymbol{\mathsf{C}}}(\boldsymbol{x},\boldsymbol{x}',t,t'),\boldsymbol{\phi}_j(\boldsymbol{x}',t') \rangle_{(\boldsymbol{x}',t')},\boldsymbol{\phi}_k(\boldsymbol{x},t)\rangle_{(\boldsymbol{x},t)}\nonumber\\ &\Rightarrow E_\zeta[a_{j,\zeta}a_{k,\zeta}^{*}] = \lambda_j \delta_{jk}, \end{align}

where $\delta _{jk}$ is the Kronecker Delta function and $()^{*}$ denotes complex conjugation.

For a zero-mean, statistically stationary flow, the correlation tensor is dependent only on the time interval, $\tau = t'-t$ and, thus, a Fourier transform in time can be used to compute the cross-spectral density tensor, ${\boldsymbol{\mathsf{S}}}$. The relations between the two are

(2.5)\begin{equation} \left.\begin{gathered} {\boldsymbol{\mathsf{S}}}(\boldsymbol{x},\boldsymbol{x}',f) = \int_{-\infty}^{\infty} {\boldsymbol{\mathsf{C}}}(\boldsymbol{x},\boldsymbol{x}',\tau) \exp({-2\mathrm{i}{\rm \pi} f\tau})\,{\rm d}\tau, \\ {\boldsymbol{\mathsf{C}}}(\boldsymbol{x},\boldsymbol{x}',\tau) = \int_{-\infty}^{\infty} {\boldsymbol{\mathsf{S}}}(\boldsymbol{x},\boldsymbol{x}',f) \exp({2\mathrm{i}{\rm \pi} f\tau})\,{\rm d}f, \end{gathered}\right\}\end{equation}

where $f$ is the frequency. Following the proof outlined in Towne et al. (Reference Towne, Schmidt and Colonius2018, pp. 859–860), we can relate the eigenvalue problem of ${\boldsymbol{\mathsf{C}}}$ with that of ${\boldsymbol{\mathsf{S}}}$. Assuming the inner product to be a weighted integral in space and time, and substituting for ${\boldsymbol{\mathsf{C}}}$ in (2.5) into (2.2), we get

(2.6)\begin{equation} \int_{-\infty}^{\infty}\int_\varOmega\int_{-\infty}^{\infty} {\boldsymbol{\mathsf{S}}}(\boldsymbol{x},\boldsymbol{x}',f) \exp({2\mathrm{i}{\rm \pi} f(t-t')})\,{\rm d}f \,{\boldsymbol{\mathsf{W}}}(\boldsymbol{x}') \boldsymbol{\phi}_i(\boldsymbol{x}',t')\,{{\rm d}\kern0.06em x}'\,{\rm d}t' = \lambda_i\boldsymbol{\phi}_i(\boldsymbol{x},t), \end{equation}

where ${\boldsymbol{\mathsf{W}}}$ is the weight matrix associated with the quadrature on the curvilinear grid on the spatial domain, $\varOmega$. From the definition of a Fourier transform in time of $\boldsymbol {{\phi }}_i(\boldsymbol {x}',t')$, $\boldsymbol {\hat {\phi }}_i(\boldsymbol {x}',f)$, we obtain

(2.7)\begin{align} \int_{-\infty}^{\infty}\int_\varOmega {\boldsymbol{\mathsf{S}}}(\boldsymbol{x},\boldsymbol{x}',f) \,{\boldsymbol{\mathsf{W}}}(\boldsymbol{x}') \hat{\boldsymbol{\phi}}_i(\boldsymbol{x}',f) \exp({2\mathrm{i}{\rm \pi} ft})\,{{\rm d}\kern0.06em x}'\,{\rm d}f = \int_{-\infty}^{\infty} \lambda_i\hat{\boldsymbol{\phi}}_i(\boldsymbol{x},f)\exp({2\mathrm{i}{\rm \pi} ft})\,{\rm d}f. \end{align}

To solve further, we assume the following:

(2.8)\begin{equation} \hat{\boldsymbol{\phi}}_i(\boldsymbol{x},f) = \boldsymbol{\psi}_i(\boldsymbol{x},f_0)\delta(f-f_0). \end{equation}

Substituting in (2.7) and eliminating $\exp {(2\mathrm {i}{\rm \pi} f_0t)}$ on both sides leads to the eigenvalue problem:

(2.9)\begin{equation} \int_\varOmega {\boldsymbol{\mathsf{S}}}(\boldsymbol{x},\boldsymbol{x}',f_0) \,{\boldsymbol{\mathsf{W}}}(\boldsymbol{x}') \boldsymbol{{\psi}}_i(\boldsymbol{x}',f_0) \,{{\rm d}\kern0.06em x}' = \lambda_i\boldsymbol{\psi}_i(\boldsymbol{x},f_0). \end{equation}

Here, $\boldsymbol {\psi }_i(\boldsymbol {x},f_0)$ and $\lambda _i$ are the $i$th eigenfunction and eigenvalue at a given frequency, $f_0$, respectively, with the former referred to as an SPOD mode.

Flow fields were reconstructed using the SPOD mode at a given $f_0$ by reverting back to the time domain. Using an inverse Fourier transform in time in (2.8), we have

(2.10)\begin{equation} \boldsymbol{\phi}_i(\boldsymbol{x},t) = \int_{-\infty}^{\infty} \boldsymbol{\psi}_i(\boldsymbol{x},f_0)\delta(f-f_0)\exp({2\mathrm{i}{\rm \pi} ft})\,{\rm d}f = \boldsymbol{\psi}_i(\boldsymbol{x},f_0) \exp{(2{\rm \pi}\mathrm{i} f_0 t)}. \end{equation}

In addition, although the expansion coefficients, $a_{i,\zeta }$, vary for different realisations, we are interested in one that best captures the entire ensemble. Based on (2.4), we achieve this by choosing an ideal realisation with an ‘average’ expansion coefficient

(2.11)\begin{equation} a_{i,\zeta_0} = \sqrt{E_\zeta\left[|a_{i,\zeta}|^{2}\right]} = \sqrt\lambda_i . \end{equation}

Note that SPOD is carried out after subtracting the mean flow field. Thus, the reconstruction based on the required SPOD mode is obtained by summing up the mean with the real part of the truncated basis expansion. From (2.3), (2.10) and (2.11) we obtain

(2.12)\begin{equation} \tilde{\boldsymbol{q}}(\boldsymbol{x},t) = \bar{\boldsymbol{q}}(\boldsymbol{x},t) + \mathrm{Re}\{\sqrt{\lambda_i}\boldsymbol{\psi}_i(\boldsymbol{x},f_0) \exp{(2{\rm \pi}\mathrm{i} f_0 t)}\}, \end{equation}

where $\tilde {\boldsymbol {q}}$ and $\bar {\boldsymbol {q}}$ represent reconstructed and mean quantities, respectively. This approach is sufficient for the present study, but we note that a more general, low-rank flow reconstruction using SPOD has been implemented for a compressible turbulent jet at a low Mach number in a recently published study (Nekkanti & Schmidt Reference Nekkanti and Schmidt2021).

2.3.2. Implementation

Indexing the eigenvalues such that $\lambda _1 > \lambda _2 > \lambda _3 > \cdots,$ we have the most-energetic SPOD mode as $\boldsymbol {\psi _1}$ with higher indices representing lower energies. The first two SPOD modes were computed in this study using the memory-efficient streaming algorithm and the software tool provided in Schmidt & Towne (Reference Schmidt and Towne2019). Negligible energy content was observed for $\boldsymbol {\psi _2}$ for all cases studied (comparisons of $\lambda _i$ provided as supplementary material) and thus only features of $\boldsymbol {\psi _1}$ are reported. The weight matrix, ${\boldsymbol{\mathsf{W}}}$, was chosen to represent a weighted 2-norm (see Schmidt & Colonius Reference Schmidt and Colonius2020, § B.2, p. 1027), with its elements being the area associated with each grid point related to the 2D data used for SPOD. The area was approximated as that of the quadrilateral with vertices as the centroids of the cells surrounding a given grid point. SPOD based on an alternative choice of ${\boldsymbol{\mathsf{W}}} = {\boldsymbol{\mathsf{I}}}$ (identity matrix) was also implemented and found to have no significant effect on the qualitative features of the mode shapes or the eigenvalue spectra.

To perform SPOD, two sets of data were used. The first, henceforth referred to as ‘Data-2D’, is based on the instantaneous flow field variables density, velocity vector and pressure, all extracted from the $z = 0$ plane (no span-averaging). The data based on all these variables at a given time instant were combined together into a ‘snapshot’. Snapshots were stored at regular time intervals of $\Delta t = 0.08$ (sampling frequency, $F_s = 12.5$ implying ${\approx }125$ snapshots per buffet cycle). These were then divided into blocks of total time, $T_\zeta$, each block representing a realisation of the stochastic process. Their number is further increased by using a 50 % overlap and this ensemble is used to compute ${\boldsymbol{\mathsf{S}}}$ using the Welch's method. As noted previously, at least 10 buffet cycles were simulated for all cases where buffet occurs. To examine the low-frequency buffet which occurs at $St \approx 0.1$, it was ensured that at least four cycles occur in a block by choosing $T_\zeta = 44$ (frequency resolution, $\Delta F_\zeta \approx 0.02$). In addition to buffet, we also observed high-frequency coherent structures that we refer to as ‘wake’ modes ($St \sim O(1$), see figure 23a). To obtain a better spectral estimate of the energy associated with these modes, the number of blocks for SPOD was increased, which reduces the duration associated with each individual block (and, hence, the frequency resolution). For these modes, $T_\zeta = 5$ was chosen, which allows for at least five cycles to be captured while increasing the total number of blocks. As the expected value increases with an increase in the number of realisations in an ensemble, the increase in blocks allows for a better spectral estimate. A higher sampling frequency of $F_s = 125$ was also examined, but no new features were observed in the spectrum.

The second set of data used for SPOD, Data-SpAv, is based on the variation of instantaneous span-averaged wall-pressure and skin-friction coefficients and the local Mach number on the monitor curve C5 shown in figure 1(b). This was only used for flow reconstruction related to a simple one-dimensional representation (in space) of the flow field (e.g., figure 27). The reconstructed $C_p$ allows to examine the pressure variations on the aerofoil surface, and the reconstructed $C_f$ shows where the flow separates and the reconstructed Mach number giving the approximate position of the shock wave. We found that because the coefficients are span-averaged, this further reduces the non-coherent noise that is present in the SPOD modes. We note that the coherent features extracted are not significantly different between the two approaches and, thus, we have used the second approach exclusively for examining the reconstructed flow field only on the aerofoil surface and the curve C5.

Flow reconstruction based on the desired SPOD mode was carried out based on (2.12). However, because the reconstruction is carried out only at a specified frequency $f_0$, it is convenient to represent the reconstructed flow based on the phase of the sinusoidal cycle instead of time, $t$. Thus, we use $\phi = 2{\rm \pi} f_0 t$ to compute $\tilde {\boldsymbol {q}}(\boldsymbol {x},\phi )$, with $\phi = 0^{\circ }$ chosen as when the lift coefficient is the maximum.

3.2.1. Buffet onset and deep buffet

The temporal evolution of $C_L$ for a few select cases of $M$ in the range $0.5 \leqslant M\leqslant 0.775$ is shown in figure 5(a). For all cases simulated in the range $0.5 \leqslant M \leqslant 0.68$ (see table 1), the variation of $C_L$ is similar to that shown for $M =0.6$ in the figure, with the equilibrium values attained past transients being approximately the same. This occurs even though the flow is entirely subsonic for $M = 0.5$ and becomes transonic at and above the approximate critical Mach number of $M = 0.6$. Interestingly, at $M = 0.69$, sinusoidal oscillations can be discerned during the initial evolution ($t < 50$), but they slowly dampen with time, with the $C_L$ stabilising at a value close to that of $M = 0.6$. By contrast, a small increase in $M$ to 0.7 (reference case) leads to relatively large-amplitude sustained oscillations implying an abrupt onset of buffet with $M$, which is similar to the results reported for turbulent buffet in Giannelis, Levinski & Vio (Reference Giannelis, Levinski and Vio2018). This trend suggests that a marginally stable eigenmode becomes unstable as the parameter, $M$, is increased from $M=0.69$ to the onset value of $M = 0.7$. This is characteristic of a supercritical Hopf bifurcation and supports the global linear instability model proposed for buffet (Crouch et al. Reference Crouch, Garbaruk and Magidov2007). Interestingly, the maximum of $C_L(t)$ at $M= 0.7$ coincides with that of the approximately steady value attained for pre-onset conditions examined. In addition, a sharp drop in mean lift coefficient past onset is apparent, with $\bar {C}_L$ (the time average of $C_L$) dropping almost linearly with $M$ in a small range of $0.7 \leqslant M \leqslant 0.775$ from approximately 0.85 to $0.3$.

Figure 5. Effect of varying freestream Mach number: (a) temporal variation of lift coefficient and (b) PSD of its fluctuating component as a function of the Strouhal number (for cases for which buffet occurs). The reference case of $M = 0.7$ is shown using a black curve. Dashed vertical lines and circles highlight the buffet Strouhal number, $St_b$, and peaks in the PSD.

The PSD estimate of the fluctuating component of the lift coefficient, $C_L'$ is shown in figure 5(b) for cases where buffet is observed. Here, $C_L'(t) = C_L(t)-\bar {C}_L$, where $\bar {C}_L$ is the time-averaged lift coefficient and the PSD is computed using a periodogram of this signal for times after buffet is established ($t > 18$ for the reference case, see figure 2 and § 2.1.5 for details). The dominant peaks and their corresponding $St$ are highlighted using circles and broken horizontal lines, respectively, with the latter indicating the buffet frequencies, $St_b$. It is evident from the figure that the buffet frequency increases monotonically in the range shown. Such a monotonic increase in frequency with $M$ for turbulent buffet has been previously reported for aerofoils (Dor et al. Reference Dor, Mignosi, Seraudie and Benoit1989; Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009; Giannelis et al. Reference Giannelis, Levinski and Vio2018; Brion et al. Reference Brion, Dandois, Mayer, Reijasse, Lutz and Jacquin2020) and swept wings (Dandois Reference Dandois2016). Note that the buffet frequency seen here is somewhat higher ($St_b \approx 0.1$, see also table 2) compared with those reported for other commonly used aerofoils (OALT25, OAT15A, NACA0012, for which $0.06 \leqslant St_b \leqslant 0.08$ (Dandois Reference Dandois2016)). However, similar results have been reported for the V2C in simulations (Szubert et al. Reference Szubert, Grossi, Jimenez Garcia, Hoarau, Hunt and Braza2015) and preliminary experiments (Zauner et al. Reference Zauner, Renou, Dandois, Brion, Moise and Sandham2021). For a different supercritical laminar aerofoil, DRA-2303, Hartmann et al. (Reference Hartmann, Klaas and Schröder2012) report a case for which the Strouhal number can be computed to be $St_b = 0.1082$, which is closer to the frequency observed here.

The flow features of pre-onset cases of $M$ are shown in figure 6 using streamwise density gradient contours. At $M = 0.65$, a small elongated pocket of supersonic region is present. Although there is no shock wave terminating the supersonic region, multiple weak pressure waves can be discerned in the flow at $x \approx 0.5$. Similar ‘wavelets’ were reported in Hilton & Fowler (Reference Hilton and Fowler1947) for low $M$ at which no shock waves were present. In contrast to $M=0.65$, a shock system is observed at $M = 0.69$, with the pressure waves now strengthening sufficiently to form oblique shock waves. As noted previously, we use the term ‘shock wave’ to refer to regions where there are steep adverse pressure gradients (few grid points) within which the local Mach number, $M_{loc}$, becomes unity. By contrast, the term ‘pressure wave’ is used when there are isolated regions of weaker pressure gradients for which $M_{loc}$ is everywhere above or below unity. Thus, shock waves have contours lines of $M_{loc} = 1$ passing through them, whereas pressure waves do not. Note that the time instant shown for $M = 0.69$ is well past the initial transient evolution and after the buffet mode is fully damped (cf. figure 5a). Interestingly, these flow features of $M = 0.69$ resemble those seen for the reference case ($M = 0.7$) in the high-lift phase (see figure 3a), which is also indicated by the evolution of $C_L$ for these cases.

Figure 6. Streamwise density gradient contours on the $x\unicode{x2013} y$ plane shown for times well past initial transients for freestream Mach numbers below buffet onset: (a) $M = 0.65$ and (b) $M = 0.69$. The sonic line is highlighted using a grey curve.

Flow features above the onset $M$ are shown for $M = 0.735$ (top) and $M = 0.75$ (bottom) in figure 7 at high- (left) and low-lift (right) phases of the buffet cycle (see also supplementary movie 1). Note that the maximum PSD($C_L'$) at buffet frequency is achieved for $M = 0.735$ (cf. figure 5b). For both cases, the flow features far from the aerofoil surface change drastically in a given cycle indicating that buffet can affect the entire supersonic region. Contours of spanwise vorticity magnitude, $|\varOmega _z|$ at these instants are shown in figure 8 for $M = 0.735$. Strongly coherent, concentrated vortices that shed into the wake can be seen in the low-lift phase. Indeed, for all cases where buffet is observed in this study, the vortices were found to be strongest when the sonic edge is at its most upstream position (see also figures 12, 18 and 21). This is also supported by the results based on a time-frequency analysis reported later in § 4.1.1. This contradicts the model proposed in Hartmann et al. (Reference Hartmann, Feldhusen and Schröder2013), which requires that the vortices that reach the TE are most intense when the shock wave is at its most downstream position and vice versa (see their figure 15, p. 14).

Figure 7. Streamwise density gradient contours on the $x\unicode{x2013} y$ plane shown for the high (a,c) and low (b,d) lift phases of the buffet cycle for the cases $M = 0.735$ (a,b) and 0.75 (c,d). The sonic line is highlighted using a grey curve.

Figure 8. Contours of spanwise vorticity magnitude on the $x\unicode{x2013} y$ plane shown for (a) high- and (b) low-lift phases of the buffet cycle for the case $M = 0.735$. The sonic line is highlighted using a grey curve.

The spatiotemporal variation of the pressure gradient along the C5 monitor curve is shown in figure 9. Multiple transitions between locally supersonic and subsonic regions that were observed for $M =0.7$ as the sonic edge moves downstream (figure 4) are reduced here, as indicated by the reduction in the number and extent of the black curves. The approximate times associated with the local extrema of $C_L$ and $(C_D)_f$ are highlighted in the figure using dashed vertical lines (see figure 4 for the colour scheme). As expected, the maximum (minimum) $C_L$ occurs when the sonic edge is at its most downstream (upstream) position. Interestingly, similar to the reference case, the minimum (maximum) $(C_D)_f$ occurs when the sonic edge is approximately in the mid position of its upstream (downstream) motion. Downstream propagating structures in the BL, similar to those seen for the reference case, can also be observed. The convection velocities of these were observed to be slightly lower than in the reference case, with the approximate range of their phase speed being $0.55\unicode{x2013} 0.65U_\infty$, and the time required for them to reach the TE is in the range $0.5 \leqslant t_{down} \leqslant 1.5$ (which is used to comment on Lee's model in § 4.1.2).

Figure 9. Spatiotemporal variation of streamwise pressure gradient on the suction side of C5 for (a) $M = 0.735$ and (b) $M = 0.75$. The sonic line is highlighted using black curves. The approximate times associated with different phases of interest, high lift (red), low skin-friction drag (green), low lift (blue) and high skin-friction drag (brown), are highlighted using dashed lines.

Mean pressure and skin friction coefficients $\bar {C}_p$ and $\bar {C}_f$ (computed by time- and span-averaging the flow field) are compared for various $M$ in figure 10. As $M$ increases, we see an increase of $\bar {C}_p$ on the suction surface, and a reduction on the pressure surface. This accounts for the monotonic reduction in mean lift with $M$ seen in figure 5(a). The presence of large-amplitude buffet causes what has been referred to as shock smearing (Giannelis et al. Reference Giannelis, Vio and Levinski2017), where the mean pressure increases along the aerofoil surface in a gradual fashion in contrast to the abrupt increase present when a stationary shock wave occurs (Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009). This can be inferred by comparing the pre-onset cases of $M \leqslant 0.69$ with the rest. From figure 10(b), we see that the mean skin-friction remains approximately a constant on the pressure surface for all $M$, but changes significantly on the suction surface when $M$ is increased above onset. For $M = 0.6$, the presence of a short separation bubble ($C_f < 0$) can be inferred in $0.64 \leqslant x' \leqslant 0.66$, whereas for $M = 0.69$ and 0.7, there is another region of flow separation seen slightly upstream to this bubble, with the two merging for $M\geqslant 0.735$ to form a separation zone that extends from upstream of the shock foot to close to the TE. Note that although this merged separation zone is present in the mean flow, there are phases in the buffet cycle (high-lift phase) when there is flow reattachment (see § 4).

Figure 10. Variation of mean (time- and span-averaged) (a) pressure and (b) skin friction coefficients along suction (——) and pressure (— $\cdot$ —) surfaces for various Mach numbers.

3.2.2. High-Mach-number features and buffet offset

For $M \geqslant 0.8$, shock wave structures are present on both aerofoil surfaces. The temporal variation of $C_L$ for these cases and the PSD of the corresponding fluctuating component are shown in figure 11. The case of $M = 0.775$, for which shock wave structures are observed only on the suction side, is additionally shown for reference. Buffet offset is seen to occur at $M = 0.9$, with $C_L$ approximately a constant, although, as noted in § 2.2.1, the domain extent might not be sufficient and the results must be interpreted with caution. At this $M$ the mean lift is seen to be significantly higher than that of $M = 0.85$. The PSD at the buffet frequency (highlighted using circles) reduces by two orders of magnitude as $M$ is increased from $0.775$ to $0.8$. For $M = 0.8$, low-frequency oscillations are still present, but are accompanied intermittently by high-frequency oscillations of significant amplitude.

Figure 11. (a) Temporal variation of lift coefficient for various freestream Mach number and (b) the PSD of its fluctuating component past initial transient. Dashed vertical lines and circles highlight the buffet Strouhal number and peaks in the PSD. The high-frequency peak ($St \approx 1.7$) for $M = 0.85$ is highlighted using the square symbol.

The instantaneous flow features at the high- and low-lift phases based on the low-frequency oscillations for $M = 0.8$ are shown in figure 12. Both the suction and pressure sides have substantial regions of supersonic flow. This development of a supersonic region on the pressure side was observed with increasing $M$ first for $M = 0.775$, with a small (${<}0.1c$) pocket of supersonic flow occurring on the pressure side close to mid-chord in part of the buffet cycle, but without any shock wave present. From figures 11 and 12, we can see that although the temporal fluctuations in $C_L$ are low for this $M$, there are significant changes in the sonic edge's positions and the wall-normal extent of the supersonic region. Indeed, the qualitative flow features on the suction side resemble those seen for the other cases where buffet occurs (cf. figure 7).

Figure 12. Streamwise density gradient contours on the $x\unicode{x2013} y$ plane shown for the (a) high- and (b) low-lift phases of the buffet cycle for $M = 0.8$. The sonic line is highlighted using a grey curve.

For $M = 0.85$, as shown in figure 11(b), a peak at $St \approx 1.7$ in the PSD is present (highlighted by the square symbol) in addition to that at the buffet frequency. The former arises from a regular vortex shedding and is shown in figure 13(a) at an arbitrary instant. Note that only part of the domain extent is shown. For this $M$, a single shock wave exists at the TE on each aerofoil surface. The shock wave positions on both sides were observed not to change significantly from that shown, although buffet could still be discerned. This was also confirmed using SPOD (not shown). Strongly coherent vortices in the wake, indicating a Kármán vortex street, are apparent from the figure. By contrasting with the reference case, where large-amplitude shock wave motion is observed but the vortical structures are not dominant, we can make an important conclusion that vortex shedding does not directly influence buffet. That is, an increase in vortex shedding strength does not necessarily translate to an increase in buffet amplitude.

Figure 13. Instantaneous streamwise density gradient contours on the $x\unicode{x2013} y$ plane shown at (a) $M = 0.85$ and (b) $M = 0.9$. The sonic line is highlighted using a grey curve.

The flow features for the offset freestream Mach number, $M = 0.9$, are shown in figure 13(b). Unlike all the other cases simulated, the BL is now laminar up to the TE, although BL separation still occurs just upstream of the TE on the suction surface. Although quantitative differences exist when the domain is extended (see 2.2.1), these qualitative features were found to be the same. The increase in $C_L$ observed for $M = 0.9$ as compared with $M = 0.85$ (figure 11a) is likely due to this separation being relatively aft for the former. Vortices in the wake are also smaller. Similar results of subdued vortex shedding have been observed when $M$ is increased close to unity for transonic flows over circular cylinders (Murthy & Rose Reference Murthy and Rose1978) whereas the shock wave moving to the TE has been seen for the NACA 0012 aerofoil at zero incidence (Bouhadji & Braza Reference Bouhadji and Braza2003). A fish-tail structure was reported in the latter flow, similar to the one that can be discerned here in the aft of the aerofoil, starting from $x \approx 0.6$.

3.2.3. Buffet at zero incidence on a supercritical aerofoil

Possibly based on distinctions made previously (e.g. Lee (Reference Lee2001, p. 149), ‘There is some difference in the mechanisms of periodic shock motion between a lifting airfoil at incidence and a symmetrical one at zero incidence’ and Iovnovich & Raveh (Reference Iovnovich and Raveh2012, § III, p. 884)), buffet has been classified in Giannelis et al. (Reference Giannelis, Vio and Levinski2017) as either Type I or Type II. The former is typically associated with buffet on biconvex or symmetric aerofoils at zero incidences and is characterised by the presence of shock waves on both aerofoil surfaces. The latter is observed on supercritical aerofoils at relatively high incidence, with shock waves present only on the suction side. Based on the differences in the models for Type I buffet in Gibb (Reference Gibb1988) and Type II buffet in Lee (Reference Lee1990), it was proposed in Lee (Reference Lee2001) that the two types are sustained by distinct mechanisms. This distinction was also noted in Iovnovich & Raveh (Reference Iovnovich and Raveh2012), where the two were referred to as type one and two.

Here, for $0.7 \leqslant M \leqslant 0.775$, the buffet observed is clearly Type II. However, for $0.8 \leqslant M \leqslant 0.85$, we have buffet with shock waves appearing on both surfaces, albeit at a high incidence angle of $\alpha = 4^{\circ }$. To check if similar buffet features can be observed at zero incidence, we carried out simulations for an additional case of $M = 0.8$ and $\alpha = 0^{\circ }$. The temporal variation of $C_L$ past transients is shown in figure 14, with the reference case result also provided for comparison. The frequency of the oscillations seen was computed as $St \approx 0.14$ for the zero-incidence case, which is slightly higher than that of the reference case ($St \approx 0.12$), but lower than that of $M = 0.8$ case at $\alpha = 4^{\circ }$ ($St \approx 0.16$). The mean lift is approximately zero for the zero-incidence case despite the aerofoil not being symmetric (see table 2). The density gradient contours at the high- and low-lift phases are shown in figure 15 (see also supplementary movie 2). The spatiotemporal variations of the streamwise pressure gradient on both the suction and pressure sides of the curve C5 are shown in figure 16. Pockets of supersonic flow and shock wave structures are seen on both surfaces for parts of the buffet cycle, although in other parts of the buffet cycle the flow becomes mostly or entirely subsonic on one of the surfaces. The high-lift (red dashed line) and low-lift (blue dashed line) phases are characterised by the shock wave structures being at their most downstream position on the suction and pressure sides, respectively. Furthermore, the extent of the supersonic regions on either surface varies approximately $180^{\circ }$ out of phase with the other surface which can be inferred by comparing the sonic edge's positions at the high- and low-lift phases in figure 16. This suggests that the present case is of a Type I buffet, albeit on a supercritical aerofoil. This is further explored using SPOD in § 4.3.

Figure 14. Temporal variation of lift coefficient for buffet at zero and non-zero incidence angles.

Figure 15. Streamwise density gradient contours on the $x\unicode{x2013} y$ plane shown for the (a) high- and (b) low-lift phases of the buffet cycle for $M = 0.8$ at $\alpha = 0^{\circ }$. The sonic line is highlighted using a grey curve.

Figure 16. Spatiotemporal variation of streamwise pressure gradient on the (a) suction and (b) pressure sides of curve C5. The sonic line is highlighted using black curves. The approximate times associated with different phases of interest, high lift (red), low skin-friction drag (green), low lift (blue) and high skin-friction drag (brown), are highlighted using dashed lines.

4.1.1. Temporal features

A comparison of the spectra for the dominant eigenvalue, $\lambda _1$, obtained through SPOD for the reference case and one case selected from each parameter variation is shown in figure 23(a). For all cases, we see peaks at a fundamental frequency of $St \approx 0.1$, i.e. the buffet frequency and, in addition, peaks at its harmonics. Significant energy content is also present in a medium-frequency range of $1 \leqslant St \leqslant 4$ associated with what we refer to as wake modes. These are shown in the following to be associated with vortex shedding (e.g. see figure 24c). A bump in the spectrum similar to that seen here and associated with vortex shedding was reported for turbulent buffet in the DES of the OAT15A aerofoil in Grossi et al. (Reference Grossi, Braza and Hoarau2014) and organised eddy simulations of Szubert et al. (Reference Szubert, Grossi, Jimenez Garcia, Hoarau, Hunt and Braza2015). We have also extracted a similar mode at $St \approx 1.87$ for the reference case using DMD in our previous study (Zauner & Sandham Reference Zauner and Sandham2020a). Note that modes similar to the wake modes, and accompanying turbulent buffet, were also reported in Sartor et al. (Reference Sartor, Mettot and Sipp2015). This was shown for the OAT15A using resolvent analysis (see their figures 16 and 18) with a frequency range $0.95 \leqslant St \leqslant 4$ which is the same as the range of $St$ seen here. Although not frequently reported in experiments, Szubert et al. (Reference Szubert, Grossi, Jimenez Garcia, Hoarau, Hunt and Braza2015) have shown that phase-averaging implemented in experiments to capture buffet features can prevent the effects of such modes from being observed. From the figure, it is evident that a change in any parameter ($M/\alpha /\textit {Re}$) leads to an increase in the energy content of the wake modes, indicating a stronger vortex shedding beyond buffet onset. With increasing $\alpha$, the frequency associated with these modes seems to be shifted to lower values.

Figure 23. (a) Eigenvalue spectra (logarithmic scale) based on SPOD of the dominant eigenvalue for the reference case and a typical case for each parameter varied. Circles highlight the buffet peaks. (b) Scalogram based on the lift coefficient ($\log _{10} |W(C_L')|$) with the temporal variation of the same overlaid (black curve) for the case of $\alpha = 6^{\circ }$.

Figure 24. SPOD modes for the reference case shown using contour plots of real part of the density field: (a) buffet mode ($St = 0.11$), (b) its first harmonic ($St = 0.23$) and (c) a wake mode ($St = 2.1$). The red, black and blue curves represent the sonic lines based on the high-lift phase, mean flow and low-lift phase, respectively.

To check for temporal variations in the intensities and frequencies of these modes, a time-frequency analysis using the continuous wavelet transform was also carried out (see § 2.1.5 for details). The scalogram based on the transform of the lift coefficient's fluctuating component, $W(C_L')$, is shown for the case of $\alpha = 6^{\circ }$ in figure 23(b). Note that the contours are based on the logarithm of the magnitude, i.e. $\log _{10} |W(C_L')|$. The temporal variation of $C_L$ is also overlaid to identify high- and low-lift phases. A horizontal band in the contour at the buffet frequency ($St \approx 0.1$) confirms that it does not vary with time. However, it is also evident that there are periodic changes of $|W(C_L')|$ in the frequency range $1 \leqslant St \leqslant 3$. As seen from figure 23(a), this is associated with the wake modes for this case. It is evident from the figure that the intensity of these modes is approximately a maximum in the low-lift phase and vice versa. This is expected based on the vortex shedding behaviour reported in § 3.3 (see figure 18). Similar results have been reported for turbulent buffet in Szubert et al. (Reference Szubert, Grossi, Jimenez Garcia, Hoarau, Hunt and Braza2015).

4.1.2. Spatial features

The spatial structure of the SPOD modes of importance are shown for the reference case in figure 24 using the real part of the density field. From left to right, these are the buffet mode ($St = 0.11$), its first harmonic ($St = 0.23$) and a typical wake mode ($St = 2.1$), where the bump in the spectra attains a local maximum). The sonic line based on the mean local Mach number (black curve) is shown for reference along with the sonic lines at high- (red) and low-lift (blue) phases (reconstructions, see § 4.2.1). For the buffet mode we see that strong density fluctuations are present in the vicinity of the TE and near-wake which are out of phase (i.e. opposite sign) with those near the shock foot. Thus, the density increases at the shock foot when it reduces at the TE and in the wake and vice versa. The harmonic mode is of a shorter wavelength relative to the buffet mode. The mode shapes of the buffet mode and its harmonic reported here are similar to those in Poplingher et al. (Reference Poplingher, Raveh and Dowell2019) for turbulent buffet (see their figure 9).

The wake mode has a structure similar to that of the von Kármán vortex street. Similar structures were also observed in the resolvent analysis of steady RANS solutions in Sartor et al. (Reference Sartor, Mettot and Sipp2015) where they are referred to arise as Kelvin–Helmholtz instabilities. Note that the flow field contains both a free shear layer due to a separated BL and a wake. The former can lead to a Kelvin–Helmholtz inflectional instability. On the other hand, a von Kármán instability can arise in the wake and studies such as Triantafyllou & Karniadakis (Reference Triantafyllou and Karniadakis1990) and Pier & Huerre (Reference Pier and Huerre2001) have shown that this instability can arise even when only the wake is considered in isolation (i.e. without the inclusion of the solid body and upstream regions). Although the wake mode's structure resembles a von Kármán vortex street (see also figure 8), a stability analysis is required to understand the origin of the instability that leads to this mode which is beyond the scope of this study.

In addition to this, waves in the flow field that seem to be generated at the TE and propagate upstream outside the BL on both sides of the aerofoil can be identified. These appear to be the Kutta waves suggested in Lee (Reference Lee1990). Based on visualisations of the time evolution of the wake mode, the upper limit on the time required for these waves to reach the shock wave starting from the TE was estimated as $t_{up} < 1.5$. Thus, for the reference case, Lee's model predicts the buffet time period as $\tau _{Lee} = t_{up}+t_{down} < 2.1$, which is substantially lower than the actual buffet period of $\tau _{LES} =8$ observed in the LES (see § 3.1). Similar results were observed for other wake modes associated with the spectral bump and for all other cases where buffet was observed, suggesting that Lee's model is invalid. This is further considered in § 5.3.

The modification to Lee's model proposed in Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009) assumes that $t_{up}$ should be replaced with $t^{J}_{up}$, the time required for the waves to travel upstream along the pressure surface and turn around the leading edge to reach the shock foot. Such upward propagating waves are also seen in figure 24(c) on the pressure side of the aerofoil. However, we observed that they reduce in intensity by more than an order of magnitude upstream of mid chord and they are unidentifiable past the leading edge in the supersonic region, likely because their intensity is reduced to levels similar to that of noise in the SPOD modes. By using the same approximation as Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009) that the time spent in the supersonic region is negligible, we observed $t^{J}_{up} < 3$ implying that the predicted time period is $\tau _J = t^{J}_{up}+t_{down} < 3.6$ which still remains significantly lower than $\tau _{LES}$.

Features of the corresponding SPOD modes for $\alpha = 6^{\circ }$ are shown in figure 25. We see that the wave number associated with each mode is substantially lower with all modes showing a larger spatial structure. A relatively strong separation in the vicinity of the shock foot is apparent, whereas the wake mode's energy is high at the TE. These results match with the flow features observed in the simulations (cf. figures 3 and 18). Similar to the reference case, upstream propagating waves emanating from the TE can be observed for the wake mode. Based on this, it was estimated that $t_{up} < 2$ and $t^{J}_{up} < 3$ implying that $\tau _{Lee} < 3.1$ and $\tau _J < 4.1$ compared with $\tau _{LES} = 9.9$.

Figure 25. SPOD modes for $\alpha = 6^{\circ }$ shown using contour plots of real part of the density field: (a) buffet mode ($St = 0.10$), (b) its first harmonic ($St = 0.21$) and (c) a wake mode ($St = 1.3$). The red, black and blue curves represent the sonic lines based on the high-lift phase, mean flow and low-lift phase, respectively.

4.2.1. Buffet mode

Contours of axial velocity at different phases are shown in figure 26 for $\alpha = 6^{\circ }$. For this case, the phases of $\phi = 73^{\circ }$ and $253^{\circ }$ represent the minimum and maximum skin friction drag, respectively. In contrast to the actual flow field in which multiple shock structures are present, the reconstructed flow field contains only a single supersonic region enclosed by the sonic line. This is expected, because these multiple shock structures arise at a frequency higher than the buffet frequency (see, e.g., figure 4). For convenience, the approximate aft part of the sonic line that terminates the supersonic region is referred to as a shock wave. The BL remains attached over most of the aerofoil when the shock wave is at its most downstream position ($\phi = 0^{\circ }$). As the shock wave moves upstream, we see that at $\phi = 73^{\circ }$, the flow is fully separated beyond the foot of the shock wave ($\tilde {u}_x < 0$). As the shock wave reaches its most upstream position ($\phi = 180^{\circ }$), the BL is attached until $x \approx 0.7$. We found the qualitative features described previously to be common to all cases considered.

Figure 26. Reconstructed flow field based on the buffet mode for $\alpha = 6^{\circ }$ shown using axial velocity contour at (a) high-lift, (b) low-skin-friction-drag, (c) low-lift and (d) high-skin-friction-drag phases. The sonic line is highlighted using a black curve.

To better understand the dynamical behaviour of buffet, it is useful to examine the variation of flow features on the aerofoil surface at different phases. This spatiotemporal variation is reconstructed based on the data set Data-SpAv (see § 2.3.2) and shown using $x'-\phi$ diagrams for the different cases in figure 27 (right-hand side contours). Unlike the preceding sections where the streamwise pressure gradient on the curve C5 was used to scrutinise shock wave structures, we show here the variation of the pressure coefficient on the aerofoil surface. In addition, we overlay the contour line, $\tilde {C}_f = 0$ (green curve), which delineates flow reversal ($\tilde {C}_f \leqslant 0$). To include the shock wave position, we also overlay the sonic line based on the local Mach number (black curve). We emphasise that only this quantity is computed based on the curve C5 whereas the rest are based on the aerofoil surface. As before, the phases at which minimum and maximum $(\tilde {C}_D)_f$ occur are highlighted using dashed horizontal lines.

Figure 27. Temporal variation of the shock wave's strength shown using $\tilde {M}_{eff}$ (left-hand side subplot) and spatiotemporal contours of $\tilde {C}_p$ on the aerofoil suction surface (right-hand side contours) for (a) reference, (b) $M = 0.735$, (c) $\alpha = 6^{\circ }$ and (d) $\textit {Re} = 1.5\times 10^{6}$ cases. The isolines for $\tilde {M}_{loc} = 1$ (solid black curve) and $\tilde {C}_f = 0$ (solid green curve) and the phases associated with maximum (dashed brown line) and minimum (dashed green line) $(\tilde {C}_D)_f$ are also shown for reference.

The shock wave's most upstream and downstream positions (black curve) occur when the lift is close to its minimum ($\phi = 180^{\circ }$) and maximum ($\phi = 0^{\circ }$), although there seems to be a small phase lag between the two. As alluded to in § 3.1, the maximum and minimum $(\tilde {C}_D)_f$ occur when the flow is ‘least’ and ‘most’ separated, as can be inferred from the chordwise extent of separation delineated by the $\tilde {C}_f = 0$ isoline (green curve). Indeed, the BL remains separated up to the TE ($\tilde {C}_f < 0$) at the phase when the minimum $(\tilde {C}_D)_f$ is attained (dashed horizontal green line) for all cases. By contrast, during the phase of maximum $(\tilde {C}_D)_f$ (dashed horizontal brown line), the BL is fully attached from the leading edge to TE for the cases $\textit {Re} = 1.5\times 10^{6}$ and $\alpha = 6^{\circ }$, whereas the flow reattaches ($\tilde {C}_f > 0$) downstream to the shock wave (black curve) and remains attached up to the TE for the other cases. Thus, we can conclude that the surface coefficients $\tilde {C}_L$ and $(\tilde {C}_D)_f$ are approximate indicators of the shock wave position and the BL separation, respectively. The reasons for these relations can be explained by observing that a negative $\tilde {C}_p$ exists upstream of the shock wave, and an increase in the chordwise extent of this low-pressure region leads to an increased lift. Similarly, as mentioned previously, an attached BL leads to larger skin friction due to large positive velocity gradients on the aerofoil surface, whereas a separated BL and flow reversal above the surface contribute to a low $(\tilde {C}_D)_f$.

Another aspect of interest is the temporal variation of the shock wave's strength as it is expected to play an important role in BL separation and, thus, on buffet. One estimate of this strength is the ratio of the instantaneous pressure immediately downstream and upstream of the sonic line in the $x-\phi$ diagram, $\tilde {p}(x^{+}_S,\phi )/\tilde {p}(x^{-}_S,\phi )$, where $x_S$ represents the shock wave's streamwise position (i.e. $x_S = x(\tilde {M}_{loc} = 1)$) and the symbols ‘$+$’ and ‘$-$’ indicate downstream and upstream, respectively. However, this choice was found to be sensitive to the choice of $x^{+}_S$ and $x^{-}_S$. This is possibly because the actual flow field has multiple shock wave structures and the single shock wave seen in the modal reconstruction is only an approximation, implying that choosing points both upstream and downstream of the shock introduces relatively large errors in the estimation of the pressure ratio. Instead, we found the upstream effective Mach number, which is based on the instantaneous upstream velocity in the inertial frame of reference moving at a speed equal to the instantaneous speed of the shock wave (i.e. $\textrm {d}\kern0.06em x_S/\textrm {d}\phi$), to be more robust with regard to the choice of $x^{-}_S$. This is given by

(4.1)\begin{equation} \tilde{M}_{eff}(x^{-}_S,\phi) = \tilde{M}_{loc}(x^{-}_S,\phi) - \left(\frac{{\rm d}\kern0.06em x_S}{{\rm d}t}\right)\frac{1}{a(x^{-}_S)}, \end{equation}

where $a$ represents the local speed of sound. The effective Mach number can then be used to represent the shock strength based on the Rankine–Hugoniot conditions (Gibb Reference Gibb1988). However, note that this approximation neglects the effect of shock wave acceleration. The left-hand side subplots in figure 27 show the effective upstream Mach number's temporal variation for $x_S^{-} = x_S - 0.1$. It is seen from the plots that $\tilde {M}_{eff}$ reaches a maximum or minimum when the shock is approximately at the mid-point of its upstream or downstream excursion, respectively. This implies that the shock wave is always strongest/weakest close to its mean streamwise position during its upstream/downstream motion. A similar result based on Mach number was shown in Fukushima & Kawai (Reference Fukushima and Kawai2018) (their figure 18), and other studies have also noted the same trend (Tijdeman & Seebass Reference Tijdeman and Seebass1980; Lee Reference Lee2001; Iovnovich & Raveh Reference Iovnovich and Raveh2012; Hartmann et al. Reference Hartmann, Feldhusen and Schröder2013). The plots show that there is a small interval in which $\tilde {M}_{eff} < 1$. This is possibly not of physical significance, because there are shock wave structures present at all times suggesting that it arises due to the approximations in the estimate of shock strength.

From figure 27, the following sequence of events can be inferred for all cases.

1. (i) When the shock wave is at its most downstream position, the BL is separated upstream of the shock wave, but reattaches downstream, forming a separation bubble.

2. (ii) As the shock moves towards the leading edge, it strengthens. The upstream separation point moves in the same direction and the reattachment point moves further downstream (or vanishes), increasing the chordwise extent of the separation region.

3. (iii) The BL is ‘most’ separated when the shock wave is approximately at the mid-point of its upstream excursion, where its strength is close to its maximum.

4. (iv) Past the mid-point, the BL separation weakens leading to the formation of a separation bubble that reattaches at the shock foot for all cases except $M = 0.735$, where this occurs later.

5. (v) Once the shock wave reaches its most upstream position and starts to move downstream, the separation point upstream of the shock wave abruptly moves downstream by almost $\delta x' = 0.1$ for the $M = 0.735$ and reference case, whereas it completely vanishes for the other cases.

6. (vi) As the shock wave reaches the midpoint of its downstream excursion it strength is relatively weakened and the flow is ‘least’ separated.

7. (vii) The reattached flow is associated with a pressure rise in the TE region, with this back pressure exceeding the freestream pressure (red contours, $\tilde {C}_p > 0 \Rightarrow \tilde {p} > \tilde {p}_\infty$). The chordwise extent of this region of increased pressure attains a maximum before the shock wave reaches its most downstream position.

The above description indicates that the laminar buffet observed in this study is essentially related to a moving shock wave of temporally varying strength whose motion is accompanied by BL separation and reattachment. Importantly, the following two features are common to all buffet cases studied: a phase lag between the shock wave position and the separation extent and the build up of back pressure. In addition, it is interesting to note that, for all cases, there are phases in the buffet cycle where the flow separation extends up to the TE, indicating that the aerofoil is intermittently stalled. However, we emphasise that it is difficult to disentangle cause, correlation and effect. The possible implications of these observations are discussed in § 5.

4.2.2. Wake mode

The contours of reconstructed axial velocity field based on the wake mode for the case of $\alpha = 6^{\circ }$ are shown in figure 28 at phases corresponding to the highest ($\phi = 0^{\circ }$) and lowest ($\phi = 180^{\circ }$) values of $\tilde {C}_L$. Note that the phase, $\phi$, considered here is based on the time period of the wake mode and not the buffet mode. Thus, the high-lift phase signifies the phase at which the wake mode interaction with the mean flow induces the maximum lift obtained for the reconstructed flow field and is not related to the buffet mode. In the plots, the sonic line (black curve) indicates that the shock foot extends and retracts at different phases, similar to the shock foot motion observed in Dandois et al. (Reference Dandois, Mary and Brion2018) for the OALT25 at $M = 0.735$, $\alpha = 4^{\circ }$ and $Re = 3\times 10^{6}$, implications of which are discussed in § 5.1. The streamwise positions of the sonic edge at $\phi = 0^{\circ }$ and $180^{\circ }$ are $x \approx 0.6$ and 0.5, respectively, implying that the shock foot moves approximately a streamwise distance of 10 % chord. This indicates that the wake mode can cause significant variations in the flow field. However, these are limited to the shock foot as the approximate wall-normal height up to which differences in the sonic line at these two phases are significant is only around 10 % (contrast with figures 26(a) and 26(c) for the reconstruction based on the buffet mode, where the streamwise distance moved is 30 % chord with the sonic line changing drastically in the wall-normal direction also).

Figure 28. Reconstructed flow field based on the wake mode for $\alpha = 6^{\circ }$ shown using axial velocity contour at (a) high- and (b) low-lift phases. The sonic line is highlighted using a black curve.