1. Introduction
Dynamic stall is a topic of great interest in unsteady aerodynamics since it can lead to aerodynamic forces and moments severe enough to cause catastrophic structural failure (McCroskey Reference McCroskey1981; Corke & Thomas Reference Corke and Thomas2015). Stall control efforts are most effective before the formation of the dynamic stall vortex (Chandrasekhara Reference Chandrasekhara2007), a characteristic feature of ‘deep’ dynamic stall. Therefore, characterizing stall onset is of crucial importance for control efforts to be deployed in a timely manner. Various criteria for dynamic stall onset based on unsteady aerodynamic coefficients have been explored to formulate first-order, semi-empirical, dynamic stall models (Leishman & Beddoes Reference Leishman and Beddoes1989; Sheng, Galbraith & Coton Reference Sheng, Galbraith and Coton2005). Several stall criteria have been proposed to narrow down the identification of stall onset to a finer degree in time. These include the leading edge suction parameter or the 
$LESP$ (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014), which is pressure based, and the boundary enstrophy flux or the 
$BEF$ (Sudharsan, Ganapathysubramanian & Sharma Reference Sudharsan, Ganapathysubramanian and Sharma2022), which is vorticity based. In the present work, we analyse these criteria applied to dynamic stall at low Reynolds number.
Low-Reynolds-number, 
$Re \sim O(10^4)$, unsteady flows are dominated by a rich variety of coherent vortical structures such as shear layer vortices, a dynamic stall vortex (DSV) system comprising multiple vortices and induced secondary vortex flow structures. These vortices generally start out as laminar with strong spanwise coherence and undergo transition as the airfoil angle of attack (
$\alpha$) increases. In their large-eddy simulations (LES) of a plunging airfoil at an 
$Re$ of 60 000, Visbal (Reference Visbal2011) found that the DSV system and the shear layer vortices independently undergo flow transition. In the same study, they observed that the flow remained laminar through the entire plunging cycle at an 
$Re$ of 1000.
For flow over a manoeuvring airfoil at low 
$Re$ and initially small 
$\alpha$, the laminar boundary layer is susceptible to separation when subject to an adverse pressure gradient (APG), resulting in the formation of a laminar separation bubble (LSB). If the APG remains low, the free-shear layer above the bubble does not transition and instead reattaches laminarly. This leads to a long LSB extending over most of the airfoil with fully laminar flow (Shyy et al. Reference Shyy, Lian, Tang, Viieru and Liu2007). The structure of the LSB is highly sensitive to 
$Re$, as determined from experiments on an NACA 
$66_3$-
$018$ airfoil by O'Meara & Mueller (Reference O'Meara and Mueller1987), where the LSB more than doubled in length when the 
$Re$ was decreased from 140 000 to 50 000. Ol et al. (Reference Ol, Bernal, Kang and Shyy2009) reported a similar increase in LSB size with reducing 
$Re$ in their experiments in a free-surface water tunnel where shallow and deep stall cases were investigated for a pitch-plunge motion of a rigid airfoil in the range 
$10\,000 \le Re \le 60\,000$. As 
$\alpha$ or 
$Re$ is increased, the APG downstream of the suction peak increases, encouraging spanwise instabilities that lead to earlier transition and reattachment, and shorter LSBs. Therefore, during a pitch-up motion at moderate pitch rates, the DSV system could remain laminar through part of the motion depending on the 
$Re$. Unsteady flow at low 
$Re$ can consist of a mix of laminar, transitional and turbulent flow regimes in space and time.
In contrast, transitional effects are spatially localized at high 
$Re$ since the shear layer transitions close to the leading edge, with turbulent flow prevailing over most of the airfoil. It has also been observed that a single coherent DSV is more typical of high 
$Re$, while a DSV system consisting of one or more laminar/transitional vortices is likely at low 
$Re$ (Galbraith & Visbal Reference Galbraith and Visbal2010). At high 
$Re$, there is a pronounced effect on the moment coefficient, 
$C_m$, which undergoes severe divergence (Visbal Reference Visbal2011). In contrast, multiple, relatively less severe moment stall occurrences are typical at low 
$Re$. DSVs are shed faster and stall occurs earlier. The early stall is attributed to the stronger viscous response of the boundary layer at low 
$Re$, leading to the earlier formation of secondary structures having strong circulation that cut off the primary vortices from the feeding shear layer (Widmann & Tropea Reference Widmann and Tropea2017).
The current state of the art in characterizing unsteady stall is the 
$LESP$ (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014; Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020), which represents the camber-wise suction force at the leading edge. 
$LESP$ has been used to trigger leading edge vortex shedding in reduced-order models based on a pre-determined critical value; 
$\max (LESP)$ is a standalone criterion that has been used as a proxy for critical 
$LESP$ in some cases, e.g. Deparday & Mulleners (Reference Deparday and Mulleners2018, Reference Deparday and Mulleners2019), to identify the onset of dynamic stall. 
$BEF$ has previously been explored for 
$Re \sim {O}(10^5\unicode{x2013}10^6)$ and has been found to reach its maximum magnitude in advance of DSV formation. More generally, the 
$\max (|BEF|)$ criterion signifies the instance of maximum wall shear and indicates imminent vortex formation. Sudharsan, Narsipur & Sharma (Reference Sudharsan, Narsipur and Sharma2023) demonstrated that both the 
$\max (LESP)$ and 
$\max (|BEF|)$ criteria effectively signal stall onset, irrespective of the stall type (leading-edge, mixed or trailing-edge), at 
$Re \sim {O}(10^5\unicode{x2013}10^6)$ when the DSV significantly influences the stall process.
Low-
$Re$ unsteady flows pose further demands on these criteria. In addition to detecting DSV formation, the criteria should also identify strong vortex-shedding events, which may occur far downstream of the airfoil leading edge. The strong vorticity associated with laminar vortices (Visbal Reference Visbal2009) and their formation is accompanied by larger streamwise gradients and edge velocities (Shyy et al. Reference Shyy, Lian, Tang, Viieru and Liu2007). These observations suggest that the 
$BEF$, which is sensitive to changes in pressure gradient and vorticity, would be effective at signalling instances of vortex formation.
While some studies, e.g. Visbal (Reference Visbal2009, Reference Visbal2011), have provided descriptions of the unsteady flow field at low 
$Re$, we focus on the development of transitional instabilities in the flow and how they affect leading edge flow, where stall criteria are typically evaluated. We assess the applicability of the 
$\max (|BEF|)$ and 
$\max (LESP)$ criteria to 
$Re \sim {O}(10^4)$, both for signalling imminent DSV formation and localized vortex shedding. Expanding their applicability would enable the use of a fundamental, standalone parameter in reduced-order dynamic stall models for the prediction of vortex shedding and DSV formation over a wider range of 
$Re$. It also holds significant promise for stall control efforts.
2. Methods and datasets
Our analysis is based on wall-resolved LES carried out using the compressible flow solver FDL3DI (Gaitonde & Visbal Reference Gaitonde and Visbal1998). The flow over an SD7003 airfoil undergoing a constant-rate, pitch-up motion is simulated at low-to-moderate pitch rates (
$\varOmega _0^+ = {\varOmega _0 c}/{U_{\infty }}$, with freestream velocity 
$U_{\infty } = 1$) of 
$0.05$ and 
$0.25$ at an 
$Re$ of 10 000 and 60 000. Spanwise-periodic boundary conditions applied to the ends of a span of length 
$0.2c$ are used to simulate an infinite wing geometry. Extensive studies have been carried out on the effect of grid resolution and spanwise extent for dynamic stall simulations using the FDL3DI solver in the same 
$Re$ range by Visbal (Reference Visbal2009, Reference Visbal2011) and Garmann & Visbal (Reference Garmann and Visbal2009). The O-grid mesh used in the present study consists of 
$554 \times 380 \times 101$ points in the circumferential, radial and spanwise directions, respectively. The selected discretization is on par with the finest grids simulated in the cited studies. The simulated airfoil has a unit chord, and the farfield boundary is located approximately 100 chord lengths away, where freestream conditions are specified. The airfoil surface is modelled as an adiabatic, no-slip wall. The spanwise extent, while limited by computational cost, must be large enough to prevent the artificial suppression of large-scale spanwise instabilities and avoid numerical artefacts due to the imposed span periodicity. A uniform span-wise spacing of 
$0.002c$ between grid points is used over the spanwise extent of 
$0.2c$, which is sufficiently large for the 
$Re$ and pitch rates under consideration, based on the studies cited above. Even for the largest 
$Re$ considered in this study, the 
$y^+$ values in the static simulations remain well below 
$1$ over the entire airfoil surface. Appendix A presents some of the static simulation results. Figure 1 shows images of the grid used in the present study. A non-dimensional time step size (
$\Delta t^* = \Delta t U_{\infty }/c$) of 
$1 \times 10^{-4}$ is used for time integration, as is typical in LES simulations using FDL3DI in the literature for similar 
$Re$. Additional details on the solver are available from Sharma & Visbal (Reference Sharma and Visbal2019) and Visbal & Gaitonde (Reference Visbal and Gaitonde2002).
Grid used in the present study: (a) full view; (b) zoomed-in view and (c) trailing-edge region. Every third point in the radial and circumferential directions is shown for clarity in panels (a,b).
Table 1 lists the datasets used in the current analysis. In all cases, an SD7003 airfoil having a maximum thickness of 8.5 % chord at 
$x= 0.24$ and a camber of 1.4 % chord undergoes a constant-rate, pitch-up motion about its quarter-chord point at a freestream Mach number, 
$M_{\infty } = 0.1$. A smooth hyperbolic tangent function is used to reach the final non-dimensional pitch rate, 
$\varOmega _0^+$, as described by Sharma & Visbal (Reference Sharma and Visbal2019). The results reported herein are obtained by averaging the three-dimensional solutions in the spanwise direction. The datasets consist of two different 
$Re$ (10 000 and 60 000) at two different 
$\varOmega _0^+$ (0.05 and 0.25), for a total of four cases. The acronyms in table 1 are used to refer to the different cases in the rest of the paper. The numbers following ‘
$R$’ represent the 
$Re$ in thousands and the numbers following ‘
$p$’ refer to the non-dimensional pitch rate in hundredths.
Datasets used in the present work. In all cases, an SD7003 airfoil is pitched up at a constant rate about the quarter-chord point, at 
$M_{\infty } = 0.1$.
3. Results and discussion
3.1. Definitions of 
$BEF$ and $LESP$
We first provide the definitions of the 
$BEF$ and 
$LESP$ parameters for reference. For a two-dimensional flow field, the 
$BEF$ is the flux of the squared spanwise vorticity at the wall scaled by 
$Re$, as shown in (3.1). That is, 
$\partial (\omega ^2/2) / \partial n$, which is written as a product of vorticity (
$\omega$) and vorticity flux (
$\partial \omega / \partial n$), and scaled by 
$Re$.
\begin{equation} BEF = \frac{1}{Re}\int_{x_p}^{x_s} \omega \frac{\partial \omega}{\partial n} \,\rm{d}s,\end{equation}
where 
$\omega$ is normalized by 
$U_{\infty }/c$. The chord-wise coordinate, 
$x$, and the normal and tangential coordinates to the airfoil surface, 
$n$ and 
$s$, respectively, are normalized by 
$c$. The integral is carried out between some value of 
$x$ on the pressure side to that on the suction side (shown in figure 2), with all quantities calculated in the airfoil frame of reference. The factor 
$Re$ can be combined with the vorticity flux to write an equivalence relation with the favourable streamwise pressure gradient for small tangential accelerations, i.e. 
$(1 / Re) (\partial \omega / \partial n) = - (1 / \rho ) (\partial p / \partial s)$. Therefore, large contributions to the 
$BEF$ arise from regions where high vorticity coexists with large pressure gradients.
Schematic showing the coordinate directions and integration region for calculating 
$BEF$ and 
$LESP$.
$LESP$, given by (3.2), is a measure of the chord-wise or camber-wise suction force, 
$F_{suction}$, near the leading edge, obtained by integrating surface pressure. The 
$\hat {e}_n$ and 
$\hat {e}_x$ in the definition of 
$F_{suction}$ are unit vectors in the normal direction to the airfoil boundary and along the chord direction, respectively, as shown in figure 2. In the forward part of the airfoil (upstream of the max thickness point, located at 
$x = 0.24$ for the SD7003 airfoil), 
$F_{suction}$ points in the 
$-x$ direction, indicating a ‘suction’ force acting on the airfoil, which explains the nomenclature. The integral to obtain 
$F_{suction}$ is conventionally carried out from the maximum thickness point on the pressure side to that on the suction side. Since the SD7003 airfoil has a very small camber (
${\sim }0.014c$), the chord-wise direction has been used to compute 
$LESP$, which has negligibly small differences compared with when using the camber direction.
\begin{align}
\left.\begin{array}{@{}c@{}} LESP = \sqrt{|C_{suction}|/(2
\pi)}, \quad {\rm where} \ C_{suction} =
F_{suction}/(q_{\infty}c), \quad q_\infty = \rho_\infty
U^2_\infty\\ {\rm and}\quad F_{suction}
= q_{\infty} \displaystyle \int_{x_p}^{x_s} C_p
\hat{e}_n \cdot \hat{e}_{x}\, \rm{d}s.
\end{array}\right\}
\end{align} In reduced-order models based on unsteady thin airfoil theory, the instance of vortex shedding from the leading edge is determined by the 
$LESP$ reaching a critical value (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). The idea is based on the suction at the leading edge being limited (via vortex shedding) to a certain critical value that can be supported. The 
$LESP$ is a measure of the 
$A_0$ coefficient in the Fourier series expansion of vorticity from classical thin airfoil theory (Katz & Plotkin Reference Katz and Plotkin2001); 
$A_0$ represents the vorticity at the leading edge. A leading-edge vortex is, therefore, shed any time the 
$LESP$ reaches the critical value, thus limiting the maximum 
$LESP$ in the model to the critical value. This approach requires the critical 
$LESP$ value to be determined a priori (via simulations or measurements) instead of real-time parameter tracking; critical 
$LESP$ depends on airfoil geometry and operating conditions. Here, 
$\max (LESP)$ has been used as a proxy for critical 
$LESP$ to circumvent this impediment (Deparday & Mulleners Reference Deparday and Mulleners2018, Reference Deparday and Mulleners2019), and will be evaluated in the present work for comparison.
The critical value of 
$LESP$ is its value at the time instant when the profile of the skin friction coefficient (
$C_f$) over the suction surface of the airfoil near the leading edge first develops an inflection point (Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020). This ‘
$C_f$-signature’ criterion signifies the development of instabilities within the LSB, leading to the bursting of the LSB, which is followed by DSV formation. Note that the 
$C_f$-signature criterion per se is hard to use as a stall indicator because it is spatially localized, and the location of the inflection point is not fixed or known a priori. Spatially integrated quantities such as the 
$LESP$ and 
$BEF$ are, therefore, preferred.
The results from the numerical simulations are described in detail in the following sections.
3.2. Case R60-p05
We begin by describing the results obtained at an 
$Re$ of 60 000 and pitch rate of 
$0.05$. Figure 3 shows the space–time contours of 
$-C_p$ and 
$C_f$ with the normalized chord-wise distance along the airfoil suction surface (
$x$) as the abscissa and angle of attack (
$\alpha$) as the ordinate. The shear layer that separates from the airfoil leading edge develops inviscid, Kelvin–Helmholtz (K-H)-type instabilities downstream, resulting in the shedding of laminar vortices at 60 % chord. The shed vortices subsequently transition to turbulence, while the instability and transition locations on the surface move upstream as the airfoil pitches up. These are pointed out in the 
$C_f$ contours in figure 3(b).
Space–time contours for Case R60-p05. (a) 
$-C_p$ and (b) 
$C_f$.
At 
$\alpha \sim 11^{\circ }$, the initial K-H instabilities give way to viscous instabilities as an LSB is established close to the leading edge. A couple of strong laminar vortices are shed at this time, as pointed out in figure 3(b). The separated shear layer transitions to turbulence and reattaches as the increasing APG amplifies the instabilities. Figure 4 shows isosurfaces of the Q-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988) coloured by spanwise vorticity values at various instances. Stall occurs primarily due to flow breakdown at the leading edge. Contours of spanwise vorticity, 
$\omega$, as the unsteady motion progresses are included in supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.753.
Isosurfaces of Q-criterion (value 100) coloured by spanwise vorticity contours (inlaid legend) showing the shear layer undergoing (a) a K-H instability and subsequent transition to turbulence, (b) upstream propagation of the K-H instability and (c) development of viscous instabilities close to the leading edge, for Case R60-p05. (a) 
$\alpha = 6.9^{\circ }$, (b) 
$\alpha = 9.3^{\circ }$ and (c) 
$\alpha = 10.7^{\circ }$.
A DSV is formed from the shed vorticity near the leading edge. The DSV grows and convects downstream, as seen from the imprint on the airfoil suction surface in the contours shown in figure 3. The downstream convection of the DSV results in a moment stall, and its subsequent shedding results in a lift stall, as observed from the aerodynamic coefficients (
$C_l$, 
$C_d$ and 
$C_m$) plotted against 
$\alpha$ in figure 5(a–c). The non-dimensional time, 
$t^*$, is plotted in panel (a) for reference. Figure 5(d) shows the maximum magnitude of 
$C_p$ within the first 
$5\,\%$ chord from the leading edge. Note that 
$C_p$ is a point quantity, in contrast to the aerodynamic coefficients plotted in panels (a)–(c), which are integrated quantities. As the airfoil pitches up, the stagnation point moves downstream on the pressure side, even as the shear layer slowly moves away from the surface on the suction side. This leads to a net increase in the curvature of the shear layer near the leading edge, leading to a rise in 
$C_p$. However, when strong vortices are generated and shed at 
$\alpha \sim 10^{\circ }$, there is a slight reduction in the slope of the 
$\max (|C_p|)$ curve. This reduction occurs because the growth of these vortices induces strong secondary vorticity beneath them, resulting in the shear layer being pushed away from the surface at a faster rate, reducing its curvature. Once the vortices are shed, the shear layer movement away from the surface again slows down, leading to a further increase in 
$C_p$ magnitude. After the leading edge flow breaks down at 
$13.5^{\circ }$, the magnitude of 
$C_p$ begins to reduce.
Variation with 
$\alpha$ of (a–c) aerodynamic coefficients and (d) 
$\max (|C_p|)$ near the first 5 % of airfoil chord, for Case R60-p05.
Figure 6 shows 
$C_p$ profiles around the time when strong laminar vortices are shed. At 
$\alpha = 10.1^{\circ }$, the strengths of the shed vortices markedly increase, as pointed out in the 
$C_f$ contours in figure 3(b). In addition to the instance of DSV formation, we are also interested in the instances of strong laminar vortex shedding. At low 
$Re$, these vortex cores have 
$C_p$ magnitudes comparable to peak 
$|C_p|$ at the leading edge, and they induce a large viscous response at the wall, significantly impacting the onset of dynamic stall and its progression. Therefore, effective criteria for vortex shedding need to demonstrate critical behaviour around these events.
$C_p$ profiles on the airfoil suction surface for different 
$\alpha$ around the time of strong laminar vortex shedding. The profiles are shifted along the ordinate by 
$\Delta C_p = 2.25$ units to avoid clutter. The profile corresponding to relatively strong vortex shedding is highlighted in red.
Our prior studies at moderate-to-high 
$Re$ have shown that 
$BEF$ calculated by integrating over 1 % chord is sufficient to indicate imminent DSV formation. The small region of integration is justified because the effect of DSV formation is strongly felt close to the airfoil leading edge through the collapse of leading-edge 
$C_p$. In contrast, laminar vortex shedding, which occurs at low 
$Re$, is a spatially localized event not limited to the leading-edge region. Therefore, we calculate 
$BEF$ by integrating over different chord lengths to capture localized vortex-shedding events. The lowest integration length is set to 
$x=0.05$ since the spatial scales are larger at these low 
$Re$. 
$LESP$ is typically calculated by integrating up to the point of maximum thickness, according to the definition provided by Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020). For comparison, we integrate both parameters over a range of fractional chord lengths. The suction force used in the definition of 
$LESP$ given by (3.2) has a length scale 
$c$ in the denominator. We retain the same length scale 
$c$ while varying the integration length used. Setting this length scale to the integration length would simply scale the 
$LESP$ curves by a constant. We also use a low-pass filter to remove the high-frequency fluctuations in time that are captured in the LES for 
$LESP$, 
$BEF$ and 
$\max (|C_p|)$ near the leading edge. A Gaussian filter with a half-width equal to the inverse of a non-dimensional frequency 
$f^+ = f c / U_{\infty }$ of 20 is used.
Figure 7 plots the variation of 
$LESP$ and 
$|BEF|$ obtained by integrating over different chord lengths. Both parameters reach their respective global maxima between 
$13^\circ$ and 
$14^{\circ }$, which is ahead of DSV formation. This point reflects the leading edge 
$|C_p|$ reaching a maximum, following which DSV formation begins. In addition, two secondary peaks, corresponding to strong laminar vortices that are shed, are observed at 
$\alpha \sim 10^{\circ }$ for the 
$|BEF|$ curves with integration length 
$\ge 0.25$. Note that there is no change in 
$BEF$ across integration lengths unless there is an additional downstream source of vorticity at a given time. That is, the 
$BEF$ curves all coincide at initial times and vary only as the flow develops. This is in contrast to the 
$LESP$, where the values vary with integration length even at the initial times (note the vertical displacement of the 
$LESP$ curves at low 
$\alpha$). Another aspect of 
$BEF$ behaviour is the contribution from vortex shedding being in the same sense as the contribution from the leading edge (see figure 9 and the related description by Sudharsan et al. (Reference Sudharsan, Ganapathysubramanian and Sharma2022) for a detailed explanation on the sense of 
$BEF$ contributions over the airfoil).
(a) 
$LESP$ and (b) 
$|BEF|$ integrated over different chord lengths, plotted against 
$\alpha$, for Case R60-p05.
Since the effect of the laminar vortices is felt locally, they are indicated by the 
$BEF$ when the pertinent area is included in the integration region. Calculating 
$BEF$ over different integration lengths and comparing the difference can also identify the region over the airfoil where vortex shedding occurs. This is demonstrated in detail in § 3.6. The 
$LESP$ curves do not show a peak around the instance of laminar vortex shedding since only the chord/camber-wise component of the suction force features in the 
$LESP$ definition (see (3.2)). A vortex shed away from the leading edge has a more pronounced effect on the normal component (as opposed to the camber-wise component) of the suction force. The curves corresponding to 
$LESP$ and 
$|BEF|$ integrated up to less than 
$x=0.25$ show a slight decrease in slope around the instance of laminar vortex shedding. This is attributed to a temporary reduction in the slope of leading-edge 
$C_p$ magnitude due to the change in the curvature of the shear layer caused by the downstream vortex shedding.
These results demonstrate that the 
$BEF$ can be used to identify imminent vortex shedding occurring away from the airfoil leading edge.
3.3. Case R60-p25
Next, we discuss the results obtained at the same 
$Re$ of 60 000, albeit with a higher pitch rate of 0.25. The unsteady lag effects due to the higher pitch rate serve to delay the angle of attack at which stall occurs, in comparison with Case R60-p05. Moreover, the higher pitch rate also promotes a more conventional leading-edge stall, characterized by the breakdown of the LSB and a drop in leading-edge 
$|C_p|$. Space–time contours of 
$-C_p$ and 
$C_f$ are shown in figure 8. As in the previous case, the inflectional velocity profiles in the shear layer subjected to an APG develop an inviscid K-H instability, leading to vortex roll-up and transition to turbulence. This occurs close to 
$x=0.6$ initially, as pointed out in the 
$C_f$ contours (figure 8b).
Space–time contours for Case R60-p25. (a) 
$-C_p$ and (b) 
$C_f$.
Between approximately 
$15^{\circ }$ and 
$20^{\circ }$, a sequence of events transpires in rapid succession. Spanwise vorticity contours in figure 9 show the dynamic flow field during this time. In panel (a), the shedding of vortices due to K-H instabilities in the shear layer and transition to turbulence on the aft section of the airfoil are evident. Panel (b) shows a time instance after the establishment of the LSB with laminar reattachment. In the downstream region of this reattached flow, a second separated shear layer develops, experiencing K-H instabilities that lead to the roll-up of laminar vortices and downstream transition to turbulence. After its initial establishment, the LSB's existence is short-lived, as it quickly succumbs to viscous instabilities magnified by the increasing APG. During this time, strong vortices are shed from the rear of the LSB. Almost simultaneously, the K-H instabilities also undergo substantial amplification, giving rise to a series of shear layer vortices that roll up downstream, as shown in panel (c). This case serves as an illustration of a transition mechanism wherein both viscous instabilities within the LSB and the inviscid K-H instabilities are simultaneously amplified. This phenomenon is clearly discernible from the 
$C_p$ and 
$C_f$ contours in figure 8, where a horizontal row of vortices is observed shortly after the formation of the LSB. Similar behaviour has previously been observed in water tunnel experiments for a pitching-plunging airfoil undergoing deep stall at 
$Re=20\,000$ (see figure 9 of Ol et al. Reference Ol, Bernal, Kang and Shyy2009).
Case R60-p25. Vorticity contours showing (a–c) the simultaneous amplification of the K-H instability in the shear layer and the viscous instabilities within the LSB leading to the roll-up of vortices; (d–f) growth of the leading-edge DSV and shear layer vortices. Purple represents clockwise vorticity, while green represents counter-clockwise vorticity. (a) 
$\alpha = 15.0^{\circ }$, (b) 
$\alpha = 15.5^{\circ }$, (c) 
$\alpha = 16.7^{\circ }$, (d) 
$\alpha = 18.4^{\circ }$, (e) 
$\alpha = 19.6^{\circ }$ and (f) 
$\alpha = 20.7^{\circ }$.
Following the onset of instabilities, the LSB collapses and a coherent leading-edge DSV is formed. The subsequent growth of the DSV and the shear layer vortices are showcased in figure 9(d–f). The DSV begins to convect downstream as it grows, with a higher convection speed compared with that of the downstream shear layer vortices. It entrains a pair of previously shed shear layer vortices further downstream (marked in figure 8b). During its downstream propagation, the DSV maintains its proximity to the airfoil surface until 
$\alpha \sim 45^{\circ }$. A more comprehensive visual representation of this sequence is included in the accompanying video (see supplementary movie 2).
Figures 10(a)–10(c) show the variation of the unsteady aerodynamic coefficients with 
$\alpha$. Moment stall occurs at 
$\alpha \sim 25^{\circ }$ when the DSV convects past the quarter-chord point. Maximum lift occurs at 
$\alpha \sim 33^{\circ }$, which is marked as the lift stall point in the figure. However, the lift remains elevated until 
$45^{\circ }$, attributed to the proximity of the DSV to the airfoil surface. As the DSV progressively convects away from the airfoil surface, the reduction in vortex-induced contribution to the bound circulation leads to a decline in 
$C_l$. Figure 10(d) shows the variation of peak leading-edge 
$|C_p|$. The collapse of the LSB is accompanied by a reduction 
$\max (|C_p|)$ near the leading edge at 
$\alpha \sim 18^{\circ }$, as is typical in a bubble-bursting, leading-edge stall.
Variation with 
$\alpha$ of (a–c) aerodynamic coefficients and (d) 
$\max (|C_p|)$ near the first 5 % of airfoil chord, for Case R60-p25.
The variation with 
$\alpha$ of 
$LESP$ and 
$BEF$ is shown in figure 11. Both parameters exhibit distinct peaks for all integration lengths, attributed to the decrease in leading-edge 
$|C_p|$ as the LSB breaks down. The larger values of 
$|BEF|$ for integration lengths greater than 
$0.05$ are due to the large vorticity and vorticity flux associated with the vortices shed downstream of the LSB as it becomes unstable. Due to the greater spatial distance of the LSB from the geometric leading edge and its extended length, DSV formation occurs at 
$x=0.1$ rather than very close to the leading edge. This is consistent with the expected behaviour at low 
$Re$ (Gaster Reference Gaster1967).
(a) 
$LESP$ and (b) 
$|BEF|$ integrated over different chord lengths, plotted against 
$\alpha$, for Case R60-p25.
Figure 11 shows that the instance of the peak in 
$LESP$ moves aft (to higher 
$\alpha$) with increasing integration length. This is due to the increasing contribution from the suction induced by the DSV as it grows and influences a larger portion of the airfoil between the leading edge and the maximum thickness point. Hence, even though the leading-edge 
$C_p$ magnitude begins to drop in advance of 
$20^{\circ }$, the 
$LESP$ curve obtained by integrating up to 
$x=0.25$ continues to rise until approximately 
$22^\circ$. In contrast, the peak 
$|BEF|$ location occurs before 
$20^{\circ }$ for all integration lengths. The leading edge 
$C_p$ peak representing stall onset is captured by the 
$BEF$ curve corresponding to integration up to 
$0.05$. For integration lengths larger than 
$x=0.05$, the peak 
$BEF$ value is much higher due to the strong vortices shed from the rear of the LSB around the same time. The vortex shedding event is captured by the 
$BEF$ for this case as well when appropriate regions of the airfoil are included. There is another peak at approximately 
$25^{\circ }$ for integration lengths 
$\geq 0.25$, which corresponds to an increase in induced secondary vorticity as the DSV gets stronger as it grows. While 
$LESP$ captures the stall onset point, the localized vortex shedding event is missed.
3.4. Case R10-p05
For the lower 
$Re$ cases, the DSV remains laminar through nearly the entire manoeuvre. These cases are also characterized by a DSV system consisting of several large-scale laminar vortices.
Figure 12 shows the vorticity contours at specific instances. A region of reversed flow at the trailing edge propagates upstream as the airfoil starts to pitch up. At 
$\alpha \sim 10^{\circ }$, the separated shear layer from the leading edge becomes susceptible to instabilities and sheds spanwise vortices. Panel (a) (
$\alpha \sim 11^{\circ }$) shows the initial stages of the unstable, separated shear layer. The instability location moves upstream as the airfoil continues to pitch up. The large-scale shear layer vortices induce secondary vorticity beneath them. The induced secondary vorticity develops into coherent, counter-clockwise vortices, which act to cut off the primary shear layer vortex that is formed at approximately mid chord.
Vorticity contours showing the formation of a DSV system through shear layer vortex and wall interactions for Case R10-p05. (a) 
$\alpha = 11.0^{\circ }$, (b) 
$\alpha = 13.3^{\circ }$, (c) 
$\alpha = 14.8^{\circ }$, (d) 
$\alpha = 15.7^{\circ }$, (e) 
$\alpha = 18.3^{\circ }$ and (f) 
$\alpha = 20.3^{\circ }$.
Panel (b) (
$\alpha \sim 13.3^{\circ }$) shows the first instance of secondary vorticity (SVa) cutting off the primary shear layer vortex (PVa), and the roll-up of the shear layer upstream into another vortex (PVb). The shedding of shear layer vortices and induced secondary vortices result in multiple instances of vortex entrainment. Therefore, the DSV system in the present case comprises multiple shear layer vortices entrained together as 
$\alpha$ increases. Note that the DSV system is centred downstream of the quarter-chord point of the airfoil during its incipience. Panels (c) (
$\alpha \sim 14.8^{\circ }$) and (d) (
$\alpha \sim 15.7^{\circ }$) show the merging of PVa and PVb into a single coherent vortex and its downstream convection, respectively. During the later part of the unsteady motion, as shown in panels (e,f), the DSV, continuing to entrain more leading-edge vortices, transitions to turbulence as it propagates downstream. The DSV weakens due to the entrained counter-clockwise vorticity (shown in green contours) and viscous dissipation. The dominant, large-scale vortices located away from the surface of the airfoil generate and influence the motion of smaller-scale secondary vortices and, over time, tertiary vortices close to the wall, that affect surface quantities and make it challenging to interpret the space–time plots of 
$C_p$ and 
$C_f$. Supplementary movie 3 provides a clearer picture of the sequence of events as the airfoil pitches up.
Figures 13(a)–13(c) show the variation with 
$\alpha$ of the unsteady aerodynamic coefficients. Large wiggles are observed in 
$C_l$, 
$C_d$ and 
$C_m$ before stall due to the strong laminar vortices shed from the shear layer. Here, 
$C_l$ continues to increase until 
$\alpha \sim 18.8^{\circ }$ due to vortex-induced suction/lift from multiple shear layer vortices. However, large fluctuations in 
$\max (|C_p|)$ near the leading edge (shown in figure 13d) are observed much earlier as large variations in shed vorticity occur downstream. The peaks and the valley in 
$\max (|C_p|)$, identified by red x markers in figure 13, are due to the vortices highlighted in figure 12(b–d). This undulatory behaviour of 
$\max (|C_p|)$ is attributed to the changes in the curvature of the shear layer at the leading edge due to the vortices shed downstream. As 
$\alpha$ increases beyond 
$12.5^\circ$ and vortex shedding occurs around mid-chord, the induced secondary vorticity pushes the shear layer away from the airfoil suction surface, reducing its curvature and causing 
$\max (|C_p|)$ to dip. After 
$15^\circ$, the (PVa + PVb) DSV system detaches from the shear layer, arresting the movement of the shear layer away from the surface, thereby temporarily but sharply increasing 
$\max (|C_p|)$ near the leading edge.
Variation with 
$\alpha$ of (a–c) aerodynamic coefficients and (d) 
$\max (|C_p|)$ near the first 5 % of airfoil chord, for Case R10-p05. The red x markers correspond to the flow fields highlighted in panels (b–d) in figure 12.
The variation with 
$\alpha$ of 
$LESP$ and 
$BEF$ integrated over various chord lengths is shown in figure 14. The shedding of the downstream shear layer vortices that affects leading-edge 
$C_p$ also impacts the variation in 
$LESP$ and 
$BEF$. For the 
$BEF$, integrated lengths up to 
$0.25$ follow the same trend as 
$\max (|C_p|)$ near the leading edge. The larger integrated lengths show larger 
$|BEF|$ values at the peak locations, corresponding to the vorticity shed downstream. The increase in 
$BEF$ around 
$13^{\circ }$ corresponds to increased induced secondary vorticity from the roll-up of PVa. The increase to the second peak around 
$15.0^{\circ }$ corresponds to the entrainment of PVa by PVb during which PVa moves closer to the wall, causing increased secondary vorticity. The small increase to another peak around 
$18.3^{\circ }$ corresponds to another instance of vortex merging with the DSV system, which again induces increased secondary vorticity. The 
$LESP$ curves for all integration lengths show nearly the same trend as 
$\max (|C_p|)$ near the leading edge. Note that there is no additional local (in space and time) contribution to 
$LESP$ from the vortices formed around mid-chord; the curves show the same trend and are almost simply displaced along the ordinate. The largest integration length shows a deviation in trend at high 
$\alpha$ since the chord-wise force contribution from the shed vortices is directed in the opposite direction over the rear of the airfoil. It is the local contributions around vortex formation that distinguish the 
$BEF$ from the 
$LESP$. This is elucidated in § 3.6 by delineating specific contributions from different airfoil sections.
(a) 
$LESP$ and (b) 
$|BEF|$ integrated over different chord lengths, plotted against 
$\alpha$, for Case R10-p05.
3.5. Case R10-p25
The higher pitch rate case at lower 
$Re$ is also characterized by the shedding of multiple shear layer vortices with the flow over the airfoil remaining laminar throughout. The difference from the lower pitch rate case (R10-p05) is the formation of a stronger DSV system further upstream with more pronounced lag effects due to unsteadiness. Supplementary movie 4 clearly illustrates the sequence of events during the unsteady manoeuvre using spanwise vorticity contours.
The sequence of events remains similar to Case R10-p05 but is postponed to higher 
$\alpha$ due to the pronounced unsteady lag. Figure 15 shows span-averaged vorticity contours at a few instances during the pitch-up manoeuvre. The upstream propagation of trailing edge reversed flow and the development of instabilities in the shear layer leading to vortex formation are observed in this case as well. However, due to the larger APG encountered by the flow, stronger leading-edge vortices are shed farther upstream (at approximately 
$x=0.12$), beginning at 
$\alpha \sim 20^{\circ }$ (figure 15a). Induced secondary (counter-clockwise) vorticity acts to cut off the shear layer, which forms clockwise vortices near the airfoil leading edge (figure 15b–d). When the clockwise vortex (marked as ‘PVa’) rolls up, it cuts off the downstream shear layer from the leading edge, leading to its roll-up into vortices downstream (as pointed out in figure 15b). The secondary vorticity itself lifts up due to the induction from the leading edge vortices, rolls up and is pinched off by the clockwise vorticity induced by it (figure 15c). This process repeats a few times, and multiple vortices are shed from the leading edge. Similar to Case R10-p05, the DSV system continues to entrain these shed vortices and grows in size. DSV-induced suction on the airfoil surface is strongest at 
$\alpha \sim 33^{\circ }$ (figure 15e). The effect of the DSV system on the surface drops as it moves farther away and becomes more diffuse due to viscosity and annihilation from the entrained secondary vorticity (figure 15f). The DSV system also entrains the downstream shear layer vortex around this time.
Vorticity contours at different instances during the pitch-up manoeuvre for Case R10-p25. Similar behaviour as in Case R10-p05, but is delayed (in 
$t^*$ and 
$\alpha$) due to the higher pitch rate, and vortex shedding occurs closer to the leading edge. (a) 
$\alpha = 20.2^{\circ }$, (b) 
$\alpha = 22.8^{\circ }$, (c) 
$\alpha = 25.7^{\circ }$, (d) 
$\alpha = 29.7^{\circ }$, (e) 
$\alpha = 33.2^{\circ }$ and (f) 
$\alpha = 49.7^{\circ }$.
Figures 16(a)–16(c) show the variation of the unsteady aerodynamic coefficients. Maximum lift occurs around 
$\alpha \sim 33^{\circ }$, when the DSV system induces the strongest suction over the airfoil surface. The lift continues to remain high until approximately 
$45^{\circ }$, owing to the DSV system remaining relatively close to the airfoil surface and continuing to interact with the newly shed leading-edge vortices and induced secondary vorticity. The magnitude of leading edge 
$C_p$ shows undulations (see figure 16d) due to the effect of vortices shed downstream, similar to the previous case.
Variation with 
$\alpha$ of (a–c) aerodynamic coefficients and (d) 
$\max (|C_p|)$ near the first 5 % of airfoil chord, for Case R10-p25.
Figure 17 shows the variation of 
$LESP$ and 
$|BEF|$ for the present case for different integration lengths. Both sets of curves are characterized by a series of peaks for all integration lengths. The peak corresponding to the maximum leading edge 
$|C_p|$ is clearly observed from 
$BEF$ at the lowest integration length. In contrast to the previous case, the newly shed shear layer vortices remain very close to the leading edge (for example, see PVb in figure 15b). Therefore, the contribution to 
$LESP$ from the suction induced by these vortices is significant, leading to multiple peaks. A series of significant spikes in 
$|BEF|$ is observed for larger integration lengths. The first instance of such a peak around 
$22^{\circ }$ corresponds to the roll-up of vortex PVa (shown in figure 15a). The subsequent instances correspond to the DSV system entraining a new leading edge vortex, during which it moves closer to the airfoil surface, inducing strong secondary vorticity. The trend of the 
$BEF$ localizing the points corresponding to peak 
$|C_p|$ and DSV-induced secondary vorticity peaks holds for the present case as well.
(a) 
$LESP$ and (b) 
$|BEF|$ integrated over different chord lengths, plotted against 
$\alpha$, for Case R10-p25.
3.6. Comparison of results
We next distinguish between the 
$\max (LESP)$ and 
$\max (|BEF|)$ criteria. We divide the airfoil into multiple sections, as shown in figure 18, and identify contributions from different regions of the airfoil to the two parameters. Since the definition of 
$LESP$ involves a square root, we instead compare the suction force coefficient, 
$C_{suction}$ (see (3.2)). The contributions to both 
$C_{suction}$ and 
$BEF$ are primarily negative, so we plot the contributions to the negative of these quantities.
Different sections of the airfoil referred to in subsequent figures.
Figure 19 shows the contributions to 
$C_{suction}$ and 
$BEF$ from the airfoil sections marked in figure 18 for Case R60-p05. The contribution from the pressure side (PS) to either parameter is negligible. The largest contributions are from the leading edge region (marked ‘LE’) for both 
$C_{suction}$ and 
$BEF$. The Post-LE and THK regions contribute to 
$C_{suction}$ throughout the unsteady manoeuvre. The effect of increasing (or decreasing) leading edge 
$|C_p|$ is felt at these regions since there is a modest chord-wise component of the surface normal. In contrast, the contribution to 
$BEF$ from variation in leading-edge 
$|C_p|$ is limited to the LE region. The Post-LE region does not contribute to 
$BEF$ because, even though 
$|C_p|$ varies with 
$\alpha$ in this region, the pressure gradient (similar behaviour to vorticity flux) drops to zero. The local peaks in the 
$BEF$ curves that occur before the peak in the LE curve correspond to local vortex shedding, while those following it correspond to the DSV convecting downstream and growing stronger. As noted in § 3.2, strong laminar vortices are shed from the THK and MID regions at 
$\alpha \sim 10^{\circ }$. These are reflected as peaks in the corresponding 
$BEF$ curves in figure 19(b). A strong vortex shed from the rear of the LSB is also captured in the 
$BEF$ curve corresponding to the THK region.
Contributions to (a) 
$C_{suction}$ and (b) 
$BEF$ from different airfoil segments, for Case R60-p05.
At higher 
$\alpha$, the MID and AFT regions show an opposite contribution to 
$C_{suction}$ relative to other regions. This arises from the convection of the DSV centre downstream, which induces high 
$|C_p|$ in the MID and AFT regions where the chord-wise component of the surface normal points in the 
$+x$ direction. The peaks in 
$BEF$ from the MID and AFT regions at higher angles of attack (
$\alpha \sim 16.5^{\circ }$ and 
$20^{\circ }$) occur due to strong secondary vorticity induced by the DSV as it grows stronger while convecting downstream. These peaks occur later in time for 
$LESP$ compared with 
$BEF$ since the chordwise component of the surface normal is more significant towards the rear of the airfoil.
Case R10-p25 is also briefly presented in a similar way. Figure 20 shows the contributions to 
$C_{suction}$ and 
$BEF$ from different airfoil sections shown in figure 18. The contribution to 
$C_{suction}$ from LE includes multiple peaks, due to the leading edge vortices remaining farther upstream for this case. In contrast, the LE contribution to 
$BEF$ shows a clear peak, following the same trend as 
$\max (|C_p|)$ near the leading edge (see figure 16d). Specific instances of vortex roll-up and vortex entrainment, which induce increased secondary vorticity into the DSV system are observed from the other regions. These include the roll-up of PVa and PVb, the merging of PVa and PVb, and the entrainment of further leading edge vortices and the downstream shear layer vortex into the DSV system. The contributions to 
$C_{suction}$ from the Post-LE and THK regions are influenced by the variation of leading edge 
$C_p$ and leading edge vortices. Opposite contributions to 
$C_{suction}$ due to DSV convection are captured in the MID and AFT regions. A similar breakdown of contributions to 
$C_{suction}$ and 
$BEF$ for the remaining cases, namely R60-p25 and R10-p05, is presented in Appendix B.
Contributions to (a) 
$C_{suction}$ and (b) 
$BEF$ from different airfoil segments, for Case R10-p25.
Figure 21 illustrates the overall picture that emerges from the results of the four simulations, including the performance of 
$LESP$ and 
$BEF$ in identifying critical flow events. At 
$Re = 10\,000$ (cases R10-p05 and R10-p25), the flow remains laminar almost through the entire manoeuvre. There is a large region of reversed flow, and the 
$C_f$ is negative over nearly the entire suction surface of the airfoil. Due to the larger viscous response, the DSV system comprises multiple large-scale laminar vortices. These vortices are shed further upstream (close to the airfoil leading edge) for the high pitch-rate case.
Illustration of observations from all four cases considered in the present study.
Cases with 
$Re = 60\,000$ (R60-p05 and R60-p25) demonstrate a more typical LSB-bursting, leading-edge stall. During the establishment of the LSB, strong, small-scale vortices are shed. For the higher pitch rate case, the DSV forms farther upstream and is stronger, as reported in the literature (Acharya & Metwally Reference Acharya and Metwally1992). The DSV is fed by the vorticity shed from the leading edge, but due to the higher 
$Re$, it does not organize into large coherent vortices before merging with the DSV.
Of the two stall criteria under investigation, 
$\max (LESP)$ mainly captures variations in leading edge 
$C_p$. The 
$\max (|BEF|)$ criterion captures the variation in leading edge 
$C_p$ and also enables the localization of vortex shedding events occurring anywhere on the airfoil surface in space and time. This includes strong, small-scale vortex shedding that is not reflected in the leading edge 
$C_p$ as well as large-scale vortex shedding events located away from the leading edge. These events are directly captured in the spatial region where they occur. Specific instances of vortex roll-up and entrainment, which trigger an increase in induced secondary vorticity, are captured by the 
$BEF$. The versatility of the 
$BEF$ parameter makes it a suitable candidate for flow control at low 
$Re$.
4. Conclusion
We investigate airfoil dynamic stall at low Reynolds numbers where laminar, transitional and turbulent regimes can coexist, giving rise to rich fluid dynamics. The dynamic stall process is strongly influenced by multiple vortices that shed from different regions of the airfoil not limited to the leading edge. Current state-of-the-art stall onset criteria, based on the 
$LESP$ and 
$BEF$ parameters, are calculated by integrating these quantities around the leading edge and hence may not directly capture vortex shedding occurring away from it.
We evaluate the 
$\max (LESP)$ and 
$\max (|BEF|)$ criteria over an extended integration region for low-
$Re$ (
${\sim }{O}(10^4)$) dynamic stall for an LES dataset consisting of the SD7003 airfoil undergoing a constant-rate, pitch-up manoeuvre at two 
$Re$ values and two pitch rates. The highly resolved LES results provide insights into the unsteady flow phenomena (instabilities, shear layer dynamics, vortex formation, pairing, shedding, dissipation, etc.), which are commented on. The ability of the 
$LESP$ and 
$BEF$ parameters to capture vortex-shedding events is analysed using these results.
Our analyses indicate that the 
$\max (|BEF|)$ criterion captures strong vortex-shedding events when the integration region includes the locations of these events. These events are reflected as local maxima in the 
$BEF$ curves. The 
$LESP$, being based on the camber-wise component of pressure force, captures only the effect of vortex shedding events that significantly affect 
$C_p$ near the leading edge. The 
$BEF$ parameter differs from 
$LESP$ in that it directly captures the flow events such as (a) strong, small-scale vortex shedding that do not significantly impact leading edge 
$C_p$, and (b) large-scale shed vortices, localizing them in space and time. It captures instances of vortex roll-up and entrainment, which induce secondary vorticity. The use of 
$BEF$ to localize stall onset and vortex-shedding events in space and time holds promise for flow control.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.753.
Acknowledgements
We acknowledge the computational resources provided by Argonne Leadership Computing Facility and Iowa State University.
Funding
This material is based upon work supported by the National Science Foundation (Grants CBET-1935255 and 1554196) and the US Air Force Office of Scientific Research (Award # FA9550-23-1-0016).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Static simulations
Static simulations were performed using FDL3DI at 
$\alpha = 4^{\circ }$ for both 
$Re$; only the results for 
$Re=60\,000$ are shown for brevity. The simulations were run for approximately 
$40$ characteristic convective times (
${c}/{U_{\infty }}$) to ensure that the forces and moments reached statistical stationarity. Comparisons of surface pressure and skin friction coefficient distributions, 
$C_p$ and 
$C_f$, respectively, are made with XFOIL (Drela Reference Drela1989). XFOIL results are obtained with the Ncrit parameter set to 
$9$, corresponding to a low freestream-turbulence intensity. Figure 22 shows good agreement between LES and XFOIL results; the transition on the suction surface occurs around 
$60\,\%$ chord in both.
Comparison of 
$C_p$ and 
$C_f$ distributions between static LES and XFOIL for an 
$Re$ of 60 000 at 
$\alpha = 4^{\circ }$. (a) Surface pressure coeff., 
$C_p$ and (b) skin friction coeff., 
$C_f$.
Appendix B. Contributions to $C_{suction}$ and $BEF$ from different airfoil regions
The contributions to 
$C_{suction}$ and 
$BEF$ from different sections of the airfoil (see figure 18) for Cases R60-p25 and R10-p05 are shown in figure 23. Critical flow events identified by the two parameters are annotated in the figure.
Contributions to 
$C_{suction}$ and 
$BEF$ from different airfoil segments for cases (a,b) R60-p25 and (c,d) R10-p05. (a) 
$C_{suction}$, R60-p25, (b) 
$BEF$, R60-p25, (c) 
$C_{suction}$, R10-p05 and (d) 
$BEF$, R10-p05.
Open access
