2.1. Governing equations

In this study, we consider the non-dimensional compressible N–S equations with external forces as the governing equations. The equation is non-dimensionalized using the velocity, density and temperature of free stream flow and chord length of NACA0012 airfoil. The compressible N–S equations consist of the equations of mass, momentum and pressure as

(2.1) \begin{align} \frac {\partial \rho }{\partial t} &= - \frac {\partial }{\partial x_{\!j}} (\rho u_{\!j}), \\[-12pt]\nonumber \end{align}

(2.2) \begin{align} \frac {\partial \rho u_i}{\partial t} &= - \frac {\partial }{\partial x_{\!j}} (\rho u_i u_{\!j}) - \frac {\partial p}{\partial x_i} + \frac {\partial \tau _{\textit{ij}}}{\partial x_{\!j}} + \rho h_i, \\[-12pt]\nonumber \end{align}

(2.3) \begin{align} \frac {\partial p}{\partial t} &= - \frac {\partial }{\partial x_{\!j}} (p u_{\!j}) - (\gamma -1) \left \{ p \frac {\partial u_{\!j}}{\partial x_{\!j}} - \tau _{\textit{ij}} \frac {\partial u_i}{\partial x_{\!j}} - \frac {\partial \mathcal{Q}_{\!j}}{\partial x_{\!j}} \right \}\!, \end{align}

where $\rho , u_x, u_y, p, \gamma$ represent the density, velocity in $x$ -direction, velocity in $y$ -direction, pressure and specific heat ratio, respectively. Here, $h_x, h_y$ are the external forces, and the viscous stress tensor $\tau _{\textit{ij}}$ and heat flux $\mathcal{Q}_i$ are expressed as

(2.4) \begin{equation} \tau _{\textit{ij}} = \frac {\mu }{\textit{Re}} \left ( \frac {\partial u_i}{\partial x_{\!j}} + \frac {\partial u_{\!j}}{\partial x_i} - \frac {2}{3} \delta _{\textit{ij}} \frac {\partial u_k}{\partial x_k} \right )\!, \end{equation}

and

(2.5) \begin{equation} \mathcal{Q}_i = \frac {\kappa }{(\gamma -1) M^2 \textit{Re} \textit{Pr}} \frac {\partial \varTheta }{\partial x_i}, \end{equation}

where $M, \textit{Re}, \textit{Pr}$ are the Mach number, Reynolds number and Prandtl number, respectively. In this study, we fixed $\textit{Pr} = 0.7$ and $\gamma = 1.4$ for the air flow around an airfoil. In the above equation, $\varTheta$ represents the temperature, and the state equation is given by

(2.6) \begin{equation} p = \frac {\rho \varTheta }{\gamma M^2}. \end{equation}

For the coefficients of viscosity $\mu$ and thermal conductivity $\kappa$ . For the viscosity, we applied Sutherland law of viscosity as

(2.7) \begin{equation} \mu ^* = 1.458\times 10^{-5} \frac {\varTheta ^{* 3/2}}{\varTheta ^* + 110.4},\quad\varTheta ^* \geqslant 110.4 \,\mathrm{K} \end{equation}

where $\varTheta ^* [\mathrm{K}]$ is the dimensional temperature. We assumed the free stream temperature was $273.15 \mathrm{K}$ for all the computation conditions.

2.2. Problem settings

2.2.1. Definition of base flow and external forcing

In this section, we define the problem setting of our optimal forcing analysis. We consider a solution trajectory of the governing equation of compressible N–S equations as the base flow, $\boldsymbol{Q} (t, \boldsymbol{x})$ , over a range of $0 \leqslant t \leqslant t_{\!f}$ , where $t_{\!f}$ is a given evaluation time. Here, $\boldsymbol{Q}$ is the state vector of the flow field, specifically $\boldsymbol{Q} = (\bar {\rho }, \bar {\rho } \bar {u}_x, \bar {\rho } \bar {u}_y, \bar {p})^T$ where $\bar {\rho }, \bar {u}_x. \bar {u}_y, \bar {p}$ are the density, $x$ -velocity, $y$ -velocity and pressure of the base flow, respectively. We assume that the base flow $\boldsymbol{Q} (t, \boldsymbol{x})$ is dependent on both time and space, as the assumption of a steady base flow is not required in the subsequent discussion.

We consider the application of external forcing, $\boldsymbol{h} (t, \boldsymbol{x})$ , to the flow field, assuming that the flow changes from $\boldsymbol{Q} (t, \boldsymbol{x})$ to $\boldsymbol{q} (t, \boldsymbol{x})$ due to the forcing. We define the disturbance flow field as $\boldsymbol{q}' (t, \boldsymbol{x}) = \boldsymbol{q} (t, \boldsymbol{x}) - \boldsymbol{Q} (t, \boldsymbol{x})$ . Here, $\boldsymbol{q} (t, \boldsymbol{x})$ represents the flow field perturbed by the external forces. For these external forces, we consider $\boldsymbol{h} = (0, h_x, h_y, 0)^T$ where $h_x$ and $h_y$ are the body forces in the $x$ - and $y$ -directions, respectively. Based on the above descriptions, the perturbed flow field $\boldsymbol{q} (t, \boldsymbol{x})$ satisfies the compressible N–S equations with external forces as

(2.8) \begin{equation} \partial _t \boldsymbol{q} (t, \boldsymbol{x}) = \mathcal{F} [\boldsymbol{q} (t, \boldsymbol{x})] + \rho \boldsymbol{h} (t, \boldsymbol{x}), \end{equation}

where $\mathcal{F}$ formally represents the right-hand side of the compressible N–S equations. Furthermore, we assume the separation of variables in the external forces as

(2.9) \begin{equation} \boldsymbol{h} (t, \boldsymbol{x}) = f(t) \boldsymbol{g} (\boldsymbol{x}), \end{equation}

where $f(t)$ is a scalar function representing the temporal distribution and $\boldsymbol{g}(\boldsymbol{x})$ is a vector function representing the spatial distribution. In flow control applications, this assumption is considered natural because the separation of variables into spatial and temporal distributions corresponds to fixing the location of control devices while modulating their actuation over time.

2.2.2. Variational formulation

In this setting, we define the problem as determining the optimal forcing field that changes the base flow under appropriate constraint conditions. We introduce an evaluation function, $\mathcal{G}[t_{\!f}, \boldsymbol{q}, \boldsymbol{Q}]$ to quantify these flow field changes. For selecting the evaluation function, several approaches have been proposed. These include the incompressible disturbance energy (Huang & Hack Reference Huang and Hack2020), compressible extensions of disturbance energy (Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996) and the time-average of total entropy production (Monokrousos et al. Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011; Duguet et al. Reference Duguet, Monokrousos, Brandt and Henningson2013; Taniguchi et al. Reference Taniguchi, Ohmichi and Suzuki2024a ). Specifically, Taniguchi et al. (Reference Taniguchi, Ohmichi and Suzuki2024a ) discussed the differences in nonlinear optimal perturbation fields for subsonic flow around a NACA0012 airfoil when using the first and third metrics. They showed that the latter exhibits high sensitivity to changes related to wall friction. We define the evaluation function as the time-averaged entropy production, defined as

