The structure and dynamics of the laminar separation bubble

Abstract

A novel selective mode decomposition, proper orthogonal decomposition and dynamic mode decomposition methods are used to analyse large-eddy simulation data of the flow field about a NACA0012 airfoil at low Reynolds numbers of $5\times 10^4$ and $9\times 10^4$, and at near-stall conditions. The objective of the analysis is to investigate the structure of the laminar separation bubble (LSB) and its associated low-frequency flow oscillation (LFO). It is shown that the flow field can be decomposed into three dominant flow modes: two low-frequency modes (LFO-Mode-1 and LFO-Mode-2) that govern an interplay of a triad of vortices and sustain the LFO phenomenon, and a high-frequency oscillating (HFO) mode featuring travelling Kelvin–Helmholtz waves along the wake of the airfoil. The structure and dynamics of the LSB depend on the energy content of these three dominant flow modes. At angles of attack lower than the stall angle of attack and above the angle of a full stall, the flow is dominated by the HFO mode. At angles of attack above the stall angle of attack the LFO-Mode-2 overtakes the HFO mode, triggers instability in the LSB and initiates the LFO phenomenon. Previous studies peg the structure, stability and bursting conditions of the separation bubble to local flow parameters. However, the amplitude of these local flow parameters is dependent on the energy content of the three dominant flow modes. Thus, the present work proposes a more robust bursting criterion that is based on global eigenmodes.

1. Introduction

Airfoils operating at low Reynolds numbers tend to induce a laminar separation bubble (LSB) on their suction surface. A laminar boundary layer develops near the airfoil leading edge, detaches, travels away from the airfoil surface due to the adverse pressure gradient (APG) and creates a region of a separated flow. The flow then undergoes a transition to turbulence and reattaches to the airfoil surface, and the LSB forms in the separated flow region as seen in figure 1 (Melvill Jones 1934). The LSB has attracted much attention because it plays an important part in determining the behaviour of the boundary layer on the airfoil surface and consequently the stalling characteristics of the airfoil. Furthermore, the LSB was also observed on circular cylinders at critical Reynolds number and the interaction of the shock wave and boundary layer on a flat plate at supersonic speeds (Tani 1964).

Figure 1. (a,b) Streamline patterns of the time-averaged flow field plotted around a NACA0012 airfoil at the stall angle of attack.

The averaged shape of the LSB on the $x$–$y$ plane was first described by Horton (1968). The flow separates at point ${\boldsymbol {S}}$ and reattaches at point $\boldsymbol {R}$, and the length of the LSB is measured by the distance between these two points along the airfoil chord. Beneath the LSB lies a small counter-rotating bubble called a secondary bubble. The turbulent boundary layer downstream of the bubble has more momentum near the surface to resist the APG and avoid a new separation. The length of the LSB depends on the angle of attack, the Reynolds number and the type of airfoil section. If the Reynolds number based on the displacement thickness is greater than $400$–$500$, the bubble is termed short. Whereas, if the Reynolds number is less than this band of values, the bubble is termed long (Owen & Klanfer 1953).

Statistics of the flow field and the shape of the LSB in the spanwise direction in numerical simulation depend on the extent of the computational domain in this direction. Several numerical studies examined the effect of the spanwise width on the statistics of the flow field (Mary & Sagaut 2002; Almutairi, Jones & Sandham 2010; Alferez, Mary & Lamballais 2013; Zhang & Samtaney 2016). Alferez et al. (2013) estimated the two-point correlation coefficient of the fluctuating streamwise velocity component along the spanwise direction and observed that the correlation coefficient oscillates around zero when the separation between the two points is a half-chord. The authors concluded that their one-chord-width computational domain is sufficient to produce a flow that is independent of the extent of the flow in this direction. Zhang & Samtaney (2016) investigated the flow field around a NACA0012 airfoil at Reynolds number of $5\times 10^4$ using computational domains of spanwise extent of $0.1$–$0.8$ chords. The different spanwise widths of the computational domain result in close prediction of the time-averaged aerodynamic quantities, and the instantaneous lift and drag coefficients at the spanwise width of $0.4$ and $0.8$ are similar. However, flow separation in the vicinity of the leading edge is notably three-dimensional for the spanwise width of $0.8$, while remaining two-dimensional in the remaining aspect ratios.

Mccullough & Gault (1951) classified the airfoil stall into three main categories: (1) leading-edge, (2) thin-airfoil and (3) trailing-edge stall. The leading-edge stall results from the flow separation near the leading edge without any subsequent reattachment downstream of the separation. In the thin-airfoil stall, the flow reattaches downstream the separation, and then the reattachment point moves towards the trailing edge as the angle of attack increases. The trailing-edge stall initiates at the trailing edge where the flow separates, and the separation point moves towards the leading edge as the angle of attack increases. Gaster (1967) who was the first to investigate the stability of the LSB found that bubble bursting occurs either by a gradual increment in the bubble length or by a suddenly discontinuous event.

Observations at pre-stall conditions have shown that the LSB is stable, and a short bubble forms on the suction surface of the airfoil and on a flat plate with an adjustable pressure distribution. The structure and the length of the LSB are altered by changing the conditions of the incoming free stream or by smoothly changing the angle of attack (Gaster 1967; Rist & Maucher 2002; Marxen & Rist 2010; Marxen & Henningson 2011; Alferez et al. 2013; Yarusevych & Kotsonis 2017; Istvan & Yarusevych 2018). At near-stall conditions, the LSB becomes unstable and bursts to form a long bubble. Consequently, the flow field switches between an attached phase and a separated phase, and triggers a low-frequency flow oscillation (LFO) (Zaman, Bar-Sever & Mangalam 1987; Zaman, McKinzie & Rumsey 1989; Bragg et al. 1996; Broeren & Bragg 1998, 1999, 2001; Rinoie & Takemura 2004; Tanaka 2004; Almutairi & Alqadi 2013; Eljack 2017; Elawad & Eljack 2019; Eljack & Soria 2020; Aniffa & Mandal 2023a,b; Dellacasagrande et al. 2023).

Eljack et al. (2021) used the compressible flow datasets generated by Eljack (2017) at ${\textit {Re}}_c = 5\times 10^4$ and by Elawad & Eljack (2019) at ${\textit {Re}}_c = 9\times 10^4$ to characterize the flow field around a NACA0012 airfoil, and demarcated the angle of attack range into three regimes: a pre-stall attached flow, a post-stall intermittently separating flow and a post-stall fully separated flow. The authors used a conditional time-averaging to describe the flow field in detail and investigate the effects of the angle of attack on the characteristics of the flow field. The authors reported that at angles of attack lower than the stall angle of attack, the LSB remains intact, and the mean flow field is attached. The instantaneous flow exhibits an LFO due to shear-layer flapping and vortex shedding; however, the amplitude of oscillations is considerably small compared with the mean flow. Thus, the mean flow remains attached downstream of the LSB. At angles of attack higher than the stall angle of attack, a trailing-edge bubble (TEB) forms in the vicinity of the trailing edge on the suction surface of the airfoil. Thus, the mean flow field contains two bubbles, the LSB and a TEB, bounded by four half-saddles. The first half-saddle point is the separation point in the vicinity of the leading edge ($\boldsymbol {S_1}$). The second half-saddle point is the reattachment point located just downstream the LSB ($\boldsymbol {R_1}$). The third half-saddle is the separation point near the trailing edge ($\boldsymbol {S_2}$). The fourth half-saddle point is the reattachment point located at the trailing edge ($\boldsymbol {R_2}$). The first two half-saddle points constitute the LSB, and the latter two half-saddle points constitute the TEB as seen in figure 2. At relatively high angles of attack, when the two bubbles start to merge, the two half-saddle points $\boldsymbol {R_1}$ and $\boldsymbol {S_2}$ move away from the wall and form a full-saddle point around mid-chord at the angle of attack of maximum LFO. The instantaneous shape of the LSB is altered between a short bubble and a long bubble and sometimes merges with the TEB depending on the angle of attack. Consequently, the flow field undergoes a transition regime in which the LFO develops until it reaches a quasi-periodic switching between separated and attached flow. At angles of attack higher than the angle of attack of maximum LFO, the flow field and the aerodynamic characteristics are overwhelmed by quasi-periodic and self-sustained LFO until the airfoil approaches the angle of a full stall, the LSB and the TEB merge to form an open bubble in the mean sense and the mean flow remains massively separated. The authors concluded that most of the observations reported in the literature about the LSB and its associated LFO are neither thresholds nor indicators for the inception of the instability, but rather are consequences of it.

Figure 2. (a,b) Streamline patterns of the time-averaged flow field plotted around a NACA0012 airfoil at the angle of attack of maximum amplitude of the LFO (SP, saddle point).

The dynamics of the flow plays a crucial role in the underlying mechanism behind the bursting and reformation of the LSB and its associated LFO. The dynamics of any turbulent or transitional flow is dominated by an organized motion known as coherent structures. These are recognizable flow structures that survive dissipation by viscosity for a relatively long time and have a typical lifecycle. The origin of such coherent structures is an unstable flow mode that feeds on the mean flow. The temporal and spatial evolution of such organized flow motion is of major importance because the instantaneous solution of the flow is dependent on it. If the dynamics of the flow is well understood and described mathematically, researchers will develop better models that predict the phenomenon, thus reducing the huge cost of modelling the phenomenon numerically using large-eddy simulation (LES) or experimentally in expensive test rigs. Furthermore, such coherent motion contains most of the kinetic energy and acts as a primary mechanism that dissipates energy. If the dynamics is better understood, more efficient and effective aerodynamic shapes can be engineered. Most importantly, such organized flow motion induces oscillations in the aerodynamic forces, vibrations, noise and drag. Thus, suppressing or energizing such coherent structures greatly improves the aerodynamic performance of airfoils. Hence, researchers and engineers could develop smart control means that remove the undesired effects of the phenomenon, utilize the phenomenon when it presents a control opportunity and improve the aerodynamic performance of airfoils.

The aforementioned organized motion is embedded into stochastic spatiotemporal data, and its extraction is not as straightforward as it might seem. The low-order statistics characterize the flow field in the mean sense. However, valuable information is lost in the averaging process. For instance, the spatial evolution of coherent structures in the flow field cannot be captured using the low-order statistical moments. Furthermore, the dynamics of the flow and the evolution of the various flow modes in time cannot be described using low-order statistics. High-order statistical methods like proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) are used to extract deterministic organized flow motion from stochastic spatiotemporal data. The power of the POD method lies in the fact that the decomposition of the flow field in the POD eigenfunctions converges optimally fast in the energy sense ($L^2$-norm). Whereas, the power of the DMD method is that it provides growth rates in addition to the shape of the dynamic modes of the flow at various frequencies.

Most previous works on the structure and dynamics of LSBs induced on flat plates and airfoils at near-stall conditions were experimental. The measurements were mostly acquired at a single point or a line, and two-dimensional simultaneous measurements are rare. Since the POD and the DMD methods are applied primarily to plane data or three-dimensional data, application of the POD and the DMD methods to investigate the LSB and its associated LFO phenomenon in the flow field around an airfoil at near-stall conditions is not so common in the literature. Classic POD, spectral POD, extended POD and DMD were applied to particle image velocimetry data of LSBs formed over airfoils, flow over flat plates with adjustable pressure distribution featuring LSBs and over and behind a two-dimensional forward–backward-facing step. The extracted flow modes feature vortex shedding, counter-rotating vortices downstream of the bubble and flapping motion (Lengani et al. 2014; Istvan & Yarusevych 2018; Alessandri et al. 2019; Fang & Tachie 2019; Verdoya et al. 2021; Toppings & Yarusevych 2022; Weiss et al. 2022; Dellacasagrande et al. 2023; Aniffa & Mandal 2023b; Braud, Podvin & Deparday 2024).

