2.1. Biglobal resolvent analysis

We consider the spatially discretized nonlinear governing equation,

(2.1) \begin{equation} \frac {{\rm d} \boldsymbol {q}}{{\rm d} t}=\mathcal{N}(\boldsymbol {q}), \end{equation}

where $\boldsymbol {q}=(\rho ,\rho \boldsymbol {u},\rho e)\in \mathbb{R}^{n}$ represents the state vector, $\rho$ is the density, $\boldsymbol {u}=(u_x,u_y,u_z)$ is the velocity vector with components along the streamwise ( $x$ ), cross-stream ( $y$ ) and spanwise ( $z$ ) directions, and $e$ is the energy. Here, $\mathcal{N}\in \mathbb{R}^{n\times n}$ is the nonlinear Navier–Stokes evolution operator, and $n=N\times 5$ , where $N$ is the number of cells in the spatial discretization and the number of state variables is five. We now consider the flow field $\boldsymbol {q}=\boldsymbol {q}_b+\boldsymbol {q}^\prime$ to be composed of the sum of a stationary base flow $\boldsymbol {q}_b$ and a fluctuating component $\boldsymbol {q}^\prime$ of small amplitude. The base flow is considered to be the time- and spanwise-averaged flow. By substituting the decomposition in (2.1), and performing a Taylor expansion, we obtain the spatially discretized linearized governing equation for the fluctuating perturbation component, $\boldsymbol {q}^\prime$ :

(2.2) \begin{equation} \frac {{\rm d} \boldsymbol {q}^\prime }{{\rm d} t}=\mathcal{L}\boldsymbol {q}^\prime +\boldsymbol {f}^\prime . \end{equation}

Here, $\mathcal{L}\equiv \boldsymbol{\nabla }_{\boldsymbol {q}} \mathcal{N} |_{\boldsymbol {q}_b} \in \mathbb{R}^{n\times n}$ is the linearized Navier–Stokes operator about the base flow, $\boldsymbol {q}^\prime \in \mathbb{R}^{n}$ is the perturbation and $\boldsymbol {f}^\prime \in \mathbb{R}^{n}$ collects the nonlinear terms (McKeon & Sharma Reference McKeon and Sharma2010; Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024). Due to the temporal and spanwise homogeneity of the base flow, the response ( $\boldsymbol {q}^\prime$ ) and forcing ( $\boldsymbol {f}^\prime$ ) can be decomposed through a spatiotemporal Fourier transform as follows:

(2.3) \begin{align} \boldsymbol {q}^\prime (x,y,z,t)=\int ^\infty _{-\infty }\int ^\infty _{-\infty }\hat {\boldsymbol {q}}_{\omega ,\beta }(x,y)e^{-i\omega t}e^{i\beta z} {\rm d}\omega {\rm d}\beta ,\nonumber\\ \boldsymbol {f}^\prime (x,y,z,t)=\int ^\infty _{-\infty }\int ^\infty _{-\infty }\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }(x,y)e^{-i\omega t}e^{i\beta z }{\rm d}\omega {\rm d}\beta . \end{align}

Here, $\beta$ and $\omega$ indicate the spanwise wavenumber and temporal frequency, respectively. By substituting these expressions into (2.2), we find the following input–output relationship:

(2.4) \begin{equation} \hat {\boldsymbol {q}}_{\omega ,\beta }(\boldsymbol {x})= (-i\omega \boldsymbol{I}-\mathcal{L}_\beta )^{-1}\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }(\boldsymbol {x})=\boldsymbol{H}_{\! \omega ,\beta }\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }(\boldsymbol {x}), \end{equation}

where $\boldsymbol{H}_{\! \omega ,\beta }\in \mathbb{C}^{n\times n}$ is the resolvent operator that acts as the transfer function between the forcing $\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }$ and the response $\hat {\boldsymbol {q}}_{\omega ,\beta }$ at frequency $\omega$ and spanwise wavenumber $\beta$ .

Among all the possible forcing and response pairs that satisfy (2.4), resolvent analysis finds the optimal pair, i.e. those that maximize the energy at given values of $\omega$ and $\beta$ ,

(2.5) \begin{equation} \sigma ^2=\max _{\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }}\frac {\big \langle \hat {\boldsymbol {q}}_{\omega ,\beta },\hat {\boldsymbol {q}}_{\omega ,\beta }\big \rangle _{E}}{\big \langle \kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta },\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }\big \rangle _{E}}=\max _{\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }}\frac {\big|\big|\hat {\boldsymbol {q}}_{\omega ,\beta }\big|\big|^2_{E}}{\big|\big|\kern2pt \hat{\boldsymbol{\kern-2pt f}}_{\! \omega ,\beta }\big|\big|^2_{E}}, \end{equation}

where $||\boldsymbol{\cdot }||_{E}$ is a suitable energy norm. In this work we consider the Chu norm (Chu Reference Chu1965; George & Sujith Reference George and Sujith2011), which is expressed as

(2.6) \begin{equation} E_{{Chu}}=\frac {1}{2}\int _V\left (\bar {\rho }|\boldsymbol {u}|^2+\frac {a^2\rho ^2}{\gamma \bar {\rho }}+\frac {\bar {\rho } c_v T^2}{\bar {T}}\right ){\rm d}V, \end{equation}