(2.10) \begin{equation} \mathcal{G} [t_{\!f}, \boldsymbol{q} (t, \boldsymbol{x}); \boldsymbol{Q} (t, \boldsymbol{x})] = \int _0^{t_{\!f}} P[\boldsymbol{q} (t, \boldsymbol{x})] {\rm d}t, \end{equation}

and $P[\boldsymbol{q} (t, \boldsymbol{x})]$ is the integral of the entropy production of the instantaneous flow field over the entire computational domain, $\varOmega$ . This term is given by

(2.11) \begin{equation} P[\boldsymbol{q}(t, \boldsymbol{x})] = \int _{\varOmega } \frac {\kappa }{\varTheta ^2} |\boldsymbol{\nabla }\varTheta |^2 + \frac {\mu }{2 \varTheta } S:S - \frac {2}{3} \frac {\mu }{\varTheta } \left ( \frac {\partial u_{\!j}}{\partial x_{\!j}} \right )^2 {\rm d}V, \end{equation}

with

(2.12) \begin{equation} S_{\textit{ij}} = \frac {\partial u_i}{\partial x_{\!j}} + \frac {\partial u_{\!j}}{\partial x_i}, \end{equation}

where $\mu$ and $\kappa$ are the coefficients of viscosity and thermal conductivity, respectively. In this study, we employed the evaluation function in (2.10) because it provides a globally defined measure that is closely correlated with the recovery of wall friction, thereby serving as a promising indicator of flow reattachment.

As a constraint condition, we consider a constant energy input for the external forcing. To determine its spatial distribution, we extended the approach of nonlinear optimal perturbation analysis by Pringle & Kerswell (Reference Pringle and Kerswell2010), Cherubini et al. (Reference Cherubini, De Palma, Robinet and Bottaro2011) and Huang & Hack (Reference Huang and Hack2020) to external forces. We define the constraints for optimal forcing as

(2.13) \begin{equation} E[\boldsymbol{g}; \boldsymbol{Q}] := \int _{\varOmega } \frac {1}{2} ( g_x^2 + g_y^2 ) {\rm d}V = \delta _0^2, \end{equation}

where $\delta _0$ is the given magnitude of the external forces. Without loss of generality, we can assume the $L^2$ norm of the temporal distribution is unity, given by

(2.14) \begin{equation} || f(t) ||_2^2 = 1. \end{equation}

In this study, we use the variational formulation of optimal forcing analysis based on the nonlinear optimal perturbation method (Pringle & Kerswell Reference Pringle and Kerswell2010; Cherubini & De Palma Reference Cherubini and De Palma2013). In this formulation, we introduce the Lagrangian functional

(2.15) \begin{align} \begin{split} \mathcal{L} &= \mathcal{G} [t_{\!f}; \boldsymbol{q} (t, \boldsymbol{x}); \boldsymbol{Q} (t, \boldsymbol{x})] - \chi (||f (t)||_2^2 - 1) - \mu (E[\boldsymbol{g} (\boldsymbol{x}); \boldsymbol{Q} (t_0, \boldsymbol{x})] - \delta _0^2) \\ &- \int _{t_0}^{t_{\!f}} \boldsymbol{\lambda }^* (t, \boldsymbol{x}) \left \{ \partial _t \boldsymbol{q} (t, \boldsymbol{x}) - \mathcal{F} [\boldsymbol{q} (t, \boldsymbol{x})] - \rho f(t) \boldsymbol{g} (\boldsymbol{x}) \right \} {\rm d}t \end{split} \end{align}

where $\mu , \chi $ and $\boldsymbol{\lambda } = \boldsymbol{\lambda } (t, \boldsymbol{x})$ are the Lagrange multipliers.

2.2.3. Spatial distribution of forcing

In this study, we use a DBD plasma actuator as the flow control device. Physically, a DBD plasma actuator consists of two thin electrodes mounted on the surface and a high-voltage power supply. Flow is induced by the plasma discharge between the two electrodes (Post & Corke Reference Post and Corke2004; Jukes et al. Reference Jukes, Choi, Johnson and Scott2006). Given this set-up, it is reasonable to assume that the installation location of the plasma actuator can be freely determined. The wall-normal distribution, determined by the plasma’s extension from the electrodes, typically exhibits a Gaussian-like profile (see figure 10 in Post & Corke (Reference Post and Corke2004) for PIV measurements of the actual flow induced by a DBD actuator). Thus, modelling the DBD plasma actuator as a body force with separated wall-tangential and wall-normal distributions is a reasonable approach. However, a potential drawback of this formulation is that it permits the external forcing to be distributed over the entire airfoil surface, creating an unphysical scenario for a localized actuator. Consequently, we will revisit this issue in § 4, specifically discussing the applicability of extracting and applying the forcing only in locally concentrated regions to ensure a physically realizable control strategy.

Based on the above discussion, we further assume a separation of variables in the spatial distribution of forcing, $\boldsymbol{g} (\boldsymbol{x})$ . For the wall-tangential direction $\xi$ and wall-normal direction $\eta$ , we assume that

(2.16) \begin{equation} \boldsymbol{g} (\boldsymbol{x}) = g_{\eta } (\eta ) \boldsymbol{g}_{\xi } (\boldsymbol{x}), \end{equation}

where $\boldsymbol{g}_{\xi } (\boldsymbol{x})$ is the $\xi$ -directional vector function, and $g_{\eta } (\eta )$ is the $\eta$ -directional scalar function. Under this assumption, we consider $\boldsymbol{g}_{\xi } (\boldsymbol{x})$ as the target function for optimization, while $g_{\eta } (\eta )$ is a function determined by the source of external forcing. In the optimization, we do not impose any further assumptions on the wall-tangential distribution of the optimal forcing to ensure the flexibility of the optimization results.

For wall-normal distribution, we model the external force source as a plasma actuator (Post & Corke Reference Post and Corke2004), where the assumption in (2.16) is a much simpler model compared with the previous studies (Asada et al. Reference Asada, Nonomura, Aono, Sato, Okada and Fujii2015). This simplified assumption is incorporated to facilitate the optimization of the forcing, whereas rather detailed analysis would be required for practical application in control systems. For $g_{\boldsymbol{x}} (\eta )$ , we adopt a Gaussian distribution, with a modification to enforce the boundary condition $\boldsymbol{h} (\boldsymbol{x}) = 0$ at $\eta =0$ as