Recently, Eljack & Soria (2020) utilized a conditional time-averaging, the POD method and a conditional phase-averaging to reveal the underlying mechanism that generates and sustains the LFO phenomenon. Those authors analysed LES data at ${\textit {Re}}_c = 5\times 10^4$ and showed that a triad of three vortices, two co-rotating vortices (TCV) and a secondary vortex counter-rotating with them, is behind the quasi-periodic self-sustained bursting and reformation of the LSB and its associated LFO phenomenon. They reported that a global oscillation in the flow field around the airfoil is observed in all of the investigated angles of attack, including at an angle of attack of zero. They have shown that the global oscillation in the flow field is self-sustained due to an interplay of the triad of vortices. The authors found that when the direction of the oscillating flow is clockwise, in the streamwise direction on the suction surface of the airfoil, it adds momentum to the boundary layer and helps it to remain attached against the APG and vice versa. However, the authors did not discuss the bursting conditions of the LSB in their model. Also, the dynamics of the flow and how various flow modes interact with each other are lost in the averaging processes. Furthermore, the short LSB formed over the airfoil at ${\textit {Re}}_c = 5\times 10^4$ spans more than $40\,\%$ of the airfoil chord; thus, an investigation at a higher Reynolds number is merited. The model based on the existence of a triad of vortices in the region of the LSB presented by Eljack & Soria (2020) agrees very well with the hypothesis that there are sub-bubbles within the LSB (Dallmann 1983; Dallmann, Vollmers & Su 1997; Theofilis, Sherwin & Abdessemed 2004).

Various bursting conditions for the LSB over airfoils at pre-stall and near-stall conditions and on a flat plate were set by Owen & Klanfer (1953), Gaster (1967) and Diwan, Chetan & Ramesh (2006). Gaster (1967) found that the structure of the LSB depends on two parameters. The first parameter is the Reynolds number of the separated boundary layer and the second parameter is a function of the pressure rise over the region occupied by the bubble. Then, he determined conditions for the bursting of short bubbles by a unique relationship between these two parameters. Diwan et al. (2006) hypothesized that the height of the LSB has effects on the bursting of the LSB. They suggested a bursting parameter given by the velocity gradient across the bubble and the height of the bubble. The bursting parameter is plotted versus the Reynolds number based on the height of the bubble and the velocity at the maximum height of the LSB. The bursting condition is set such that the bubble is short if the bursting parameter is greater than $-28$. However, the bursting criteria of Gaster (1967) and Diwan et al. (2006) use changes in local flow parameters in the vicinity of the leading edge and across the bubble. Thus, the bursting criteria are not universal in classifying and distinguishing short bubbles from long bubbles as pointed out by Diwan et al. (2006) and confirmed by Alferez et al. (2013).

The objective of the present paper is to carry out a detailed dissection of the flow field and shed some light on the dynamics of the flow about a NACA0012 airfoil at ${\textit {Re}}_c = 5\times 10^4$ and $9\times 10^4$, and near-stall conditions. The POD and the DMD methods were applied to datasets sampled on the $x$–$y$ plane including the velocity components, the pressure and the aerodynamic coefficients. The datasets span four low-frequency cycles and were locally time-averaged every $50$ time steps and ensemble-averaged in the spanwise direction on the fly before they were recorded. The conditional time-averaging introduced by Eljack & Soria (2020) is further developed to implement a selective mode decomposition (SMD) method to the flow field. The SMD method utilizes time histories of aerodynamic forces to recover the flow modes that induce oscillations in these forces. The SMD method recovers flow modes based on time-averaging of the data; therefore, it estimates the amplitude of each of the flow modes accurately. Whereas, the POD method recovers flow modes under the condition of optimizing the energy content of each mode; thus, it overestimates the energy content of some of the flow modes. The amplitudes of the POD and the SMD coefficients are investigated to examine the effects of the percentage of energy content of the most dominant eigenmodes on the bursting of the LSB. The following section introduces the SMD method. The computational set-up and validations of the results are presented in § 3. Section 4 presents the application of the POD, the DMD and the SMD methods. The dynamics of the flow and the bursting criterion are discussed in § 5 and 6, respectively.

2. Selective mode decomposition

Applied research focuses on suppressing or enhancing phenomena like mixing, aerodynamic forces, skin-friction drag, noise, vibrations, induced acoustic waves and shockwave instability. The POD method recovers flow modes based on their energy content. However, the most dominant mode is not necessarily the most important flow mode that induces the phenomenon of interest. Also, there is reasonable doubt that the POD modes are mathematical flow structures that satisfy Lumely's condition and objectively recover the most energetic flow modes but do not necessarily physically exist in the flow. Furthermore, one of the major shortcomings of the POD method is that the POD eigenfunctions, used as input parameters in flow control, depend on the flow field from which they were extracted. The dynamics of the flow is naturally integrated on the solid surface and induces aerodynamic forces. That is, the pressure field makes the flow feel the wall, and the dynamics of the flow creates fluctuations in the aerodynamic forces. Thus, the time history of an aerodynamic force can be used to recover the flow mode that induces fluctuations in it. Directly identifying the flow mode associated with the phenomenon of interest would provide a direct answer to which flow mode induces the undesired effect; consequently, it can be directly suppressed by applying suitable control. Furthermore, the fluctuating aerodynamic forces change automatically when the flow field changes. Thus, the fluctuating aerodynamic forces can be used as an alternative to POD eigenfunctions to extract selective flow modes. In the present work, the conditional time-averaging suggested by Eljack & Soria (2020) was used to recover flow modes corresponding to each aerodynamic force.

The time history of the lift coefficient ($\textit{C}_{\!_\textit{L}}(t)$) can be used as a reference to conditionally time-average the flow field (Eljack & Soria 2020; Eljack et al. 2021). The time average of an instantaneous variable of the flow field like the streamwise velocity component ($\boldsymbol{\mathsf{u}}_1(x,y,t)$) or its discrete form ($\boldsymbol{\mathsf{u}}_1(x,y,n)$) can be defined on three different levels: a mean-lift ($\overline {\textit{C}_{\!_\textit{L}}(t)}$), a high-lift ($\widehat {\textit{C}_{\!_{\textit{L}}}(t)}$) and a low-lift () time average. Figure 3 illustrates the concept of conditional time-averaging that is based on the lift coefficient. The mean-lift time average of the streamwise velocity ($\overline {\boldsymbol{\mathsf{U}}_1}(x,y)$) is simply the time average of all data points of the streamwise velocity component. The high-lift average of the streamwise velocity ($\widehat {\boldsymbol{\mathsf{U}}_1}(x,y)$) is the time average of data points of $\boldsymbol{\mathsf{u}}_1$ that correspond to $\textit{C}_{\!_\textit{L}}(t) > \overline {\textit{C}_{\!_\textit{L}}}$. The low-lift average of the streamwise velocity () is the time average of data points of $\boldsymbol{\mathsf{u}}_1$ that correspond to $\textit{C}_{\!_\textit{L}}(t) < \overline {\textit{C}_{\!_\textit{L}}}$. It is implemented by taking the mean of the lift coefficient at each angle of attack. At the beginning, $\widehat {\boldsymbol{\mathsf{U}}_1}(x,y)$ and are initialized with zeros. If $\textit{C}_{\!_\textit{L}}(n) > \overline {\textit{C}_{\!_\textit{L}}}$, then $\boldsymbol{\mathsf{u}}_1(x,y, n)$ is added to $\widehat {\boldsymbol{\mathsf{U}}_1}(x,y)$ and a counter $n_a$ is increased by one. If $\textit{C}_{\!_\textit{L}}(n) < \overline {\textit{C}_{\!_\textit{L}}}$, then $\boldsymbol{\mathsf{u}}_1(x,y, n)$ is added to and a counter $n_b$ is increased by one. The if condition statement is run for all of the snapshots of the data, and the final matrices $\widehat {\boldsymbol{\mathsf{U}}_1}(x,y)$ and are divided by the final integers $n_a$ and $n_b$, respectively. The conditional time-averaging is synchronized for all flow variables so that at a given data point ($n$), all flow variables are allocated to either the high-lift or low-lift regimes. In this sense, one can define a time-averaged surplus ($\Delta {\overline {\boldsymbol{\mathsf{U}}_1}}^ +$) and deficit ($\Delta {\overline {\boldsymbol{\mathsf{U}}_1}}^{-}$) of the streamwise velocity above and below that of the mean streamwise velocity, respectively. Here $\Delta {\overline {\boldsymbol{\mathsf{U}}_1}}^+ = \widehat {\boldsymbol{\mathsf{U}}_1} - \overline {\boldsymbol{\mathsf{U}}_1}$ and :

(2.1)

where $M$, $N$ and $L$ are the number of data points of $\textit{C}_{\!_\textit{L}}$, data points of $\textit{C}_{\!_\textit{L}}$ that have values larger than $\overline {\textit{C}_{\!_\textit{L}}}$ and data points of $\textit{C}_{\!_\textit{L}}$ that have values less than $\overline {\textit{C}_{\!_\textit{L}}}$, respectively. Also, $\mathbb {S}=\{n|n\in \mathbb {N}, n = 1, 2, 3,\ldots, M\}$; $\mathbb {S}_a=\{n|n\in \mathbb {N}, n = 1, 2, 3,\ldots, N; \textit{C}_{\!_\textit{L}}(n) > \overline {\textit{C}_{\!_\textit{L}}}\}$; $\mathbb {S}_b=\{n|n\in \mathbb {N}, n = 1, 2, 3,\ldots, L; \textit{C}_{\!_\textit{L}}(n) < \overline {\textit{C}_{\!_\textit{L}}}\}$. The surplus and deficit of the streamwise velocity are mirror symmetric. Thus, the spatial streamwise velocity mode that is causing fluctuations in the lift coefficient can be represented by the average of the surplus and the deficit. That is, $\overline {\boldsymbol{\mathsf{U}}_1}^+ = \frac {1}{2}(\Delta {\overline {\boldsymbol{\mathsf{U}}_1}}^+ - \Delta {\overline {\boldsymbol{\mathsf{U}}_1}}^{-})$. In addition to the streamwise velocity, this process can be done to the wall-normal velocity and the pressure. Consequently, the spatial SMD mode that is causing fluctuations in the lift coefficient $\boldsymbol {\tilde {\varPsi }}_{\textit{C}_{\!_\textit{L}}}$ is given by

(2.2)\begin{equation} \boldsymbol{\tilde{\varPsi}}_{\textit{C}_{\!_\textit{L}}}=\left[ \begin{array}{c} \overline{\boldsymbol{\mathsf{U}}_1}^+ \\ \overline{\boldsymbol{\mathsf{U}}_2}^+ \\ \overline{\boldsymbol{\mathsf{P}}}^+ \\ \end{array} \right]. \end{equation}

The spatial SMD mode corresponding to the lift coefficient is orthogonal but not normalized and can be normalized to obtain the orthonormal spatial SMD mode $\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{L}}}$ using the following formula:

(2.3)\begin{equation} \boldsymbol{\varPsi}_{\textit{C}_{\!_\textit{L}}}(x,y) = \frac{\boldsymbol{\tilde{\varPsi}}_{\textit{C}_{\!_\textit{L}}}(x,y) }{\left[ \left(\boldsymbol{\tilde{\varPsi}}_{\textit{C}_{\!_\textit{L}}}(x,y)\right) \left(\boldsymbol{\tilde{\varPsi}}_{\textit{C}_{\!_\textit{L}}}(x,y)\right)^{\intercal} \right]^{{1}/{2}}}. \end{equation}

The SMD coefficient corresponding to the lift coefficient ($\textit {a}_{\textit{C}_{\!_\textit{L}}}(t)$) can be determined using