where the variables with $(\,\bar {\boldsymbol{\cdot }}\,)$ represent the base flow quantities. To take into account the energy norm, a similarity transform is applied to the resolvent operator. Therefore, we consider instead $\boldsymbol{H}_{\! \omega ,\beta }^W= W \kern-1pt \boldsymbol{H}_{\! \omega ,\beta } W^{-1}$ . Here, $W$ is a volume-weighted matrix that allows us to express the energy norm in terms of the $L_2$ norm for the singular value decomposition.

A singular value decomposition can now be performed on the weighted matrix. By retaining only the first $m\ll n$ singular values and right/left singular vectors, we find a low-rank approximation of $\boldsymbol{H}^W_{\omega ,\beta }$ :

(2.7) \begin{equation} \boldsymbol{H}^W_{\omega ,\beta }\approx \boldsymbol{U}\varSigma \boldsymbol{V^*}. \end{equation}

The columns of $\boldsymbol{U}=[\boldsymbol {u}_1, \boldsymbol {u}_2, \ldots , \boldsymbol {u}_m]\in \mathbb{C}^{n\times m}$ and $\boldsymbol{V}=[\boldsymbol{v}_1, \boldsymbol{v}_2, \ldots , \boldsymbol{v}_m]\in \mathbb{C}^{n\times m}$ hold the response and forcing modes, respectively, while $\varSigma =\text{diag}(\sigma _1,\sigma _2,\ldots ,\sigma _m)\in \mathbb{R}^{m}$ retains the gains of the corresponding forcing–response pairs.

In the present work, some cases are characterized by a linear dynamics that present eigenvalues with positive growth rate. For this reason, we use discounted resolvent analysis (Jovanovic Reference Jovanovic2004; Jovanović & Bamieh Reference Jovanović and Bamieh2005), that considers a Laplace transform instead of a Fourier transform. This modification is equivalent to temporally damping the forcing and response by $e^{-\gamma t}$ , which translates to considering the dynamic over finite time scales. A zero temporal damping $\gamma =0$ , when there are no unstable poles, corresponds to investigating the asymptotic (infinite-time) dynamics. Introducing $\gamma \neq 0$ , we consider the dynamics over a finite-time horizon $t_\gamma =1/\gamma$ . To apply discounting, the integration line of the inverse Laplace transform is taken above all positive real parts of the eigenvalues of the linear operator. The discounted parameter $\gamma$ is thus introduced, which must satisfy $\gamma \gt \text{Re}\{\lambda\}$ , with $\lambda$ being the eigenvalue of the linear operator with the largest real part. With discounting, the resolvent operator now reads

(2.8) \begin{equation} \boldsymbol{H}_{\! \omega ,\beta }=[( \gamma -i\omega )\boldsymbol{I}-\mathcal{L}_\beta ]^{-1}. \end{equation}

The singular value decomposition of the resolvent operator is approximated using the Krylov subspace projection method, with a subspace dimension for the reduced-order problem set to $m=24$ . Both the simulation of the base flow and resolvent analysis are performed within the compressible flow solver CharLES (Khalighi et al. Reference Khalighi, Ham, Nichols, Lele and Moin2011), coupled with the PETSc and SLEPc libraries (Roman et al. Reference Roman, Campos, Romero and Tomás2016; Balay et al. Reference Balay2020) for the computation of the singular value decomposition. Further details on the formulation and implementation of resolvent analysis used herein can be found in Rolandi et al. (Reference Rolandi, Ribeiro, Yeh and Taira2024).

Figure 2. Computational set-up used for the base flow computation and resolvent analysis.

2.2. Computational set-up

We simulate the spanwise periodic flow around a two-dimensional NACA0012 airfoil at angles of attack $14^\circ$ and $9^\circ$ for a Reynolds number of ${\textit{Re}}=U_\infty c/\nu =1000,\;2500,\;5000$ and $10\,000$ and a Mach number of $M_\infty =0.1$ . Here, $U_\infty$ indicates the free stream velocity, $c$ the chord of the airfoil and $\nu$ the kinematic viscosity. At these Reynolds numbers, the flow is three-dimensional, so a spanwise domain length $L_z$ of one chord is considered for the computational domain. This choice is based on observations indicating that the flow transitions to a three-dimensional state with a spanwise wavelength of approximately $\lambda _z\approx c/3$ (Gupta et al. Reference Gupta, Zhao, Sharma, Agrawal, Hourigan and Thompson2023). The base flows are simulated with direct numerical simulations (DNS) for ${\textit{Re}}=1000$ and $2500$ and wall-resolved large eddy simulations (LES) for ${\textit{Re}}=5000$ and $10\,000$ . For the wall-resolved large eddy simulation, we use the Vreman subgrid-scale model (Vreman Reference Vreman2004). The boundary conditions are comprised of a no-slip adiabatic condition on the airfoil surface, uniform constant velocity, pressure, and temperature at the inlet and on far-field boundaries, a zero-pressure gradient at the outlet and periodic boundary conditions on the lateral sides. For each case, the time step is selected such that the local Courant–Friedrichs–Lewy number is less than unity.

The base flow is computed by time- and spanwise-averaging the instantaneous flow field over ${\sim} 100$ convective times ( $t_c=tU_\infty /c$ ) after the transient is flushed out. Convergence of the base flow is checked by ensuring that the relative variation on the velocity flow field over the last five convective times satisfies