(2.17) \begin{equation} g_{\eta } (\eta ) = \frac {\eta }{\sqrt {2 \pi \sigma ^2}} \exp \left ( -\frac {(\eta -\mu )^2}{2 \sigma ^2} \right )\!, \end{equation}

where $\sigma$ and $\mu$ are constant free parameters that should be determined to suit the targeted control device. In this study, we used the values $\sigma = 0.001$ and $\mu = 0.005$ to align the forcing distribution with the height of the laminar separation bubble in the base flow.

2.2.4. Euler–Lagrange equation

To solve the maximization problem of the Lagrangian in (2.15), we derive the Euler–Lagrange equations. Taking the first variation of (2.15) yields

(2.18) \begin{align} \frac {\partial \mathcal{L}}{\partial \boldsymbol{q} (t_{\!f})} = 0 \Rightarrow & \boldsymbol{\lambda } (t_{\!f}) = 0, \\[-12pt]\nonumber \end{align}

(2.19) \begin{align} \frac {\partial \mathcal{L}}{\partial \boldsymbol{g}_{\xi }} = 0 \Rightarrow & - \mu\! \int g_{\eta }^2{\rm d}\eta \boldsymbol{g}_{\xi } + \int _0^{t_{\!f}} f(t) \int \rho \boldsymbol{\lambda } g_{\eta } {\rm d}\eta {\rm d}t = 0, \\[-12pt]\nonumber \end{align}

(2.20) \begin{align} \frac {\partial \mathcal{L}}{\partial \boldsymbol{q} (t)} = 0 \Rightarrow & \partial _t \boldsymbol{\lambda } (t) + \left ( \cfrac {\partial \left ( \mathcal{F}[\boldsymbol{q} (t)] + \rho (t) f(t) \boldsymbol{g} (t) \right )}{\partial \boldsymbol{q} (t)} \right )^{\dagger } \boldsymbol{\lambda } (t) = \left ( \frac {\partial P[\boldsymbol{q} (t)]}{\partial \boldsymbol{q} (t)} \right )^{\dagger }\!, \\[-12pt]\nonumber \end{align}

(2.21) \begin{align} \frac {\partial \mathcal{L}}{\partial f (t)} = 0 \Rightarrow & - 2 \chi f(t) + \langle \boldsymbol{\lambda } (t), \rho (t) \boldsymbol{g} (t) \rangle = 0. \end{align}

Here, the dagger $(\dagger )$ denotes the adjoint of the linear operator. For the derivation of the compressible adjoint equations, we followed the formulation by Ohmichi & Yamada (Reference Ohmichi and Yamada2021). In the above equations, we omitted the explicit spatial variable $\boldsymbol{x}$ from the arguments of functions for brevity. By taking the first variation of $\mathcal{L}$ with respect to the Lagrange multipliers, we recover the constraint conditions in (2.8), (2.13) and (2.14).

To update the spatial and temporal distribution of the optimal forcing, the standard steepest ascent method (Pringle & Kerswell Reference Pringle and Kerswell2010; Cherubini & De Palma Reference Cherubini and De Palma2013) can be applied. Algorithm 1 outlines the optimization process for the NLOF analysis. In this approach, the step size is independently determined by adopting the approach by Pringle et al. (Reference Pringle, Willis and Kerswell2012). The remaining Lagrange adjoint multipliers, $\mu$ and $\chi$ , are determined from the constraint conditions. The convergence tolerances for the residuals, $\epsilon _{\textit{spa}}, \epsilon _{\textit{tem}}$ , will be discussed in § 3.1.

2.3. Numerical methods

2.3.1. Flow condition

In this study, we investigated two-dimensional flow around a NACA0012 airfoil in the high-subsonic, low-Reynolds-number regime. Following the study by Bouhadji & Braza (Reference Bouhadji and Braza2003), we selected the conditions of $M=0.75$ and $\textit{Re}=2.0\times 10^4$ . The angle of attack ( $ \textit{AoA} $ ) was set to $2^{\circ }$ , a condition where low-frequency oscillation is not observed.

Algorithm 1 Optimization process of NLOF analysis

2.3.2. Computational domain and grid

We present the computational domain and grid used for the NLOF analysis. We set the origin of the coordinates at the midpoint of the airfoil’s chord line, and the chord length $c$ was taken as the reference length. The computational domain was a circular region with a radius of $60 c$ , including a sponge layer at its outer edge, as will be detailed in § 2.3.3.

Figure 1 shows the numerical grid near the airfoil. We considered a rounded trailing edge with a radius $0.0015 c$ , following Wolf & Lele (Reference Wolf and Lele2012). A body-fitted O-type computational grid was employed. This topology allows us to apply a simple method to define the normal vector, avoiding the singular point at the trailing edge. The grid was clustered near the airfoil surface to resolve the separated flow. The minimal grid spacing was $1.8\times 10^{-5}$ in the radial direction and $1.7\times 10^{-4}$ in the angular direction. The number of the grid points was $256$ in the radial direction and $1024$ in the angular direction.

Figure 1. Computational grid in every five points near the airfoil at $\textit{Ao}A=2^{\circ }$ .

2.3.4. Numerical validation

This section presents the numerical validation of our computational methods. Figure 2 compares the time-averaged lift coefficients at $M=0.2$ with previous studies, specifically the experimental values from Ohtake et al. (Reference Ohtake, Nakae and Motohashi2007) and two-dimensional laminar calculation results from Lee et al. (Reference Lee, Nonomura, Oyama and Fujii2017). At $\textit{Re}=10^4$ , our calculated values agree well with both numerical and experimental data at low angles of attack, although they deviate from experimental values at high angles of attack. For the higher Reynolds number case, shown in figure 2(b), our current calculations also show good agreement with the numerical values reported by Lee et al. (Reference Lee, Nonomura, Oyama and Fujii2017), with the exception of a larger deviation at an $\textit{AoA}=3^{\circ }$ . In comparison, Kim, Chang & Chung (Reference Kim, Chang and Chung2011) reported $C_l = 0.40$ at $\textit{Re}=2.3\times 10^4$ and $C_l = 0.49$ at $\textit{Re}=3.3\times 10^4$ for a NACA0012 airfoil in a wind tunnel experiment at this $ \textit{AoA} $ . Our present calculated values fall within this expected range. Thus, we conclude that the current numerical scheme accurately captures the dependency of the lift coefficients on the $ \textit{AoA} $ within the flow conditions considered in this study.

Second, we examine the mesh convergence of the calculation results. Table 1 presents the mesh dependency of the lift and drag coefficients. From this table, we conclude that the current numerical solutions is sufficiently converged with respect to the computational grid.

Table 1. Mesh dependency of time-averaged lift and drag coefficients at $M=0.75, \textit{Re}=2\times 10^4, \textit{AoA}=2^{\circ }$ .

