2.1 Problem description

We consider separated flows over a NACA 0012 airfoil at two angles of attack of $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ for a moderate chord-based Reynolds number $Re_{L_{c}}\equiv v_{\infty }L_{c}/\unicode[STIX]{x1D708}_{\infty }=23\,000$ and a free-stream Mach number $M_{\infty }\equiv v_{\infty }/a_{\infty }=0.3$ , as shown in figure 2. Here, $v_{\infty }$ is the free-stream velocity, $L_{c}$ is the chord length, $a_{\infty }$ is the free-stream sonic speed and $\unicode[STIX]{x1D708}_{\infty }$ is the kinematic viscosity. To perform active flow control, a thermal actuator is placed across the span near the leading edge. This actuator introduces oscillatory heat flux at a prescribed frequency and spanwise profile as an open-loop actuation input. The details of this thermal actuator will be discussed in § 2.3.

2.2 Simulation set-up

Figure 2. The problem description. Separated flow over a NACA 0012 airfoil (shown for $\unicode[STIX]{x1D6FC}=6^{\circ }$ ) at free-stream Mach number $M_{\infty }=0.3$ and chord-based Reynolds number $Re_{L_{c}}=23\,000$ .

We perform LES to simulate spanwise-periodic flows over the airfoil using the finite-volume compressible flow solver CharLES (Khalighi et al. Reference Khalighi, Nichols, Ham, Lele and Moin2011; Brès et al. Reference Brès, Ham, Nichols and Lele2017), which is second-order accurate in space and third-order accurate in time. Vremen’s subgrid-scale model (Vreman Reference Vreman2004) is utilized in the LES. We utilize a C-shaped computational mesh, with the airfoil positioned with its leading edge at $x/L_{c}=y/L_{c}=0$ , as shown in figure 3. The extent of the computational domain is $x/L_{c}\in [-19,26]$ , $y/L_{c}\in [-20,20]$ and $z/L_{c}\in [-0.1,0.1]$ in the streamwise, transverse and spanwise directions, respectively. This domain is discretized with approximately 35 million grid cells. We have examined the grid convergence by comparing the flow field and aerodynamics forces from this mesh to two other meshes that are further refined in the near field with a total of 63 and 82 million grid cells. From each mesh, the force data are collected for the developed flow over 80 convective time units. The time-averaged wall-normal velocity profiles were converged within $1.5\,\%$ and aerodynamic forces were also observed to be insensitive to the grid resolution of the three meshes.

For the fluid properties, we use the specific heat ratio $\unicode[STIX]{x1D6FE}=1.4$ and the Prandtl number $Pr=0.7$ , which are representative for standard air. The temperature-varying dynamic viscosity, $\unicode[STIX]{x1D707}(T)$ , is evaluated with a power law as $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{\infty }(T/T_{\infty })^{0.76}$ , where $\unicode[STIX]{x1D707}_{\infty }$ and $T_{\infty }$ are the free-stream dynamic viscosity and temperature, respectively (Garnier, Adams & Sagaut Reference Garnier, Adams and Sagaut2009). The power law models the dynamic viscosity variation for standard air in the range $T/T_{\infty }\in [0.5,1.7]$ . This range is suitable for the current study with local thermal inputs, where we observe that the maximum temperature fluctuation is within $42\,\%$ of $T_{\infty }$ for all controlled flows.

Figure 3. The computational domains ( $x$ – $y$ plane, shown for $\unicode[STIX]{x1D6FC}=6^{\circ }$ ) for LES and resolvent analysis. The near-field mesh (right) is shown along with instantaneous spanwise vorticity from LES and streamwise velocity mode from resolvent analysis. For both meshes, uniform $\unicode[STIX]{x0394}x$ is adopted in $x/L_{c}\in [1.5,6]$ to resolve the wake structures.

The simulations are performed with the Dirichlet boundary condition specified at the far-field boundary as $[\unicode[STIX]{x1D70C},v_{x},v_{y},v_{z},T]=[\unicode[STIX]{x1D70C}_{\infty },v_{\infty },0,0,T_{\infty }]$ , where $\unicode[STIX]{x1D70C}$ is the density, $v_{x}$ , $v_{y}$ and $v_{z}$ are respectively the streamwise, transverse and spanwise velocities and $T$ is the temperature. Over the airfoil, the no-slip adiabatic boundary condition is prescribed, except for where the actuator is placed for controlled cases. Along the outlet boundary, a sponge layer (Freund Reference Freund1997) is applied over $x/L_{c}\in [15,25]$ with the target state set to the running-averaged flow over 10 acoustic time units. Time integration is performed at a constant time step of $\unicode[STIX]{x0394}tv_{\infty }/L_{c}=4.14\times 10^{-5}$ , corresponding to a maximum Courant–Friedrichs–Lewy number of $0.84$ . Further details regarding the meshing strategy and computational set-up are reported in Yeh, Munday & Taira (Reference Yeh, Munday and Taira2017a ).

2.3 Actuator model

The thermal actuator is implemented as an oscillatory energy-flux boundary condition to model the fundamental effects of thermoacoustic and plasma-based actuators in the LES. It is prescribed in the energy equation as an unsteady Neumann boundary condition, along with no-slip boundary condition for the momentum equation in the compressible Navier–Stokes equations. The actuator model is expressed as

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}^{+}(\unicode[STIX]{x1D714}^{+},k_{z}^{+})=\frac{1}{4}\hat{\unicode[STIX]{x1D719}}\sin (\unicode[STIX]{x1D714}^{+}t)[1+\cos (k_{z}^{+}z)]\left\{1+\cos \left[\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D70E}_{a}}(x-x_{a})\right]\right\}, & & \displaystyle\end{eqnarray}$$

where $(x-x_{a})/\unicode[STIX]{x1D70E}_{a}\in [-0.5,0.5]$ . This expression provides the boundary heat flux input with a compact spatial support in the form of a Hanning window centred at $x_{a}/L_{c}=0.03$ on the suction surface with width of $\unicode[STIX]{x1D70E}_{a}/L_{c}=0.04$ , as illustrated in figure 2. The actuator introduces the open-loop control input at the prescribed actuation frequency, $\unicode[STIX]{x1D714}^{+}$ , and spanwise wavenumber, $k_{z}^{+}$ . They are parameterized in the LES of controlled cases and reported in terms of the actuation Strouhal number $St^{+}=\unicode[STIX]{x1D714}^{+}L_{c}/(2\unicode[STIX]{x03C0}v_{\infty })$ and the normalized wavenumber $k_{z}^{+}L_{c}$ throughout this paper. Due to the choice of the spanwise extent for the computational domain ( $z/L_{c}\in [-0.1,0.1]$ ), the actuation wavenumbers $k_{z}^{+}L_{c}$ are restricted to integer multiples of $10\unicode[STIX]{x03C0}$ . Hence, we only consider the use of $k_{z}^{+}L_{c}=0$ , $10\unicode[STIX]{x03C0}$ , $20\unicode[STIX]{x03C0}$ and $40\unicode[STIX]{x03C0}$ in the LES of controlled flows. In the actuator model (2.1), the actuation amplitude $\hat{\unicode[STIX]{x1D719}}$ is selected such that the normalized total actuation power

(2.2) $$\begin{eqnarray}\displaystyle E^{+}=\frac{\frac{1}{4}\hat{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D70E}_{a}}{\frac{1}{2}\unicode[STIX]{x1D70C}_{\infty }v_{\infty }^{3}(L_{c}\sin \unicode[STIX]{x1D6FC})}=0.0902 & & \displaystyle\end{eqnarray}$$