(2.9) \begin{equation} \frac {||\bar {\boldsymbol {u}}(t_{\textit{end}})-\bar {\boldsymbol {u}}(t_{\textit{end}}-5t_c)||_2}{||\bar {\boldsymbol {u}}(t_{\textit{end}})||_2}\lt 0.01, \end{equation}

where $\bar {\boldsymbol {u}}$ is the time- and spanwise-averaged velocity flow field.

Details of the mesh used in this study are shown in figure 2. The red-framed domain indicates the computational domain used for computing the base flow and is extended to $x/c\in [-20,25]$ , $y/c\in [-20,20]$ and $z/c\in [0, 1]$ . The origin is positioned at the leading edge of the airfoil. For the biglobal resolvent analysis, the time- and spanwise-averaged base flow is interpolated onto a smaller two-dimensional grid, the yellow-framed grid, whose extent is $(x/c,y/c)\in [-14,14]\times [-4.5,4.5]$ . This reduction is possible because the domain and grid resolution requirements for computing resolvent modes differ significantly from those used in unsteady simulations (Rolandi et al. Reference Rolandi, Ribeiro, Yeh and Taira2024). Specifically, the modal structures are concentrated near the airfoil, eliminating the necessity for an extended domain, which is instead essential for the base flow simulations. Additionally, the grid refinement near the airfoil focuses on both the downstream and upstream regions, as the forcing modes develop upstream due to the convective nature of the amplification mechanisms.

We performed a grid convergence study using three meshes for the unsteady simulation to compute the base flow. In addition to these meshes, two meshes were tested for the resolvent analysis. In table 1, we report the results of the grid convergence study performed on the highest Reynolds number, ${\textit{Re}}=10\,000$ . The combinations of the meshes used for the base flow computation and resolvent analysis are summarized in Case 1, Case 2 and Case 3. Cells $_{\textit{BF}}$ and Cells $_{\textit{R}}$ refer to the number of cells of the base flow mesh and resolvent analysis mesh, respectively. Number of cells on the $xy$ -plane ( $N_{\textit{xy}}$ ) and along the $z$ -direction ( $N_{z}$ ) are reported. The computational domain of Case 2 and Case 3 for the unsteady simulation of this Reynolds number extends over $x/c\in [-20,25]$ , $y/c\in [-20,20]$ and $z/c\in [0, 0.4]$ , given that the spanwise structures are smaller compared with the lower Reynolds number cases. The values of mean drag and lift coefficients from the unsteady simulation are shown, together with the frequency $St^*$ of maximum gain from the resolvent analysis. Resolvent analysis for the grid study was performed considering $\beta =0$ and $\gamma =1.25$ . Overall, we can see convergence in the mean drag and lift coefficients and a relative difference between the results of $St^*$ within $5\,\%$ . In figure 3, we present the $0$ -contour of the streamwise velocity for the time- and spanwise-averaged base flow for the three different cases. Additionally, we display the streamwise velocity contours of the response mode at $\omega /2\pi =2.7$ . The results of both the base flow and resolvent analysis with the different meshes are in good agreement. The meshes of Case 2 for both the unsteady simulation and resolvent analysis were used in this work.

Table 1. Mesh convergence check for the ${\textit{Re}}=10\,000$ case. Values of time-averaged drag ( $\bar {C}_D$ ) and lift ( $\bar {C}_L$ ) coefficients from the unsteady simulation together with frequency of maximum amplification from resolvent analysis at $\beta =0$ and $\gamma =1.25$ for the different meshes tested.

Figure 3. Mesh convergence check for the ${\textit{Re}}=10\,000$ case. (a) Contour of time- and spanwise-averaged streamwise velocity $\bar {u}_x=0$ . (b) Contours of the streamwise velocity of the first response mode at frequency $\omega /2\pi =2.7$ , $\beta =0$ and $\gamma =1.25$ .

In table 2, we provide details on the computational set-up in terms of number of grid points on the $xy$ -plane $N_{\textit{xy}}$ , the number of grid points along the suction side of the airfoil section $N_{\textit{airfoil}}$ (symmetric with respect to the pressure side), the number of grid points in the spanwise direction $N_z$ and its spanwise extent $L_z$ , initial cell spacing on the wake $\Delta x$ and vertical off-wall spacing $\Delta y$ . The discretization details are reported for the unsteady simulation at the different Reynolds numbers, and for the mesh used for the resolvent analysis.

Table 2. Mesh details for the unsteady simulation of the different Reynolds numbers and resolvent analysis.

Figure 4. Instantaneous flow fields around a NACA0012 airfoil at $\alpha =14^\circ$ and different Reynolds numbers. Visualization of isosurfaces of Q-criterion $Q=0.05$ , coloured by streamwise velocity, and $Q=0.005$ in translucent.

Figure 5. Lift coefficient $C_L$ , drag coefficient $C_D$ and lift spectra $\hat {C}_L$ at the different Reynolds numbers.

3.1. Unsteady flows