Figure 2. Variation of time-averaged lift coefficient with $ \textit{AoA} $ at $M=0.2$ : (a) $\textit{Re}=10^4$ ; (b) $\textit{Re}=3 \times 10^4$ .

2.4. Data analysis methods

The application of optimal forcing induces spatiotemporal transient changes in the flow field. To analyse the unsteady behaviour of the flow changes, we adopted the variational mode decomposition-based non-stationary coherent structure (VMD-NCS) method proposed by Ohmichi (Reference Ohmichi2024). The effectiveness of VMD-NCS has been previously demonstrated in revealing non-stationary fluid phenomena, for example, in the analysis of buffet cells (Ohmichi & Sugioka Reference Ohmichi and Sugioka2026 ).

Let $\boldsymbol{q} (t, \boldsymbol{x})$ represent flow data sampled at discrete time steps $t=t_1, t_2, \ldots , t_N$ . First, we apply proper orthogonal decomposition (POD). We assume that the data can be projected onto an $M$ -dimensional subspace spanned by POD basis vectors $\tilde {\boldsymbol{q}}_{\!j} (\boldsymbol{x})$ for $j=1, 2, \ldots , M$ . The flow field can then be represented as

(2.22) \begin{equation} \boldsymbol{q} (t, \boldsymbol{x}) = \sum _{j=1}^M a_{\!j}(t) \tilde {\boldsymbol{q}}_{\!j} (\boldsymbol{x}) \end{equation}

where $a_{\!j} (t) = \langle \boldsymbol{q} (t, \boldsymbol{x}), \tilde {\boldsymbol{q}} \rangle$ , and $\langle , \rangle$ denotes the inner product. Following Hanifi et al. (Reference Hanifi, Schmid and Henningson1996), we define the inner product for compressible flow fields as

(2.23) \begin{equation} \langle \boldsymbol{q}_1, \boldsymbol{q}_2 \rangle = \frac {1}{2} \int _\varOmega \frac {\bar {\varTheta } \rho _1 \rho _2}{\gamma \bar {\rho } M^2} + \bar {\rho } ({u_x}_1 {u_x}_2 + {u_y}_1 {u_y}_2) + \frac {\bar {\rho } \varTheta _1 \varTheta _2}{\gamma \tilde {\gamma } \bar {\varTheta } M^2} {\rm d}V, \end{equation}

where $\boldsymbol{q}_1, \boldsymbol{q}_2$ are arbitrary state vectors of the flow field. Consequently, the flow data can be entirely reconstructed as a point $\boldsymbol{a} (t) = (a_1(t), a_2(t), \ldots , a_M(t))^T$ in the POD basis coordinate system.

The core idea of VMD-NCS (Ohmichi Reference Ohmichi2024) is to apply multivariate VMD (Rehman & Aftab Reference Rehman and Aftab2019) to this $M$ -dimensional time series data, extracting modes termed intrinsic coherent structures (ICS) after mapping them back to the physical space. We briefly introduce this concept, starting with VMD analysis, which is originally proposed in the context of signal analysis by Dragomiretskiy & Zosso (Reference Dragomiretskiy and Zosso2013) as an alternative to empirical mode decomposition. In this analysis, the signal data is converted into an analytic signal via the Hilbert transform and decomposed into a limited number of VMD modes which possess local frequency bands around the centre frequency. The extension of VMD into multivariate signals was given by Rehman & Aftab (Reference Rehman and Aftab2019), based on the assumption that all the variables of the data possess common oscillation properties. For an application to fluid dynamics, the validity of this assumption was demonstrated by Ohmichi (Reference Ohmichi2024).

The mathematical formulation of multivariate VMD (MVMD) is thereby given as follows. The decomposition of multivariate signal $\boldsymbol{a}(t) = (a_1 (t), a_2 (t), \ldots , a_M (t))^T$ into $K$ band-limited MVMD modes $\{(v^{k,1} (t), v^{k,2}, \ldots , v^{k,M})^T\}_k$ is formulated as an inverse problem to recover the signal as a linear combination of MVMD modes. Introducing Tikhonov regularization, this is achieved by solving the optimization problem of the functional (equation (14) in Rehman & Aftab (Reference Rehman and Aftab2019))

(2.24) \begin{equation} \sum _{m=1}^M || a_m (t)-\sum _{k=1}^K v^{k, m} (t) ||_2^2 + \alpha \sum _{m=1}^M \sum _{k=1}^K || \partial _t [v_+^{k,m} (t) \mathrm{e}^{-\mathrm{i} \omega _k t}] ||_2^2, \end{equation}

where $v_+^{k, m}$ is the analytical signal of $v^{k, m}$ and $\omega _k$ is the centre frequency of the $k$ th MVMD mode. Note that this regularization term enforces a narrow-band frequency distribution around the centre frequency $\omega _k$ for the $k$ th MVMD mode, while their linear combination recovers the original signal.

The number of ICS modes, $K$ , and the regularization parameter $\alpha$ are user-defined. In the actual implementation, alternating direction method of multipliers is employed; for details, refer to Rehman & Aftab (Reference Rehman and Aftab2019) and Ohmichi (Reference Ohmichi2024). The $k$ th ICS, $\boldsymbol{\phi }_k$ is then reconstructed using the POD basis as

(2.25) \begin{equation} \boldsymbol{\phi }_k(t, \boldsymbol{x}) = \sum _{m=1}^M v^{k, m} (t) \boldsymbol{\tilde {q}}_m (\boldsymbol{x}), \end{equation}

where $\phi _k$ represents the extracted $k$ th ICS in the physical space. This ICS represents the spatiotemporal structure within a narrow temporal frequency band around the centre frequency $\omega _k$ .

2.5. Characteristic of the base flow field

This section presents the characteristics of the base flow field at $M=0.75, \textit{Re}=20\,000$ , with an $ \textit{AoA} $ of $2^{\circ }$ . Figure 3 shows the power spectral density (PSD) of the lift coefficient of the base flow field, including a comparison with the $\textit{AoA}=0^{\circ }$ case. In this analysis, we used the Hann window function to suppress spectral leakage. The (dimensionless) sampling period was $1.07\times 10^{-4}$ at $\textit{AoA}=0^{\circ }$ and $1.12 \times 10^{-4}$ at $\textit{AoA}=2^{\circ }$ . Data were obtained after sufficient time evolution from the initial condition.

Figure 3. The PSD of lift coefficients of the base flow at $\textit{AoA}=0^{\circ }$ and $2^{\circ }$ .