for all controlled cases throughout this work. This magnitude is representative of those used in thermally actuated flow control studies (Corke, Enloe & Wilkinson Reference Corke, Enloe and Wilkinson2010; Sinha et al. Reference Sinha, Alkandry, Kearney-Fischer, Samimy and Colonius2012; Akins, Singh & Little Reference Akins, Singh and Little2015; Yeh et al. Reference Yeh, Munday, Taira and Munson2015). For this thermal actuator, Yeh, Munday & Taira (Reference Yeh, Munday and Taira2017b ) have investigated its control mechanism and flow control capability in free shear layers. The thermal input from the actuator translates to vortical perturbations in the forms of oscillatory surface vorticity flux and baroclinic torque. The thermal actuation is capable of exciting fundamental and subharmonic instabilities, which is suitable for modifying the shear-layer dynamics herein.

2.4 Baseline simulations

We validate the baseline simulations at angles of attack of $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ by comparing the surface pressure distribution and the aerodynamic forces to those reported in the literature for $Re_{L_{c}}=23\,000$ . Throughout this study, the pressure coefficient $C_{p}$ , lift coefficient $C_{L}$ and drag coefficients $C_{D}$ are defined as

(2.3a-c ) $$\begin{eqnarray}\displaystyle C_{p}=\frac{p-p_{\infty }}{\frac{1}{2}\unicode[STIX]{x1D70C}_{\infty }v_{\infty }^{2}},\quad C_{L}=\frac{F_{L}}{\frac{1}{2}\unicode[STIX]{x1D70C}_{\infty }v_{\infty }^{2}A},\quad C_{D}=\frac{F_{D}}{\frac{1}{2}\unicode[STIX]{x1D70C}_{\infty }v_{\infty }^{2}A}, & & \displaystyle\end{eqnarray}$$

where $F_{L}$ and $F_{D}$ are the total lift and drag forces on the airfoil, respectively, and $A$ is the planform area of the airfoil. The time-averaged aerodynamic forces and surface pressure profile are respectively presented in table 1 and figure 4. We found reasonable agreements with those reported by Kim et al. (Reference Kim, Yang, Chang and Chung2009), Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013) and Munday & Taira (Reference Munday and Taira2018). We note that the numerical study of Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013) was conducted using implicit LES and Munday & Taira (Reference Munday and Taira2018) reported the results from incompressible LES. The discrepancy in the surface pressure with the experimental measurement by Kim et al. (Reference Kim, Yang, Chang and Chung2009) can be attributed to the different transverse blockage ratios ( $0.26\,\%$ for the present study).

Table 1. The time-averaged drag and lift coefficients on a NACA 0012 airfoil at $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ at $Re_{L_{c}}=23\,000$ . The present study performs compressible LES at a free-stream Mach number $M_{\infty }=0.3$ , in comparison with the results from incompressible LES by Munday & Taira (Reference Munday and Taira2018) and implicit LES by Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013) at $M_{\infty }=0.2$ .

Figure 4. Surface pressure profiles for (a) $\unicode[STIX]{x1D6FC}=6^{\circ }$ and (b) $\unicode[STIX]{x1D6FC}=9^{\circ }$ over the chord-wise coordinate $x_{c}/L_{c}=(x\cos \unicode[STIX]{x1D6FC}+y\sin \unicode[STIX]{x1D6FC})/L_{c}$ , in comparison with those reported by Kim et al. (Reference Kim, Yang, Chang and Chung2009), Kojima et al. (Reference Kojima, Nonomura, Oyama and Fujii2013) and Munday & Taira (Reference Munday and Taira2018).

Figure 5. (a,b) Baseline flow visualization using $Q$ -criterion (iso-surface of $QL_{c}^{2}/v_{\infty }^{2}=50$ coloured by streamwise velocity) and spanwise-average $\text{TKE}=\overline{(v_{x}^{\prime 2}+v_{y}^{\prime 2}+v_{z}^{\prime 2})}/v_{\infty }^{2}$ in the background. (c,d) The time-averaged streamlines. The contour line for $\bar{v}_{x}=0$ is shown and will be used for characterizing the extent of the separation region throughout this study.

The instantaneous flow fields and time-average streamlines for the baseline flows at $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ are shown in figure 5. The iso-surface of $Q$ -criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988) is used to visualize the vortical structures. The contour line of time- and spanwise-averaged streamwise velocity $\bar{v}_{x}=0$ is also shown to identify the flow separation and reattachment. This contour line is also shown on top of the time-average streamlines, where we see the contour line extends through the separation bubble for each case. For both angles of attack, laminar separation is observed near the leading edge and forms a shear layer. The shear layer rolls up over the suction surface and evolves into spanwise vortices. This roll-up process leads to increasing turbulent kinetic energy (TKE) within the shear layer. Farther downstream, these spanwise vortices break up and lose their spanwise coherence, resulting in the laminar–turbulent transition. Within this roll-up and transition process, one common feature in the pressure profiles in figure 4 is the ‘plateau’ observed for both angles of attack. Such a plateau in the pressure profile is also observed by Marxen, Lang & Rist (Reference Marxen, Lang and Rist2013) and Benton & Visbal (Reference Benton and Visbal2018) in the transition process that takes place over a laminar separation bubble. The transition process is accompanied by the maximum TKE over the airfoil at $x/L_{c}\approx 0.6$ for $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $x/L_{c}\approx 0.5$ for $\unicode[STIX]{x1D6FC}=9^{\circ }$ . The roll-up and break-up processes both result in momentum mixing and entraining of the free stream, leading to flow reattachment for $\unicode[STIX]{x1D6FC}=6^{\circ }$ at $x/L_{c}\approx 0.85$ . Over the airfoil at $\unicode[STIX]{x1D6FC}=9^{\circ }$ , the flow is in full stall. To quantitatively characterize the stall condition, we calculate the potential-flow lift $C_{L,0}$ using the panel method (Hess Reference Hess1990) to mark a theoretical upper bound of the lift for both angles of attack. The flow over the airfoil at $\unicode[STIX]{x1D6FC}=6^{\circ }$ reattaches and achieves $84\,\%$ of $C_{L,0}$ . Whereas for the $\unicode[STIX]{x1D6FC}=9^{\circ }$ airfoil, while experiencing deep stall, provides only $52\,\%$ of the potential flow lift. This difference in the stall condition will be reflected in the control effectiveness to be discussed in § 4.

The excitation of shear-layer instabilities serves as the key to separation control (Greenblatt & Wygnanski Reference Greenblatt and Wygnanski2000). For the laminar separation bubble that is observed in both baseline flows, Häggmark, Bakchinov & Alfredsson (Reference Häggmark, Bakchinov and Alfredsson2000) have experimentally shown that the Kelvin–Helmholtz instability dominates the laminar–turbulent transition. In order to leverage the Kelvin–Helmholtz instability for flow control, we place the thermal actuator slightly upstream of the separation point such that the perturbations can be introduced at the onset of the shear layer.

3.1 Formulation

Let us consider the compressible Navier–Stokes equations expressed as

(3.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}t}={\mathcal{N}}(\boldsymbol{q})+\boldsymbol{f}^{+}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{N}}$ is the nonlinear Navier–Stokes operator that acts on the flow state variable $\boldsymbol{q}=[\unicode[STIX]{x1D70C},v_{x},v_{y},v_{z},T]^{\text{T}}$ and $\boldsymbol{f}^{+}$ represents the external actuation input from active flow control. Note that in general the external forcing $\boldsymbol{f}^{+}$ can be absent. We perform the Reynolds decomposition of $\boldsymbol{q}=\bar{\boldsymbol{q}}+\check{\boldsymbol{q}}$ so that the flow state variable $\boldsymbol{q}$ is decomposed into a statistically stationary long-time mean component $\bar{\boldsymbol{q}}$ and a fluctuating component $\check{\boldsymbol{q}}$ . Substituting $\boldsymbol{q}$ with its Reynolds decomposition into the Navier–Stokes equations (3.1) yields