(2.4)\begin{equation} \textit{a}_{\textit{C}_{\!_\textit{L}}}(t)=\int_{Domain}\boldsymbol{\mathsf{D}}(x,y,t)\boldsymbol{\varPsi}_{\textit{C}_{\!_\textit{L}}}(x,y) \,{\rm d}\kern0.7pt x\,{\rm d}y, \end{equation}

where $\boldsymbol{\mathsf{D}}(x,y,t)$ is the spanwise ensemble-averaged fluctuating flow. The fluctuating flow field corresponding to the lift coefficient ($\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{L}}}$) can be reconstructed from

(2.5)\begin{equation} \boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{L}}}(x,y,t)=\textit{a}_{\textit{C}_{\!_\textit{L}}}(t)\boldsymbol{\varPsi}_{\textit{C}_{\!_\textit{L}}}(x,y). \end{equation}

The SMD mode causing the fluctuations in the drag coefficient $\boldsymbol {\tilde {\varPsi }}_{\textit{C}_{\!_\textit{D}}}$ can be reconstructed similarly. However, the fluctuating flow field corresponding to the lift coefficient must first be subtracted from the fluctuating flow as follows:

(2.6)\begin{equation} \dot{\boldsymbol{\mathsf{D}}}(x,y,t)=\boldsymbol{\mathsf{D}}(x,y,t)-\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{L}}}(x,y,t). \end{equation}

The orthonormal spatial SMD mode $\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{D}}}$, the SMD coefficient $\textit {a}_{\textit{C}_{\!_\textit{D}}}(t)$ and the fluctuating flow field $\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{D}}}$ corresponding to the drag coefficient can be extracted from the fluctuating flow field $\dot {\boldsymbol{\mathsf{D}}}(x,y,t)$ using (2.2), (2.4) and (2.5). This process can be repeated until all flow modes featuring instabilities that are fixed in space and grow in time (absolute instabilities) are recovered. This method cannot be used to recover flow modes that feature convecting instabilities like the travelling Kelvin–Helmholtz waves along the wake of the airfoil that represents a high-frequency oscillating (HFO) mode. However, the spectrum of the lift coefficient or time history of the wall-normal velocity at a probe downstream of the trailing edge can be used to determine the frequency of the most dominant HFO mode. Consequently, a phase average of the flow at the selected frequency would recover the most dominant HFO mode.

Figure 3. The conditional time-averaging process. Time history of the lift coefficient at ${\textit {Re}}_c = 9\times 10^4$ and angle of attack of $11.0^{\circ }$. The dashed, solid and dash-dotted horizontal lines show the high-lift, mean-lift and low-lift time average, respectively.

The objective of presenting the SMD method in this work is to confirm the findings of the POD method and identify which flow mode causes fluctuations in the aerodynamic forces. However, we ended up with a simple and easy-to-implement method to recover selective flow modes, as will be seen later. Most importantly, the SMD method is based on time-averaging of the data. The method is not general to all laminar, transitional and turbulent flows, but it can be implemented for similar flow fields where the effect of the dynamics of the flow is integrated on a solid surface and induces fluctuations in aerodynamic forces. Fluctuations in the drag coefficient of the flow field about a cylinder can be used as a basis function to recover the coherent structures that induce these fluctuations in the drag coefficient, time history of the side force on a cylinder body at high angles of attack can be used to recover the flow mode that causes the side force and the time history of the oscillation of a shockwave interacting with a boundary layer can be used to reconstruct the flow mode that causes these oscillations in the shockwave. The SMD method is not computationally demanding especially for three-dimensional configurations, can be implemented on the fly and the basis functions change automatically with their associated flow modes. Thus, the basis functions are not affected by applying flow control or changing the flow artificially.

3. Computational set-up and validations

The LES code utilized in the present simulations is an LES version (Almutairi 2010) of the direct numerical simulation code written and validated by Jones (2008). The Navier–Stokes equations were discretized using a fourth-order explicit central difference scheme for spatial discretization in the interior points. The fourth-order boundary scheme of Carpenter, Nordström & Gottlieb (1999) was used to treat points near and at the boundary. To preserve the spatial characteristics, the transformation metrics tensor $\grave {\boldsymbol {\xi }_{ij}}$ was evaluated using the same fourth-order scheme. Temporal discretization was performed using a low-storage fourth-order Runge–Kutta scheme. The solution stability was improved by implementing an entropy splitting scheme Sandham, Li & Yee (2002). The entropy splitting constant $\beta$ was set equal to $2.0$ (Jones 2008). The integral characteristic boundary condition was applied at the free-stream and the far-field boundaries (Sandhu & Sandham 1994). The zonal characteristic boundary condition was applied at the downstream exit boundary (Sandberg & Sandham 2006) which is considered as a non-reflected boundary condition to overcome the circulation effects at the boundaries. The adiabatic and no-slip conditions were applied on the airfoil surface. The LSB, its associated LFO and the dynamics of the flow are inherently statistically two-dimensional. Therefore, a periodic boundary condition was applied in the spanwise direction for each step of the Runge–Kutta time steps. The internal branch-cut boundary was updated at each step of the Runge–Kutta scheme. In terms of the airfoil chord $(c)$, the dimensions of the computational domain were set as follows: $L_{\xi }=5c$ in the wake region (from the airfoil trailing edge to the outflow boundary in the streamwise direction), $L_\eta =7.3c$ in the wall-normal direction (the C-grid radius) and $L_{\zeta }=0.5c$ in the spanwise direction, where $\xi$, $\eta$ and $\zeta$ are the curvilinear coordinates. The LSB formation, elongation and bursting processes should be well resolved; therefore, the grid resolution near the airfoil surface, especially on the suction side, is critical. Three different C-grids were constructed for each Reynolds number with various distributions in $\xi$, $\eta$ and $\zeta$ directions as shown in table 1. Here $\Delta x^+$, $y^+$ and $\Delta z^+$ are the grid resolutions in the $\xi$, $\eta$ and $\zeta$ directions, respectively; and $N_\xi$, $N_\eta$ and $N_\zeta$ are the grid points in the $\xi$, $\eta$ and $\zeta$ directions, respectively.

Table 1. Computational grid parameters.

Large-eddy simulations were carried out for the flow around a NACA0012 airfoil at Mach number of $0.4$, Reynolds number of $5\times 10^4$ and $9\times 10^4$ and several angles of attack near stall. Reducing the free-stream Mach number of the numerical simulation to simulate the flow field at a near-incompressible regime (${M}_\infty < 0.2$) would considerably increase the computational cost. The entire domain was initialized using the free-stream conditions ($\rho _{\infty }=1$, $\rho _{\infty }\textit{U}_{\infty }=1$, $\rho _{\infty }\textit{V}_{\infty }=0$, $\rho _{\infty }\textit{W}_{\infty }=0$ and $\textit{T}_{\infty }=1$). The simulations were performed with a time step of $10^{-4}$ non-dimensional time units. Samples for statistics were collected once the transition of the simulations had decayed and the flow became stationary in time after $50$ non-dimensional time units. Aerodynamic coefficients were sampled for each angle of attack at a frequency of 10 000 to generate two-and-a-half million samples over the time period of $250$ non-dimensional time units. The locally time-averaged and spanwise ensemble-averaged pressure, velocity components and Reynolds stresses were sampled every $50$ time steps on the $x$–$y$ plane. A dataset of 20 000 $x$–$y$ planes was recorded at a frequency of $204$ at each angle of attack. A grid sensitivity study was first carried out at two representative angles of attack and it was concluded that Grid-$2$ gives the best grid distributions as seen in figure 4.

Figure 4. Sensitivity of the (a,c) mean pressure and (b,d) skin-friction coefficients to grid refinement at two representative angles of attack: $9.25^{\circ }$ at a Reynolds number of $5\times 10^4$ and $11.0^{\circ }$ at a Reynolds number of $9\times 10^4$.

Figure 5 shows a comparison of the LES data of the present work with the LES data of Alferez et al. (2013) and Alferez (2014). The figure displays a comparison of the predicted variations of the pressure coefficient, the skin-friction coefficient and the mean streamwise velocity profiles at sixteen $x_s/c$ locations starting at a location upstream of the separation point (S) to a location downstream of the reattachment point (R). As seen in the figure, the LES data at ${\textit {Re}}_c = 9\times 10^4$, ${M}_{\infty }=0.4$ and stall angle of attack of $10.25^{\circ }$ compare very well with the LES data of Alferez (2014) at ${\textit {Re}}_c = 10^5$, ${M}_{\infty }=0.16$ and angle of attack of $10.55^{\circ }$ which is slightly lower than the stall angle of attack ($10.55^{\circ } \leqslant \alpha _s \leqslant 10.8^{\circ }$). The discrepancy between the LES data of the present work and that of Alferez (2014) is due to the significant difference in the free-stream Mach number and consequently the effect of compressibility. The effects of compressibility are expressed by the shift in the magnitude of the pressure and skin-friction coefficients and the profiles of the streamwise velocity in the separated region of the flow. Recently, Benton & Visbal (2020) showed that increased compressibility causes an earlier inception of dynamic stall as a consequence of the bursting of the LSB. Thus, similar effects of compressibility in the case of static stall affect the transition process of the separated shear layer, the characteristics of the LSB and the development of the boundary layer. The reader is referred to Eljack (2017), Elawad & Eljack (2019), Eljack & Soria (2020) and Eljack et al. (2021) for more details on the mathematical modelling, computational set-up and validation of the data.

Figure 5. Profiles of the (a) mean pressure coefficient ($\textit{C}_{\!_\textit{P}}$), (b) skin-friction coefficient ($\textit{C}_{\!_\textit{f}}$) and (c) mean streamwise velocity scaled by the local external velocity ($\textit{U}_1/\textit{U}_e$) plotted versus the distance from the airfoil leading edge computed along the curvilinear coordinate on the airfoil suction side (${x_s}/c$) with ${x_s}/c = 0$ at the stagnation point. Here $y_n/c$ is the vertical distance measured from the airfoil surface. The separation (S), transition (T) and reattachment (R) locations are indicated by the filled black circles (a,b) and the dash-dotted vertical lines (c).

4. Application of the POD, the DMD and the SMD methods

The locally time-averaged and spanwise ensemble-averaged data were utilized to implement the snapshot POD, the DMD and the SMD methods. The locally time-averaged and spanwise ensemble-averaged fluctuating streamwise velocity, wall-normal velocity and pressure were used to estimate the two-point correlation in time using 20 000 data points that span $100$ non-dimensional time units or four low-frequency cycles. The eigenvalue problem ($\boldsymbol{\mathsf{C}}\boldsymbol {\varPhi }=\boldsymbol {\varLambda }\boldsymbol {\varPhi }$) was solved, and the POD eigenvalues and eigenvectors were obtained. Here $\boldsymbol {C}$ represents the correlation matrix in time. The locally time-averaged and spanwise ensemble-averaged data were formulated into a two-dimensional matrix to implement the DMD method. The rows contain the streamwise velocity component ($\boldsymbol {u_1}(x, y, t)$), the wall-normal velocity component ($\boldsymbol {u_2}(x, y, t)$), the pressure ($\kern 1.5pt\boldsymbol {p}(x, y, t)$), the lift coefficient ($\textit{C}_{\!_\textit{L}}(t)$) and the drag coefficient ($\textit{C}_{\!_\textit{D}}(t)$). The columns contain the temporal variation of these variables in time that span 20 000 data points, $100$ non-dimensional time units or four low-frequency cycles. The companion matrix ($\boldsymbol {S}$) was then formulated, and the eigenvalue problem was solved for the eigenvalues ($\boldsymbol {\varGamma }$) and the eigenvectors ($\boldsymbol {\varOmega }$). After that, the growth rates and phase velocities were calculated, and the dynamic modes were constructed. All the utilized flow variables have the same eigenvectors (dynamic modes shapes) and eigenvalues (growth rates and phase velocities). However, the amplitude of each dynamic mode at different frequencies depends on the flow variable. Thus, even though all the analysed flow variables share the same eigenvalues and eigenvectors, each flow variable has its respective spectrum.