In this figure, we observe that the flow field at $\textit{AoA}=0^{\circ }$ exhibits a primary peak frequency at $St=0.074$ and a secondary peak frequency at $St=2.37$ , where the first peak frequency shows a reasonable value compared with $0.06$ at $\textit{Re}=10\,000$ in Bouhadji & Braza (Reference Bouhadji and Braza2003) and $0.07$ at $\textit{Re}=50\,000$ in three-dimensional calculations by Moise et al. (Reference Moise, Zauner and Sandham2024). Each peak frequency is accompanied by its integer multiples (harmonics). In contrast, the flow field at $\textit{AoA}=2^{\circ }$ shows its primary peak frequency at $St=2.24$ . Although much weaker compared with the dominant frequency, a peak frequency at $St\approx 0.089$ is also present as a counterpart of the low-frequency oscillation at $\textit{AoA}=0^{\circ }$ . Therefore, the flow field at $\textit{AoA}=0^{\circ }$ is characterized by two dominant peak frequencies, indicating a quasiperiodic behaviour where the low-frequency oscillation is a dominant phenomenon. In contrast, the flow field at $\textit{AoA}=2^{\circ }$ exhibits a more periodic oscillation, with a negligible contribution from low-frequency oscillation.

Figure 4 presents the instantaneous flow field and time-history of the lift coefficient at $\textit{AoA}=2^{\circ }$ . This figure shows massive separation over the entire suction surface of the airfoil, accompanied by the shedding of von Kármán vortices in the wake. These features of the flow field show a qualitative agreement with the low-Reynolds number flow around a thick airfoil in smoke-wire visualizations reported by Kim et al. (Reference Kim, Chang and Chung2011). This vortex shedding is dominant, and a corresponding periodic variation in the lift coefficient is observed in figure 4(b). In this figure, the mean lift coefficient is negative, which is consistent with the negative lift slope characteristics in low Reynolds number flows around airfoils reported by Ohtake et al. (Reference Ohtake, Nakae and Motohashi2007), Yonemoto et al. (Reference Yonemoto, Takato, Ochi and Fujie2008) and Winslow et al. (Reference Winslow, Otsuka, Govindarajan and Chopra2018).

Figure 4. Instantaneous flow field at $\phi = 0$ and the time-history of $C_l$ at $\textit{AoA}=2^{\circ }$ . The black curves in the instantaneous flow field represent the separation line zero-streamwise velocity isocontours $(u_x=0)$ .

3.1. Optimization of evaluation function

We first examine the convergence of the solution obtained from the NLOF analysis solution. To evaluate the increase in entropy production, we define the percentage ratio of the evaluation function, $\mathcal{F} (t_{\!f})$ , as

(3.1) \begin{equation} \mathcal{F}(t_{\!f}) \times 10^{-2} = \frac {\mathcal{G}[t_{\!f}, \boldsymbol{q}(t, \boldsymbol{x})]}{\mathcal{G}[t_{\!f}, \boldsymbol{Q}(t, \boldsymbol{x})]} - 1. \end{equation}

Figure 5 shows a typical history of the optimization process. In this figure, the variation in the optimization is defined as $|\mathcal{F}^{(n+1)} (t_{\!f}) - \mathcal{F}^{(n)} (t_{\!f})|$ during the update from step $n$ to $n+1$ . We also present the residuals of the optimization process regarding the temporal and spatial distributions. As shown in this figure, the evaluation function was gradually updated and increased towards its maximum value. In the early optimization process, the evaluation function values improved rapidly towards their maximum values. When the evaluation function value increased and its fluctuation fell below 0.01, the convergence rate of the optimization calculation significantly decreased. This reduction in the variation of the evaluation function was proportional to the residuals of the optimization process, whose values were also influenced by the selection of the evaluation time. We set $\epsilon _{\textit{spa}} = 5 \times 10^{-3}$ and $\epsilon _{\textit{tem}} = 10^{-5}$ as the convergence tolerances for the residuals, ensuring that the variation of the evaluation function is less than $0.01$ , a level at which no significant changes of the optimal forces were observed.

The convergence properties of the optimization process were also affected by the selection of the step size. Selecting different step sizes in the temporal and spatial directions suppresses numerical oscillations during the optimization calculations. However, we confirmed that if either step size rapidly approaches zero, the likelihood of becoming trapped in a local optimum increases. As demonstrated in figure 5(b), the residuals gradually decreased, but completely eliminating numerical oscillations remained difficult. The primary factor hindering convergence is phase shifts in the time domain, which, coupled with the update of the spatial distribution, impede the optimization process. However, since these temporal shifts do not affect the changes in the separation region examined in this study, further improving of convergence efficiency is left as a future challenge.

Figure 5. History of the evaluation function in the optimization process at $t_{\!f} = 2.0, \delta =0.01, \phi =-\pi$ . (a) Evaluation function value and its variation during the optimization. (b) Residuals for spatial and temporal distributions (semi-log plot).

3.2. Characteristic of optimal forcing

In this section, we investigate the spatiotemporal characteristics of the optimal forcing. For the following discussion, we introduce the growth rate of the disturbance field, defined as

(3.2) \begin{equation} G(t_{\!f}) = \frac {1}{\delta _0^2} \int _{\varOmega } \left . \frac {1}{2 \bar {\rho }} m_i' m_i' \right |_{t=t_{\!f}} {\rm d}V. \end{equation}

Figure 6 shows the spatial and temporal distributions of the optimal forcing for different amplitudes. The growth rate of the disturbance field is found to be inversely related to the amplitude of the external forcing: $G=2.4 \times 10^{4} (\mathcal{F}=1.78)$ at $\delta _0=0.001$ whereas $G = 78.7 (\mathcal{F}=9.8)$ at $\delta _0=0.02$ . This significant reduction in gain is attributed to the nonlinear saturation of the disturbance growth. Conversely, the entropy production increases with increasing $\delta _0$ . This is because the definition in (2.11) is based on the total flow field, irrespective of its decomposition into base and disturbance components. In figure 6(a), the optimal forcing exhibits a locally concentrated spatial distribution just upstream of the laminar separation bubble of the base flow. In this linear regime ( $\delta _0=0.001$ ), the temporal distribution shows a periodic pattern. Thus, for the small $\delta _0$ case, the application of optimal forcing exploits the amplification mechanism of the disturbance field over the separation bubble, as reported in De Pando, Schmid & Sipp (Reference De Pando, Schmid and Sipp2017) and Thomareis & Papadakis (Reference Thomareis and Papadakis2018).

Increasing the amplitude of external forces significantly alters the spatiotemporal characteristics. Spatially, the forcing concentrates in two primary regions: near the leading edge and around the separation point of the base flow field. The balance of intensity between these regions is governed by $\delta _0$ and $t_{\!f}$ . In the leading-edge region (especially for large $\delta _0$ and small $t_{\!f}$ ), the forcing exhibits a streamwise-elongated spatial structure within the laminar boundary layer. This component is associated with a monotonically decreasing temporal profile, which serves to increase wall friction and generate high entropy production. Conversely, in the separation region, the forcing acts on the reverse flow. While this typically induces periodicity in linear regimes, high-amplitude forcing causes the perturbed flow to deviate significantly from the base flow, breaking the periodicity. Consequently, the temporal distribution in this region takes the form of a damped oscillation, reflecting the physical need for a strong initial impulse to shift the flow state. The superposition of these components results in a temporal distribution characterized by a damped oscillation with a positive mean value.