(3.2) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\bar{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\check{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}=\underbrace{{\mathcal{L}}(\check{\boldsymbol{q}},\bar{\boldsymbol{q}}\unicode[STIX]{x1D735}\check{\boldsymbol{q}},\check{\boldsymbol{q}}\unicode[STIX]{x1D735}\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FB}^{2}\check{\boldsymbol{q}},\cdots \,)}_{{\mathcal{L}}_{\bar{\boldsymbol{ q}}}(\check{\boldsymbol{q}})}+\underbrace{{\mathcal{N}}(\bar{\boldsymbol{q}})+\bar{f}(\check{\boldsymbol{q}}^{n})+\boldsymbol{f}^{+}}_{\check{\boldsymbol{u}}}. & & \displaystyle\end{eqnarray}$$

With the Reynolds decomposition, the linear operations for $\check{\boldsymbol{q}}$ are extracted from the operation of ${\mathcal{N}}(\bar{\boldsymbol{q}}+\check{\boldsymbol{q}})$ . We collect these terms that are linear with respect to $\check{\boldsymbol{q}}$ and denote them as ${\mathcal{L}}_{\bar{\boldsymbol{q}}}(\check{\boldsymbol{q}})$ . The term ${\mathcal{N}}(\bar{\boldsymbol{q}})$ accounts for the Navier–Stokes operation taking place only on $\bar{\boldsymbol{q}}$ , and $\bar{f}(\check{\boldsymbol{q}}^{n})$ collects the nonlinear higher-order terms for $\check{\boldsymbol{q}}$ in $O(\boldsymbol{q}^{n})$ , where $n>1$ . In particular, we note that ${\mathcal{N}}(\bar{\boldsymbol{q}})+\bar{f}(\check{\boldsymbol{q}}^{n})$ can be interpreted as the internal forcing in the turbulent flow due to the nonlinear interaction (Farrell & Ioannou Reference Farrell and Ioannou1994; McKeon & Sharma Reference McKeon and Sharma2010). This internal forcing together with the external forcing $\boldsymbol{f}^{+}$ is further denoted as $\check{\boldsymbol{u}}$ . Observing that the base flow is taken to be time-invariant, i.e.  $\unicode[STIX]{x2202}_{t}\bar{\boldsymbol{q}}=0$ , equation (3.2) can be simplified as

(3.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\check{\boldsymbol{q}}}{\unicode[STIX]{x2202}t}={\mathcal{L}}_{\bar{\boldsymbol{q}}}(\check{\boldsymbol{q}})+\check{\boldsymbol{u}}. & & \displaystyle\end{eqnarray}$$

We note that, for a statistically stationary turbulent flow where $\bar{\boldsymbol{q}}$ is known a priori and is not affected by the sustained forcing input $\check{\boldsymbol{u}}$ , the linear operator ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ is also time-invariant and can be explicitly constructed with the knowledge of $\bar{\boldsymbol{q}}$ . Therefore, equation (3.3) can be viewed as a linear system for $\check{\boldsymbol{q}}$ . The nonlinear interaction in $\check{\boldsymbol{u}}$ supports the base flow $\bar{\boldsymbol{q}}$ via its mean component and forces $\check{\boldsymbol{q}}$ with its fluctuating component (McKeon & Sharma Reference McKeon and Sharma2010). Thus far, no assumptions have been made in the formulation except for the statistical stationarity of the turbulent flow about which the Navier–Stokes equations are rewritten in the above form.

Now, we cast equation (3.3) for the spanwise-periodic flow over the airfoil. Considering the two-dimensional airfoil geometry in this study, the time- and spanwise-average flow obtained from the baseline flow simulation is used as the mean component so that the Reynolds decomposition can be written as

(3.4) $$\begin{eqnarray}\displaystyle \boldsymbol{q}(x,y,z,t)=\bar{\boldsymbol{q}}(x,y)+\check{\boldsymbol{q}}(x,y,z,t). & & \displaystyle\end{eqnarray}$$

The spanwise-periodic setup in the present study allows for the biglobal-mode representation for $\check{\boldsymbol{q}}$ and $\check{\boldsymbol{u}}$ as the sum of temporal and spanwise Fourier modes (Theofilis Reference Theofilis2003) respectively as

(3.5) $$\begin{eqnarray}\displaystyle \check{\boldsymbol{q}}(x,y,z,t)=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}(x,y)\text{e}^{\text{i}(k_{z}z-\unicode[STIX]{x1D714}t)}\,\text{d}\unicode[STIX]{x1D714}\,\text{d}k_{z} & & \displaystyle\end{eqnarray}$$

and

(3.6) $$\begin{eqnarray}\displaystyle \check{\boldsymbol{u}}(x,y,z,t)=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}(x,y)\text{e}^{\text{i}(k_{z}z-\unicode[STIX]{x1D714}t)}\,\text{d}\unicode[STIX]{x1D714}\,\text{d}k_{z}. & & \displaystyle\end{eqnarray}$$

Here, $\text{i}=\sqrt{-1}$ , $\unicode[STIX]{x1D714}$ is the complex radian frequency, $k_{z}$ is the real spanwise wavenumber and $\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}$ and $\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}$ are the biglobal modes for spanwise wavenumber $k_{z}$ and temporal frequency $\unicode[STIX]{x1D714}$ . Substituting the modal expressions (3.5) and (3.6) for $\check{\boldsymbol{q}}$ and $\check{\boldsymbol{u}}$ into (3.3), we arrive at the linearized Navier–Stokes equations in Fourier space:

(3.7) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x1D714}\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}={\mathcal{L}}_{\bar{\boldsymbol{q}}}(\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}};k_{z})+\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}. & & \displaystyle\end{eqnarray}$$

By treating $\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}$ as a known forcing term, equation (3.7) (or (3.3) equivalently) represents an inhomogeneous linear differential equation that governs the time evolution of perturbation $\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}$ , with $\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}$ being the inhomogeneous forcing term on the right-hand side. Its general solution comprises a homogeneous solution and a particular solution. The homogeneous solution can be found by solving equation (3.7) without the forcing term. That is,

(3.8) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x1D714}\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}={\mathcal{L}}_{\bar{\boldsymbol{q}}}(\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}};k_{z}), & & \displaystyle\end{eqnarray}$$

which forms an eigenvalue problem, which reminds us that the homogeneous solution is associated with the spectrum of ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ . On the other hand, the particular solution of (3.7) can be expressed as

(3.9) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}=[-\text{i}\unicode[STIX]{x1D714}-{\mathcal{L}}_{\bar{\boldsymbol{q}}}(k_{z})]^{-1}\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}, & & \displaystyle\end{eqnarray}$$

where the operator

(3.10) $$\begin{eqnarray}\displaystyle {\mathcal{H}}_{\bar{\boldsymbol{q}}}(k_{z},\unicode[STIX]{x1D714})=[-\text{i}\unicode[STIX]{x1D714}-{\mathcal{L}}_{\bar{\boldsymbol{q}}}(k_{z})]^{-1}, & & \displaystyle\end{eqnarray}$$

is referred to as the resolvent and is associated with the pseudospectrum of ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ (Trefethen & Embree Reference Trefethen and Embree2005).

Our objective is not to solve differential equation (3.7), which requires knowledge of the initial condition and the explicit forcing $\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}$ . However, we characterize its general solution by analysing the spectrum and pseudospectrum of the linear operator ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ . Moreover, we note that particular solution (3.9) describes a linear operation that takes place between a sustained input $\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}$ and the harmonic output $\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}$ through the resolvent operator ${\mathcal{H}}_{\bar{\boldsymbol{q}}}(k_{z},\unicode[STIX]{x1D714})$ . For this reason, the pseudospectrum of ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ , which captures the energy amplification through the input–output process, is the main focus of this study on active flow control.