The flow fields around a NACA0012 airfoil at an angle of attack of $14^\circ$ for the considered Reynolds numbers are shown in figure 4(a–d). At these Reynolds numbers, the flow is unsteady and exhibits the characteristic von Kármán vortex shedding in the wake region. The flow around the airfoil at $\alpha =14^\circ$ undergoes a transition from steady state to periodic at approximately ${\textit{Re}}\approx 380$ through a Hopf bifurcation (Rolandi et al. Reference Rolandi, Jardin, Fontane, Gressier and Joly2022). As the Reynolds number increases further, the periodic flow transitions to three-dimensional dynamics through a period-doubling bifurcation, also known as the Mode C instability (Sheard et al. Reference Sheard, Thompson and Hourigan2005a ; Meneghini et al. Reference Meneghini, Carmo, Tsiloufas, Gioria and Aranha2011; Rolandi Reference Rolandi2021). The Reynolds number at which the flow becomes three-dimensional at this angle of attack is between ${\textit{Re}}=750$ and $1000$ (Gupta et al. Reference Gupta, Zhao, Sharma, Agrawal, Hourigan and Thompson2023). This transition is characterized by the emergence of a subharmonic component of the vortex shedding, effectively doubling the flow periodicity. Notably, period-doubling bifurcations are also associated with the onset of chaos in fluid flows (Pulliam & Vastano Reference Pulliam and Vastano1993), through the so-called period-doubling cascade.

The three-dimensional flow resulting from this transition features spanwise structures that develop across the airfoil. The Mode C instability develops in the stretched region between two consecutive vortices, called the braid region, and results in the formation of elongated streamwise vortices. At ${\textit{Re}}=1000$ , the spanwise wavelength of the three-dimensional structures is approximately $\lambda _z\approx c/3$ , consistent with previous studies (Gupta et al. Reference Gupta, Zhao, Sharma, Agrawal, Hourigan and Thompson2023), where the stability of the periodic solution was examined through Floquet analysis. As the Reynolds number increases further, the shear layer separating from the leading edge becomes unstable. The shear layer rolls up closer to the leading edge, forming two-dimensional spanwise vortical structures, that are related to the Kelvin–Helmholtz instability. The characteristic length of the vortical spanwise elongated structures decreases with the Reynolds numbers (see figure 4 c–d). Finally, these structures are convected downstream, where they break down into smaller, three-dimensional structures, contributing to the increasing level of turbulence in the flow.

We report the time evolution of the drag and lift coefficients in figure 5, together with the frequency spectra of the lift coefficients in terms of Strouhal number $St=c f/U_\infty$ , where $f$ indicates the frequency. For the case of ${\textit{Re}}=1000$ , we observe a clear peak at $St= 0.72$ and its harmonics, due to the periodicity of the flow at this Reynolds number. A clear peak and its harmonics are also visible at ${\textit{Re}}=2500$ , with the low-frequency contents becoming increasingly energetic. At this Reynolds number, the oscillations have higher amplitude compared with ${\textit{Re}}=1000$ and the peak is at a lower frequency corresponding to $St=0.67$ . Increasing the Reynolds number further, the frequency peak of the wake oscillation for both ${\textit{Re}}=5000$ and $10\,000$ occurs at ${\textit{St}}\approx 0.67$ , as for the ${\textit{Re}}=2500$ case. However, the amplitudes of the spectral peak reduce with the Reynolds number, consistent with the decreasing oscillatory type behaviour. At these Reynolds numbers, we also observe a broader frequency spectrum and the emergence of slower dynamics, visible in the $C_L$ and $C_D$ variations, which are associated with the interaction between the separated shear-layer roll-up and the airfoil (Zaman, McKinzie & Rumsey Reference Zaman, McKinzie and Rumsey1989; Mukai, Enomoto & Aoyama Reference Mukai, Enomoto and Aoyama2006). For ${\textit{Re}} = 5000$ , these slower dynamics occur at $St \approx 0.05$ and $St \approx 0.28$ , while for ${\textit{Re}} = 10\,000$ , they appear at $St \approx 0.04$ and $St \approx 0.24$ .

Figure 6. Instantaneous lift coefficient at ${\textit{Re}}=5000$ : (a) low-pass-filtered over $St\lt 0.1$ ; (b) bandpass-filtered over $0.2\lt St\lt 0.4$ . Instantaneous spanwise vorticity fields are shown at instants indicated by the red dots.

To further investigate the slower dynamics observed in the lift spectra at ${\textit{Re}}=5000$ and $10\,000$ , we focus on the ${\textit{Re}}=5000$ case. Figure 6 presents the frequency-filtered time evolution of the lift coefficient for the ${\textit{Re}}=5000$ case. Specifically, we examine the low-frequency range $0\lt St\lt 0.1$ and the midfrequency range $0.2\lt St\lt 0.4$ , targeting the spectral peaks at ${\textit{St}}\approx 0.05$ and ${\textit{St}}\approx 0.28$ , respectively, as reported in the lift spectra in figure 5 for the ${\textit{Re}}=5000$ case. Instantaneous spanwise vorticity fields corresponding to the crests and troughs of the filtered signals are also shown. The vorticity fields at the troughs of the low-frequency signal, at $t=t_L^1,\;t_L^3$ and $t_L^5$ in figure 6(a), have similar flow characteristics. The separated shear layer is elongated and rolls up farther downstream from the leading edge, with vortex cores forming at streamwise locations $x\gt 0.7$ . This positioning keeps the vorticity detached from the airfoil’s suction surface, thereby limiting the development of induced vorticity of opposite-sign near the surface. As a result, the flow is massively separated from the airfoil.