For $t_{\!f} \gt 1.0$ , the optimal forcing exhibits a non-uniform spatial distribution, with the temporal distribution dominated by a damped oscillation. Table 2 lists the calculated evaluation function values and disturbance growth rates. Notably, at larger evaluation times $(t_{\!f} = 4.0, 5.0)$ , the disturbance growth rate exhibits a significantly stronger dependence on the phase of the base flow compared with the entropy production. This discrepancy implies that distinct flow trajectories can yield similar evaluation function values. Consequently, the presence of local optima within the current optimization framework cannot be ruled out.

Table 2. Calculated evaluation function value ( $\mathcal{F}$ ) at $\delta _0 = 0.01$ for various phases of the base flow. Values in brackets indicate the disturbance growth rate.

Figure 6. Spatial distributions of $g_x$ at $t_{\!f} =2.0, \phi = -\pi$ for (a) $\delta _0 = 0.001$ and (b) $\delta _0 = 0.02$ ; (c) the temporal distributions at $\textit{AoA}=2^{\circ }$ . The red dotted line indicates $\bar {u}_x=0$ contour of the base flow.

For longer evaluation time ( $t_{\!f}$ ), the periodic oscillation component in $f(t)$ becomes dominant; consequently, the time-averaged entropy production approaches a nearly constant value. As shown in table 2, the evaluation function value, averaged over the base flow phases, is $\mathcal{F}=9.0$ at $t_{\!f}=4.0$ and $9.5$ at $t_{\!f}=5.0$ . This trend persists with further increases in evaluation time: for example, $\mathcal{F}=9.4$ at $\phi =-\pi , t_{\!f} = 6.0$ and $\mathcal{F}=9.6$ at $\phi =-\pi , t_{\!f}=7.0$ . This consistency suggests that the flow mechanisms induced by the optimal forcing remain qualitatively consistent as the evaluation time increases. Furthermore, the table indicates that the evaluation function is largely independent of the base flow phase. The unsteadiness at $\textit{AoA}=2^{\circ }$ is primarily driven by von Kármán vortex shedding, whereas the spatial distribution of the optimal forcing is highly concentrated near the separation point of the base flow field. Given that the downstream vortex shedding exerts minimal influence on the upstream separation point, the observed independence of the control performance from the shedding phase is physically reasonable.

Figure 7 illustrates the spatial and temporal characteristics of the optimal forcing. Spatially, the forcing is primarily concentrated at the separation point of the base flow. We confirmed a high degree of similarity in the distributions across different $t_{\!f}$ , with the exception that the upstream component near the leading edge gradually fades. This suggests that for longer $t_{\!f}$ , the periodic forcing at the separation point becomes the dominant control mechanisms. Notably, for $t_{\!f}\gt 2.0$ , the forcing patterns remain largely invariant.

To quantify the temporal characteristics, we analysed the frequency of $f(t)$ using VMD (Dragomiretskiy & Zosso Reference Dragomiretskiy and Zosso2013) combined with a Hanning window. Averaging the results from the eight cases ( $t_{\!f} = 4.0, 5.0$ ) listed in table 2 yields a dominant frequency of $St=2.11$ with a standard deviation of only $0.08$ . This low standard deviation demonstrates the robustness of the optimal forcing frequency, confirming its independence from the base flow phase and minor variations in evaluation time.

Figure 7. Spatial and temporal characteristics of the optimal forcing at $\phi =-\pi , \delta _0 = 0.01, \textit{AoA}=2^{\circ }$ : spatial distributions for (a) $t_{\!f} = 2.0$ and (b) $t_{\!f}=5.0$ ; (c) temporal distributions $t_{\!f}$ . The red dotted line indicates the $\bar {u}_x=0$ contour of the base flow.

Figure 8. Changes in the flow field induced by the application of optimal forcing at $\textit{AoA}=2^{\circ }, \phi =-\pi$ with an amplitude of $\delta _0=0.01$ . The skin friction is shown for the case of $t_{\!f}=4.0$ , where the black curves represent the $\bar {u}_x=0$ isocontours of the base flow, and the black dashed line indicates the evaluation time $t=t_{\!f}$ . The forcing was applied from $t=0$ to $t=t_{\!f}$ . Here (a) skin friction; (b) lift coefficient.

3.3. Changes of the flow field

In this section, we examine the changes in the flow field induced by the application of optimal forcing. The external forcing produced clockwise vortices, which intensified as they propagated downstream from the separation point by extracting the energy from the shear layer of the base flow. Figure 8 shows the characteristics of the perturbed flow field in terms of skin friction and lift coefficients. An increase in the evaluation function values correlates with the integral value of the square of the skin friction coefficient, as defined in (2.10). In figure 8(a), a significant increase in skin friction is observed in association with vortex advection; the vortex trails appear as lines with a positive slope (indicating downstream movement). Moreover, the momentum transfer induced by the vortex advection caused the separation region to contract. Consequently, the separation point, which had been oscillating between points $-0.175$ and $-0.163$ in the base flow field, shifted downstream. This corresponds to the downstream shift of the separation point trace during the forcing application shown in figure 8(a). The contraction of the separation region accelerated the flow, resulting in a pressure decrease on the suction surface and a subsequent increase in the lift coefficient, as shown in figure 8(b). In this figure, we observe that the lift coefficient increases continuously with $t_{\!f}$ . This corresponds to the fact that the temporal distribution of the optimal forcing exhibits sustained periodicity for long $t_{\!f}$ , as shown in figure 7(c). Therefore, we conclude that the lift increase is directly driven by vortex convection and can be sustained by the periodic passage of vortices through the separation region. This mechanism will be further investigated in the next section. As shown in figure 8(a), we also observe that the transient changes to the flow field persist even after the external forcing is removed, which will be discussed in § 3.4.

Figure 9 illustrates the physical process of vortex formation and convection along the airfoil. In this figure, the external forcing applied at the separation point of the base flow produces clockwise-rotating vortices that convect downstream. These vortices evolve by extracting energy from the shear layer between the free stream and the separated flow region. This process is analogous to the Kelvin–Helmholtz instability in inviscid flow. Concurrently with this vortex evolution in the shear layer, strong vorticity formation near the wall is observed due to the interaction between the shear layer and the wall surface. This eventually leads to the generation of secondary vortices (see $t_{\!f}=1.25$ in the same figure), resulting in a local increase in skin friction. This vortex shedding process shows a qualitative similarity to the results of resolvent analysis over a separation bubble, as reported by Ricciardi, Wolf & Taira (Reference Ricciardi, Wolf and Taira2022). However, the present results demonstrate the emission of vortices with constant convective speed, as evidenced by the vortices trails in figure 8(a), which prevents complexities such as vortex pairing. This stable emission of vortices was consistently observed in the present simulations across multiple evaluation times.