With the knowledge of $\bar{\boldsymbol{q}}$ and the appropriate boundary conditions, the linear operator ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ can be explicitly constructed in its discretized form $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ for a prescribed spanwise wavenumber $k_{z}$ . We can rewrite equation (3.7) in discrete form as

(3.11) $$\begin{eqnarray}\displaystyle -\text{i}\unicode[STIX]{x1D714}\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}=\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}+\hat{\boldsymbol{u}}_{k_{z},\unicode[STIX]{x1D714}}, & & \displaystyle\end{eqnarray}$$

where the operation of ${\mathcal{L}}_{\bar{\boldsymbol{q}}}$ on $\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}$ is represented by a matrix–vector multiplication of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}$ . The modal wavenumber $k_{z}$ is embedded in $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ since it emerges from the spatial differentiation in the construction of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ . With the discrete linear operator $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ constructed, its spectrum and pseudospectrum can be found numerically. Below, we document the domain discretization and boundary conditions for constructing the discrete linear operator $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ . The numerical approach for computing its spectrum and pseudospectrum is also offered.

3.2 Numerical set-up

The discretization for equation (3.7) is performed on the computational mesh shown in figure 3 highlighted in orange on top of the LES domain. This two-dimensional rectangular domain has an extent of $x/L_{c}\in [-15,16]$ , $y/L_{c}\in [-12,12]$ and is composed of approximately $0.14$ million grid points. The mean turbulent flow $\bar{\boldsymbol{q}}=[\bar{\unicode[STIX]{x1D70C}},\bar{u},\bar{v},\bar{w},\bar{T}]^{\text{T}}$ obtained from the baseline simulation on the LES mesh is interpolated onto this coarser mesh to perform the stability (spectral) and resolvent (pseudospectral) analyses. For the far-field boundary and over the airfoil, the Dirichlet boundary condition is set for $[\unicode[STIX]{x1D70C}^{\prime },u^{\prime },v^{\prime },w^{\prime }]=[0,0,0,0]$ and the Neumann boundary condition is set for $T^{\prime }$ such that $\boldsymbol{e}_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T^{\prime }=0$ , where $\boldsymbol{e}_{n}$ is the unit normal boundary vector. At the outlet boundary, the same Neumann boundary condition is set for all flow variables. With these boundary conditions and the baseline mean flow $\bar{\boldsymbol{q}}$ , we construct the linear operator in its discrete form $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ for a chosen spanwise wavenumber $k_{z}$ .

In the current study, the size of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ is approximately $0.7$ million $\times$ $0.7$ million. Considering the large size of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ , the implicitly restarted Arnoldi method (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998) is used to handle the large-scale eigenvalue problem to solve for its spectrum and pseudospectrum. The eigenvalues and the resolvent norm (for pseudospectrum) are computed with a Krylov space of $128$ vectors and a residual tolerance of $10^{-10}$ . The domain size and mesh resolution were examined to ensure that the results converge to at least six significant digits. Further detail of grid convergence is documented in appendix A.

3.3 Spectrum and pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$

The mean-flow-based linear operator $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ can be characterized by its spectrum (eigenvalues) and pseudospectrum. They describe the dynamical response of the fluid-flow system as they arise from the general solution of the linearized Navier–Stokes equations (3.7).

3.3.1 Spectrum

The eigenvalue problem arising from the homogeneous problem (3.8) can be expressed in its discretized form

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}=-\text{i}\unicode[STIX]{x1D714}\hat{\boldsymbol{q}}_{k_{z},\unicode[STIX]{x1D714}}, & & \displaystyle\end{eqnarray}$$

where $-\text{i}\unicode[STIX]{x1D714}$ and $\hat{\boldsymbol{q}}$ are the eigenvalue and eigenmode, respectively. The eigenvalue $-\text{i}\unicode[STIX]{x1D714}=-\text{i}\unicode[STIX]{x1D714}_{r}+\unicode[STIX]{x1D714}_{i}$ determines the temporal stability with modal frequency $\unicode[STIX]{x1D714}_{r}$ and growth (or decay) rate $\unicode[STIX]{x1D714}_{i}$ . An instability is identified if the complex modal frequency $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ resides on the positive imaginary plane with $\unicode[STIX]{x1D714}_{i}>0$ . Upon prescribing a modal wavenumber $k_{z}$ for $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ , the eigenvalue problem (3.12) can be viewed as the biglobal linear stability analysis (Theofilis Reference Theofilis2011) at $k_{z}$ with the mean turbulent flow $\bar{\boldsymbol{q}}$ as the base state.

Figure 6. Spectrum (a) and dominant eigenmodes (b) of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z}=0)$ for $\unicode[STIX]{x1D6FC}=9^{\circ }$ mean flow. Three representative dominant modes shown in (b) are circled with $\circ$ (blue) or $\circ$ (red) in the spectrum (a). Spurious eigenvalues that associate with unphysical standing waves (Lesshafft Reference Lesshafft2018) are coloured in grey in the spectrum (a). An eigenmode is considered to be physical if $99.9\,\%$ of its energy (equation (3.14)) lies in the box of $y/L_{c}\in [-5,5]$ and $x/L_{c}\in [-2,16]$ , where $x/L_{c}=16$ is the computational outlet where Neumann boundary condition is imposed. The magenta dashed lines in (a) highlight the frequencies of the dominant wake modes and are also shown in the frequency spectra in figure 7(a,b). The shear layer over the separation bubble is identified by the time-averaged spanwise vorticity $\bar{\unicode[STIX]{x1D701}}_{z}$ as shown in (c) and is marked with dashed lines to highlight the shear-layer structures in the eigenmodes.

We show in figure 6 the results of the spectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z}=0)$ and three representative eigenmodes for $\unicode[STIX]{x1D6FC}=9^{\circ }$ . We note that the spectrum is symmetric about the $\unicode[STIX]{x1D714}_{i}$ axis, since the modal phase velocity does not exhibit preferential spanwise direction due to the two-dimensional geometry of the airfoil. Thus, in figure 6, we only show the spectrum on the positive frequency plane ( $\unicode[STIX]{x1D714}_{r}\geqslant 0$ ). In the spectrum, two branches can be identified: the wake-mode branch and the shear-layer-mode branch. These two branches can be characterized by the frequency bandwidth of the eigenvalues or through examination of their modal structures. Three eigenmodes are chosen in the spectrum with $\circ$ and their modal structures are visualized in figure 6(b) with the streamwise velocity profile $\hat{u}$ : (1) the dominant shear-layer mode; (2) the dominant wake mode; and (3) a coupling mode of shear layer and wake. On top of each modal structure, a dashed line is shown to mark the location of the time-averaged shear layer. This line is determined by examining the time-averaged spanwise vorticity $\bar{\unicode[STIX]{x1D701}}_{z}$ for its local maximum magnitude over the separation bubble, as shown in figure 6(c). The shear-layer mode presents distinctively strong structure along the shear layer. While the shear-layer mode gradually vanishes in the wake, the wake-mode structure extends farther downstream and resembles the pattern of the von Kármán vortex street behind a bluff body. On the wake branch, the frequencies of four dominant modes are highlighted with magenta lines. These frequencies, marked again in the frequency spectrum of lift ${\hat{C}}_{L}$ in figure 7(b), are found to be in agreement with the peaks obtained from LES. Similar agreement holds for $\unicode[STIX]{x1D6FC}=6^{\circ }$ results in figure 7(a). The agreement between the $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ spectrum and the dominant frequency identified from the baseline flow shows that the nonlinear vortex-shedding physics can be revealed by the linear analysis. Comparing the lift spectra for $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ in figure 7(a,b), we find that the frequency content of ${\hat{C}}_{L}$ scales well with the frontal-height-based Strouhal number $St_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D714}(L_{c}\sin \unicode[STIX]{x1D6FC})/2\unicode[STIX]{x03C0}v_{\infty }$ . The $St_{\unicode[STIX]{x1D6FC}}$ scaling for the lift spectra has been studied by Fage & Johansen (Reference Fage and Johansen1927), reporting the appearance of the ${\hat{C}}_{L}$ peaks near $St_{\unicode[STIX]{x1D6FC}}\approx 0.2$ .