The lift coefficient was utilized as a reference to conditionally time-average the streamwise velocity, wall-normal velocity and pressure to determine the SMD mode corresponding to the lift coefficient using (2.2), (2.4) and (2.5). The spatial SMD mode corresponding to the drag coefficient was then calculated as follows:

(4.1)\begin{equation} \boldsymbol{\tilde{\varPsi}}_{\textit{C}_{\!_\textit{D}}}=\left[ \begin{array}{c} {}\overline{\dot{\boldsymbol{\mathsf{U}}}_1}^+ \\ {}\overline{\dot{\boldsymbol{\mathsf{U}}}_2}^+ \\ {}\overline{\dot{\boldsymbol{\mathsf{P}}}}^+ \end{array} \right], \end{equation}

where $\dot {\boldsymbol{\mathsf{D}}}$ was calculated using (2.6). Mode $\boldsymbol {\tilde {\varPsi }}_{\textit{C}_{\!_\textit{D}}}$ was normalized using (2.3) to obtain the orthonormal spatial SMD mode ($\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{D}}}$), and the SMD coefficient corresponding to the drag coefficient (${a}_{\textit{C}_{\!_\textit{D}}}(t)$) was then estimated from

(4.2)\begin{equation} {a}_{\textit{C}_{\!_\textit{D}}}(t)=\int_{Domain}\dot{\boldsymbol{\mathsf{D}}}(x,y,t)\boldsymbol{\varPsi}_{\textit{C}_{\!_\textit{D}}} (x,y)\,{\rm d}\kern0.7pt x\,{\rm d}y. \end{equation}

The fluctuating flow corresponding to the drag coefficient, $\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{D}}}$, was reconstructed from

(4.3)\begin{equation} \boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{D}}}(x,y,t)=\textit{a}_{\textit{C}_{\!_\textit{D}}}(t)\boldsymbol{\varPsi}_{\textit{C}_{\!_\textit{D}}}(x,y). \end{equation}

The fluctuating flow field corresponding to the drag coefficient was then subtracted to obtain the remaining flow field ($\ddot {\boldsymbol{\mathsf{D}}}(x,y,t)$) as follows:

(4.4)\begin{equation} \ddot{\boldsymbol{\mathsf{D}}}(x,y,t)=\boldsymbol{\mathsf{D}}(x,y,t)-\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{L}}}(x,y,t)-\boldsymbol{\mathsf{D}}_{\textit{C}_{\!_\textit{D}}}(x,y,t). \end{equation}

Once the flow modes that induce oscillations in the lift and the drag are subtracted from the flow field, the remaining flow consists of high-frequency flow modes with various frequencies and amplitudes that represent the evolution of the flow field away from the wall either along the shear layer (Kelvin–Helmholtz instability) or along the wake of the airfoil (travelling Kelvin–Helmholtz waves). We are not interested in all of the flow modes; thus, the leading HFO mode can be recovered using its frequency. The spectrum of the lift coefficient was used to determine the frequency of the most dominant HFO mode. Consequently, the fluctuating flow field ($\ddot {\boldsymbol{\mathsf{D}}}(x,y,t)$) was phase-averaged at the dominant frequency to recover the leading HFO mode. Previous work on airfoils at near-stall conditions has shown that the most dominant high-frequency mode oscillates at a frequency of about $1.0\ {\rm Hz}$. Thus, a data record of $100$ non-dimensional units would constitute about $100$ cycles of oscillations and ensures the statistical accuracy of the phase-averaging process. The DMD spectra have shown that the HFO mode is driven by the wall-normal velocity. Thus, at angles of attack where the lift coefficient spectrum is not smooth enough to determine the most dominant frequency, the spectrum of the time history of the wall-normal velocity sampled at a probe located downstream of the airfoil trailing edge was used to determine the dominant frequency of the HFO mode.

4.1. The POD modes

The POD eigenvalues were utilized to estimate the cumulative energy content using (B7). Figure 6 shows the cumulative energy plotted versus the number of POD modes used to estimate it. The range of angles of attack is $\alpha =9.25^{\circ }$–$9.8^{\circ }$ for ${\textit {Re}}_c = 5\times 10^4$ and $\alpha =10.25^{\circ }$–$10.8^{\circ }$ for ${\textit {Re}}_c = 9\times 10^4$ on the top part of the figure, and $\alpha =9.8^{\circ }$–$10.5^{\circ }$ and $\alpha =10.8^{\circ }$–$11.2^{\circ }$ on the bottom part of the figure. The angles of attack of $9.8^{\circ }$ at ${\textit {Re}}_c=5\times 10^4$ and $10.8^{\circ }$ at ${\textit {Re}}_c=9\times 10^4$ are duplicated in both panels because they are the angles at which the POD modes have the fastest convergence to the total energy; that is, the minimum number of POD modes used to attain more than $99\,\%$ of the total energy. Thus, they are displayed in both panels to compare the convergence of POD modes at other angles of attack. These are the angles of attack at which the LFO has the maximum amplitude of oscillations (Eljack et al. 2021). At the stall angle of attack, the cumulative energy converges slowly towards $1$ as the number of POD modes approaches $250$, as seen in the figure. The cumulative energy converges faster as the angle of attack increases. The fastest cumulative energy convergence was achieved at the angles of attack of maximum LFO. The convergence process of the POD energy slows down as the angle of attack increases above the angle of attack of maximum LFO.

Figure 6. Cumulative energy of the POD modes ($\kappa _t(m)$) estimated by summing the POD eigenvalues over $m$ POD modes. The arrows indicate the direction in which the angle of attack $\alpha$ increases.

4.2. The DMD modes

The DMD spectra were estimated for the lift coefficient, the drag coefficient, the pressure, the streamwise velocity and the wall-normal velocity. The DMD spectra of all flow variables were used to identify the most dominant DMD flow modes. There is one low-frequency peak in the spectrum of each of the analysed flow variables featuring the low-frequency oscillating mode one (LFO-Mode-1) for all of the investigated angles of attack. However, the spectrum of the drag coefficient exhibits a low-frequency peak at a different Strouhal number from that of the LFO-Mode-1. Thus, the DMD flow mode that dominates the drag spectrum represents a second low-frequency oscillating mode (LFO-Mode-2). The LFO-Mode-1 and the LFO-Mode-2 exchange the dominance of the spectra of the lift coefficient, the drag coefficient, the pressure and the streamwise velocity. There is no significant high-frequency peak in any of the spectra of these variables. However, there is a pronounced high-frequency peak in the spectrum of the wall-normal velocity featuring the HFO mode. Hence, the DMD spectra of the lift coefficient, the drag coefficient and the wall-normal velocity were used to identify the LFO-Mode-1, the LFO-Mode-2 and the HFO mode, respectively. Once each of the three modes was identified, its frequency was used to locate its corresponding dynamic mode in the spectrum of the other flow variables. Thus, the DMD spectra of the flow variables show three dominant modes. However, each of the three dominant modes is only dominant in the DMD spectrum of the variable used to identify it.

Figure 7 shows the DMD spectra for the five variables. The left-hand side of the figure shows the spectra for angles of attack of $9.25^{\circ }$–$10.5^{\circ }$ at Reynolds number of $5\times 10^4$ and the right-hand side shows the spectra for angles of attack of $10.25^{\circ }$–$11.2^{\circ }$ at Reynolds number of $9\times 10^4$. The points in the spectra representing the LFO-Mode-1, the LFO-Mode-2 and the HFO are indicated by filled black circles, filled red circles and filled blue circles, respectively. The LFO-Mode-1 and the LFO-Mode-2 dominate the spectra of the lift coefficient and the drag coefficient, respectively, for all of the investigated angles of attack and Reynolds numbers of $5\times 10^4$ and $9\times 10^4$. At Reynolds number of $5\times 10^4$, the LFO-Mode-1 dominates the spectra of the pressure, the streamwise velocity and the wall-normal velocity at angles of attack higher than the angle of attack of $9.5^{\circ }$ and lower than the angle of attack of $10.1^{\circ }$. It is worth mentioning that the LFO phenomenon initiates at angles of attack higher than the angle of attack of $9.25^{\circ }$, becomes fully developed at the angle of attack of $9.8^{\circ }$ and loses momentum at angles of attack higher than the angle of attack of $10.0^{\circ }$. The LFO-Mode-1 dominates the DMD spectra of the pressure and streamwise velocity at all of the investigated angles of attack at Reynolds number of $9\times 10^4$.

Figure 7. DMD spectra of (a) the lift coefficient ($\textit{C}_{\!_\textit{L}}$), (b) the drag coefficient ($\textit{C}_{\!_\textit{D}}$), (c) the pressure ($\kern 1.5pt p$), (d) the streamwise velocity ($u_1$) and (e) the wall-normal velocity ($u_2$). The filled black circles denote the LFO-Mode-1, the filled red circles display the LFO-Mode-2 and the filled blue circles indicate the HFO mode.

The spectra of the wall-normal velocity are dominated by one low-frequency mode and a broad band of high-frequency modes for each of the investigated angles of attack. At Reynolds number of $5\times 10^4$ and angle of attack of $9.25^{\circ }$ the spectrum is dominated by the HFO mode as seen in the left-hand side of the figure. As the angle of attack increases above $9.25^{\circ }$, a sudden and drastic change occurs, and the dominant mode shifts to the LFO-Mode-1 and the LFO-Mode-2. The two low-frequency modes dominate the spectra of the wall-normal velocity component until the angle of attack is raised above $10.0^{\circ }$; then, the HFO mode dominates the spectra again, as seen in the figure. At Reynolds number of $9\times 10^4$, the HFO mode dominates the spectra at angles of attack lower than or equal to the angle of attack of $10.6^{\circ }$. Whereas, the LFO-Mode-1 dominates the spectra at angles of attack higher than or equal to the angle of attack of $10.8^{\circ }$. It is worth mentioning that the LFO phenomenon is pronounced in the amplitude of oscillation and uniform in its lifecycle at angles of attack higher than or equal to $\alpha = 10.8^{\circ }$. In summary, the LFO-Mode-1 and the LFO-Mode-2 (absolute instability) dominate the spectrum when the LFO phenomenon exists in the flow, and the HFO mode (convective instability) dominates the spectrum when the LFO does not exist in the flow. Thus, the LFO phenomenon occurs in the flow when there is absolute instability. The HFO mode has an insignificant magnitude in the DMD spectrum of the lift coefficient at all of the investigated angles of attack. This is indicative that the HFO mode does not contribute directly to the oscillations in the lift coefficient. It is worth noting that the DMD analysis decomposes the fluctuating flow into many low-frequency modes as seen in figure 7. Whereas, the POD method recovers only two low-frequency modes that reconstruct the fluctuating flow favourably as will be discussed later.

Figure 8 shows plots of the Strouhal number of the most dominant flow modes versus the angle of attack. The profile of the Strouhal number of the LFO-Mode-1 for the lift coefficient, the pressure and the streamwise velocity is similar to that obtained using the fast Fourier transform algorithm by Eljack et al. (2021). The wall-normal velocity component is dominated by a high-frequency flow mode at angles of attack of $\alpha \leqslant 9.25^{\circ }$ and $\alpha \geqslant 10.1^{\circ }$, and dominated by a low-frequency mode at angles of attack of $9.4^{\circ } \leqslant \alpha \leqslant 10.0^{\circ }$ for Reynolds number of $5\times 10^4$. At Reynolds number of $9\times 10^4$, the HFO dominates at angles of attack lower than or equal to $10.6^{\circ }$. The vertical axis is broken into two ranges in the plot for the wall-normal velocity component to capture both the low-frequency and the high-frequency ranges that dominate it.