In contrast, at the crests of the low-frequency signal, at $t=t_L^2,\;t_L^4$ and $t_L^6$ , the shear layer rolls up much closer to the suction surface $x\lt 0.7$ , leading to the generation of strong opposite-sign vorticity. This causes the flow to locally form a secondary recirculation region with higher and positive streamwise velocity on the airfoil surface, generating local low-pressure zones, which enhance the lift force. This suggest that the low-frequency behaviour is associated with the formation of the secondary recirculation region (Zaman et al. Reference Zaman, McKinzie and Rumsey1989; Broeren & Bragg Reference Broeren and Bragg1998; Mukai et al. Reference Mukai, Enomoto and Aoyama2006). Therefore, the flow alternates between two states: (i) a massive separated flow with a negligible secondary recirculation region ( $t=t_L^1,\;t_L^3$ and $t_L^5$ ) and (ii) a separated flow coexisting with a secondary recirculation region ( $t=t_L^2,\;t_L^4$ and $t_L^6$ ). The two states depend on the roll-up location of the separated shear layer. We also note that the low frequency dynamics ( ${\textit{St}}\approx 0.05$ ) are an order of magnitude lower than the vortex shedding ( ${\textit{St}}\approx 0.67$ ), which is in agreement with previous studies (Zaman et al. Reference Zaman, McKinzie and Rumsey1989; Broeren & Bragg Reference Broeren and Bragg1998), which associate this mechanism with the breathing of the laminar separation bubble.

We now consider the midfrequency signal, bandpass-filtered over $0.2\lt St\lt 0.4$ . In figure 6(b), we observe that at the crests of the signal, at $t=t_M^1,\;t_M^3$ and $t_M^5$ , the shear-layer vortex has not fully rolled-up. On the contrary, at the through of the signal, $t_M^2,\;t_M^4$ and $t_M^6$ , the leading-edge vortex is rolled-up, forming a defined single vortical structure detached from the downstream same sign negative vorticity. Some works link these dynamics to vortex pairing (Tang et al. Reference Tang, Wang, Wang, Hong, Yao and Tang2021).

Figure 7. Contours of energy spectra at ${\textit{Re}}=1000$ , $2500$ , $5000$ and $10\,000$ . Contours are shown at different streamwise locations $x/c$ along $y/c\in [0,0.25]$ . Black dashed lines indicate the dominant frequency peaks associated with the lift coefficient, see figure 5. Red horizontal lines indicate the vertical locations considered in figure 8.

The lift coefficient signal effectively captures the dynamics of large vortex shedding in the wake, and the slow dynamics of the separation region. However, our interest also lies in the dynamics of the separated shear layer, which manifests as the leading-edge vortex roll-up. To investigate the shear-layer dynamics, we compute the energy spectra from the shear layer at various streamwise positions: $x/c = 0.25, 0.5, 0.75$ and $1$ , along $0 \lt y/c \lt 0.25$ , as shown in figure 7.

At ${\textit{Re}} = 1000$ , the flow is periodic, with the energy spectra at each position showing a dominant peak corresponding to vortex shedding. At $x/c = 1$ , the subharmonic becomes significant, indicating a period-doubling instability of the NACA0012 periodic shedding that leads the transition to three-dimensional flow (Sheard et al. Reference Sheard, Thompson, Hourigan and Leweke2005b ; Meneghini et al. Reference Meneghini, Carmo, Tsiloufas, Gioria and Aranha2011; Rolandi Reference Rolandi2021; Gupta et al. Reference Gupta, Zhao, Sharma, Agrawal, Hourigan and Thompson2023). At ${\textit{Re}} = 2500$ , the harmonic components become more prominent, and at even higher Reynolds numbers, the spectrum broadens significantly, especially as the streamwise position approaches the wake. Also for the higher Reynolds number cases, we observe the presence of subharmonics, which, in this case, correspond to vortex pairing within the shear layer. For ${\textit{Re}} \geqslant 2500$ , the highest energy amplitudes at $x/c = 0.25$ and $0.5$ are concentrated at a specific cross-stream location, approximately $y/c \approx 0.1$ . This identifies a broad region of high-frequency, high-amplitude unsteadiness associated with the shear-layer dynamics at the subsequent streamwise positions.

For the specific locations indicated by the red dashed lines in figure 7, we also present the spectra of the vertical velocity in figure 8. As the Reynolds number increases, we observe a gradual intensification of the harmonic peaks, reflecting enhanced nonlinear interactions (Deissler Reference Deissler1969). This is particularly evident at ${\textit{Re}} = 2500$ , where higher-order harmonics becomes clearly distinguishable, and at ${\textit{Re}} = 5000$ where the first harmonic is highly energetic. At the highest Reynolds number, ${\textit{Re}} = 10\,000$ , the spectrum becomes noticeably broader. The dominant energy shifts towards higher frequencies, with the maximum peak occurring at $St \approx 2.6$ , which indicates the shear-layer roll-up, as will be discussed later.

Figure 8. Vertical velocity spectra $\hat {u}_y$ at locations indicated by the the red dashed lines in figure 7 at the different Reynolds numbers.

3.2. Oscillator and amplifier dynamics

The spectral change observed in figures 7 and 8 reflects the transition between oscillator- to amplifier-type behaviour (Towne et al. Reference Towne, Colonius, Jordan, Cavalieri and Bres2015; Rosenberg, Symon & McKeon Reference Rosenberg, Symon and McKeon2019; Symon et al. Reference Symon, Sipp and McKeon2019). Oscillator flows exhibit intrinsic global instabilities that give rise to self-sustained unsteadiness. In contrast, amplifiers are convectively unstable and amplify external disturbances.