The linear operator $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ is observed to be unstable for $k_{z}L_{c}=0$ as it possesses eigenvalues with positive growth rates. In fact, $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ is found to be unstable for $k_{z}L_{c}\lesssim 8\unicode[STIX]{x03C0}$ . The identification of the critical $k_{z}$ that yields instability in $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ is not the focus the present study. However, we leave a cautionary note here that its unstable nature for low $k_{z}$ necessitates further care when performing the resolvent analysis of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ , which will be discussed in detail in § 3.5.

3.3.2 Pseudospectrum

Figure 7. Frequency spectra of lift ${\hat{C}}_{L}$ from baseline LES (a,b) and pseudospectra of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ with $k_{z}L_{c}=0$ (c,d) for both $\unicode[STIX]{x1D6FC}=6^{\circ }$ (top) and $9^{\circ }$ (bottom). The pseudospectra are constructed over the complex $\unicode[STIX]{x1D714}$ -plane by seeking the leading singular value $\unicode[STIX]{x1D70E}_{1}$ in (3.16). Magenta dots in (c,d) depict the eigenvalues of the corresponding $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ . Over the horizontal axis, two frequency scales are provided: the Fage–Johansen Strouhal number $St_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D714}(L_{c}\sin \unicode[STIX]{x1D6FC})/2\unicode[STIX]{x03C0}v_{\infty }$ on the upper axis and the chord-based Strouhal number $St=\unicode[STIX]{x1D714}L_{c}/2\unicode[STIX]{x03C0}v_{\infty }$ on the lower axis. In each panel, the magenta dashed lines mark the frequencies of dominant wake modes from the spectra of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ .

A normal operator satisfies $\unicode[STIX]{x1D647}\unicode[STIX]{x1D647}^{\ast }=\unicode[STIX]{x1D647}^{\ast }\unicode[STIX]{x1D647}$ , where the superscript $^{\ast }$ denotes the Hermitian transpose. It has orthonormal eigenmodes with corresponding eigenvalues that govern the dynamical behaviour. For a non-normal operator (i.e.  $\unicode[STIX]{x1D647}\unicode[STIX]{x1D647}^{\ast }\neq \unicode[STIX]{x1D647}^{\ast }\unicode[STIX]{x1D647}$ ), its transient behaviour is not described simply by the eigenvalues and eigenvectors. Instead of just the spectrum, the pseudospectrum is needed to analyse the dynamics resulting from a non-normal operator. Trefethen & Embree (Reference Trefethen and Embree2005) examined pseudospectra of non-normal operators and explained how they align with the dynamical behaviours governed by these operators. In fluid-flow systems, non-normality appears in shear-dominant flows (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Schmid & Henningson Reference Schmid and Henningson2001; McKeon & Sharma Reference McKeon and Sharma2010). From the baseline flows, we readily identify the presence of strong shear particularly over the separation bubble. This region of strong shear can be recognized in the mean flow profile about which $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ is constructed.

We have mentioned that the pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ arises from the resolvent operator ${\mathcal{H}}_{\bar{\boldsymbol{q}}}(k_{z},\unicode[STIX]{x1D714})$ in the particular solution (3.9). Here, we work with the discrete resolvent operator

(3.13) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}(k_{z},\unicode[STIX]{x1D714})=[-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})]^{-1}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix. The pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ is to be mapped out over the complex $\unicode[STIX]{x1D714}$ plane by seeking a 2-norm measure through the SVD of its resolvent matrix $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}$ . An appropriate 2-norm for this study can be introduced as the weighted inner product between two state vectors:

(3.14) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{q}_{1},\boldsymbol{q}_{2}\rangle _{E}=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{q}_{1}^{\ast }\mathsf{diag}\left(\frac{R\bar{T}}{\bar{\unicode[STIX]{x1D70C}}},\bar{\unicode[STIX]{x1D70C}},\bar{\unicode[STIX]{x1D70C}},\bar{\unicode[STIX]{x1D70C}},\frac{R\bar{\unicode[STIX]{x1D70C}}}{(\unicode[STIX]{x1D6FE}-1)\bar{T}}\right)\boldsymbol{q}_{2}\,\text{d}\boldsymbol{x}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the domain of interest and $R$ is the ideal gas constant. The inner product $\langle \boldsymbol{q}_{1},\boldsymbol{q}_{2}\rangle _{E}$ is referred to as the energy norm (Schmid & Henningson Reference Schmid and Henningson2001). We adopt the compressible disturbance energy proposed by Chu (Reference Chu1965) and use this 2-norm for the computation of pseudospectra. For the discrete flow fields, the energy norm is evaluated as

(3.15) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{q}_{1},\boldsymbol{q}_{2}\rangle _{E}=\boldsymbol{q}_{1}^{\ast }\unicode[STIX]{x1D652}\boldsymbol{q}_{2}, & & \displaystyle\end{eqnarray}$$

where the weight matrix $\unicode[STIX]{x1D652}$ is the numerical quadrature that accounts for both the energy weight and spatial integration. By introducing the similarity transformation of $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}\mapsto \unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}^{1/2}\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}\unicode[STIX]{x1D652}^{-1/2}$ , the energy norm for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}$ can be handled within the 2-norm framework for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}$ (Trefethen & Embree Reference Trefethen and Embree2005). Also, the similarity transformation performed for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}$ translates to $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ and preserves its eigenvalues. The pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ with respect to the energy norm (3.14) can be evaluated through the SVD of $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}$ as

(3.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}(k_{z},\unicode[STIX]{x1D714})=\unicode[STIX]{x1D64C}_{\unicode[STIX]{x1D652}}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}^{\ast }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D64C}_{\unicode[STIX]{x1D652}}$ and $\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}$ are the left- and right-singular vector, respectively. By seeking the leading singular value $\unicode[STIX]{x1D70E}_{1}$ in $\unicode[STIX]{x1D72E}$ , the pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ is obtained at the complex $\unicode[STIX]{x1D714}$ .

Following the approach, in figure 7(c,d), we present the pseudospectra of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}(k_{z})$ with respect to the energy norm for both $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ with $k_{z}=0$ , along with the frequency spectra of the lift coefficients from LES (figure 7 a,b). For all four panels, we provide two different frequency scalings over the horizontal axes: the Fage–Johansen Strouhal number $St_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D714}(L_{c}\sin \unicode[STIX]{x1D6FC})/2\unicode[STIX]{x03C0}v_{\infty }$ on the top, and the chord-based Strouhal number $St=\unicode[STIX]{x1D714}L_{c}/2\unicode[STIX]{x03C0}v_{\infty }$ on the bottom. Comparing the results from two angles of attack, we observe that, while the lift spectra scale well with $St_{\unicode[STIX]{x1D6FC}}$ , the general behaviour of the pseudospectra agrees better with $St$ , especially in the high $\unicode[STIX]{x1D714}_{i}$ region. The pseudospectra levels spread out from the region where most of the shear-layer eigenmodes reside for both angles of attack. This observation can be explained by the highly non-normal nature of these shear-layer modes, whose structures are supported by the separation bubble above the airfoil that exhibits the strongest shear in the mean flow. The high non-normality in these shear-layer modes expands the pseudospectral radius about them such that they are centred by the roll-off in the pseudospectra levels. Therefore, instead of the $St_{\unicode[STIX]{x1D6FC}}$ scaling which emphasizes the wake physics, the shear-layer-dominated behaviour is better supported by the $St$ scaling for the pseudospectra.