Figure 8. Strouhal number (St) of the most dominant flow modes versus the angle of attack and Reynolds number for the LFO-Mode-1 (left) and the HFO mode (right).

4.3. The SMD modes

Time history of the lift coefficient ($\textit{C}_{\!_\textit{L}}(t)$) was utilized to obtain the orthonormal spatial SMD mode ($\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{L}}}$) using (2.2) and (2.3). The flow field was then projected onto $\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{L}}}$ to obtain the SMD coefficient corresponding to the lift coefficient (${a}_{\textit{C}_{\!_\textit{L}}}(t)$) using (2.4). After that, the orthonormal spatial SMD mode ($\boldsymbol {\varPsi }_{\textit{C}_{\!_\textit{D}}}$) and the SMD coefficient corresponding to the drag coefficient (${a}_{\textit{C}_{\!_\textit{D}}}(t)$) were duly estimated. The fluctuating flow along the wake of the airfoil mimics the oscillating behaviour of the wall-normal velocity. Thus, the wall-normal velocity component was probed at four locations downstream of the airfoil trailing edge, and the most dominant frequency was determined. The value of the frequency was used to estimate the rate at which the HFO mode repeats periodically; consequently, phase-averaging of the fluctuating flow at this frequency was performed to recover the spatial HFO mode. Equation (2.3) was then used to obtain the orthonormal spatial SMD mode ($\boldsymbol {\varPsi }_{v}$), where the subscript $v$ denotes the wall-normal velocity used to extract the mode. After that, the fluctuating flow field was projected onto $\boldsymbol {\varPsi }_{v}$ to obtain the SMD coefficient corresponding to the HFO mode (${a}_{v}(t)$). The SMD modes corresponding to the lift and the drag coefficients represent the LFO-Mode-1 and LFO-Mode-2, respectively. Figure 9 shows that the SMD modes corresponding to the lift coefficient, the drag coefficient and the wall-normal velocity compare very well with their corresponding POD modes. The LFO phenomenon is fully developed at an angle of attack of $9.8^{\circ }$ and Reynolds number of $5\times 10^4$, and an angle of attack of $10.8^{\circ }$ and Reynolds number of $9\times 10^4$. Thus, the LFO repeats periodically with some disturbed cycles. However, at lower angles of attack, the LFO becomes stochastic, and the POD and SMD coefficients do not compare with each other.

Figure 9. Comparison of the POD coefficients (blue line) and the SMD coefficients (red line) for (a,b) the LFO-Mode-1, (c,d) the LFO-Mode-2 and (e,f) the HFO mode.

4.4. The spatial LFO-Mode-1

Figure 10 shows the DMD spectrum of the lift coefficient and the DMD spatial mode corresponding to the LFO-Mode-1. The DMD spectrum shows the growth rates of the LFO-Mode-1 and the LFO-Mode-2 indicated by the filled black and red circles, respectively. The growth rates of the LFO-Mode-1 and the LFO-Mode-2 are both positive at this angle of attack, indicating that the two flow modes are growing. At other angles of attack, the growth rates of the LFO-Mode-1 and the LFO-Mode-2 have negative or positive signs, indicating that the LFO-Mode-1 and LFO-Mode-2 are decaying or growing, respectively. The LFO is quasi-periodic; thus, the growth rates for the leading DMD modes should be zero. However, the data from which the DMD spectrum was extracted spans four low-frequency cycles. Should a longer data record be used, a zero growth rate would have been obtained. The size of the circles denotes the relative amplitude of each dynamic mode. The relative magnitude of the LFO-Mode-1 and the LFO-Mode-2 is proportional to the angle of attack and reaches its maximum amplitude at the angle of attack of maximum LFO. At angle of attack of $9.8^{\circ }$ at Reynolds number of $5\times 10^4$ and angle of attack of $11.0^{\circ }$ at Reynolds number of $9\times 10^4$, the LFO-Mode-1 and LFO-Mode-2 contain more than $50\,\%$ of the turbulent kinetic energy and the remaining kinetic energy is almost equally distributed among higher modes as seen in figure 6. Therefore, only a few low-frequency modes are shown in figure 10, and the higher modes become insignificant. This is indicative that the LFO process is fully developed at these angles of attack, and the oscillating flow is more pronounced in magnitude and more coherent in shape. The relative amplitude of the LFO-Mode-1 and the LFO-Mode-2 becomes inversely proportional to the angle of attack when the angle increases above the angle of attack of maximum LFO. The figure shows only the DMD spectrum of the lift coefficient; however, the DMD spectra of the drag, the pressure, the streamwise velocity and the wall-normal velocity are similar to that of the lift coefficient. The LFO-Mode-1 features the triad of vortices as seen in figure 10(b). The triad of vortices is driven by the oscillating streamwise velocity component across the laminar portion of the separated shear layer. The magnitude of the oscillating velocity increases as the angle of attack increases (Eljack & Soria 2020).

Figure 10. The LFO-Mode-1 for angle of attack of $9.8^{\circ }$ and Reynolds number of $5\times 10^4$. (a) The DMD spectrum of the lift coefficient, where $\omega _r$ is the growth rate and $St$ is the Strouhal number. (b) Streamline patterns superimposed on colour maps of the spanwise vorticity $\omega _z$ for the DMD spatial mode.

Figure 11 shows the spatial LFO-Mode-1 constructed using the POD, the DMD and the SMD methods. The POD spatial mode for the LFO-Mode-1 was constructed by multiplying its orthonormal spatial POD mode by the average amplitude of its corresponding coefficient $\overline {|b^{(1)}(t)|}$. The POD method recovers the flow modes based on their energy content. Thus, the spatial LFO-Mode-1 constructed using the POD method is contaminated by other high-frequency flow modes. Hence, the shape and evolution of the triad of vortices are obscured by high-frequency flow modes. However, POD captured the general features of this flow mode. In the DMD method, the LFO-Mode-1 is the dynamic mode corresponding to the low-frequency peak in the DMD spectrum of the lift coefficient. Thus, the LFO-Mode-1 represents the instability that features globally oscillating flow around the airfoil, the oscillating pressure along the airfoil chord and the process that creates and sustains the triad of vortices. Hence, the triad of vortices and the general features of this flow mode are accurately captured by the DMD method. However, the LFO-Mode-1 is presented by several dynamic modes; whereas, it was objectively recovered by a single POD mode. Equation (2.2) was used to calculate the pressure, the streamwise velocity and the wall-normal velocity corresponding to the spatial SMD mode of the LFO-Mode-1 based on the time history of the lift coefficient. The zoom-in box shows the triad of vortices constructed using the SMD method. The lift coefficient is dominated by a low frequency corresponding to the LFO-Mode-1. Consequently, the spatial LFO-Mode-1 constructed using the SMD method is similar to that constructed using the DMD method and constitutes a coherent triad of vortices. Furthermore, the SMD method effectively recovered the LFO-Mode-1 in a single flow mode. Thus, the SMD method combines the strengths of the POD and the DMD methods by objectively recovering the flow mode and capturing the instability. As seen in the figure, the LFO-Mode-1 is a global flow mode that drives and sustains the triad of vortices. The spatial LFO-Mode-1 creates an adverse oscillating pressure gradient when the sign of its amplitude is positive, and a favourable oscillating pressure gradient when the sign of its amplitude is negative (see figure 9). The LFO-Mode-1 preserves shape over all of the investigated angles of attack; however, the magnitude of oscillation increases as the angle of attack increases.

Figure 11. Streamline patterns superimposed on colour maps of the oscillating pressure field for the spatial LFO-Mode-1 constructed using the POD, the DMD and the SMD methods: (a,c,e) $\alpha =9.8^{\circ }$ and ${\textit {Re}}_c = 5\times 10^4$; (b,d,f) $\alpha =11.0^{\circ }$ and ${\textit {Re}}_c = 9\times 10^4$.

4.5. The spatial LFO-Mode-2

The spatial LFO-Mode-2 constructed using the POD, the DMD and the SMD methods is shown in figure 12. The POD mode corresponding to the LFO-Mode-2 was constructed by multiplying its orthonormal spatial POD mode by the average amplitude of its corresponding POD coefficient $\overline {|b^{(2)}(t)|}$. The time history of the drag coefficient was used to estimate the pressure, streamwise velocity and wall-normal velocity corresponding to the SMD construction of the LFO-Mode-2. The LFO-Mode-2 originates and evolves on the suction surface of the airfoil and features the expansion and advection of the upstream vortex of the TCV (Eljack & Soria 2020). The LFO-Mode-2 constructed using the POD and the SMD methods are similar in shape and magnitude as seen in the figure. The POD and the SMD analyses extracted only two low-frequency modes. Whereas, the DMD analysis decomposed the low-frequency oscillating flow into many low-frequency modes. The low-frequency modes extracted by the POD and the SMD methods are more coherent and energetic compared with those extracted by the DMD method. Furthermore, the DMD and the SMD methods accurately captured the instability. Thus, combining the three methods provided all the information that governs and controls the dynamics of the flow.

Figure 12. Streamline patterns superimposed on colour maps of the oscillating pressure field for the spatial LFO-Mode-2 constructed using the POD, the DMD and the SMD methods: (a,c,e) $\alpha =9.8^{\circ }$ and ${\textit {Re}}_c = 5\times 10^4$; (b,d,f) $\alpha =11.0^{\circ }$ and ${\textit {Re}}_c = 9\times 10^4$.

4.6. The spatial HFO mode

Figure 13 shows the spatial HFO mode constructed using the POD, the DMD and the SMD methods. As seen in the figure, the HFO mode features large oscillations that originate at the leading edge and advect towards the trailing edge. At the trailing edge, the HFO mode interacts with the local flow instability, stretched by the strong shear, and energizes the vortex shedding. The magnitude of the oscillation is proportional to the angle of attack. It is interesting to note that there is no high-frequency peak in the spectra of the lift and the drag coefficients. The aerodynamic coefficients are direct consequences of variations in the flow field adjacent to the airfoil surface. Furthermore, the pressure field communicates the dynamics of the flow away from the wall to the flow field adjacent to the wall. However, when low-frequency modes of significant amplitude are present in the flow, the high-frequency modes do not much affect the flow variables near the wall and will not be reflected in the aerodynamic forces. Consequently, such an important flow feature will not show up as a peak in the spectra of the aerodynamic coefficients unless enough blocks of data are used and the spectra are smoothed out.

Figure 13. Streamline patterns superimposed on colour maps of the oscillating wall-normal velocity field for the spatial HFO mode constructed using the POD, the DMD and the SMD methods: (a,c,e) $\alpha =9.8^{\circ }$ and ${\textit {Re}}_c = 5\times 10^4$; (b,d,f) $\alpha =11.0^{\circ }$ and ${\textit {Re}}_c = 9\times 10^4$.

