1 Introduction
Transonic flow over an airfoil exhibits complex dynamics from the interplay among the unsteady shock motion, flow separation due to the shock/boundary layer interaction, and the generation and propagation of pressure waves. Under certain ranges of angle of attack and Mach number, the shock wave oscillates vigorously, creating a highly unsteady flow (Lee Reference Lee2001). This shock-wave oscillation phenomenon, known as transonic buffet, is caused by a self-sustaining aerodynamic instability. The large-scale motion of the shock induces violent fluctuations in aerodynamic forces and leads to structural vibrations. The level of unsteadiness induced by transonic buffet can be so large that it can compromise the integrity of aircraft structure and flight safety. For these reasons, transonic buffet is a major limiting factor for aircraft flight envelopes and has been the focus of many studies, as summarized recently by Giannelis, Vio & Levinski (Reference Giannelis, Vio and Levinski2017).
Various experimental and numerical studies have been devoted to the characterization of transonic buffet. Past studies have revealed that the typical oscillating frequency of shock waves is very low, with the chord-based Strouhal number being approximately 0.06 (Deck Reference Deck2005; Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009). This low frequency is distinguished from other unsteady phenomena around the airfoil related to the velocity fluctuations in the separation region and vortex shedding in the wake (Sartor, Mettot & Sipp Reference Sartor, Mettot and Sipp2014). More recently, some researchers have reported that a transonic buffet on a three-dimensional wing has spanwise instability with a much higher frequency, where 
$0.2<St<0.6$ (Dandois Reference Dandois2016; Ohmichi, Ishida & Hashimoto Reference Ohmichi, Ishida and Hashimoto2018).
Separation on the airfoil induced by shock-wave/boundary layer interaction (SWBLI) plays an important role in transonic buffet. Early experiments (Levy Reference Levy1978; Seegmiller, Marvin & Levy Jr. Reference Seegmiller, Marvin and Levy1978) showed that the onset of shock oscillation occurred when the separation generated at the shock foot reaches the trailing edge of the airfoil. Recent numerical studies (Iovnovich & Raveh Reference Iovnovich and Raveh2012; Grossi, Braza & Hoarau Reference Grossi, Braza and Hoarau2014; Fukushima & Kawai Reference Fukushima and Kawai2018) reported that unsteadiness of the separation is closely tied to buffet. When the shock moves downstream from its upstream position, the shock becomes weaker and the flow behind the shock remains attached. In contrast, when the shock moves upstream from its downstream location, the shock becomes stronger, and causes large-scale separation.
A number of studies have analysed the self-sustained mechanism of shock-wave oscillation. Lee (Reference Lee2001) argued that the pressure wave propagating between the shock wave and the trailing edge of the airfoil plays a critical role in feedback to maintain shock oscillation. More recently, global linear stability analysis performed by Crouch et al. (Reference Crouch, Garbaruk, Magidov and Travin2009) extracted the growth rate of the buffet fluctuation with respect to the time-averaged flow field based on the eigenvalue problem of a linearized Navier–Stokes operator. They examined the linear stability of transonic flow around a two-dimensional NACA 0012 airfoil obtained from Reynolds-averaged Navier–Stokes (RANS) simulations. Analysing the critical angle of attack and Mach number of the buffet onset, the critical values agreed reasonably with the experimental results obtained by McDevitt & Okuno (Reference McDevitt and Okuno1985). Nonetheless, the origin of the transonic buffet is still under debate, with open questions on the validity of the linear stability formalism for base flows that are not the exact solution of the Navier–Stokes equations.
Recently, Sartor et al. (Reference Sartor, Mettot and Sipp2014) employed global stability, adjoint, and resolvent analysis to examine transonic buffet flow around a NACA 0012 airfoil at 
$Re=3\times 10^{6}$. At the buffet frequency, they used resolvent analysis to find that the essential origin of shock unsteadiness is at the shock foot region on the suction side of the airfoil. There are also reports, including Nitzsche (Reference Nitzsche2009), that identified various sources that can trigger buffet, such as the motion of flaps, pitch oscillation, and translation of the airfoil. Based on these observations, various methods have been developed to control buffet, including vortex generators (McCormick Reference McCormick1993) and trailing-edge deflection techniques (Gao et al. Reference Gao, Zhang, Kou, Liu and Ye2017).
The objective of the present study is to systematically determine the origin of two-dimensional transonic buffet flow around a NACA 0012 airfoil at a low Reynolds number of 
$Re=2000$ in the absence of any modelling terms. We reveal the input–output relationship for transonic buffet with the resolvent analysis, which can extract not only the frequency response but also identify the forcing (input) and response (output) modes of flow systems (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017, Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2019; Yeh & Taira Reference Yeh and Taira2019). To validate that the identified source of buffet is indeed the driving mechanism for sustained oscillations, we perform direct numerical simulations (DNS) to show that buffet can appear even in low-Reynolds-number flows under appropriate conditions. This analysis not only provides the fundamental insights into transonic buffet but also serves as a building block for low-Reynolds-number compressible aerodynamics (
$Re\approx O(10^{3}{-}10^{5})$). This area of aerodynamics has traditionally been overlooked but is now becoming important for developing high-efficiency wings and propulsion systems for unmanned aircraft on Mars (Anyoji et al. Reference Anyoji, Numata, Nagai and Asai2015; Munday et al. Reference Munday, Taira, Suwa, Numata and Asai2015; Koning, Johnson & Grip Reference Koning, Johnson and Grip2019).
We discuss the problem description, simulation approach and resolvent analysis in § 2. The insights gained from resolvent analysis are offered in § 3 with discussions on the source of buffet. Our findings are validated by forcing the flow at the identified source to stimulate the emergence of buffet. Finally, we offer concluding remarks in § 4.
2 Problem set-up
2.1 Problem description
We consider two-dimensional laminar transonic flows over a NACA 0012 airfoil at a free-stream Mach number 
$M_{\infty }\equiv U_{\infty }/a_{\infty }=0.85$, a chord-based Reynolds number of 
$Re_{L_{c}}\equiv U_{\infty }L_{c}/\unicode[STIX]{x1D708}_{\infty }=2000$ and an angle of attack of 
$\unicode[STIX]{x1D6FC}=3^{\circ }$. Here, 
$U_{\infty }$ is the free-stream velocity, 
$a_{\infty }$ is the free-stream sonic speed, 
$L_{c}$ is the chord length and 
$\unicode[STIX]{x1D708}_{\infty }$ is the free-stream kinetic viscosity. The Prandtl number is set to 
$Pr=0.7$. While not reported, we also considered 
$\unicode[STIX]{x1D6FC}=1^{\circ }$ and 
$5^{\circ }$. As we obtained analogous results for these angles of attack, we present only the representative case of 
$\unicode[STIX]{x1D6FC}=3^{\circ }$ herein. Shock waves form around the airfoil under these conditions, making this base flow an appropriate candidate for this study.
2.2 Flow simulation
We simulate transonic flows over the airfoil by solving the compressible Navier–Stokes equations in a two-dimensional setting using the finite-volume solver CharLES, which is second-order accurate in space and third-order accurate in time (Brès et al. Reference Brès, Ham, Nichols and Lele2017). We employ a second-order fully unstructured essentially non-oscillatory (ENO) method (Shi, Hu & Shu Reference Shi, Hu and Shu2002) for shock capturing. A C-shaped hexahedral mesh is used for the simulations, as shown in figure 1. The computational domain has an extent of 
$x_{c}/L_{c}\in [-50,50]$ and 
$y_{c}/L_{c}\in [-50,50]$, where 
$x_{c}$ and 
$y_{c}$ represent the chordwise and chord-normal directions, respectively. The airfoil is positioned with its leading edge at the origin. The domain is discretized into 70 000 cells with 100 nodes on each side of the airfoil, 90 nodes along the wake and 80 nodes in the normal direction, where the minimum wall-adjacent 
$\unicode[STIX]{x0394}y$ is set to be 
$\unicode[STIX]{x0394}y/L_{c}=1.42\times 10^{-3}$. A grid convergence study has been conducted with respect to the time-averaged flow to ensure sufficient spatial resolution for this computational grid. At the far-field boundary, the free-stream condition is prescribed as 
$[\unicode[STIX]{x1D70C},v_{x_{c}},v_{y_{c}},T]=[\unicode[STIX]{x1D70C}_{\infty },U_{\infty }\cos \unicode[STIX]{x1D6FC},U_{\infty }\sin \unicode[STIX]{x1D6FC},T_{\infty }]$, where 
$\unicode[STIX]{x1D70C}$ is density and 
$T$ is temperature. The free-stream angle is changed through the prescription of the far-field flow velocities 
$v_{x_{c}}$ and 
$v_{y_{c}}$ in the 
$x_{c}$ and 
$y_{c}$ directions, respectively. A no-slip adiabatic condition is prescribed over the airfoil. Along the outlet boundary, a sponge layer (Freund Reference Freund1997) is applied over 
$x_{c}/L_{c}\in [40,50]$, with the target state being the running-averaged flow over 
$tU_{\infty }/L_{c}=1$. The time integration is performed at a constant Courant–Friedrichs–Lewy number of 1.
Figure 1. The computational grid used for numerical simulation and resolvent analysis. (a) The whole computational domain. The enlarged views around the airfoil are shown with the (b) instantaneous magnitude of normalized vorticity and (c) time-averaged velocity fields.
2.3 Resolvent analysis
The objective of the present study is to reveal the input–output relationship for two-dimensional transonic buffet and identify the origin of the energetic shock-wave oscillations. In this study, we focus our resolvent analysis on two-dimensional modes, which are essential in all buffet phenomena (Iovnovich & Raveh Reference Iovnovich and Raveh2012; Sartor et al. Reference Sartor, Mettot and Sipp2014; Ohmichi et al. Reference Ohmichi, Ishida and Hashimoto2018). Note that resolvent analysis itself is not limited to two-dimensional flow. Further resolvent analysis of transonic buffet over three-dimensional wings could extract the spanwise modes, which have been reported in previous studies (Iovnovich & Raveh Reference Iovnovich and Raveh2015; Ohmichi et al. Reference Ohmichi, Ishida and Hashimoto2018).
We consider the Reynolds decomposition of the flow variable 
$\boldsymbol{q}(\boldsymbol{x},t)=[\unicode[STIX]{x1D70C},v_{x},v_{y},T]$ into the sum of a time-averaged base state 
$\bar{\boldsymbol{q}}(\boldsymbol{x})$ and the statistically stationary fluctuating component 
$\boldsymbol{q}^{\prime }(\boldsymbol{x},t)$, which enables us to express the Navier–Stokes equation as
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{q}^{\prime }=\unicode[STIX]{x1D647}_{\bar{\boldsymbol{ q}}}\boldsymbol{q}^{\prime }+\unicode[STIX]{x1D63D}\boldsymbol{f}^{\prime },\end{eqnarray}$$ where 
$\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ is the linearized Navier–Stokes operator about the mean flow 
$\bar{\boldsymbol{q}}$ (McKeon & Sharma Reference McKeon and Sharma2010) and 
$\unicode[STIX]{x1D63D}$ is an input matrix that can spatially window the forcing input 
$\boldsymbol{f}^{\prime }$. We construct the discrete linear operator 
$\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$, incorporating the boundary conditions of 
$[\unicode[STIX]{x1D70C}^{\prime },v_{x_{c}}^{\prime },v_{y_{c}}^{\prime },\unicode[STIX]{x1D735}_{n}T^{\prime }]=0$ over the airfoil and the far field, and 
$\unicode[STIX]{x1D735}_{n}[\unicode[STIX]{x1D70C}^{\prime },v_{x_{c}}^{\prime },v_{y_{c}}^{\prime },T^{\prime }]=0$ at the computational outlet, where 
$\unicode[STIX]{x1D735}_{n}$ denotes the surface-normal gradient. The spatial discretization for 
$\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ is performed on the same mesh used in the flow simulation, as shown in figure 1. When the time-averaged flow is chosen as the base state in resolvent analysis, the finite-amplitude nonlinear terms with respect to 
$\boldsymbol{q}^{\prime }$ are incorporated in 
$\boldsymbol{f}^{\prime }$, which can be interpreted as a sustained forcing input within the natural feedback system (McKeon & Sharma Reference McKeon and Sharma2010). Moreover, we consider an observable output vector 
$\boldsymbol{y}$ given by 
$\boldsymbol{y}=\unicode[STIX]{x1D63E}\boldsymbol{q}^{\prime }$, such that spatial windowing 
$\unicode[STIX]{x1D63E}$ can be applied in general (Jeun, Nichols & Jovanović Reference Jeun, Nichols and Jovanović2016; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018).
With the Fourier representation 
$\boldsymbol{q}^{\prime }(\boldsymbol{x},t)=\int _{-\infty }^{\infty }\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D714}}(\boldsymbol{x})\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}\,\text{d}\unicode[STIX]{x1D714}$, equation (2.1) can be expressed in frequency space as
$$\begin{eqnarray}-\text{i}\unicode[STIX]{x1D714}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D714}}+\unicode[STIX]{x1D63D}\hat{\boldsymbol{f}}_{\unicode[STIX]{x1D714}},\quad \text{with }\hat{\boldsymbol{y}}_{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D63E}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D714}}.\end{eqnarray}$$ The input–output relationship between the input (forcing) 
$\hat{\boldsymbol{f}}_{\unicode[STIX]{x1D714}}$ and the output (response) 
$\hat{\boldsymbol{y}}_{\unicode[STIX]{x1D714}}$ at a specified frequency 
$\unicode[STIX]{x1D714}$ becomes
$$\begin{eqnarray}\hat{\boldsymbol{y}}_{\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})\hat{\boldsymbol{f}}_{\unicode[STIX]{x1D714}},\quad \text{where }\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})=\unicode[STIX]{x1D63E}[-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}]^{-1}\unicode[STIX]{x1D63D}.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})$ is referred to as the resolvent (operator). It serves as a transfer function that amplifies (or attenuates) the harmonic forcing input 
$\hat{\boldsymbol{f}}_{\unicode[STIX]{x1D714}}$ and maps it to the response 
$\hat{\boldsymbol{y}}_{\unicode[STIX]{x1D714}}$. In this study, the base flows are found to be stable according the eigenvalues of the linear operators 
$\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$. Therefore, we consider a real-valued frequency 
$\unicode[STIX]{x1D714}$ in constructing 
$\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})$ for our resolvent analysis (Jovanović Reference Jovanović2004; Yeh & Taira Reference Yeh and Taira2019).
Resolvent analysis identifies the dominant directions along which 
$\hat{\boldsymbol{f}}_{\unicode[STIX]{x1D714}}$ can be most amplified through 
$\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})$ to form the corresponding response 
$\hat{\boldsymbol{y}}_{\unicode[STIX]{x1D714}}$. This is accomplished by performing the singular value decomposition (SVD) of the resolvent
$$\begin{eqnarray}\unicode[STIX]{x1D64D}_{\bar{\boldsymbol{q}}}(\unicode[STIX]{x1D714})=\unicode[STIX]{x1D654}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D641}^{\ast },\end{eqnarray}$$ where 
$\unicode[STIX]{x1D641}^{\ast }$ denotes the Hermitian of 
$\unicode[STIX]{x1D641}$. Resolvent analysis interprets left and right singular vectors 
$\unicode[STIX]{x1D654}=[\,\hat{\boldsymbol{y}}_{1},\hat{\boldsymbol{y}}_{2},\ldots ,\hat{\boldsymbol{y}}_{m}]$ and 
$\unicode[STIX]{x1D641}=[\,\hat{\boldsymbol{f}}_{1},\hat{\boldsymbol{f}}_{2},\ldots ,\hat{\boldsymbol{f}}_{m}]$, respectively, as response modes and forcing modes, with the magnitude-ranked singular values 
$\unicode[STIX]{x1D72E}=\text{diag}(\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{m})$ being the amplification (gain) for the corresponding forcing–response pair at frequency 
$\unicode[STIX]{x1D714}$. We examine the modal structures of the dominant forcing and response modes 
$[\,\hat{\boldsymbol{f}}_{1},\hat{\boldsymbol{y}}_{1}]$ to study the origin of transonic buffet over this canonical airfoil.
The application of a spatial window to the global response 
$\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D714}}$ can be implemented by specifying 
$\unicode[STIX]{x1D63E}$ as a diagonal weight matrix with unit weights inside the window and zeros otherwise. This window is chosen to highlight the location over the suction surface where oscillatory motion of the shock wave appears and then limits the output to be within inside of window, as is illustrated in figure 3. While matrix 
$\unicode[STIX]{x1D63D}$ for forcing can be designed in a similar manner, we use 
$\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D644}$ such that there is no restriction on the region of forcing to identify the source of buffet. The use of window 
$\unicode[STIX]{x1D63E}$ in the resolvent analysis can help reveal the optimal energy amplification from forcing to a local response in the shock region over the suction surface of the airfoil.
3 Results
3.1 Base flow
We simulate the unsteady two-dimensional flow over a NACA 0012 airfoil at 
$\unicode[STIX]{x1D6FC}=3^{\circ }$ for 
$M_{\infty }=0.85$ and 
$Re_{L_{c}}=2000$. An instantaneous vorticity field from DNS is shown in figure 2(a). The laminar vortex sheets generated from the pressure and suction sides of the wing roll into vortices in the wake. This representative transonic flow at low Reynolds number exhibits a distinct shedding frequency of 
$St\equiv fL_{c}/U_{\infty }=1.0$. This frequency can be detected from the power spectral density (PSD) of the velocity probe at 
$(x,y)/L_{c}=(2.10,-0.11)$, along with its harmonics, in figure 2(b). For this flow, we do not directly observe buffet over the airfoil. Nitzsche (Reference Nitzsche2009) investigated the frequency response of shock waves to small perturbations and found that the external forcing to transonic flow around an airfoil can trigger shock-wave oscillations even if the base flow is stable. He has reported that the maximum gain was observed close to the buffet onset angle. The present resolvent analysis holds great potential to reveal the optimal input–output relationship with respect to the onset of buffet.
Figure 2. Unsteady transonic flow over a NACA 0012 airfoil. (a) Instantaneous vorticity field (
$\unicode[STIX]{x1D6FA}_{z}L_{c}/a_{\infty }$). (b) Power spectral density of 
$v_{y}/a_{\infty }$ probed at 
$(x,y)/L_{c}=(2.10,-0.11)$.
To perform resolvent analysis, we consider the time-averaged flow to be the base state. The temporal averaging is performed over 
$tU_{\infty }/L_{c}\in [0,350]$ to ensure statistical convergence. This long window for averaging allows us to set the minimum frequency in the resolvent analysis at 
$2.9\times 10^{-3}$, which is sufficiently small compared with the dominant buffet frequency of 
$St\approx 0.07$ (Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009; Sartor et al. Reference Sartor, Mettot and Sipp2014). The time-averaged streamwise velocity and the numerical schlieren fields are shown in figure 3. The solid and broken curves in figures 3(a) and 3(b) indicate the contour lines for 
$M=1$ and 
$v_{x}=0$, respectively. The area enclosed in the solid curve corresponds to the supersonic region, while the region surrounded by the broken curves represents the separated flow region. Weak shock waves are observed in the supersonic regions. The supersonic region on the suction side is further divided into two stages of acceleration–deceleration phases according to the 
$\unicode[STIX]{x1D706}$-type structure of the shock waves (Delery Reference Delery1985). The first stage corresponds to the initial flow acceleration near the leading edge. The flow surpasses the sonic speed and decelerates after the first shock wave. Over the airfoil, the flow undergoes the second stage of acceleration and forms the second shock structure near the trailing edge. Strong inverse pressure gradients due to 
$\unicode[STIX]{x1D706}$-type shock waves are observed aft of the first acceleration.
Figure 3. Time-averaged flow field around the NACA 0012 airfoil. (a) Streamwise velocity (
$v_{x}/a_{\infty }$). (b) Numerical schlieren with density gradient magnitude (
$\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\Vert L_{c}/\unicode[STIX]{x1D70C}_{\infty }$). Regions enclosed by the solid and broken curves represent the supersonic and recirculation regions, respectively. The dashed-dotted box depicts the spatial window for the filtered resolvent analysis.
3.2 Resolvent analysis
As transonic buffet is characterized by the energetic shock-wave oscillation over the airfoil, the cause of buffet may be uncovered by closely studying the input–output process at the relevant fluctuation frequency using the resolvent analysis described in § 2.3. Let us first present in figure 4(a) the distribution of the leading gain 
$\unicode[STIX]{x1D70E}_{1}$ over a range of frequencies from resolvent analysis without windowing. The leading gain 
$\unicode[STIX]{x1D70E}_{1}$ shows a distinct peak at 
$St=1$ corresponding to the von Kármán shedding frequency. However, there are no other distinct peaks appearing in the gain distribution for 
$\unicode[STIX]{x1D70E}_{1}$. This is expected since the dominant physics for this low-Reynolds-number base flow is the wake shedding. We do not observe noticeable peaks even for subdominant modes.
Figure 4. (a) Original and (b) windowed resolvent gains over 
$St$.
Let us present the dominant resolvent modes for 
$St=1$ without windowing. Since the gain peaks at the shedding frequency, response modes exhibit modal shapes that are representative of convective dynamics in the wake, as visualized in figure 5. The forcing modes are present in the vicinity of the airfoil where shear is large, as can be observed from figure 3. This is in line with our understanding of the formation of wake vortices from the roll up of shear layers. These observations are also in agreement with modal analysis for subsonic flows over airfoils (Zhang & Samtaney Reference Zhang and Samtaney2016; Yeh & Taira Reference Yeh and Taira2019).
Figure 5. Response and forcing modes normalized by the velocity magnitude for 
$St=1$.
Next, let us consider the gain distribution from the windowed resolvent analysis. Here, the output is windowed by operator 
$\unicode[STIX]{x1D63E}$ with a spatial profile that is restricted to the region in the vicinity of the shock, as shown in figure 3. With the output of resolvent analysis focusing on the fluctuations close to the standing shock, we find that the dominant gain 
$\unicode[STIX]{x1D70E}_{1}$ shows a peak at the buffet frequency of 
$St=0.06$, which is in agreement with past numerical and experimental findings (Deck Reference Deck2005; Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009; Sartor et al. Reference Sartor, Mettot and Sipp2014; Dandois Reference Dandois2016), which were conducted at high Reynolds numbers of 
$Re_{L_{c}}>O(10^{6})$. For subdominant gains, we observe peaks appearing at superharmonics that capture smaller-scale oscillations. We can observe a minor peak in 
$\unicode[STIX]{x1D70E}_{1}$ at 
$St=1$. However, the overall amplification process for the shock region is focused towards the low-frequency components at the buffet frequency. This suggests that the amplification physics for buffet is buried underneath the global input–output dynamics at this low Reynolds number. We also note that the conclusive assessment on the buffet physics is robust against the choice of the response window. A prominent peak at the buffet frequency 
$St=0.06$ with similar forcing structure is still revealed using a different window that covers the entire suction surface over 
$x_{c}\in [0,1]$ and 
$y_{c}\in [0,\infty )$.
We visualize the response and forcing modes that correspond to the buffet frequency of 
$St=0.06$ as identified from the gain distribution with windowed resolvent analysis. These modes in figure 6 are quite different from the modes at 
$St=1$. While the response modes for 
$St=1$ exhibit convective oscillations in the wake region, the response modes for 
$St=0.06$ are predominantly supported in the region of the standing shock. The windowed analysis is able to highlight the response in the region of the shock. For these reasons, the response modes can be seen as the oscillatory mechanism for buffet.
Figure 6. The original (top) and the windowed (bottom) resolvent modes normalized by the velocity magnitude for 
$St=0.06$. The dashed-dotted lines depicts the position of the shock on the suction side.
The input to generate the buffet unsteadiness can be identified by studying the forcing modes. We can observe from figure 6 that the forcing modes from both resolvent analysis with and without windowing show the boundary layer around the airfoil to be the source of buffet. In particular, from the windowed forcing mode, we find that the region corresponding to the shock foot is the most sensitive region to instigate the buffet phenomenon – which is in agreement with Sartor et al. (Reference Sartor, Mettot and Sipp2014), who have carried out resolvent analysis at 
$Re=3\times 10^{6}$. In a number of past studies (Deck Reference Deck2005; Grossi et al. Reference Grossi, Braza and Hoarau2014; Fukushima & Kawai Reference Fukushima and Kawai2018), the SWBLI at the shock foot exhibited turbulent boundary layer separation, which was suspected to play a role in the emergence of buffet. The current result by the windowed resolvent analysis suggests that the low-frequency buffet can be stimulated even though the boundary layer is entirely laminar at a significantly lower Reynolds number of 
$Re=1000$. Interestingly, the response modes do not show fluctuations about the weaker pressure-side shock wave observed in the numerical schlieren of figure 3(b).
Let us also bring attention to the vicinity of the trailing edge for the original forcing modes in figure 6. We can observe that there is a compact region at the trailing edge that appears in the original forcing modes. Since the dominant response mode covers the wing surface, its fluctuation alters the circulation of the airfoil, making the trailing edge act as sensitive location due to its singular nature from the cusp geometry. For this reason, we see the trailing edge also supporting the fluctuations in the original forcing mode. This observation is in tune with the report of Nitzsche (Reference Nitzsche2009) and past control techniques using a trailing-edge deflector to effectively attenuate buffet (Gao et al. Reference Gao, Zhang, Kou, Liu and Ye2017). The blockage of flow around the trailing edge by the deflector can obstruct the change in circulation in a kinematic manner. The current resolvent analysis with windowing is able to identify these sources of buffet in a clear manner. While the present discussion is concerned only with 
$\unicode[STIX]{x1D6FC}=3^{\circ }$, we note that similar observations are made for other angles of attack of 
$\unicode[STIX]{x1D6FC}=1^{\circ }$ and 
$5^{\circ }$.
Figure 7. Transonic flows over the airfoil with perturbations added in the form of the dominant forcing modes from original (a,c) and windowed resolvent analysis (b,d). The motions of the shock wave are visualized in (a,b) with the blue contour lines in four phases 
$\unicode[STIX]{x1D719}_{f}$ with respect to the forced buffet cycle, on top of the red lines for the unforced case. Both contour lines display levels of 
$\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\Vert L_{c}/\unicode[STIX]{x1D70C}_{\infty }$. Visualized in (c,d) are differences of the r.m.s. of streamwise velocity fields 
$\unicode[STIX]{x0394}v_{rms,x}/a_{\infty }$, where 
$\unicode[STIX]{x0394}v_{rms,x}\equiv v_{rms,x,forced}-v_{rms,x,unforced}$.
To corroborate our findings from resolvent analysis, we perform companion simulations of transonic flows over the airfoil with perturbations added in the shape of forcing modes. We add the forcing to the flow as a body force, 
$\boldsymbol{f}_{body}$, with the spatial profiles given by the dominant forcing modes at 
$St=0.06$ (see figure 6). The amplitude of forcing is chosen such that the momentum coefficient 
$C_{\unicode[STIX]{x1D707}}\equiv \Vert \,\boldsymbol{f}_{body}\Vert /({\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{\infty }U_{\infty }^{2}L_{c})=0.003$, which is low for this viscous flow problem (Munday & Taira Reference Munday and Taira2018). These perturbations with the spatial profiles of the primary forcing modes without and with windowing do indeed stimulate the emergence of buffet over the airfoil, as shown in figure 7(a,b). Even at this low Reynolds number, the standing shock is found to oscillate violently with large-amplitude motion over the airfoil. The amplification in the oscillation from the forcing input to the flow response (including nonlinear effects) is assessed using 
$\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D705}\unicode[STIX]{x0394}\Vert \boldsymbol{u}_{rms}\Vert /\Vert \,\boldsymbol{f}_{body}\Vert$, where 
$\unicode[STIX]{x0394}\Vert \,\boldsymbol{u}_{rms}\Vert \equiv \Vert \boldsymbol{u}_{rms}\Vert _{forced}-\Vert \boldsymbol{u}_{rms}\Vert _{unforced}$ and 
$\unicode[STIX]{x1D705}\equiv \unicode[STIX]{x1D70C}_{\infty }\,f^{+}$, showing higher amplification of 
$\unicode[STIX]{x1D6E4}=9.4$ using the forcing mode from the original resolvent analysis, compared to 
$\unicode[STIX]{x1D6E4}=7.0$ using that from the windowed analysis. The higher amplification agrees with the higher gain at 
$St=0.06$ obtained from the original resolvent analysis (see figure 4), which measures the global response to the forcing input. Although the use of the windowed resolvent forcing profile shows lower amplification, the shock oscillation coupled with the shear layer fluctuation at the shock foot is more energetic, as observed in figure 7(d). These results validate the insights gained from resolvent analysis. Perturbations at the foot of the standing shock becomes amplified through the linear dynamics about the base flow to oscillate the shock in a violent manner, even at a Reynolds number that is commonly not associated with buffet.
4 Conclusions
We examined the origin of two-dimensional transonic buffet over a NACA 0012 airfoil at a Reynolds number of 2000 through resolvent analysis. While the transonic base state at this Reynolds number exhibits unsteadiness only due to the von Kármán shedding, we show that the amplification mechanism for buffet is present in the global dynamics. The response and forcing modes from the windowed resolvent analysis revealed that the source of transonic buffet lies at the shock foot. Perturbations within the boundary layer at the shock foot can be amplified to produce oscillations about the standing shock over the airfoil with a low frequency of 
$St=0.06$. We also noted that such perturbations are closely tied to the change in the flow around the trailing edge, suggesting the effectiveness of trailing-edge buffet control devices. The identified buffet mechanism was validated by companion DNS that generated buffet by perturbing the flow through the forcing modes. The results show that even at a Reynolds number much lower than what is traditionally associated with transonic buffet, we are able to instigate buffet through the hidden amplification mechanisms. This mechanism does not require the flow to be at high Reynolds number. The findings from the present study offers fundamental insights not only into the origin of buffet but also for low-Reynolds-number compressible aerodynamics with the growing interest in the development of Martian aircraft.
Acknowledgements
Y.K. and M.K. gratefully acknowledge the Overseas Travel Assistance Program supported by the TUAT President’s Office. C.-A.Y. and K.T. gratefully acknowledge support from the US Office of Naval Research (Program Manager: Dr B. Holm-Hansen, grant no. N00014-19-1-2460) and the US Air Force Office of Scientific Research (Program Managers: Drs G. Abate and D. Smith, grant no. FA9550-18-1-0040).
Declaration of interests
The authors report no conflict of interest.