Figure 9. Spatial distribution of vorticity in the instantaneous flow field at $\textit{AoA}=2^{\circ }, \phi =-\pi$ . The optimal forcing was calculated for $t_{\!f} = 5.0, \delta _0 = 0.01$ .

3.4. Transient changes of the flow field

Next, we investigate the characteristics of the flow field changes driven by the advection of vortices induced by the optimal forcing. As shown in figure 8(a), this flow response exhibits a low-frequency oscillation that persists as transient dynamics even after the forcing is removed. Since this behaviour shares characteristics with transonic buffet in high-subsonic flows, analysing these transient dynamics can shed light on the underlying oscillation mechanisms. In the remainder of this section, we present the results of VMD-NCS analysis for the case of $t_{\!f} = 2.0$ . The flow evolution was computed from $t=-4.12$ to $40.4$ , with optimal forcing applied during the interval $0 \leqslant t \leqslant 2.0$ . We extracted $512$ snapshots at a sampling interval of $\Delta t = 0.087$ . The VMD-NCS parameters were set to $K=7$ and $\alpha = 5000$ . The dependence of the analysis results on these parameters is examined in Appendix.

Figure 10 shows the time history of the magnitude of each ICS during its evolution, based on the definition in (2.23). In this figure, the dominant ICS at $St=1.991$ is closely related to the primary peak frequency at $St=2.24$ of the base flow field (figure 3), exhibiting a slight frequency shift under transient conditions. This dominant ICS at $St=1.991$ exhibits a temporal decrease in magnitude between $0.0 \leqslant t \leqslant 25.0$ , accompanied by the evolution of intermediate modes in the flow field. During this transient decay, two adjacent ICS at $St=1.794$ and $St=2.091$ emerge, governing the subsequent flow changes. As shown in figure 11 at $t=8.08$ , both of these ICS are characterized by the convection of vortices through the massive separation region on the suction surface. This observation strongly supports that the hypothesis that the breathing mechanism of the separation bubble is primarily driven by the momentum transfer associated with these convecting vortices.

Figure 10. Time evolution of the norm of each ICS in the flow field perturbed by optimal forcing ( $t_{\!f}=2.0, \delta _0=0.01$ ). The forcing is activated during the interval $t=0$ to $t=2.0$ .

Following the development of these transient modes, the magnitude of the ICS at $St=0.092$ increases significantly. The frequency of this mode corresponds to the low-frequency oscillation of the base flow field in figure 3, and it also corresponds to the low-frequency oscillation observed after the removal of external forcing, as shown in figure 8(a). This low-frequency oscillation induced by external forces persists for a much longer duration ( $t \leqslant 30$ ) than the original excitation period $0 \leqslant t \leqslant 2$ . Recently, the connection of this low-frequency oscillation and transonic buffet is suggested by Moise et al. (Reference Moise, Zauner and Sandham2024), and clarifying the self-sustaining mechanism of this low-frequency oscillation may be critical to shed light the mechanism of low-frequency oscillations and hence the transonic buffet.

Figure 11. Instantaneous distribution of $u_x$ in ICS at $t=8.08$ : (a) St = 1.794; (b) St = 2.091.

Figure 12. Spatiotemporal evolution of the ICS mode at $St=0.092$ , visualized over a half-period.

To explain the mechanism maintaining these low-frequency oscillations, a feedback loop mechanism similar to that observed in trailing-edge vortex noise (Desquesnes, Terracol & Sagaut Reference Desquesnes, Terracol and Sagaut2007; Ricciardi et al. Reference Ricciardi, Wolf and Taira2022) may be applicable. Figure 12 presents the instantaneous flow field of the ICS at $St=0.092$ . This figure visualizes a half-period of the flow field associated with the centre frequency of $St=0.092$ , showing the expansion and contraction of the separation region. The emission of acoustic waves from the trailing edge, accompanied by this low-frequency oscillation, is clearly observed; these waves interact with the airfoil surface, inducing velocity disturbances, particularly near the leading edge. At $T=21.5$ and $22.9$ in the same figure, velocity fluctuations are generated near the leading edge on the suction surface. These disturbances evolve into vortices that convect downstream, modulating the boundary layer state from deceleration to acceleration. As discussed regarding figure 8(a), the expansion and contraction of the separated boundary layer is attributed to momentum transfer caused by vortex convection along the airfoil. Therefore, for oscillations to be self-sustained, a mechanism that continuously supplies velocity fluctuations to the flow field is required. The velocity disturbances shown in figure 12 indicate that the interaction between acoustic waves and the boundary layer serves as this source.

The decay of the amplitude for the ICS at $St=0.092$ after $t \geqslant 20$ indicates that the vortex generation induced by acoustic waves is insufficient to maintain the oscillation in the present flow field. Conversely, when the acoustic wave intensity is sufficiently high and the boundary layer is close to neutral stability, this low-frequency oscillation is expected to be sustained. Indeed, at a slightly higher Reynolds number of $\textit{Re}=5\times 10^4$ , the onset of the buffet phenomenon has been reported by Moise et al. (Reference Moise, Zauner and Sandham2024).

4.1. Simplification of temporal and spatial distribution for control application

In this section, we examine the application of optimal forcing to suppress flow separation around an airfoil. As demonstrated in figure 8, applying external forcing over an appropriate time interval has the potential to increase the lift. This section focuses on using such external forcing to control the separated flow. Before implementing the control, however, the spatiotemporal distribution of the optimal forcing requires several simplifications to match the constraints of a realistic flow control system.

Based on the characteristics of the optimal forcing over long evaluation times (specifically $t_{\!f} = 4.0, 5.0$ ), we focus on forcing applied at the separation point of the base flow. The external forcing exhibits a concentrated distribution in the range $-0.35 \leqslant x \leqslant 0$ on the suction surface. We obtained a simplified spatial distribution by setting all values outside this range to zero and applying a cosine-tapered window function to ensure continuity at the domain boundaries. Specifically, the cosine-tapered function was treated as a one-dimensional function in the $x$ -direction; it was set to unity within $[-0.35, 0]$ and included transition regions over $[-0.4, -0.35]$ and $[0, 0.1]$ . This operation was applied only to the $x$ -component, while the $y$ -component was set to zero. The spatial profile was derived specifically from the spatial distribution of the optimal forcing at $\phi =-\pi , t_{\!f}=5.0, \delta _0=0.01$ , as shown in figure 7(b).

For the temporal distribution, we adopted a sinusoidal pattern with a frequency of $St=2.11$ , consistent with the results in § 3.2. We utilized a pulsed application of the optimal forcing with a repetition interval $\Delta T$ , defined as