At low Reynolds number, ${\textit{Re}} \lt 2500$ , the flow exhibits narrow-band spectral peaks with clear harmonics, characteristic of an oscillator-type dynamics. As the Reynolds number increases, ${\textit{Re}} \gt 2500$ , the frequency spectrum broadens significantly and energy shifts towards higher frequencies. This broadening reflects the transition to an amplifier-type regime dominated by shear-layer instabilities.

3.3. Time-averaged base flows

The unsteady flow field is averaged in time and along the spanwise direction to obtain the base flow for the resolvent analysis. The resulting streamwise velocity fields are shown in figure 9 for the considered Reynolds numbers. As the Reynolds number increases, the shear layer separating from the leading edge becomes progressively thinner and curves due to stronger reverse flow in the recirculation region.

Figure 9. Time- and spanwise-averaged (base flow) streamwise velocity around a NACA0012 wing at $\alpha =14^\circ$ and (a) ${\textit{Re}}=1000$ , (b) $2500$ , (c) $5000$ and (d) $10\,000$ . (e) Contour of time- and spanwise-averaged streamwise velocity $\bar {u}_x=0$ .

The corresponding zero streamwise velocity contours of the time- and span-averaged base flows are visualized in figure 9(e). At ${\textit{Re}}=1000$ , the base flow features a single recirculation region, which is elongated compared with those at higher Reynolds numbers. When the Reynolds number is increased to ${\textit{Re}}=2500$ , the recirculation region shortens before being stretched again at higher Reynolds numbers. For ${\textit{Re}}\geqslant 2500$ , the base flow exhibits a secondary recirculation region on the airfoil’s suction side, which shifts upstream and becomes thinner as the Reynolds number increases.

The secondary recirculation region is caused by the induced vorticity generated in the separated region on the suction side, as previously discussed. Additionally, the separation point of the primary recirculation region shifts upstream with increasing Reynolds number (Counsil & Goni Boulama Reference Counsil and Goni Boulama2013; Brunner et al. Reference Brunner, Kiefer, Hansen and Hultmark2021). This behaviour is due to a higher adverse pressure gradient near the leading edge of the airfoil due to the thinning of the boundary layer over the airfoil upstream of the separation point with increasing Reynolds number.

4.1. Eigenvalue decomposition of the linear operator

Let us examine the eigenvalues of the linearized Navier–Stokes operator at the different Reynolds numbers by solving the eigenvalue problem

(4.1) \begin{equation} \mathcal{L}\boldsymbol{\phi }=-i\lambda \boldsymbol{\phi }, \end{equation}

where $-i\lambda =-i(\lambda _r+i\lambda _i)$ is the complex eigenvalue and $\boldsymbol{\phi }$ is the corresponding eigenvector. Here, we consider $\mathcal{L}=\mathcal{L}_{\beta =0}$ , while the effect of spanwise wavenumber $\beta \neq 0$ will be explored in § 4.3. This allows us to consider appropriate ranges for the discount parameter depending on whether we consider two-dimensional $\beta =0$ or three-dimensional $\beta \neq 0$ perturbations.

In figure 10(a), the eigenvalue with the largest real part for each ${\textit{Re}}$ is shown in the complex plane. Both the growth rate and frequency increase with increasing Reynolds number. The linear dynamics has an eigenvalue with a positive real part $\text{Re}\{-i\lambda \}=\lambda _i$ only at ${\textit{Re}}=10\,000$ , with $\lambda _i=0.0963$ . For the lower Reynolds number cases, all the eigenvalues have negative real parts, despite the unsteady nature of the nonlinear dynamics. In this regard, it should be considered that the eigenvalue analysis depends on the choice of the base flow and that there are conditions for the validity of mean flow stability analysis (Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016).

In fact, it is important to consider that stability analysis based on the time-averaged flow (not an exact solution of the Navier–Stokes equations) needs to be performed with care. Moreover, it has been shown that flows exhibiting quasimonochromatic oscillations can yield different linear responses when analysed about the time-averaged flow versus the underlying unstable fixed point (Barkley Reference Barkley2006; Turton, Tuckerman & Barkley Reference Turton, Tuckerman and Barkley2015; Rolandi Reference Rolandi2021). The results of $\text{Re}\{-i\lambda \}\lt 0$ even at ${\textit{Re}}=1000$ are somewhat counterintuitive, as oscillator-type flows are expected to exhibit a global instability at the vortex shedding frequency. However, for this type of flow, this expectation applies to stability analysis about the fixed point, not the time-averaged flow. In fact, previous studies have shown that the stability analysis about the fixed-point solution at ${\textit{Re}}=1000$ and $\alpha =14^\circ$ exhibits a positive growth rate (Rolandi et al. Reference Rolandi, Jardin, Fontane, Gressier and Joly2022). Nevertheless, we present here the eigenvalue analysis of the linear operator linearized around the time-averaged flow, as this constitutes a necessary step for the subsequent discounted resolvent analysis.

From the result shown in figure 10(a), we observe a monotonic increase of $\lambda _i$ with respect to the Reynolds number, with a linear increase for ${\textit{Re}}\geqslant 2500$ . We can then infer that the eigenvalue based on the time-averaged base flow crosses the imaginary axis at ${\textit{Re}}\approx 5900$ . The corresponding eigenvectors are shown in figure 10(b). At ${\textit{Re}}=1000$ the modal structure is mainly concentrated in the wake region, while at higher Reynolds numbers they also exhibit structures in the shear-layer region. As the Reynolds number increases, the position of maximum modal structure amplitude, indicated by the dashed circle, shifts from the far wake (at ${\textit{Re}}=1000$ ) to the near wake (at ${\textit{Re}}=10\,000$ ). Furthermore, for ${\textit{Re}}\geqslant 2500$ , the modal structure appears not only in the wake but also in the shear layer. This reflects the fact that for these Reynolds numbers the shear layer rolls up and interacts with the wake dynamic.