Figure 7 shows that the spectra of the wall-normal velocity peak significantly at high frequency. Whereas, the spectra of the streamwise velocity and the pressure show no significant peak at high frequency. This is indicative that the HFO is primarily driven by the wall-normal velocity component. However, the streamwise velocity and the pressure indirectly affect the behaviour of the HFO mode. Furthermore, the spectra of the wall-normal velocity are interchangeably dominated by the HFO mode and the two low-frequency modes. Thus, the two low-frequency modes and the HFO mode are interlinked and play a profound role in the dynamics of the flow. Copious amounts of research work have concentrated on the characteristics of the HFO mode downstream of the airfoil trailing edge and the evolution of the Kelvin–Helmholtz instability in the vicinity of the leading edge (Jones, Sandberg & Sandham 2008, 2010; Almutairi, Eljack & Alqadi 2017). Most of those studies found that there are many aspects of similarities and concluded that the HFO mode generates acoustic waves that travel upstream, energizes the Kelvin–Helmholtz instability in the vicinity of the leading edge and forces early transition. However, the HFO mode originates in the vicinity of the airfoil leading edge, then travels downstream towards the airfoil trailing edge and evolves along the wake of the airfoil as seen in figure 13. Hence the similarity between the characteristics of the HFO mode and the evolution of the Kelvin–Helmholtz instability. An experiment in a water-tunnel or solving the incompressible Navier–Stokes equations numerically would answer the question of whether or not the acoustic waves generated by the HFO mode affect the transition process and the LFO phenomenon.

It is worth noting that the HFO flow along the wake of the airfoil is sinusoidal. Thus, the summation of all Fourier modes would approximate this flow accurately. However, each one of the Fourier modes does not necessarily exist in the flow. On the contrary, the POD algorithm looks for the most energetic structures in the flow convected at the same frequency and collects them in a single flow mode. The POD method continues this optimization process for all frequencies until there is not much energy left in the remaining flow. The SMD method selectively phase-averages at the chosen frequency. Consequently, the SMD method recovers the time-averaged shape of flow structures that physically exist in the flow field. The flow modes recovered by the POD and the SMD methods converge towards similar flow modes when the fluctuating flow field is closely approximated by Fourier modes.

5. The dynamics of the flow

Figures 11–13 show an ‘average’ spatial distribution of the LFO-Mode-1, the LFO-Mode-2 and the HFO mode without any description of the temporal evolution of these modes. In this section, the temporal evolution of the flow modes is investigated, and the flow field is reconstructed using the most dominant flow modes to examine the dynamics of the flow.

5.1. Temporal evolution of the flow modes

The POD eigenfunctions ($\boldsymbol {\phi }^{(n)}$) are used to reconstruct the flow modes and as basis functions for low-dimensional modelling without much attention to their amplitude or how they evolve in time. In the present study, the POD eigenfunctions were scaled with their corresponding eigenvalues ($\boldsymbol {\lambda }^{(n)}$), and the percentage of the scaled amplitude of each of the temporal POD modes is estimated at each instant in time using (B9). Figure 14 shows time histories of the POD eigenfunctions, the percentage of energy content in each POD mode and the fluctuating lift and the fluctuating drag coefficients at $\alpha =9.8^{\circ }$ and ${\textit {Re}}_c = 5\times 10^4$, and angles of attack of $10.8^{\circ }$ and $11.0^{\circ }$ at ${\textit {Re}}_c = 9\times 10^4$. The angles of attack of $9.8^{\circ }$ and $11.0^{\circ }$ are the angles of attack at which the LFO has the maximum amplitude of oscillations. The LFO at the angle of attack of $10.8^{\circ }$ exhibits three uniform cycles followed by a disturbed cycle, thus representing a special case that is worth presenting. As seen in the figure, the scaled POD eigenfunctions make more sense since they display the temporal evolution of the flow modes with their actual amplitude of energy content compared with other POD modes, and most importantly, they mimic time histories of the fluctuating lift and the fluctuating drag coefficients. At angles of attack where the flow field exhibits an LFO, the LFO-Mode-1 peaks at more than 75 %, and the LFO-Mode-2 peaks at more than 25 %. However, additional data at various Reynolds numbers and angles of attack are required to verify these thresholds.

Figure 14. Time histories of (a–c) the POD eigenfunctions ($\boldsymbol {\phi }^{(n)}(t)$), (d–f) percentage of the energy content in POD mode $n$ ($\boldsymbol {\beta }^{(n)}(t)$) and (g–i) the fluctuating lift and the fluctuating drag coefficients. The black, red and blue lines indicate the LFO-Mode-1, the LFO-Mode-2 and the HFO mode, respectively. The black and red lines display the fluctuating lift and the fluctuating drag coefficients, respectively.

The LFO-Mode-1 and the LFO-Mode-2 signals are correlated, and their magnitudes are interlinked. That is, if the LFO-Mode-1 peaks at a relatively high amplitude, so will the LFO-Mode-2, and vice versa. When the periodic cycle of the LFO-Mode-1 is disturbed, so is that of the LFO-Mode-2 as seen in the figure. The two low-frequency modes drive the dynamics of the LFO and play a contradicting role. That is, the LFO-Mode-1 decays to a minimum, during which the LFO-Mode-2 grows to a maximum. Thus, when the LFO-Mode-1 vanishes and the flow reaches a temporary equilibrium, the LFO-Mode-2 amplifies and starts a new disequilibrium. The process continues periodically as discussed by Eljack & Soria (2020). Therefore, these two low-frequency modes are interlinked and govern the dynamics of the LFO. The LFO-Mode-1 generates the triad of vortices and feeds energy into the upstream vortex of the TCV. When the upstream vortex of the triad of vortices is saturated with energy, it expands, and the LFO-Mode-2 dominates the flow until the oscillating flow switches its direction, then the LFO-Mode-1 takes over again. The HFO mode fluctuates more energetically when the flow is separated as seen in the figure. When the LFO-Mode-2 peaks during an attached flow phase, the upstream vortex of the TCV expands and separates the flow. Vortex shedding from the separated flow energizes the shedding along the wake. Consequently, the HFO mode is energized. That is, the interaction takes place, and the HFO mode becomes more energetic whenever the LFO-Mode-2 peaks to a maximum value. When the upstream vortex expands and advects downstream, it gets attracted by the low-pressure region at the trailing edge and consequently energizes the trailing-edge shedding and the HFO mode as seen in the spatial HFO mode in figure 13.

The time history of the percentage of the energy content in LFO-Mode-1 is not Gaussian as seen in figure 14. The signal exhibits sharp peaks at higher amplitudes and troughs at smaller amplitudes. When the percentage is positive, the oscillating pressure of the LFO-Mode-1 has a negative value in the vicinity of the LSB and a positive value in the vicinity of the trailing edge as seen in figure 11. Consequently, an adverse oscillating pressure works against the LFO-Mode-2 and delays the ejection of the upstream vortex. Thus, the percentage of LFO-Mode-1 peaks sharply at higher amplitudes. On the contrary, the oscillating pressure is favourable when the percentage of the LFO-Mode-1 is negative. Thus, the oscillating pressure helps the LFO-Mode-2 in ejecting the upstream vortex, and the LFO-Mode-1 exhibits troughs at smaller amplitudes.

Another interesting observation is the phase shift in the time histories of the LFO-Mode-1, the LFO-Mode-2, the fluctuating lift and the fluctuating drag. The LFO-Mode-1 and the fluctuating lift are leading the LFO-Mode-2 and the fluctuating drag by a phase shift of ${\rm \pi} /2$. The LFO-Mode-1 has a peak when the flow is attached and a trough when the flow is fully separated. Consequently, the fluctuating lift coefficient mimics the LFO-Mode-1. The upstream vortex ejects and separates (attaches) the flow whenever the LFO-Mode-2 has a peak (trough). Consequently, the oscillating pressure on the suction surface of the airfoil, the pressure difference across the airfoil in the streamwise direction and the fluctuating drag coefficient oscillate accordingly. Thus, the fluctuating drag coefficient mimics the LFO-Mode-2. The pressure field makes the flow feel the wall, and the dynamics of the flow creates oscillations in the aerodynamic forces. However, there is a time lag between the aerodynamic forces and the POD coefficients. Thus, the LFO-Mode-1 and the LFO-Mode-2 are systematically leading the fluctuating lift and the fluctuating drag with a few degrees in phase shift as seen in the figure.

5.2. Reconstruction of the flow field

The orthonormal POD spatial modes and the POD coefficients can be used to reconstruct the original flow field using any subset of the POD modes. The flow field was reconstructed and probed at selected locations using the LFO-Mode-1; the LFO-Mode-1 and the LFO-Mode-2; and the LFO-Mode-1, the LFO-Mode-2 and the HFO mode. A comparison of the reconstructed signals probed at different locations shows that the streamwise velocity, the wall-normal velocity and the pressure have interesting behaviour in the vicinity of the leading edge, mid-chord and downstream of the trailing edge on the suction surface of the airfoil. Figure 15 shows comparisons of the original LES data and the reconstructed POD data for the streamwise velocity, the wall-normal velocity and the pressure at angle of attack of $9.8^{\circ }$ and Reynolds number of $5\times 10^4$. The leading-edge, the mid-chord and the trailing-edge probe were chosen to be in a location inside the upstream vortex of the TCV, at the mid-chord point close to the wall and at about $0.4$ chords downstream of the trailing edge, respectively. The grey solid line indicates the original LES data. The black, red and blue solid lines display the reconstructed data using the LFO-Mode-1; the LFO-Mode-1 and the LFO-Mode-2; and the LFO-Mode-1, the LFO-Mode-2 and the HFO mode, respectively. Figures 15(a,d,g), 15(b,e,h) and 15(c,f,i) show signals of the flow variables in the vicinity of the airfoil leading edge, mid-chord and downstream of the trailing edge, respectively. The three most dominant POD modes (the LFO-Mode-1, the LFO-Mode-2 and the HFO mode) reconstructed the oscillating flow favourably. However, these three dominant flow modes contain a maximum of $65\,\%$ of the energy of the flow in all of the investigated angles of attack. Thus, at least $35\,\%$ of the energy content is not accounted for by these three dominant modes. It is noted that these three dominant flow modes represent the periodic motion of the flow ‘total fluctuations’ and not the random motion of the flow ‘turbulent fluctuations’. Thus, the remaining $35\,\%$ of energy content is ‘turbulent fluctuation’; therefore, it is not reflected by these three modes.

Figure 15. Time histories of the LES data of (a–c) the streamwise velocity, (d–f) the wall-normal velocity and (g–i) the pressure probed at the vicinity of the leading edge, mid-chord and downstream of the trailing edge on the suction surface of the airfoil compared with their corresponding POD reconstruction of the data at angle of attack of $9.8^{\circ }$ and Reynolds number of $5\times 10^4$. Grey solid line, LES data; black solid line, reconstructed data using the LFO-Mode-1; red solid line, reconstructed data using the LFO-Mode-1 and the LFO-Mode-2; blue solid line, reconstructed data using the LFO-Mode-1, the LFO-Mode-2 and the HFO mode.

The leading-edge probe shows that the velocity components and the pressure are reconstructed mostly by the LFO-Mode-1 and the LFO-Mode-2; thus, the HFO mode does not contribute much to any of the flow variables. This is indicative that the HFO mode does not directly influence the flow at this location. Furthermore, the original data of the velocity components are mostly recovered by the LFO-Mode-2. While the LFO-Mode-1 has a small and uniform effect, it is the LFO-Mode-2 that shaped the signal to its original LES form. This is indicative that the velocity components in this location are mostly influenced by the expansion and advection of the upstream vortex of the TCV rather than being influenced by the LFO-Mode-1. The pressure field is mostly recovered by the LFO-Mode-1 in this location with a small and limited effect of the LFO-Mode-2 and the HFO mode. The trailing-edge probe shows that the flow variables are overwhelmed by the HFO mode and its subharmonics. However, the low-frequency pattern exists in all of the flow variables at this location. Thus, the velocity components are influenced exclusively by the LFO-Mode-1 and the HFO mode. Whereas, the pressure at this location is influenced by the LFO-Mode-1 and the HFO mode in addition to a limited effect of the LFO-Mode-2.

5.3. Spatial evolution of the flow field