(4.1) \begin{equation} \hat {f} (t) = \left \{ \begin{array}{ll} - A\sin (2 \pi \omega s), & s \leqslant 1/\omega, \\ 0, & s \gt 1/\omega, \end{array} \right . \end{equation}

where $A$ is the forcing amplitude, $\omega (=2.11)$ is the forcing frequency, and the time is given by $t = m \Delta T + s$ with $m \in \mathrm{N}$ and $0 \leqslant s \lt \Delta T$ .

Figure 13 shows the spatiotemporal distribution of the external forcing for the impulsive control of the flow field. Compared with the original distribution in figure 7(b), the present distribution eliminates low-amplitude regions while retaining the key structural features. Regarding the temporal statistics, the time-averaged power of the external force applied over one period $\Delta T$ can be calculated as

(4.2) \begin{equation} \mathcal{P}^2 = \frac {1}{\Delta T} \int _0^{\Delta T} | \hat {f} (t) |^2 {\rm d}t = \frac {A^2}{2 \omega } \frac {1}{\Delta T}. \end{equation}

This equation indicates that the power is inversely proportional to the pulse interval $\Delta T$ . Since the present optimal control analysis does not inherently determine the optimal pulse interval $\Delta T$ for long-term control, the appropriate interval will be determined through numerical experiments in the remainder of this section.

Figure 13. Spatiotemporal characteristics of the simplified external forcing used for impulsive flow control: (a) spatial distribution and (b) temporal distribution.

Figure 14. Time evolution of $C_l$ during the impulsive application of external forcing.

4.2. Characteristics of forced flow field

Figure 14 shows the time-evolution of $C_l$ under the impulsive application of external forcing. In this figure, the $C_l$ trajectory gradually deviates from the base flow and transitions into a new state sustained by the external force, exhibiting a significant increase in $C_l$ compared with the base flow value (except for the case $A=3.0, \Delta T=0.65$ ). We found that the strong interactions between vortices caused flow unsteadiness when the pulse interval was short $(\Delta T\lt 0.75)$ . The intensity and pulse interval of the external force significantly affect the characteristics of the final flow field. Since the disturbance field cumulatively modifies the mean flow, we simulated the flow evolution for a sufficiently long duration ( $t \geqslant 200$ ) to confirm the convergence of the flow statistics. In this figure, the early evolution of $C_l$ exhibits similarities to the numerical results of separation control using plasma actuators reported by Asada et al. (Reference Asada, Nonomura, Aono, Sato, Okada and Fujii2015) (see figure 6 therein). They noted that the impulsive actuation of a plasma actuator initially causes a decrease in the lift coefficient; this tendency was also observed in the present calculations.

Figure 15 presents the converged aerodynamic coefficients under impulsive forcing. We define the control efficiency $E_{{l, d}}$ as the gain in lift (drag) coefficient per unit input power, given by

(4.3) \begin{equation} E_{l, d} = \frac {\bar {C}_{l, d} [\boldsymbol{q}] - \bar {C}_{l, d} [\boldsymbol{Q}]}{\mathcal{P}}, \end{equation}

where $\bar {C}_{l, d} [\boldsymbol{q}]$ and $\bar {C}_{l, d} [\boldsymbol{Q}]$ represent the time-averaged lift (drag) coefficients for the forced ( $\boldsymbol{q}$ ) and base ( $\boldsymbol{Q}$ ) flows, respectively. A single pulse is applied per interval $\Delta T$ , with amplitude $A$ and a spatial distribution based on figure 7. When the pulse interval is short ( $\Delta T \leqslant 0.75$ ), close to the forcing period, complex vortex interactions prevent the flow from stabilizing. For larger intervals, the lift enhancement generally diminishes as $\Delta T$ increases. This occurs because the lift gain relies on vortex-induced momentum exchange, which scales with the pulse frequency $(1/\Delta T)$ . However, at high amplitudes, nonlinear vortex interactions alter this trend, producing a local maximum in lift at $\Delta T=1.2$ . Compared with the base flow values of $\bar {C}_l = -0.8572$ and $\bar {C}_d = 0.05183$ listed in table 1, the forced flow yields a higher lift coefficient than the base flow in all cases, except when the amplitude is low ( $A=1$ ) and the interval is long $\Delta T \geqslant 1.5$ .

Figure 15. Converged aerodynamic coefficient under impulsive forcing: (a) time-averaged lift coefficient $\bar {C}_l$ and (b) drag coefficient $\bar {C}_d$ . Error bars denote the standard deviation of the fluctuations.

In terms of power efficiency ( $E_l$ ), the optimal strategy depends on the pulse interval. At small intervals ( $\Delta T=0.75$ ), low-amplitude forcing is more efficient as it leverages the linear amplification of disturbances. At longer intervals ( $\Delta T=1.5$ ), however, high-amplitude forcing achieves better efficiency through nonlinear vortex interactions. Notably, the reduced frequency $F^+ = 1.25$ corresponding to the peak lift ( $\Delta T=1.2, A=4$ ) agrees well with the value of $F^+=1.2$ reported for stalled airfoils at low Reynolds numbers by Greenblatt et al. (Reference Greenblatt, Goksel, Rechenberg, Schule, Romann and Paschereit2008). This reduced frequency is defined as $F^+ = X_{te}/ (\Delta T U_{\infty })$ , where $X_{te}$ is the distance from the forcing location to the trailing edge and $U_{\infty }$ is the free stream velocity.

Next, we investigate the physical mechanism underlying the increase in $C_l$ in the forced flow field. Figure 16 shows the spatial distribution of streamwise velocity for both the base flow and the forced flow. In the time-averaged base flow, a separation region extends beyond the trailing edge of the airfoil. This separation induces flow acceleration on the lower surface, which explains the negative lift observed in the base flow (Ohtake et al. Reference Ohtake, Nakae and Motohashi2007; Winslow et al. Reference Winslow, Otsuka, Govindarajan and Chopra2018). The application of external forcing facilitates momentum transfer between the vortices and the mean flow, leading to the reattachment of the separated flow on the upper surface, as shown in figures 16(b) and 16(c). This modification of the mean flow field causes acceleration on the upper surface, resulting in an increase in $C_l$ due to the pressure difference between the upper and lower surfaces. We note that this acceleration region exhibits characteristics of shock waves, indicated by a local Mach number greater than unity. However, the generation of these shock waves does not induce significant changes in the flow field due to their weak intensity.

Figure 16. Spatial distribution of time-averaged streamwise velocity ( $u_x$ ) for the base and forced flows. The red dashed line indicates the $u_x = 0$ contour: (a) base flow; (b) $A=2.0, \Delta T=0.75$ ; (c) $A=4.0, \Delta T=1.2$ .