4.2. Temporal discounting

Once the eigenvalues are found, the discount parameter for the resolvent analysis is chosen such that $\gamma \gt \max \{\text{Re}\{-i\lambda \},0\}$ . This allows us to consider the overall forced dynamics over a time scale shorter than the time scale associated with $\text{Re}\{-i\lambda \}$ . Longer time scales should not be considered, because the implication of having a positive real part $\text{Re}\{-i\lambda \}\gt 0$ would make the response seemingly unbounded as $t\to \infty$ , masking the effect of forcing.

Here, a value of $\gamma \gt 0.0963$ for ${\textit{Re}}=10\,000$ corresponds to dynamics within a finite time horizon $t_\gamma \lt 10.4$ . The effect of discounting on the energy amplification is shown in figure 11. The variations of the first singular value as a function of the frequency are plotted for $\gamma =\{0.15,0.20,0.25,0.30,0.40,0.625,1.25\}$ corresponding to $0.8\leqslant t_\gamma =1/\gamma \leqslant 6.67$ . The frequency is considered in terms of Strouhal number $St=c\omega /(U_\infty 2\pi )$ . Considering this range of $t_\gamma$ , in what follows we will refer to $t_\gamma =0.8$ as the short time scale and $t_\gamma =6.67$ as the long time scale.

Figure 11. Variation of the first singular value $\sigma _1$ with respect to frequency for different finite-time horizon $t_\gamma \in [0.8;6.67]$ at ${\textit{Re}}=1000$ , $2500$ , $5000$ and $10\,000$ . Dashed grey lines indicate the frequencies of maximum gain at short and long time scales.

Figure 12. Contours of streamwise velocity component of the first response mode at the frequencies of maximum gain at short and long time scales. Here ( ) line frame indicates the mode at the lower frequency peak, ${\textit{St}}_W$ , and ( ) line frame indicates the mode at the higher frequency peak, ${\textit{St}}_S$ .

At ${\textit{Re}}=1000$ , only one peak emerges as $t_\gamma$ varies, while for ${\textit{Re}}\geqslant 2500$ two distinct peaks are observed, as indicated by the dashed lines in figure 11. The first peak appears at a high frequency on the short time scale, while on the long time scale, another peak at a lower frequency dominates. For all Reynolds numbers, the long time scale peak occurs at a frequency corresponding to the eigenvalue with the maximum real part. This is because, with lower $\gamma$ (higher $t_\gamma$ ), the inverse Laplace integration is closer to such an eigenvalue, and the norm of the resolvent operator increases at that frequency.

By comparing figure 11 with figure 8, we observe a clear correspondence between the gain peaks and the spectral content of the vertical velocity. In addition to the low-frequency peaks ${\textit{St}}_W$ associated with vortex shedding, the gain shifts towards higher frequencies as the Reynolds number increases. At ${\textit{Re}} = 2500$ , the second peak ${\textit{St}}_S^{2500}$ aligns with the first harmonic observed in figure 8, while at higher Reynolds numbers, the influence of higher-order harmonics remains evident in the gain distribution over longer time scales. Notably, at ${\textit{Re}} = 10\,000$ , the high-frequency peak in the gain ${\textit{St}}_S^{10\,000}$ matches the highest frequency observed in the vertical velocity spectra in figure 8.

In figure 12, we show the streamwise velocity components of the first response mode at short and long time scale and at the frequencies ${\textit{St}}_W$ and ${\textit{St}}_S$ , indicated in figure 11. Frequencies ${\textit{St}}_W$ and ${\textit{St}}_S$ correspond to the short and long time scale peaks, respectively. At the lower frequencies, ${\textit{St}}_W$ , we observe that the structures emerge in the wake and highlight the coupling between the leading and trailing edge, we will therefore refer to this mode as the wake mode. We also note that the modal structures at the lowest Reynolds number are similar to the structures revealed from non-modal stability analysis of low-Reynolds number flow (He et al. Reference He, Gioria, Pérez and Theofilis2017) when increasing the time horizon. On the other hand, at the higher frequency peak, ${\textit{St}}_S$ , the structures are present in the shear-layer region detaching from the leading edge. We therefore refer to this mode as the shear-layer mode. For both the wake and shear-layer modes, we observe that the structures at $t_{\gamma ,{\textit{short}}}$ are closer to the airfoil, while they develop downstream when increasing $t_\gamma$ . This corresponds to the fact that the perturbation has more time to grow, and translates into the higher energy gain shown in figure 11. Further discussion on the effects of discounting can be found in Appendix A.

Figure 13. The eigenvalues with the largest real components for different $\beta$ at ${\textit{Re}}=10\,000$ .

4.3. Spanwise wavenumber effects