3.4 Resolvent analysis for active flow control

To provide physical interpretation for the right- and left-singular vectors ( $\unicode[STIX]{x1D64C}_{\unicode[STIX]{x1D652}}$ and $\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}$ ) of the SVD (3.16), let us recall the resolvent operator as part of the particular solution:

(3.17) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}=\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}\hat{\boldsymbol{u}}. & & \displaystyle\end{eqnarray}$$

Here, we have dropped the subscripts $k_{z}$ and $\unicode[STIX]{x1D714}$ for simplicity. The similarity transformation for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}$ can be brought into the particular solution as $\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{q}}=(\unicode[STIX]{x1D652}^{1/2}\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}\unicode[STIX]{x1D652}^{-1/2})\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{u}}$ . With the SVD for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}$ in (3.16), the particular solution can be rewritten considering the energy norm as

(3.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{q}}=(\unicode[STIX]{x1D64C}_{\unicode[STIX]{x1D652}}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}^{\ast })\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{u}}. & & \displaystyle\end{eqnarray}$$

Starting from the right-hand side of this equation, we see the projection of the weighted forcing $\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{u}}$ onto the vector space spanned by the right-singular vectors $\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}$ . This projection takes the inner product with respect to the energy norm and decomposes $\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{u}}$ into the vector components in $\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}$ with a series of projection coefficients. Each forcing component is amplified by the corresponding singular value in $\unicode[STIX]{x1D72E}$ , producing a set of scaled coefficients for the corresponding left-singular vectors. The output $\unicode[STIX]{x1D652}^{1/2}\hat{\boldsymbol{q}}$ is generated through the linear combination of the left-singular vectors using this set of scaled coefficients. Thus, in the SVD of $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D652}}$ , the left-singular vectors $\unicode[STIX]{x1D64C}_{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}^{1/2}[\hat{\boldsymbol{q}}_{1},\hat{\boldsymbol{q}}_{2},\ldots ,\hat{\boldsymbol{q}}_{n}]$ can be interpreted as response modes, whereas the right-singular vectors $\unicode[STIX]{x1D650}_{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}^{1/2}[\hat{\boldsymbol{u}}_{1},\hat{\boldsymbol{u}}_{2},\ldots ,\hat{\boldsymbol{u}}_{n}]$ can be interpreted as forcing modes. Each forcing–response pair is subjected to the corresponding amplification in $\unicode[STIX]{x1D72E}=\mathsf{diag}(\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{n})$ , where $\unicode[STIX]{x1D70E}_{k}$ can be arranged in a descending order. If $\unicode[STIX]{x1D70E}_{1}\gg \unicode[STIX]{x1D70E}_{2}$ , the rank-1 assumption (McKeon & Sharma Reference McKeon and Sharma2010; Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016; Gómez et al. Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016) can be appropriately made, expecting that the input–output process is dominated by the leading forcing–response pair, i.e.  $\hat{\boldsymbol{q}}\approx \hat{\boldsymbol{q}}_{1}\unicode[STIX]{x1D70E}_{1}\langle \hat{\boldsymbol{u}}_{1},\hat{\boldsymbol{u}}\rangle _{E}$ , as long as $\langle \hat{\boldsymbol{u}}_{1},\hat{\boldsymbol{u}}\rangle _{E}$ is reasonably large. This assumption will be justified shortly with the results presented in the next section.

Figure 8. Schematic demonstration of resolvent analysis: each SVD provides an optimal forcing–response pair with the associated amplification (gain) while sweeping through frequency $\unicode[STIX]{x1D714}$ and wavenumber $k_{z}$ .

Recognizing that the SVD is performed for $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}(k_{z},\unicode[STIX]{x1D714})$ for prescribed $k_{z}$ and $\unicode[STIX]{x1D714}$ , a concept of ‘Bode plot’ can be realized by sweeping through the frequency $\unicode[STIX]{x1D714}$ for each $k_{z}$ , seeking for the leading amplification (as the ‘gain’) from each SVD (Jovanović & Bamieh Reference Jovanović and Bamieh2005). Such an approach is illustrated in figure 8, where each SVD gives a leading forcing–response pair along with the associated gain. With the Bode plot constructed based on the pseudospectral analysis of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ , efficient ways of forcing may be predicted by searching for the combination of $k_{z}$ and $\unicode[STIX]{x1D714}$ that produces high gain. Such a forcing input will be highly amplified by $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}}}$ to produce perturbation $\hat{\boldsymbol{q}}$ about $\bar{\boldsymbol{q}}$ . The amplitude of perturbation may grow beyond the validity of linear regime governed by $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ . Through nonlinearity, the highly amplified perturbation can modify the mean flow $\bar{\boldsymbol{q}}$ , which is the objective of flow control. For this reason, resolvent analysis, arising from the input–output process in the particular solution (3.17), provides insightful information for flow control design. While following this approach, we provide a couple of cautionary comments on the use of resolvent analysis in the present context:

1. (i) Even though the effective forcing input predicted by the resolvent analysis may have a good chance to modify $\bar{\boldsymbol{q}}$ , the direction of the change (e.g. increase or decrease in lift) may be beyond the insights that can be provided by the amplification. The achievement of an aerodynamically favourable change may require further knowledge, such as the structure of the harmonic response rather than just the knowledge of amplifications.

2. (ii) Once the base flow $\bar{\boldsymbol{q}}$ is modified with control, the results from the analysis performed with respect to the operator for uncontrolled base state $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ may no longer be valid. However, resolvent analysis can still provide valuable insights for the effective forcing before the system departs from the linear regime about the uncontrolled $\bar{\boldsymbol{q}}$ .

We have presented a control-oriented interpretation of the results from resolvent analysis. Traditionally, resolvent analysis used in fluid mechanics deals with asymptotically stable base flows (the Lyapunov stability). With asymptotic stability, the gain obtained from the sustained forcing is bounded over the infinite-time horizon. However, the linear operators $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ for the present flows are unstable, as pointed out in figure 6. To address this matter for the present flow control effort, we discuss an extension to the standard resolvent analysis in the following section.

3.5 Finite-time horizon resolvent analysis

While the analysis of asymptotic stability requires an infinite-time horizon, the dynamical behaviour of a non-normal system within a finite-time horizon is also relevant. For an asymptotically stable system, the perturbation energy can undergo transient growth due to non-normality of the operator. Such dynamics is not described by the asymptotic behaviour of the operator, but can be characterized through an initial-value problem by specifying a finite-time horizon (Schmid & Brandt Reference Schmid and Brandt2014). Even if the system is characterized as unstable (unbounded) asymptotically, a bounded amplification can be found when a finite-time horizon is specified. For the present problem, some nonlinear dynamic processes, such as the shear-layer roll-up, the break-up of spanwise vortical structures and the vortex merging process, can all take place within a short time window. Therefore, we do not concern ourselves with the concept of asymptotic stability, but rather focus on the short-term dynamics by considering a finite-time horizon for the input–output analysis, following the approach proposed by Jovanović (Reference Jovanović2004).