Figure 16 shows the reconstructed fluctuating flow using the LFO-Mode-1, the LFO-Mode-2 and the HFO mode at angle of attack of $10.8^{\circ }$ and Reynolds number of $9\times 10^4$. The main panel displays streamlines of the fluctuating flow superimposed on colour maps of the spanwise vorticity, and the top-left inset shows streamline patterns of the reconstructed instantaneous flow superimposed on colour maps of the instantaneous streamwise velocity. Time histories of the lift and drag coefficients are also shown in the bottom-right inset. The LSB and the triad of vortices are visualized by the streamline patterns of the instantaneous and fluctuating flow, respectively. The LFO-Mode-1 extracts energy from the mean flow and stores it in the triad of vortices. When the sign of the LFO-Mode-1 is positive, the oscillating flow rotates around the airfoil in the clockwise direction and flows in the streamwise direction on the suction surface of the airfoil. Consequently, the oscillating flow adds momentum to the boundary layer and helps it to remain attached against the APG. At a specific energy threshold of the LFO-Mode-2, the upstream vortex expands and separates the flow. However, the process reverses its direction as the LFO-Mode-1 switches its sign, and the process continues in a periodic manner (Eljack & Soria 2020). The HFO mode oscillates at a relatively small amplitude, but it tends to be more energetic when the LFO-Mode-2 peaks. As discussed before, the upstream vortex expands and advects downstream whenever the LFO-Mode-2 peaks. Consequently, the HFO mode fluctuates at a relatively higher amplitude during the separated phase of the flow field.

Figure 16. Streamline patterns of the POD reconstruction of the fluctuating flow superimposed on colour maps of the spanwise vorticity (main panel), streamline patterns of the POD reconstruction of the instantaneous flow superimposed on the instantaneous streamwise velocity (top-left inset) and plot of time histories of the fluctuating lift and the drag coefficients (bottom-right inset). (a) Attached phase of the flow. (b) Separated phase of the flow (multimedia view).

The evolution of the flow shows how the triad of vortices switches the direction of the fluctuating flow and separates/attaches the instantaneous flow field as seen in supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.855. The reconstructed oscillating flow using the LFO-Mode-1, the LFO-Mode-2 and the HFO mode is added to the mean flow to obtain the reconstructed instantaneous flow field. The flow is attached, and the LFO-Mode-1 is at its maximum amplitude at the flow time of $258.25$. The flow is attached, the oscillating flow direction is clockwise and the triad of vortices is present and in their most coherent state. At the flow time of $258.5$, the upstream vortex of the triad of vortices pops up above the downstream vortex of the triad of vortices. At the flow time of $259.0$, the upstream vortex pops a little more above the upstream vortex of the TCV. At the flow time of $259.5$, the upstream vortex starts to slide above the downstream vortex of the TCV. The downstream vortex starts to merge with the upstream vortex and energizes it at the flow time of $260.0$. While the downstream vortex merges with the upstream vortex, the latter continue to expand until the two vortices form one vortex at the flow time of $260.5$. At the flow time of $261.0$, the recently formed vortex expands abruptly. It is interesting to note that the LFO-Mode-2 peaks after this instant in time, as shown in figure 14 at angle of attack of $10.8^{\circ }$. After that, the recently formed vortex continues to expand and starts to change the direction of the oscillating flow. The secondary vortex is present and intact during this process. The oscillating flow then changes its direction of rotation from the clockwise to the counter-clockwise direction. At the flow time of $264.0$, the secondary vortex vanishes, and a new downstream vortex forms. As the oscillating flow continues to change its direction of rotation, a new upstream vortex starts to form at the flow time of $266.0$. It is interesting to note how the vorticity changes sign after this instant in time and how it creates the new upstream vortex of the TCV that is rotating in the counter-clockwise direction. The newly formed upstream vortex grows in size and strength at flow times of $267.5$ and $268.5$. The flow is fully separated, and the LFO-Mode-1 is at its minimum amplitude at the flow time of $269.0$. The oscillating flow direction is counter-clockwise, and the triad of vortices is present and in their most coherent state. Consequently, the lift coefficient drops to its minimum value as seen in the figure. Reconstruction of the flow field provided a detailed description of how the flow separates. When the cycle of the LFO proceeds to a flow time larger than $270.0$, the process reverses its direction, and the fully separated flow starts to attach. Thus, a similar sequence of events takes place, as shown in supplementary movie 1. The evolution of the oscillating flow is animated to reveal how the flow separates and attaches. Thus, a detailed description of the LFO phenomenon is given, and the underlying mechanism is discussed in detail. It is important to mention that the POD spatial mode corresponding to LFO-Mode-1 does not exhibit the triad of vortices. However, the reconstruction of the flow field using the three dominant flow modes clearly shows the evolution and dynamics of the triad of vortices. None of the POD modes recover the triad of vortices as discussed in § 4.4, but the summation of two or more of the most dominant POD modes reconstructs a flow field that exhibits the triad of vortices and their evolution.

Figure 17 shows the reconstructed fluctuating flow using the SMD flow modes corresponding to the LFO-Mode-1, the LFO-Mode-2 and the HFO mode. The flow is attached at the starting point in time at $t=206.0$ as seen in figure 14. The evolution of the oscillating flow is visualized in step-by-step snapshots that show how the triad of vortices separates the flow. When the cycle of the LFO proceeds to flow time ${>}217.5$, the process reverses its direction, and the fully separated instantaneous flow starts to attach as seen in figure 14. Thus, a similar sequence of events takes place, and the SMD method accurately reconstructs the flow field using the most dominant flow modes.

Figure 17. (a–r) Streamline patterns superimposed on colour maps of the spanwise vorticity ($\omega _z$) for the SMD reconstruction of the fluctuating flow at $\alpha =9.8^{\circ }$ and Reynolds number of $5\times 10^4$.

5.4. The effect of the spanwise extent on the LSB and the LFO phenomenon

Figure 5 shows that the LES data of the present work compare very well with the LES data of Alferez et al. (2013) and Alferez (2014) representing the flow field about a NACA0012 airfoil of aspect ratio of one. However, we would like to examine the effect that the narrow span of $0.5$ chords used in the present work can have on the obtained results. For instance, stall cells would appear at relatively high angles of attack; however, they might be suppressed by the narrowness of the span, consequently questioning the accuracy of the data and the generality of the conclusions based on them.

Previous studies on two-dimensional airfoil stalling characteristics have shown a low-frequency and highly unsteady flow or a steady large-scale three-dimensional structure. The latter is termed stall cells, and there has been a considerable amount of research focused on their structure (Winkelman & Barlow 1980; Dallmann 1983; Winkelman 1990; Yon & Katz 1998; Rodríguez & Theofilis 2010, 2011; He et al. 2017). Winkelman & Barlow (1980) carried out oil-flow visualizations of the flow field on the suction surface of stalled wings of a Clark Y airfoil section at chord Reynolds numbers of $24.5\times 10^4$, $26\times 10^4$ and $38.5\times 10^4$. The ends of the wing are flush with the tunnel side walls or splitter plates and on plane rectangular wings of infinite aspect ratio. They observed that mushroom-shaped cells started to form at the onset of stall on the two-dimensional model. The three-dimensional cells coexisted with trailing-edge stall cells on the wing's surface several degrees above the stalling angle of attack. They tentatively sketched a flow field model showing the general features of a leading-edge separation bubble and trailing-edge separation. Winkelman (1990) took up the work of Winkelman & Barlow (1980) and measured the spectra of the velocity in the wake of a rectangular wing having the same airfoil section. The spectra of the wake measurements did not show any indications of a dominant low-frequency mode. Yon & Katz (1998) used fine-thread tuft-flow visualization and high-frequency response pressure transducer measurements to investigate the unsteady features of the stall cells. They studied the flow field around a wing of aspect ratio ranging from $2$ to $6$ with a NACA0015 airfoil section. They imposed the two-dimensionality using end plates that effectively eliminated the tip flow. The authors reported the existence of a low-frequency mode in their pressure spectrum but with a considerably small amplitude of oscillation.

Broeren & Bragg (2001) studied five different airfoil configurations (NACA2414, NACA64A010, LRN-1007, E374 and Ultra-Sport) by measuring the wake velocity across the spanwise direction and using mini-tufts for flow visualization. They found that all the stall types were dependent on the type of airfoil. They concluded that the LFO phenomenon always occurs in airfoils that exhibit a thin-airfoil stall or a combination of thin-airfoil and trailing-edge stall. Most importantly, they found that the low-frequency unsteadiness is essentially two-dimensional. Their conclusions were in good agreement with that of Zaman et al. (1989), who observed that the LFO takes place with airfoils exhibiting either trailing-edge- or thin-airfoil-type stalls but does not take place with the leading-edge-type stall. To this end, it is evident that the LFO phenomenon is inherently two-dimensional by nature, and neither the imposed periodic boundary condition nor the spanwise width of the computational domain affects the accuracy of the data and the conclusions based on them.

To further examine the effect of spanwise width on the shape of the LSB and the LFO phenomenon, an LES was carried out at an angle of attack of $9.4^{\circ }$, Reynolds number of $5\times 10^4$ and a spanwise width of two chords. The time history of the lift coefficient at each location in the spanwise direction $\textit{C}_{\!_\textit{L}}(x_3, t)$ was used to estimate the time-averaged surplus of the oscillating flow on the fly to obtain the three-dimensional LFO-Mode-1. Figure 18 shows the time-averaged shape of the LSB and the triad of vortices. A uniform seeding of the streamlines was done by drawing a line parallel to the spanwise direction and specifying a certain number of points along this line. Each of the streamlines that pass by any of these points is automatically completed and allowed to extend in both directions; hence, the visualization of the streamlines is user-independent. Despite the irregularities encountered in the LSB and the triad of vortices, the flow field is essentially two-dimensional, and stall cells are not present in any form as seen in the figure.

Figure 18. Streamline patterns of the mean flow and the SMD mode corresponding to the lift coefficient plotted around a NACA0012 airfoil of spanwise width of two chords at an angle of attack of $9.4^{\circ }$ and Reynolds number of $5\times 10^4$ showing the time-averaged shape of the LSB (a) and the triad of vortices (b).

5.5. The effect of compressibility on the LFO phenomenon

Classically, the bursting of the LSB and its associated LFO used to be investigated in compressible flow settings. Experiments are carried out in wind-tunnels rather than in water-tunnels and compressible Navier–Stokes equations are solved in numerical simulations rather than solving the incompressible set of equations. The vast majority of investigations being compressible is mainly due to the widely accepted hypothesis that there is an acoustic wave feedback mechanism involved. It is assumed that the acoustic waves generated at the trailing edge of the airfoil travel upstream and interact with some receptivity mechanism in the vicinity of the leading edge (Jones et al. 2008, 2010; Almutairi et al. 2017; Kurelek, Kotsonis & Yarusevych 2018; Pröbsting & Yarusevych 2021). Such feedback would force the shear layer to undergo early transition and reattach the flow. The movement of the transition location along the shear layer, causing late transition/early transition, is synchronized with the flow separation and reattachment, respectively. Thus, the acoustic wave feedback mechanism generates and sustains the LFO phenomenon. On the contrary, flow forcing using a selected set of frequencies could remove the LFO as reported by Zaman et al. (1987) and Eljack, Alqadi & Almutairi (2018). Such a conflicting role of the acoustic excitation was also noted by Zaman et al. (1989).