In this subsection, we consider the effects of the spanwise wavenumber $\beta$ . Firstly, we need to compute the eigenvalues of the linear operator $\mathcal{L}_\beta$ at varying $\beta$ , as we performed for the $\beta =0$ case. In figure 13 the unstable eigenvalues are shown in the complex plane for $\beta \in [0,8\pi ]$ . The eigenvalue with the largest real part corresponds to $\beta =4\pi$ at $St=\lambda _r/(2\pi )=0.36$ , and the real part decreases at the same frequency for increasing wavenumber. In this case, the value of the largest real part $\lambda _i=0.358$ suggests that we should consider dynamics within a time scale of $t_\gamma \approx 2.8$ , shorter compared with the $\beta =0$ investigated in the previous subsection. The effects of spanwise wavenumber $\beta$ at short ( $t_\gamma =0.8$ ) and medium ( $t_\gamma =2.5$ ) time scales are shown in figures 14 and 15, respectively. The gain distributions for the first three singular values ( $\sigma _1,\;\sigma _2$ and $\sigma _3$ ) are shown over the $\beta -St$ plane for the different Reynolds numbers. The peaks of the spectral content of the lift coefficient, shown in figure 5, are also reported for comparison.

Figure 14. Gain distributions of the first three singular values over the $\beta -St$ plane at ${\textit{Re}}=1000$ , $2500$ , $5000$ and $10\,000$ at $t_\gamma =0.8$ . Black dashed lines indicate the dominant frequency peaks associated with lift coefficients (see figure 5).

Figure 15. Gain distributions of the first three singular values over the $\beta -St$ plane at ${\textit{Re}}=1000$ , $2500$ , $5000$ and $10\,000$ at $t_\gamma =2.5$ . Black dashed lines indicate the dominant frequency peaks associated with lift coefficients, see figure 5.

For the short time scale, $t_\gamma =0.8$ , we observe that the overall distributions of $\sigma _1,\sigma _2$ and $\sigma _3$ show some similarities across different Reynolds numbers. For $\sigma _1$ , the maximum gain is achieved at $\beta =0$ . The singular values $\sigma _2$ and $\sigma _3$ are instead more sensitive to the spanwise variation. In particular, the variation of $\sigma _2$ and $\sigma _3$ across $\beta$ and $St$ are similar, with maximum values achieved at low frequency and spanwise wavenumbers that increase with the Reynolds number. Overall, the second and third singular modes, even at short time scales, seem to reflect the development of smaller structures in the flow when increasing the Reynolds number, while the first singular value mostly reflects the two-dimensional dynamics.

The matter changes when we consider longer time scales, as shown in figure 15. Increasing $t_\gamma$ , several mode switchings are observed. These appear at low frequencies, particularly close to the characteristic frequencies at the highest Reynolds numbers and evident from a change in the gain variation. This indicates that higher-order modes, which at short time scale reflect the relevance of finer spanwise structures, need more time to grow and overcome the energy of two-dimensional mechanisms that prevail at short time scale. At ${\textit{Re}}=10\,000$ , we show a zoomed-in view of the low- $St$ and high- $\beta$ parametric space, showing the emergence of a local maximum at ${\textit{St}}\approx 0.36$ and $\beta \approx 4\pi$ , which reflects the large real part eigenvalue presented in figure 13.

Spanwise effects are thus seen to affect the low-frequency dynamics over long time scale. At the highest Reynolds number, we observe that the modal structures at low frequencies and high spanwise wavenumber ( $\beta$ ) predominantly affect the recirculation region, which presents elliptical streamlines. At lower $\beta$ , the dynamics are concentrated in the wake, which was also observed for the $\beta = 0$ case in the previous section. In our analysis, elliptic instabilities emerge clearly over long time horizons, suggesting that while their growth is weak, they represent persistent, spatially localized structures in the recirculation zone. This is in contrast to higher-frequency modes, which show maximum gain near $\beta = 0$ and remain largely unaffected by the time scale. These modes are linked to shear-layer dynamics and correspond to the quasi-two-dimensional roll-up of the shear layer separating from the leading edge, as seen in figure 4(c,d).

The different frequency–wavenumber responses are summarized in figure 16. An analogous division of mechanisms was proposed by Pickering et al. (Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020) for jet flows, associating low-frequency/low-wavenumber modes with the Orr mechanism, high-frequency modes with the Kelvin–Helmholtz instability, as in our case, and low-frequency/high-wavenumber modes with streak-like features.

High- $\beta$ , low-frequency phenomena might also be linked to bursting. Indeed, numerical simulations have shown that bursting occurs when three-dimensional disturbances in the aft part of the recirculation region grow to levels sufficient to break up the shear-layer roll-up vortices, reducing their spanwise coherence (Marxen & Henningson Reference Marxen and Henningson2011; Toppings & Yarusevych Reference Toppings and Yarusevych2024).

Figure 16. Response modal structures in the frequency–wavenumber space for the ${\textit{Re}} = 10\,000$ case at $t_\gamma =2.5$ , shown alongside the corresponding primary gain $\sigma _1$ contour map. Shown are the vertical velocity components of selected response modes marked by the red dots, superimposed on base flow streamlines within the recirculation region.

Previous studies at lower Reynolds numbers also report that the two-dimensional vortex shedding mechanism is predominant (He et al. Reference He, Gioria, Pérez and Theofilis2017; Nastro et al. Reference Nastro, Robinet, Loiseau, Passaggia and Mazellier2023), in agreement with the present results at ${\textit{Re}}=1000$ . By increasing the Reynolds number, this is true for the shear layer, which present structures with high energy gains at $\beta =0$ , while less energetic three-dimensional mechanisms appear at lower frequencies.