Jovanović (Reference Jovanović2004) introduced an input–output analysis of an unstable system with an exponential discount. This analysis starts with the introduction of a temporal filter performed on both response and forcing such that $\check{\boldsymbol{q}}_{\unicode[STIX]{x1D6FD}}=\check{\boldsymbol{q}}\text{e}^{-t/t_{\unicode[STIX]{x1D6FD}}}$ and $\check{\boldsymbol{u}}_{\unicode[STIX]{x1D6FD}}=\check{\boldsymbol{u}}\text{e}^{-t/t_{\unicode[STIX]{x1D6FD}}}$ . The time constant $t_{\unicode[STIX]{x1D6FD}}>0$ is chosen such that the decay rate $\unicode[STIX]{x1D6FD}=1/t_{\unicode[STIX]{x1D6FD}}$ in the temporal filter $\text{e}^{-\unicode[STIX]{x1D6FD}t}$ overtakes the growth rate of the dominant unstable eigenvalue of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ . That is, $\unicode[STIX]{x1D6FD}>\max (\unicode[STIX]{x1D714}_{i})$ . The use of such temporal filter ensures that we examine the dominant transient growth that takes place over a time window characterized by $t_{\unicode[STIX]{x1D6FD}}$ . Therefore, the value of $t_{\unicode[STIX]{x1D6FD}}$ can be chosen according to physical interests. Upon substituting these growth-discounted modes of $\check{\boldsymbol{q}}_{\unicode[STIX]{x1D6FD}}$ and $\check{\boldsymbol{u}}_{\unicode[STIX]{x1D6FD}}$ into the Navier–Stokes equation (3.3), we have

(3.19) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D6FD}-\text{i}\unicode[STIX]{x1D714})\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}\hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6FD}}+\hat{\boldsymbol{u}}_{\unicode[STIX]{x1D6FD}}. & & \displaystyle\end{eqnarray}$$

Thus, we can express the discounted resolvent analysis as

(3.20) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}_{\unicode[STIX]{x1D6FD}}=[-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-(\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D644})]^{-1}\hat{\boldsymbol{u}}_{\unicode[STIX]{x1D6FD}}, & & \displaystyle\end{eqnarray}$$

with the discounted resolvent operator $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}$

(3.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}=[-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D644}-(\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D644})]^{-1}. & & \displaystyle\end{eqnarray}$$

This expression constructs the discounted resolvent operator $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}$ using the shifted linear operator $(\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D644})$ . The eigenvalues of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ are now shifted by $-\unicode[STIX]{x1D6FD}$ and all reside on the stable complex plane so that the standard resolvent analysis can be performed with $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}$ along the real axis of $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}$ . Note that $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}$ can also be expressed as $\unicode[STIX]{x1D643}_{\bar{\boldsymbol{q}},\unicode[STIX]{x1D6FD}}=[-\text{i}(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}]^{-1}$ , suggesting that an equivalent exercise can be performed by directly evaluating the pseudospectrum of $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ on a raised frequency axis of $(\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D6FD})$ . The traditional approach is recovered by setting $\unicode[STIX]{x1D6FD}=0$ (i.e.  $t_{\unicode[STIX]{x1D6FD}}\rightarrow \infty$ for infinite-time horizon), which was undertaken by Thomareis & Papadakis (Reference Thomareis and Papadakis2018) on the same airfoil at $Re_{L_{c}}=50\,000$ .

Figure 9. Finite-time horizon resolvent analysis with different choices of $t_{\unicode[STIX]{x1D6FD}}$ , considering $\unicode[STIX]{x1D6FC}=9^{\circ }$ mean flow with $k_{z}L_{c}=0$ . (a) Gain distribution over frequency in $St$ ; the leading resolvent (b) response and (c) forcing modes at $St=0.833$ . The lowest magnitude marked by contour lines is $1\,\%$ of the modal maximum. The streamwise extent of the modal structures shortens with decreasing $t_{\unicode[STIX]{x1D6FD}}$ .

We demonstrate this finite-time horizon resolvent analysis in figure 9 by showing representative results over various choices of $t_{\unicode[STIX]{x1D6FD}}$ . Here, we use the operator $\unicode[STIX]{x1D647}_{\bar{\boldsymbol{q}}}$ constructed with $k_{z}=0$ about the $\unicode[STIX]{x1D6FC}=9^{\circ }$ mean baseline flow and choose $t_{\unicode[STIX]{x1D6FD}}$ such that $t_{\unicode[STIX]{x1D6FD}}v_{\infty }/L_{c}=3$ , $5$ and $7$ . The results from these choices of $t_{\unicode[STIX]{x1D6FD}}$ will be compared with those from the infinite-time horizon analysis ( $t_{\unicode[STIX]{x1D6FD}}\rightarrow \infty$ ).

Let us analyse the gain distribution over frequency shown in figure 9(a). By decreasing $t_{\unicode[STIX]{x1D6FD}}$ from $7$ to $3$ , we observe that the gain over $St$ decreases with $t_{\unicode[STIX]{x1D6FD}}$ . The decrease in gain can be explained by the shorter time horizon over which the growth in perturbation energy is evaluated. It can also be understood as the decreasing pseudospectral level with increasing $\unicode[STIX]{x1D714}_{i}$ (moving away from the neutral stability axis) as we can observe in figure 7. The finite-time horizon analysis removes the spikes appearing in the gain distribution evaluated with the infinite-time horizon. The spikiness is attributed to the response of pseudospectral level to subdominant and spurious eigenmodes populating densely near the frequency ( $\unicode[STIX]{x1D714}_{r}$ ) axis, which can be seen in the spectrum in figure 6(a).

The leading response modes and forcing modes are respectively shown in figure 9(b,c) for the corresponding $t_{\unicode[STIX]{x1D6FD}}$ using $St=0.833$ as a representative case. From the response modes in figure 9(b), we observe that all choices of $t_{\unicode[STIX]{x1D6FD}}$ reveal the flow responses in the shear layer over the airfoil and in the wake. In figure 9(c), the forcing modes exhibit advective structures near the airfoil and its upstream. Note that the time scale, $t_{\unicode[STIX]{x1D6FD}}v_{\infty }/L_{c}$ , can also be interpreted as the advective length scale over the finite-time window. The streamwise coverage of the structures in both response and forcing modes is well characterized by each time constant $t_{\unicode[STIX]{x1D6FD}}$ used in the temporal filter.

The advective feature of the forcing mode motivates the use of local actuation, since the locally introduced perturbation that advects with the flow can leverage this feature as long as the forcing mode structures extend farther downstream of the actuator. Moreover, we observe that the forcing modes exhibit high level of fluctuation near the leading edge for all values of $t_{\unicode[STIX]{x1D6FD}}$ examined. This observation is in agreement with Thomareis & Papadakis (Reference Thomareis and Papadakis2018), who performed the resolvent analysis for infinite horizon. The forcing mode shape suggests that the amplification from the input–output process can be efficiently leveraged if actuation is introduced near the leading edge (Gómez & Blackburn Reference Gómez and Blackburn2017). Our choice of the actuator location ( $x_{a}/L_{c}=0.03$ ) is hence supported by the observation of the forcing mode structure.