Almutairi, Alqadi & Eljack (2015) and Almutairi et al. (2017) applied the DMD method to the pressure field of the flow field around a NACA0012 airfoil at ${\textit {Re}}_c = 1.3\times 10^5$, $M = 0.4$ and angle of attack of $11.5^{\circ }$. The authors sampled the instantaneous pressure field on the $x$–$y$ plane at a non-dimensional frequency of $658$, and the data span about one low-frequency cycle ($25$ non-dimensional time units). The DMD identified two dominant flow modes: a low-frequency mode at Strouhal number of $0.008$ featuring the bursting and reformation cycle of the LSB and a high-frequency mode featuring the trailing-edge shedding frequency. They concluded that the trailing-edge shedding induces acoustic waves that travel upstream and excite the separated shear layer via some receptivity mechanism and forces it to undergo early transition and reattach the separated flow. Once the flow is attached, the strong vortex shedding downstream of the trailing edge dies down, and the acoustic feedback is cut off; consequently, the flow separates again, and the separation/reattachment process repeats periodically. However, the authors applied the DMD method to the pressure field only; therefore, the authors did not show streamlines corresponding to the low-frequency mode. Consequently, the triad of vortices discovered and explained by Eljack & Soria (2020) was not presented in their work.

The model of the LFO phenomenon presented by Eljack & Soria (2020) and the discussion in § 5.3 above show that the circumstances that govern the stability of the LSB and its associated LFO phenomenon can be met in both compressible and incompressible flows. Figure 19 shows the variation of the mean velocity above and below that of the mean flow due to attachment and separation of the flow, respectively. The production of kinetic energy increases and an early transition occurs when the flow is attached due to the increase in the velocity gradient across the separated shear layer and vice versa. Thus, the movement of the transition location is interlinked with the oscillating flow induced by the LFO-Mode-1 and the LFO-Mode-2. Consequently, the oscillating flow affects the velocity gradient across the separated shear layer and the development of the Kelvin–Helmholtz instability. In conclusion, previous work on the LFO suggests that the acoustic wave feedback mechanism is crucial for the LFO phenomenon. Whereas, the present work shows that an interplay of three vortices induces the LFO phenomenon without excluding the acoustic waves scenario.

Figure 19. (a,b) The scaled mean streamwise velocity (${U}_1/{U}_{\infty }$) at the location of maximum production of turbulent kinetic energy along the separated shear layer at the stall angle of attack and the angle of attack of maximum LFO. Black line, time-averaged flow; red line, time-average attached flow; blue line, time-average separated flow.

Rodríguez et al. (2013, 2015) carried out several direct numerical simulations in which they solved the incompressible Navier–Stokes equations. They investigated the flow field around a NACA0012 airfoil at Reynolds number of $5\times 10^4$ and angles of attack of $5^{\circ }$, $8^{\circ }$, $9.25^{\circ }$ and $12^{\circ }$. The authors reported that they started the simulations from an initially homogeneous flow field that introduced some numerical disturbances. However, the authors did not clearly state whether or not they captured a global LFO. They have not shown time histories of the aerodynamic coefficients; thus, the reader cannot determine if the LFO phenomenon was captured. The only low-frequency flow phenomenon reported by the authors is the shear-layer flapping on the suction surface of the airfoil at Strouhal number of $0.021$. To the best of the author's knowledge, these are the only incompressible studies that investigated the LFO in the flow field around an airfoil at near-stall conditions. However, further investigations are required before a conclusion on whether or not an LFO phenomenon can be captured in an incompressible flow. Therefore, experimental measurements in a water-tunnel or solving the incompressible Navier–Stokes equations in a direct numerical simulation are merited. In the case of the numerical simulation, starting the simulation using clean flow might help capture the phenomenon. Whereas, in the water-tunnel experiment, adjusting the free-stream turbulence to the right level is crucial for capturing the phenomenon.

6. The bursting criterion of the LSB

Diwan et al. (2006) refined the bursting criterion developed by Gaster (1967) to consider the effect of the height of the LSB in addition to its length. The bursting parameter is given by $P_h=({h^2}/{\nu })({\Delta U}/{\Delta x})$, where $h$ is the height of the LSB, $\nu$ is the kinematic viscosity, $\Delta U$ is the velocity difference across the bubble and $\Delta x$ is the length of the LSB. The Reynolds number based on the height of the bubble is given by $Re_h = ({{U}_h h}/{\nu })$, where ${U}_h$ is the velocity at the maximum height of the LSB ($h$). The bubble is termed short if $P_h > -28$. Figure 20 shows a plot of the bursting parameter $P_h$ versus the Reynolds number $Re_h$. As seen in the figure, the bursting criterion suggests that the LSBs at Reynolds number of $5\times 10^4$ and angles of attack lower than the stall angle of attack are short bubbles. This is in total agreement with the definition of the short bubble and the above discussion which shows that the LSBs are long bubbles at angles of attack higher than the stall angle. The bursting parameter $P_h$ is proportional to $h^2$, and the LSBs are thinner at Reynolds number of $9\times 10^4$ than their ${\textit {Re}}_c=5\times 10^4$ counterparts. Thus, the criterion suggests that all the LSBs at Reynolds number of $9\times 10^4$ are short bubbles, as seen in the figure.

Figure 20. The bursting parameter $P_h$ plotted versus the Reynolds number $Re_h$. The arrows indicate the direction in which the angle of attack, $\alpha$, increases.

The percentage of the energy content in each POD and SMD mode is shown in figure 21. The energy percentage for the POD modes is estimated from the ratio of each POD eigenvalue ($\boldsymbol {\lambda }^{(n)}$) to the sum of all POD eigenvalues. Whereas, the energy percentage for the SMD modes is estimated from the variance of each SMD coefficient divided by the sum of the energy in all of the flow modes. Classification of the most dominant POD mode as a low-frequency or high-frequency mode was identified from the corresponding temporal POD mode. The SMD modes are selectively recovered; thus, the modes are identified accordingly. The POD mode corresponding to the LFO-Mode-1 dominates the flow field in all of the investigated angles of attack and at all of the investigated Reynolds numbers, as seen in the figure. This is because the POD method recovers flow modes based on their energy content, thus overestimating the amplitude of some flow modes. On the contrary, the SMD method recovers physical flow modes; therefore, the amplitude of the SMD mode corresponding to the lift coefficient (LFO-Mode-1) decays gradually as the angle of attack decreases until it is overtaken by the LFO-Mode-2 and the HFO mode at angles of attack lower than the stall angle of attack as seen in the figure. Consequently, the SMD mode corresponding to the HFO mode dominates the flow field at angles of attack lower than the stall angle, including at zero angle of attack. It is well known that the flow field around airfoils at zero or relatively small angles of attack is dominated by the HFO mode along the wake. Thus, the SMD method adjusts the amplitude of the flow modes to the correct magnitude and overcomes the shortcomings of the POD method.

Figure 21. The percentage of the energy content in POD mode $n$ ($\kappa ^{(n)}$) and each of the SMD modes at (a,b) angles of attack of $8.8^{\circ }$–$10.5^{\circ }$ and Reynolds number of $5\times 10^4$ and (c) angles of attack of $10.0^{\circ }$–$11.2^{\circ }$ and Reynolds number of $9\times 10^4$. The black, red and blue bars indicate the energy percentage of the LFO-Mode-1, the LFO-Mode-2 and the HFO mode, respectively.

At angles of attack lower than the stall angle, the LFO-Mode-2 is less energetic than the HFO mode. Thus, the flow remains attached with occasional separations at small amplitudes of oscillations, and the LSB remains intact. The flow at these angles of attack oscillates at a low frequency as a consequence of vortex shedding and shear-layer flapping, but the LFO-Mode-2 is not energetic enough to eject the upstream vortex, separate the flow and induce instability of the LSB. At angles of attack above the stall angle of attack and below the angle of a full stall, the LFO-Mode-2 overtakes the HFO mode and becomes the most dominant flow mode. The LFO-Mode-2 becomes energetic enough to separate the flow, and the LSB bursts and reforms intermittently. As the angle of attack increases above the angle of attack of a full stall, the HFO mode overtakes the LFO-Mode-2 again, an open bubble forms on the suction surface of the airfoil and the LFO-Mode-2 is not energetic enough to attach the flow. The LFO-Mode-1 represents an absolute instability, and the HFO features a convective instability. Thus, the instability changes from convective to absolute when the LFO-Mode-2 overtakes the HFO at the inception of the stall, triggers the instability of the bubble and initiates the LFO phenomenon. The above discussion shows that the bursting of the LSB depends on the percentage of energy content in the LFO-Mode-1, the LFO-Mode-2 and the HFO mode. Whereas, previous studies associate bursting conditions of the LSB with local flow parameters such as the value of the Reynolds number of the separated boundary layer and the pressure rise over the region occupied by the bubble (Owen & Klanfer 1953; Gaster 1967). However, these local flow parameters depend on the energy content of the most dominant flow modes. Thus, previous bursting criteria remain valid, albeit are not universal, as pointed out by Diwan et al. (2006) and confirmed by Alferez et al. (2013).

7. Conclusion

The objective of the present work was to carry out a detailed flow dissection and shed some light on the structure and dynamics of the flow around a NACA0012 airfoil at low Reynolds number and near-stall conditions. Time histories of the fluctuating aerodynamic forces have been utilized to reconstruct selective flow modes that cause oscillations in these forces, and a novel SMD method has been developed. Three distinct dominant flow modes have been identified by the POD, the DMD and the SMD methods: a globally oscillating flow mode at a low frequency (LFO-Mode-1) featuring an absolute instability that creates and sustains the triad of vortices; a locally oscillating flow mode on the suction surface of the airfoil at a low frequency (LFO-Mode-2) presenting the ejection of the upstream vortex of the TCV; and a locally oscillating flow mode along the wake of the airfoil at a high frequency (HFO mode) featuring the travelling Kelvin–Helmholtz waves along the airfoil wake (convective instability). The analysis showed that the dynamics of the LSB is associated with the percentage of energy content in each of the three dominant flow modes. At angles of attack below the stall angle of attack, the HFO mode dominates the flow, and the LSB remains intact. At angles of attack above the stall angle of attack and below the angle of attack of a full stall, the LFO-Mode-2 becomes more energetic than the HFO mode. Consequently, the LSB bursts to form a long bubble and triggers low-frequency oscillations in the flow field.

The underlying mechanism and the bursting criterion presented in the current study are for flow fields about airfoils at moderate subsonic Mach number ($0.4$). However, the change in the Mach number from $0.2$ (near incompressible) to $0.4$ is not expected to dramatically change the shape of the identified dominant flow modes, the proposed underlying mechanism and the suggested bursting criterion. The structure of the LSB and the dynamics of the flow are statistically two-dimensional. This was confirmed by analysing a dataset generated in a computational domain of a spanwise width of two chords. It is evident that neither the imposed periodic boundary condition nor the spanwise width of the computational domain affects the accuracy of the data and the conclusions based on them. The amplitudes of flow modes and their energy content depend primarily on the accuracy of the statistical method used to recover them. The SMD method is based on time-averaging of the data; thus, it accurately estimates the energy content in each flow mode within the statistical error.

The dynamics of the flow discussed in the present work shows that the physics of the LSB formed naturally on the suction surface of an airfoil exhibits some features that do not evolve and sustain in the LSB induced on a flat plate with an adjustable pressure distribution (Gaster 1967). Therefore, a comparative study of the dominating flow modes in the two configurations is merited. Another interesting model that has been used to tackle this problem is solving linearized Navier–Stokes equations to investigate instability mechanisms of the flow. However, this study showed that the fluctuating flow variables are well above $1\,\%$ of their corresponding mean values. Therefore, a thorough study to examine the percentage of the fluctuating flow at the instant of bursting is also merited. The present investigation opens the door for the optimum design of airfoils, and has profound implications for the constantly increasing applications that operate at low-Reynolds-number conditions. Finally, the SMD is a robust method that utilizes the time history of an aerodynamic force to selectively recover a flow mode that physically exists in the flow and causes fluctuations in the aerodynamic force.