We now move our attention to the choice of $t_{\unicode[STIX]{x1D6FD}}v_{\infty }/L_{c}=5$ in the rest of this work. However, we will also show that the conclusive assessment still stands for other choices in appendix B. We present the gain distribution over the $\unicode[STIX]{x1D714}$ – $k_{z}$ plane in figure 10. For each $\unicode[STIX]{x1D6FC}$ , the gain constructed from the second singular value $\unicode[STIX]{x1D70E}_{2}$ is also presented in comparison with that from $\unicode[STIX]{x1D70E}_{1}$ over the same frequency–wavenumber plane. The difference between $\unicode[STIX]{x1D70E}_{1}$ and $\unicode[STIX]{x1D70E}_{2}$ is typically greater than an order of magnitude. This gap between the leading and second singular value justifies the rank-1 assumption discussed in the previous section. Comparing the results from both angles of attack, we find that leading gain over the entire $\unicode[STIX]{x1D714}$ – $k_{z}$ plane is well scaled with the chord-based Strouhal number $St=\unicode[STIX]{x1D714}L_{c}/2\unicode[STIX]{x03C0}v_{\infty }$ and wavenumber $k_{z}L_{c}$ . The resemblance stems from the highly non-normal shear-layer modes residing near $St\approx 5$ for both angles of attack, which are observed from their pseudospectra in figure 7. Moreover, the gain exhibits a general decreasing trend with increasing $k_{z}L_{c}$ . This behaviour can be attributed to the attenuation of three-dimensional instability, which has been studied by Pierrehumbert & Widnall (Reference Pierrehumbert and Widnall1982) and Hwang, Kim & Choi (Reference Hwang, Kim and Choi2013) for free shear layer and wake, respectively.

Figure 10. Gain distribution over the $\unicode[STIX]{x1D714}$ – $k_{z}$ space for (a) $\unicode[STIX]{x1D6FC}=6^{\circ }$ and (b) $\unicode[STIX]{x1D6FC}=9^{\circ }$ , computed using $t_{\unicode[STIX]{x1D6FD}}v_{\infty }/L_{c}=5$ . Approximately $O(10)$ difference from the leading to second singular value is observed.

Figure 11. Streamwise velocity mode $\hat{v}_{x}$ , transverse velocity mode $\hat{v}_{y}$ and spanwise modal Reynolds stress $\hat{R}_{z}$ of representative $k_{z}$ – $St$ combinations for $\unicode[STIX]{x1D6FC}=9^{\circ }$ mean baseline flow. The response modes are obtained with $t_{\unicode[STIX]{x1D6FD}}v_{\infty }/L_{c}=5$ and are visualized by the contour lines of $\hat{q}/|\hat{q}|_{\infty }\in \pm [0.01,0.9]$ .

The structure of the response mode can also provide knowledge for identifying the appropriate actuation $k_{z}^{+}$ and $\unicode[STIX]{x1D714}^{+}$ that result in aerodynamically favourable control effects. Given a response mode $\hat{\boldsymbol{q}}\equiv [\hat{\unicode[STIX]{x1D70C}},\hat{v}_{x},\hat{v}_{y},\hat{v}_{z},\hat{T}]^{\text{T}}$ at specified $k_{z}$ and $\unicode[STIX]{x1D714}$ , we also evaluate the associated streamwise, transverse and spanwise Reynolds stress respectively by

(3.22a-c ) $$\begin{eqnarray}\displaystyle \hat{R}_{x}(k_{z},\unicode[STIX]{x1D714})=\text{Re}(\hat{v}_{y}^{\ast }\hat{v}_{z}),\quad \hat{R}_{y}(k_{z},\unicode[STIX]{x1D714})=\text{Re}(\hat{v}_{z}^{\ast }\hat{v}_{x}),\quad \hat{R}_{z}(k_{z},\unicode[STIX]{x1D714})=\text{Re}(\hat{v}_{x}^{\ast }\hat{v}_{y}),\quad & & \displaystyle\end{eqnarray}$$

where $\text{Re}(\cdot )$ denotes the real component of the argument. We visualize the response modes in figure 11 in their streamwise velocity $\hat{v}_{x}$ , transverse velocity $\hat{v}_{y}$ and the associated spanwise Reynolds stress $\hat{R}_{z}$ with the representative $k_{z}L_{c}$ – $St$ combinations for the mean baseline flow at $\unicode[STIX]{x1D6FC}=9^{\circ }$ . For modes of $St=1.5$ and $2.5$ , response structure develops from the shear layer above the suction surface and extends farther into the wake. Particularly for $(k_{z}L_{c},St)=(0,1.5)$ , we observe an extended wake structure in the velocity modes as well as the resolvent Reynolds stress. The Reynolds stress exhibits a pattern of von Kármán vortex shedding with negative correlation developing in the shear layer above the airfoil and positive correlation extending from the trailing edge over the bottom. By increasing either $St$ or $k_{z}L_{c}$ , the streamwise extent of the modal structure reduces to the shear layer. Further increase of frequency moves the response structure towards the leading edge where the shear layer remains thin and is capable of supporting small-scale structures from high-frequency perturbations. In § 5, we will further explore the resolvent Reynolds stress and show that it provides insight into momentum mixing that is critical for suppressing stall. We will discuss a metric that incorporates the gain and the spatial integration of the resolvent Reynolds stresses to provide a quantitative guidance to separation control.

We have performed resolvent analysis for the mean baseline flows of $\unicode[STIX]{x1D6FC}=6^{\circ }$ and $9^{\circ }$ and discussed an extension to the standard approach for the two unstable linear operators. From the gain distribution over frequency and wavenumber, we have seen the shear-layer-dominated feature for the baseline flows at both angles of attack. In the next section, we discuss the results of controlled flows and investigate important flow physics that achieves suppression of separation. The enhancement in aerodynamic performances will later be compared to the prediction of resolvent analysis.

The nonlinear physics beyond resolvent analysis

With resolvent analysis being a linear technique for the present nonlinear fluid-flow problem, we also observed some limitations of the interpretation in the aerodynamic performance and the prediction of $M(k_{z}^{+},\unicode[STIX]{x1D714}^{+})$ . Below, we comment on these limitations and identify the associated nonlinear physics.

For controlled cases with $k_{z}^{+}L_{c}=0$ , drag increases over $4\lesssim St^{+}\lesssim 10$ for both angles of attack. Such increase in drag is not captured by the value of $M$ . As discussed in the previous section, this drag increase is due to the vortex merging process that causes trailing-edge separation. Therefore, the difference between the controlled flow results and resolvent analysis can be attributed to the nonlinear nature of the merging process that transfers energy from a fundamental frequency to its subharmonics. With the energy transfer across frequency space, this nonlinear process is not captured by the linear resolvent analysis that deals with a harmonic input–output process.

Another nonlinear process that leads to differences between the LES findings and the results of resolvent analysis is the laminar–turbulent transition following the break-up of spanwise vortices. In the previous section, the transition process has been shown to be a mechanism responsible for the suppression of separation, in addition to the shear-layer excitation. Such a mechanism is particularly important for suppressing stall in the control cases with high frequency near $St^{+}\approx 10$ and $k_{z}^{+}L_{c}>0$ , leading to peak drag reduction at $St^{+}\approx 12$ for $\unicode[STIX]{x1D6FC}=9^{\circ }$ and comparable level of force enhancement across three choices of three-dimensional actuation profiles ( $k_{z}^{+}L_{c}>0$ ). However, the level of $M(k_{z},\unicode[STIX]{x1D714})$ evaluated from resolvent analysis suggests degraded mixing for $k_{z}L_{c}>0$ and high frequencies. Therefore, while the aerodynamic forces benefit from the laminar–turbulent transition, this nonlinear process is also beyond the capability of resolvent analysis to predict the force enhancement through transition by using the quantitative level of $M(k_{z},\unicode[STIX]{x1D714})$ .

Resolvent analysis as a guiding tool for separation control

We have presented a design guideline that leverages the knowledge obtained from resolvent analyses performed on mean baseline flows for suppressing stall. We evaluate the level of mode-based mixing by combining the knowledge of amplification, modal structure and a shear-layer window over the airfoil, providing a scalar function over the frequency–wavenumber space. In spite of slight deviations due to the nonlinear physics beyond the present linear modal, the control effect well correlates with lift enhancement and drag reduction for open-loop controlled flows. Such a guideline provides quantitative assessment towards selecting actuation frequency and wavenumber for effective unsteady separation control.