1. Introduction
Dynamic stall over wind turbine and helicopter blades has been a problem of great interest to fluid dynamicists owing to its severity compared to static stall, with the possibility of consequent structural damage. Dynamic stall models such as the Beddoes/Leishman model (Leishman & Beddoes Reference Leishman and Beddoes1989) provide a first-order, flow-physics-based framework for modelling unsteady flow using semi-empirical formulations. To formulate these models, various criteria for dynamic stall onset based on the unsteady aerodynamic coefficients have been explored, such as deviation of the normal force coefficient, drop in pitching moment coefficient, maximum value of chord-wise force coefficient, suction collapse at the leading edge (Sheng, Galbraith & Coton Reference Sheng, Galbraith and Coton2005). Currently, the most suitable standalone criterion that could be utilised for characterising dynamic stall onset is the maximum of the leading edge suction parameter (
$LESP$) (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014), even though its intended use was to initiate leading edge vortex shedding in a reduced-order model when 
$LESP$ reaches a critical value (details in § 3.3). However, we find that 
$LESP$ calculated from computational fluid dynamics (CFD) reaches its maximum value well after stall onset is indicated by the onset of instabilities in the laminar separation bubble (LSB). This offers scope for exploring more conservative criteria, as it is well-known that control efforts applied to dynamic stall (blowing, suction, etc.) are most effective in mitigating stall before the formation of the dynamic stall vortex (DSV) (Chandrasekhara Reference Chandrasekhara2007). Further, a consistent definition of ‘stall onset’ is currently not found in the literature. Therefore, our first step in the current work is to define a stall onset regime based on control effectiveness. Next, we propose a more physically-intuitive, vorticity-based parameter, namely, the boundary enstrophy flux (
$BEF$), to characterise stall onset. We present its advantages and demonstrate that it can serve as a more conservative criterion compared to 
$LESP$.
When the angle of attack of an aerofoil undergoing a pitch-up motion exceeds its static stall angle of attack, 
$\alpha _{ss}$, a combination of inviscid (increase in effective camber and Magnus effect at the leading edge) and viscous (LSB formation) effects allows for increasing lift beyond 
$\alpha _{ss}$ (Corke & Thomas Reference Corke and Thomas2015). The eventual stall could be the result of LSB bursting or breakdown, flow separation that begins at the trailing edge, or boundary-layer shock interactions in high-Mach-number flows. In a classic leading-edge-type stall, leading edge suction collapses when the LSB suddenly breaks down. The ejected vorticity coalesces into a coherent DSV, whose convection away from the aerofoil surface results in severe stall. Widmann & Tropea (Reference Widmann and Tropea2015) found that the DSV can be cut off from the leading edge shear layer by either a boundary layer eruption mechanism, where a sharply focused concentration of vorticity erupts into the main flow (Doligalski Reference Doligalski1994), or a bluff body detachment mechanism, where a region of reverse flow grows from the trailing edge upstream and cuts off the DSV. The reviews by Corke & Thomas (Reference Corke and Thomas2015) and Carr (Reference Carr1988) provide a detailed breakdown of the stages involved in dynamic stall.
Though the sequence of events leading to dynamic stall could be different depending on the operating conditions and aerofoil geometry, Benton & Visbal (Reference Benton and Visbal2018, Reference Benton and Visbal2019) found that the LSB was crucial in the development of the DSV for chord-based Reynolds numbers as high as 
$10^6$. Based on the myriad of flow control studies focusing on the leading edge that have been found to be effective (Beahan et al. Reference Beahan, Shih, Krothapalli, Kumar and Chandrasekhara2014; Zhao & Zhao Reference Zhao and Zhao2015; Choudhry, Arjomandi & Kelso Reference Choudhry, Arjomandi and Kelso2016), understanding flow physics near the leading edge is crucial for characterising stall onset. This necessitates first providing a consistent definition for stall onset.
McCroskey (Reference McCroskey1981) delineated four regimes of flow in the context of dynamic stall, out of which stall onset is the regime where maximum lift is produced, without excessive drag or pitching moment. On the other hand, Mulleners & Raffel (Reference Mulleners and Raffel2012) considered stall onset to occur when the DSV detaches from the leading edge. In the present work, we consider stall onset to be the period of time between the first appearance of instabilities in the coefficient of friction (
$C_f$) profile of the LSB and the formation of the DSV. The task of characterising stall onset is made non-trivial by the wide range of parameters (
$Re_c$, Mach number, pitch rate, etc.) affecting the type (leading edge, trailing edge or mixed) of stall.
Some researchers have used modal decomposition to analyse stall onset. For example, Mulleners & Raffel (Reference Mulleners and Raffel2012) used proper orthogonal decomposition (POD) on a sinusoidally pitching aerofoil and found that the temporal coefficient corresponding to the DSV growth mode reaches a maximum when the DSV detaches from the leading edge (which correlates with their definition of stall onset). Lagrangian approaches (Peacock & Haller Reference Peacock and Haller2013) – such as using finite-time Lyapunov exponent (FTLE) fields – have also been used to study dynamic stall (Wang et al. Reference Wang, Cao, Dang, Zhang and Deguchi2021). However, most of these studies used the entire pitching cycle for their calculations and were not focused on resolving localised instabilities around stall onset, as in the present work.
Since flow dynamics near the leading edge plays an important part in determining stall onset, it is reasonable to analyse parameters reaching criticality at the leading edge as candidates for characterising stall onset. Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) used 
$LESP$ to modulate leading edge vortex shedding in their reduced-order model based on unsteady thin aerofoil theory. The critical value of 
$LESP$, predetermined from CFD based on observed events, was cited to be independent of aerofoil kinematics for a given geometry and Reynolds number. Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016) found that 
$LESP$ continues to increase even after the onset of instabilities in the 
$C_f$ profile near the leading edge at low Reynolds numbers. They also found that the maximum value of 
$LESP$ varies with pitch rate for chord-based Reynolds numbers, 
$Re_c$, of 30 000 and 
$3\times 10^6$. The same conclusion was reached by Deparday & Mulleners (Reference Deparday and Mulleners2019) based on experiments at 
$Re_c = 10^6$. The invariance of the maximum value of 
$LESP$ with motion kinematics for a given aerofoil geometry and 
$Re_c$ was found to hold for low Reynolds numbers and higher pitch rates (Ramesh et al. Reference Ramesh, Granlund, Ol, Gopalarathnam and Edwards2018; Hirato et al. Reference Hirato, Shen, Gopalarathnam and Edwards2021). Deparday & Mulleners (Reference Deparday and Mulleners2018) identified a primary instability stage and vortex formation stage, by a maximum in the leading edge suction between the two stages, also delineated by a change in shear layer growth rate.
Vorticity-based approaches have yielded promising physical insights in several fundamental fluid mechanical problems. Lighthill (Reference Lighthill1986) used the counterplay of the convection and diffusion of vorticity within the boundary layer to provide an intuitive explanation for flow separation from a surface. Reynolds & Carr (Reference Reynolds and Carr1985) explained dynamic stall as the result of a disturbance in the balance between the outer potential flow and the accumulating boundary layer vorticity of the LSB. Acharya & Metwally (Reference Acharya and Metwally1992) found that the amount of vorticity fed into the DSV through the surface vorticity flux was positively correlated with increasing pitch rate and resulted in a stronger, more compact DSV. Akkala (Reference Akkala2015) calculated the vorticity budgets for flow over a plunging flat plate and concluded that the diffusive flux of vorticity from the surface correlated well with flow evolution, increasing as long as the DSV remained attached. Gehlert & Babinsky (Reference Gehlert and Babinsky2021) found that the value of vortex sheet strength on an accelerating and rotating cylinder at the unsteady separation point collapses during cylinder motion when the effect of rotation, instantaneous velocity and far-field shed vorticity is removed. In the current paper, we take a vorticity-based approach to characterising the onset of incompressible, leading-edge-type stall.
2. Methods
The focus of the present work is on the analysis of numerical results of dynamic stall simulations performed using the flow solver FDL3DI. Initial solutions from simulations previously carried out at the US Air Force Research Laboratory were repeated with much higher temporal resolution around the stall onset regime.
The governing equations and numerical methods used are available in Sharma & Visbal (Reference Sharma and Visbal2019) and Visbal & Gaitonde (Reference Visbal and Gaitonde2002). Briefly, the implicit large-eddy simulation (LES) solver, FDL3DI, is used to solve the compressible Navier–Stokes equations. The effect of subgrid scale stresses is modelled by spatial filtering to remove energy at unresolved scales. Spatial discretisation uses a sixth-order compact finite difference scheme, while time integration is carried out through an implicit, approximate factorisation technique. All the quantities presented in the paper are appropriately non-dimensionalised using the aerofoil chord, 
$c$, and free stream velocity, 
$U_{\infty }$, as the reference length and velocity scales, respectively. Aerofoils undergoing a pitch-up motion at a constant, non-dimensional pitch rate of 
$-0.05$ are simulated at chord-based Reynolds numbers (
$Re_c$) 200 000–500 000 and Mach number 0.1. We study two different aerofoils at two Reynolds numbers for a total of four cases, all of which undergo leading-edge-type stall. The symmetric NACA0012 aerofoil has a maximum thickness of 12 % chord at 
$x/c=0.3$. The SD7003 aerofoil has a maximum thickness of 8.5 % chord at 
$x/c=0.24$ and a camber of 1.4 % chord. The presented results correspond to the ‘Fine’ mesh described in the appendix of Sharma & Visbal (Reference Sharma and Visbal2019), having 410, 1341 and 134 grid points, respectively, in the surface-normal, tangential and spanwise directions. Ten per cent of the aerofoil chord length is simulated in the span direction, with periodic boundary conditions. The current analysis is based on two-dimensional (2-D) flow fields obtained by averaging along the span.
Table 1 shows the datasets used in the current analysis, along with the acronyms that will be used to refer to them. The letters refer to the aerofoil (either ‘SD’ for SD7003 or ‘NC’ for NACA0012), and the suffixed numbers refer to the Reynolds number in thousands. For example, the SD7003 aerofoil case at 
$Re_c$ 200 000 is shortened to ‘SD200’, and the NACA0012 aerofoil at 
$Re_c$ 500 000 is shortened to ‘NC500’. In all cases, the aerofoil undergoes a constant-rate pitch-up motion, at a non-dimensional pitch rate, 
$\varOmega _0^{+} = \varOmega _0 c/ U_{\infty }$, of 
$-0.05$, pivoted about the quarter-chord point. The time integration was advanced by a non-dimensional time step size of 
$2 \times 10^{-5}$. The span-averaged solution was recorded around the time of stall onset every 100 time steps, providing highly sampled datasets to resolve the events within the stall onset regime. This covers about 
$3^\circ$ of the total pitch-up motion, with data available every 
$0.006^\circ$. For one of the cases (NC200), the span-averaged solution was recorded every 2000 time steps for the entire pitch-up manoeuvre. This is a different realisation of case NC200 from the one where the solution was frequently recorded around stall onset.
Table 1. CFD datasets used in present work. All cases undergo a pitch-up motion pivoted at 
$x/c=0.25$, at a non-dimensional pitch rate of 
$-0.05$ and Mach number 0.1.
3. Results and discussion
Space–time contours of coefficients of pressure and friction (
$C_p$ and 
$C_f$, respectively) on the surface of the aerofoil are illustrative of the different flow phenomena that occur as the aerofoil undergoes an unsteady manoeuvre. Figure 1 shows space–time diagrams on the suction surface of the aerofoil for the four cases investigated here. Figures 1(a) and 1(b) show contours of 
$C_p$ and 
$C_f$, respectively, over the entire pitch-up manoeuvre for NC200. Relevant flow events are marked in the figures. The point of dynamic stall onset can be estimated roughly with such plots; it is associated with LSB bursting in these cases. Post stall onset, the flow characteristics are dominated by the DSV. The static stall angle, 
$\alpha _{ss}$, for this case was found to be 
$10.9^{\circ }$ from XFOIL (see table 3). The log of the amplification factor of the most-amplified frequency that triggers transition (
$N_{crit}$) was set to 11, which is suitable for clean or very low inflow turbulence wind tunnel flow. This is consistent with our simulations with zero inflow turbulence. At low angles of attack, the upstream propagation of the transition location is accompanied by the occurrence of 2-D instability modes, until the formation of the LSB around 
$11.2^{\circ }$. The recirculation region within the LSB grows stronger with increasing 
$\alpha$ (see large values of 
$C_p$ in figure 1a), inducing oppositely-signed vorticity beneath it, leading to its eventual burst, around 
$17^{\circ }$. The LSB coalesces into a coherent vortex structure, namely, the DSV, which grows and convects downstream (blue region moving downstream in figure 1b). There is some flow reversal at the trailing edge moving upstream, at the right edge of figure 1(b). However, this region does not interact with the LSB for this case, and the formation of the DSV and subsequent dynamic stall are clearly dependent only on flow behaviour at the aerofoil leading edge.
Figure 1. Space–time diagrams of pressure and friction coefficients (
$C_p$ and 
$C_f$, respectively) on the suction surface of aerofoils undergoing a constant-rate pitch-up manoeuvre. Results from the entire manoeuvre are shown in (a,b) for NC200, while (c–e) show space–time contours of 
$C_f$ focusing on the stall onset regimes of SD200, SD500 and NC500, respectively.
Since the focus of this paper is on investigating stall onset, panels (c), (d) and (e) of figure 1 show the simulation results over shorter time periods that correspond to the stall onset regimes for SD200, SD500 and NC500, respectively (see table 1 for the nomenclature used here). The widening of the blue region (negative 
$C_f$) corresponds to DSV formation, growth and downstream advection. The thickening of the shear layer and subsequent formation and convection of the DSV can be relatively gradual, as in SD200 (panel (c)) and NC500 (panel (e)), or abrupt, as in SD500 (panel (d)). From the 
$C_f$ contours, we also observe a stronger vortex formation for SD500 compared to NC500, though both have the same 
$Re_c$. We attribute this to the difference in the leading edge shape between the two aerofoils.
Figure 2 shows contours of the pressure coefficient, 
$C_p$, and vorticity, 
$\omega$, at two different angles of attack, one before DSV formation (panel a) and one afterwards (panel b), for NC200. The region of peak suction (corresponding to large magnitudes of 
$C_p$), which is initially confined to the leading edge, moves downstream after the LSB coalesces into a coherent vortex as shown. This corresponds to the negative 
$C_f$ (blue) region in the space–time diagram (figure 1b) moving progressively downstream with increasing angle of attack. In the bottom row of figure 2(b), oppositely-signed vorticity (red) induced by the DSV is observed, which eventually cuts off the DSV. Additional flow images are presented in Appendix A.
Figure 2. Contours of pressure coefficient, 
$C_p$, and vorticity, 
$\omega$, (a) before and (b) after DSV formation for NC200.
3.1. Stall onset regime definition
In order to characterise dynamic stall onset, it is necessary to define what it constitutes, based on observed flow features. Our definition of stall onset is based on control effectiveness, which decreases significantly after the formation of the DSV (Chandrasekhara Reference Chandrasekhara2007). This sets the upper bound for our stall onset regime as the instant when the DSV is formed. However, it is not possible to find this point exactly; in fact, DSV formation might not happen instantaneously in all but the most extreme sudden stall cases. Therefore, the upper bound is found by visual analysis of streamlines and 
$C_p$ distributions.
For the lower bound of stall onset, we consider the first time at which spikes reaching zero appear within the negative 
$C_f$ distribution of the LSB. This is similar to the 
$C_f$ signature criterion used by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) to determine the ‘critical’ 
$LESP$ point from CFD. Note that unsteady Reynolds-averaged Navier–Stokes (uRANS) CFD was used in Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) to identify the 
$C_f$ signature criterion, whereas we use wall-resolved LES. The physical rationale for considering this point to signal imminent stall is as follows. The value of 
$C_f$ is initially negative within the LSB due to the recirculating flow. Stall onset occurs from instabilities to the bubble, leading to generation of counter-clockwise vorticity below the LSB, culminating in its bursting and coalescence into a coherent DSV. These instabilities manifest as spikes in the 
$C_f$ distribution, with values that are initially negative reaching zero, and then becoming positive, until the bubble bursts. Note that using the 
$C_f$ signature criterion by itself to characterise stall onset is not feasible with the current approach, since it involves capturing a small-scale instability that appears within the noisy 
$C_f$ distribution of the transitional/turbulent LSB region at every time step. Therefore, it is used only as a rough estimate for initiation of stall. It serves as a benchmark to evaluate the critical behaviour of a criterion based on an integrated parameter, without the same drawbacks.
To sum up, our definition of the stall onset regime encompasses the time when spikes reaching zero first occur in the LSB up until the time when the DSV is seen to have formed. This is shown schematically in figure 3. The width of this region depends on the type of stall (gradual or sudden), which in turn depends on aerofoil geometry, 
$Re_c$, etc. These bounds are approximate estimates to help verify if our proposed criterion reaches its critical value between them. This will ensure that 
$BEF$ is correlated with observed flow events. Appendix B shows how these bounds were identified for one of the cases.
Figure 3. Schematic showing the lower bound (onset of instabilities within the LSB) and upper bound (DSV formation) of the stall onset regime. LE stands for leading edge.
All the time series plots presented here have been low-pass-filtered using a Gaussian filter having half-width equal to the inverse of a non-dimensional cut-off frequency 
$f^+ = fc/U_{\infty }$ of 30, unless specified otherwise. The cut-off frequency was determined based on non-stationary spectral analysis carried out using 
$C_p$ time series from points near the leading edge. While filtering shifts the peak values of the parameters studied slightly aft in time, the results remain the same qualitatively.
3.2. Boundary enstrophy flux
We focus on vorticity and vorticity-based quantities in incompressible, near-wall viscous flows since the wall is both source and sink for all vorticity. Since viewing flow separation in terms of generation of surface vorticity in response to the prevailing pressure gradient (Lighthill Reference Lighthill1986) provides a useful physical picture, we first look at the integrated boundary vorticity flux, defined by
\begin{equation} BVF = \frac{1}{Re_c}\int_{(x/c)_p}^{(x/c)_s} \frac{\partial \omega}{\partial n} \, {\rm d}s, \end{equation}
where 
$\omega$ is the spanwise component of the vorticity vector 
$\boldsymbol {\omega }$, and 
$n$ and 
$s$ are normal and tangential directions to the aerofoil surface, respectively. Note that all the variables (
$\omega, n, s$) are non-dimensional. The integral is carried out from 
$x/c$ on the pressure side up to 
$x/c$ on the suction side, as shown in figure 4. Assuming small tangential surface acceleration, it is useful to think of the boundary vorticity flux as equivalent to the favourable streamwise pressure gradient (i.e. 
$(1/Re_c)\,\partial \omega /\partial n \sim -(1/\rho )\,\partial p/\partial s$). The behaviour of 
$BVF$ is not consistent across cases and is dependent on the region of integration selected. That is, the curves do not collapse when the region of integration is varied (see figure 5a).
Figure 4. Region of integration to calculate 
$BVF$, 
$BEF$ and 
$LESP$.
Figure 5. Variation with 
$\alpha$ of (a) 
$BVF$ and (b) 
$BEF$ integrated over different 
$x/c$, for SD500. 
$BEF$ curves all collapse irrespective of region of integration, in contrast to 
$BVF$.
Next, we look at the integrated boundary enstrophy flux (
$BEF$), defined by
\begin{equation} BEF = \frac{1}{Re_c}\int_{(x/c)_p}^{(x/c)_s} \omega\,\frac{\partial \omega}{\partial n}\, {\rm d}s . \end{equation}
Enstrophy is the square of the vorticity magnitude 
$|\boldsymbol {\omega }|^2$, the analogue of kinetic energy for velocity. For a span-averaged flow field as in the present case, 
$|\boldsymbol {\omega }|^2=\omega ^2$. We look at the integral of the quantity 
$\partial {(\omega ^2 / 2)}/{\partial n} = \omega \,(\partial \omega / \partial n)$; the factor 
$1/2$ is taken for convenience. It is divided by 
$Re_c$ so that it can be viewed as a product of vorticity and 
$BVF$ (compare (3.1) and (3.2)). The 
$BEF$ can be thought of roughly as an integral of the product of the vorticity and streamwise pressure gradient. Since it involves the pressure gradient rather than pressure, it is more useful and suitable for physical interpretation. Using 
$\omega = - \partial u / \partial n$ at the wall, 
$BEF$ can be rewritten as 
$(1/\rho )\int (\partial u/\partial n) \, (\partial p/\partial s) \, {\rm d}s$.
If the boundary layer is thought to be a vortex sheet, then the wall generates new vorticity depending on the edge velocity required to maintain the no-slip boundary condition at the wall (Lighthill Reference Lighthill1986). If the normal vorticity flux at the wall is of opposite sign compared to the vorticity at the wall, then viscous diffusion 
$\nu \,{\nabla }^2 \omega$ tries to smear out the gradient by transporting vorticity away from the wall, thereby weakening the wall vorticity. If, on the other hand, 
$\omega$ and 
$\partial \omega /\partial n$ are of the same sign, then wall vorticity is strengthened by viscous diffusion of vorticity of same sign towards the wall. A positive boundary enstrophy flux, which is the product of these two quantities, therefore represents a strengthening of the wall vorticity, while a negative sign represents its weakening (Wu, Ma & Zhou Reference Wu, Ma and Zhou2015), both quantified by its magnitude.
Figure 5(b) shows 
$BEF$ integrated over different chord lengths for SD500. In contrast to the behaviour of 
$BVF$ (figure 5a), the 
$BEF$ curves all collapse. This behaviour is explained in detail in § 3.4. For SD500, 
$BEF$ gradually drops, reaches a minimum at 
$\alpha = 19.0^{\circ }$ and begins to rise again sharply around 
$\alpha = 19.4^{\circ }$. Its magnitude, 
$|BEF|$, therefore reaches its peak value at 
$\alpha = 19.0^{\circ }$. The curves corresponding to integration over different 
$x/c$ diverge only after the formation of the DSV (
$\alpha > 20^{\circ }$). These trends are consistent across the cases investigated.
3.3. 
$BEF$ in comparison with $LESP$
The leading edge suction parameter 
$LESP$ was used by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) to determine when to trigger leading edge vortex shedding in their reduced-order model (based on unsteady thin aerofoil theory) for predicting dynamic stall. They argue that the leading edge limits the maximum suction to a critical value that can be supported by the aerofoil; this is achieved by vortex shedding at the leading edge. The parameter is a measure of the chord-wise suction force at the leading edge experienced by the aerofoil. The 
$LESP$ value is set equal to the 
$A_0$ term of the Fourier series describing the bound vorticity in thin aerofoil theory. As explained in Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014), the critical value can be predetermined from CFD using the 
$C_f$ signature criterion. Note that this corresponds to the lower bound of our definition of the stall onset regime. Therefore, the instantaneous 
$LESP$ value (set to 
$A_0$), reaching the critical value (predetermined using the 
$C_f$ signature criterion from CFD), determines the initiation of leading edge vortex shedding in the reduced-order model.
Based on Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016), we have calculated 
$LESP$ from CFD as a scaled coefficient of the force in the chord-wise direction near the leading edge, in a frame fixed to the aerofoil:
\begin{equation} LESP = \sqrt{\frac{|F_s|}{2 {\rm \pi}q_{\infty} c}}, \quad {\rm where} \ F_s = q_{\infty} \int_{(x/c)_p}^{(x/c)_s} C_{p} \hat{\boldsymbol{n}} \boldsymbol{\cdot} \hat{\boldsymbol{e}}_x\,{\rm d}s, \end{equation}
where 
$\hat {\boldsymbol {n}}$ is a unit vector normal to the aerofoil surface, 
$\hat {\boldsymbol {e}}_x$ is a unit vector in the chord-wise direction, and 
$q_{\infty }$ is the dynamic head; 
$F_s$ is the chord-wise suction force on the aerofoil. We integrate the pressure coefficient over the first 25 % chord (
$x/c=0.25$) to obtain 
$LESP$, since the main contribution to the chord-wise force originates from the leading edge. The definition of 
$LESP$ was revised in Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020) but the revision is not relevant for our study (see Appendix D).
A note on using 
$LESP$ as a criterion. The critical value of 
$LESP$ in all these studies is determined from a separate analysis based on observed flow events near the leading edge (occurrence of a spike in 
$C_f$ reaching zero in Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014) and Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016), and occurrence of an inflection point in the 
$C_f$ profile in Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020). The critical value is not determined from the behaviour of the parameter itself, e.g. 
$LESP$ reaching a maximum. However, our goal in this paper is to find a standalone criterion whose behaviour can signal stall onset without requiring a separate set of analyses. Since 
$LESP$ reaches its maximum value in the vicinity of stall onset, we find that 
$\max (LESP)$ is a good proxy for its critical value and therefore have used it in that manner, as have some others in the literature (Deparday & Mulleners Reference Deparday and Mulleners2018, Reference Deparday and Mulleners2019). Also, the maximum of the chord-wise force coefficient has been used previously as a dynamic stall onset criterion in semi-empirical models for data fitting (Sheng et al. Reference Sheng, Galbraith and Coton2005). Therefore, we have compared the two standalone criteria, namely, 
$\max (BEF)$ and 
$\max (LESP)$.
Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016) found that a positive recirculation region appears near the leading edge before the maximum value of 
$LESP$ is reached. While their results were obtained using unsteady RANS for Reynolds numbers of 30 000 and 
$3 \times 10^6$, we can corroborate the same findings at moderate Reynolds numbers of 200 000 and 500 000 from wall-resolved LES calculations with more accurate wall pressure distributions. We verify this using one realisation of case NC200, where the solution was recorded for the entire pitch-up motion. The time of occurrence of these events and the corresponding 
$LESP$ values are tabulated in the last two columns of table 2. Figure 6 shows the trend of 
$LESP$ (black curve) with 
$\alpha$, with the occurrence of the two relevant events marked. Clearly, 
$LESP$ reaches its maximum value following the occurrence of the 
$C_f$ signature criterion. The same figure also shows the variation of 
$|BEF|$, with its maximum occurring earlier, closer to the 
$C_f$ signature. For consistency with the calculation of 
$LESP$, we integrate over 25 % chord to calculate 
$BEF$ here. (Based on our analysis of the contributions to 
$BEF$ in § 3.4, we recommend an integration length of 1 % chord as appropriate and sufficient.)
Table 2. Flow events studied by Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016), shown along with 
$\alpha$ and 
$LESP$ found for NC200 from present work.
Figure 6. 
$LESP$ versus 
$|BEF|$ (unfiltered) for NC200 over the entire pitch-up motion. The 
$|BEF|$ peak is closer to the 
$C_f$ signature point than the 
$LESP$ peak.
The 
$\max (LESP)$ and 
$C_f$ signature points in figure 6 are separated by a change in 
$\alpha$ of more than 
$1^{\circ }$. Stall control efforts deployed without this delay after the onset of instabilities in the LSB could improve the outcome considerably. An advantage of using 
$LESP$ is that it requires only surface pressure data, which are obtained readily using pressure probes. However, there are two difficulties associated with the interpretation of the 
$LESP$ maximum, since it is calculated by integrating some part of the leading edge pressure distribution. One is in the use of pressure, as opposed to pressure gradient, which appears in the dynamical equations. The other is the arbitrariness in choosing the chord-wise or camber-wise component of 
$C_p \hat {\boldsymbol {n}}$.
We carry out similar comparisons for the SD7003 aerofoil cases next.
Figure 7(a) shows the time variation of 
$BEF$ for SD200, relative to the stall onset regime identified. The horizontal axis is the angle of attack, 
$\alpha$, which extends over the region around stall onset. The variation in 
$\max (|C_p|)$ near the leading edge up to 5 % chord is shown in the top panel, for an indication of when the leading edge suction collapses. The middle panel shows the magnitude of 
$BEF$, along with the 
$LESP$ variation plotted for reference; 
$LESP$ lags behind 
$BEF$ by almost a degree. The bottom panel shows the timeline of events occurring within the stall onset regime. The region in purple identifies the stall onset regime, extending between the time when spikes reaching zero first occur in the negative 
$C_f$ profile within the LSB, and the time at which the DSV can be identified as having formed. The magnitude of 
$BEF$ reaches a maximum first. Next, the peak magnitude of 
$C_p$ near the leading edge, which is a point quantity, reaches its maximum. Finally, 
$LESP$ reaches its maximum, approximately when 
$|C_p|$ has dropped by 2 % of its maximum value. Further, while 
$\max (|BEF|)$ occurs earlier within the stall onset regime as defined here, 
$\max (LESP)$ lies at its farther edge.
Figure 7. Sequence of events within the stall onset regime for (a) SD200 and (b) SD500. Top: 
$\max (|C_p|)$ near the leading edge. Middle: 
$\lvert BEF \rvert$ (green) shown along with 
$LESP$ (black). Bottom: timeline of all events. Extent of the stall onset regime (purple) varies inversely with 
$Re_c$.
To understand the effect of Reynolds number, we look at the same plot for SD500, shown in figure 7(b). In this case, both 
$\max (|C_p|)$ and 
$\max (|BEF|)$ occur at the same time. Following this, 
$LESP$ reaches its peak; at this higher Reynolds number, there is a smaller delay (about 
$0.3^{\circ }$) between the magnitude of 
$BEF$ and 
$LESP$ reaching their peaks. The stall onset regime also encompasses a shorter period of time. We find that the trend of 
$BEF$ reaching a peak before 
$LESP$ is maintained for all the cases analysed. The timeline for case NC500 is available in Appendix A.
Table 3 shows the angles corresponding to the events shown in the timelines for all cases. The static stall angles of attack, 
$\alpha _{ss}$, calculated from XFOIL are also included. The last column shows the delay between 
$|BEF|$ and 
$LESP$ reaching their respective peaks for all the cases, which ranges between 
$0.3^{\circ }$ and 
$0.8^{\circ }$ of rotation.
Table 3. Angles of attack (deg.) at which different events occur for the four cases investigated, along with that of static stall, 
$\alpha _{ss}$, found from XFOIL. The last column shows the delay (in deg.) between the maximum values of 
$|BEF|$ and 
$LESP$.
The 
$LESP$ and 
$BEF$ values capture the critical instability point in different ways. The 
$LESP$ value is a measure of the chord-wise force on the aerofoil. With the coalescence of the DSV, the force vector is rotated towards the chord-normal direction, leading to a drop in its chord-wise component. The 
$BEF$ value is a measure of net wall shear, with its critical (
$=\max (|BEF|)$) value corresponding to the instant of maximum net shear prevailing at the wall (see § 3.5). While both parameters are correlated through the dynamical equations, 
$BEF$ identifies the onset of instabilities earlier. We argue that a criterion based on the accumulation of net wall shear is more physically intuitive than the chord-wise component of suction force.
3.4. Physical significance of $BEF$
Next, we present an explanation for the behaviour of 
$BEF$ and the reason why some other vorticity-based parameters fail. For simplicity, we speak in terms of the pressure gradient, which grows in proportion to the vorticity flux. Since 
$BEF$ is found by integrating the product of vorticity and surface vorticity flux (or favourable streamwise pressure gradient), it has two significant sources: one due to the large clockwise (CW) vorticity at the leading edge of the aerofoil as flow accelerates around it, and the other from the counter-clockwise (CCW) vorticity induced by the LSB between the leading edge and the maximum thickness of the aerofoil. Figure 8(a) shows the space–time contours of 
$C_f$ for NC200 through the entire pitch-up motion. The 
$x$-axis goes from the pressure surface at mid-chord up to the suction surface at mid-chord, for which we have used the variable 
$\hat {x}$ (
$\hat {x}=x$ on the suction surface and 
$\hat {x}=-x$ on the pressure surface). The 
$y$-axis is the angle of attack in degrees. We identify three regions of analysis, which are delineated by the four green lines in figure 8(a). Region 1 extends from the maximum thickness location on the pressure side to the stagnation point. Region 2 extends from the stagnation point up to the transition point on the suction surface. Both of these regions are laminar. Region 3 extends from the transition point up to the point of maximum thickness on the suction surface. Since all of our cases have insignificant trailing edge separation, it is not necessary to consider regions beyond the maximum thickness point.
Figure 8. Region 1 (max. thickness on pressure side to stagnation point), region 2 (stagnation point to transition location on suction side) and region 3 (transition location to max. thickness on suction side) (a) identified from 
$C_f$ space–time contours, and (b) corresponding 
$BEF$ (unfiltered) contributions, for NC200 (full pitch-up).
We next look at the 
$BEF$ contribution from each of these regions, with increasing angle of attack. This is shown in figure 8(b). The contribution from region 1 is insignificant, since the pressure gradient is small within that region. The only significant contributions arise from regions 2 and 3, where there is large vorticity accompanied by large pressure gradients. Region 2 is fully laminar, and the 
$BEF$ trend mimics the rise and fall in the magnitude of the CW vorticity as well as the favourable pressure gradient (FPG) at the leading edge, with increasing angle of attack. The change in 
$BEF$ slope can be gradual or abrupt, since, as with the pressure gradient, it is dependent on the stalling characteristics of the aerofoil under the given operating conditions. The contribution from region 3 is due to CCW vorticity induced beneath the LSB over the suction surface. It grows from around the time of formation of the LSB. The recirculation region of the LSB, having CW vorticity, grows stronger with increasing 
$\alpha$, inducing a region of increasingly stronger adverse pressure gradient (APG). This APG is strongest to the left of the LSB centre (see discussion in Doligalski Reference Doligalski1994). The flow images in Appendix A show the relatively large CCW vorticity (red) region at the wall to the left of the LSB (blue) region. This leads to an increase in 
$BEF$ magnitude. As the DSV coalesces and moves away from the wall, the APG region becomes weaker and more smeared out, leading to a drop in 
$BEF$. This region, being transitional or turbulent, is noisy. However, this noisy contribution makes up a small part of the total 
$BEF$.
We observe that both the CW vorticity region (at the leading edge) and the CCW vorticity region (below the LSB) contribute to 
$BEF$ in the same sense, i.e. negative. Note that positive vorticity flux corresponds to the negative of the pressure gradient. So substituting pressure gradient for vorticity flux in the definition of 
$BEF$ includes a negative sign. Therefore, FPG is designated as ‘
$+$’ and APG as ‘
$-$’. The CW vorticity region (
$-$) at the leading edge is accompanied by FPG (
$+$), leading to their product being negative. The CCW vorticity region (
$+$) at the LSB region is accompanied by APG (
$-$), also leading to a negative product. This is shown schematically in figure 9. If we simply integrated vorticity or pressure gradient over each of these regions, the contributions from different regions would partially cancel, as they have opposite signs. To circumvent this problem, if the square of either quantity is integrated, i.e. 
$\omega ^2$ or 
$(\partial p / \partial s)^2$, the contributions from the noisy LSB region are not necessarily insignificant. The contribution from the LSB region remains small only when the interaction of pressure gradient and wall vorticity is considered through their product in the calculation of 
$BEF$. Note that 
$BEF$ is calculated by integrating the product of vorticity and surface vorticity flux and not from the product of their integrals.
Figure 9. Schematic showing contributions to 
$BEF$ from the laminar leading edge and transitional/turbulent LSB regions; FPG is designated as 
$+$, since substituting pressure gradient for vorticity flux in the definition of 
$BEF$ includes a negative sign, i.e. 
$BEF\sim (1/\rho ) \int \omega \,(- \partial p / \partial s)\,{\rm {ds}}$.
The above explanation helps us to address how the 
$|BEF|$ maximum can be used. All the noise in the 
$BEF$ time signal arises from the LSB region. A low-pass filter applied to the signal yields the laminar portion, which forms the bulk of the total 
$BEF$, along with a smoothed contribution from the LSB region. The laminar 
$BEF$ contribution by itself signals the onset of stall earlier than considering 
$C_p$ or pressure-based parameters, since 
$BEF$ is based on pressure gradient. The addition of any smoothed contribution from other regions shifts the location of the peak only slightly. A question might arise here as to why the integral of the pressure gradient from the laminar region alone is not sufficient to signal stall, since it is the region that contributes the most to 
$BEF$. However, this would require knowledge of the dynamically varying spatial boundaries of the LSB region, since its inclusion would lead to a drop in the magnitude of the favourable pressure gradient. Also, say, for a plunging motion, what constitutes the ‘suction’ surface is time-varying, and it is difficult to separate the laminar region contribution from the transitional region. Whereas, as long as the region very close to the leading edge on both sides of the geometry is captured, the laminar contribution to 
$BEF$ can be identified easily. As this region is very small, integration up to 1 % chord is sufficient.
Figure 10 shows the filtered contributions from the laminar region 2 (referred to as 
$BEF_{LE}$), and transitional/turbulent region 3 (referred to as 
$BEF_{LSB}$) for the other three cases. In each case, the contribution from the LSB region is less than a fifth of that from the leading edge region. The locations of the maximum values of 
$|BEF_{LE}|$ and 
$|BEF_{LE} + BEF_{LSB}|$ are marked with an ‘
$\times$’ in the figure and are observed to be very close.
Figure 10. 
$BEF$ contribution from the leading edge (
$BEF_{LE}$), from the LSB region (
$BEF_{LSB}$), and from both (
$BEF_{LE} + BEF_{LSB}$): (a) SD200, (b) SD500, (c) NC500.
3.5. Analytical relation between $BEF$ and total pressure
Next, we relate 
$BEF$ to other flow parameters analytically using the dynamical equations. We can relate 
$BEF$ to total pressure by writing the non-dimensionalised momentum equation for an incompressible flow using the Lamb vector, 
$\boldsymbol {\omega } \boldsymbol {\times } \boldsymbol {u}$:
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\omega} \boldsymbol{\times} \boldsymbol{u} ={-} \boldsymbol{\nabla} \left(\frac{p}{\rho} + \frac{|\boldsymbol{u}|^2}{2}\right) + \frac{1}{Re_c}\nabla^2 \boldsymbol{u}. \end{equation}
By taking the divergence of (3.4) and using 
$\boldsymbol {\nabla } \boldsymbol {\cdot } (\boldsymbol {\omega } \boldsymbol {\times } \boldsymbol {u}) = \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \times \boldsymbol {\omega } - \boldsymbol {\omega } \boldsymbol {\cdot } \boldsymbol {\omega }$, we get a relation between enstrophy and total pressure. For a 2-D flow, this enstrophy–total pressure relation is given by
\begin{equation} \omega^2 = \nabla^2 \left( \frac{p}{\rho} + \frac{|\boldsymbol{u}|^2}{2}\right) - \boldsymbol{u} \boldsymbol{\cdot} \nabla^2 \boldsymbol{u}. \end{equation}
If we apply the viscous wall boundary condition at this stage, then the last term drops out. Therefore, the enstrophy at the wall is given by the Laplacian of the total pressure. However, if we take the normal derivative of (3.5) and then apply the viscous wall boundary condition and neglect curvature, we obtain (because the normal derivative of the last term from (3.5) evaluated at the wall reduces to 
$- \omega \,\partial \omega /\partial n$)
\begin{equation} 3 \omega\,\frac{\partial \omega}{\partial n} = \frac{\partial}{\partial n} \nabla^2 \left( \frac{p}{\rho} + \frac{|\boldsymbol{u}|^2}{2}\right) . \end{equation} Using the above relation in our definition of 
$BEF$ gives
\begin{equation} BEF = \frac{1}{Re_c}\int_{(x/c)_p}^{(x/c)_s} \omega\,\frac{\partial \omega}{\partial n}\, {\rm d}s = \frac{1}{3 Re_c} \int_{(x/c)_p}^{(x/c)_s} \frac{\partial}{\partial n} \left[\nabla^2 \left( \frac{p}{\rho} + q\right)\right] {\rm d}s, \end{equation}
where 
$q = |\boldsymbol {u}|^2/2$ is the kinetic energy. We have related the third derivative of total pressure to 
$BEF$. We observe that the contribution from the integral of 
$\partial ^3 q / \partial n^3$ is the most dominant, and it exhibits the same trend as 
$BEF$. Further, we find that the integral of the first derivative term 
$\partial q / \partial n$ also exhibits a similar trend to 
$|BEF|$. It reaches a maximum at the same time as 
$|BEF|$. This is shown for SD200, SD500 and NC500 in figure 11. Therefore, a good approximation to the behaviour of 
$BEF$ is given by the integral of 
$\partial q / \partial n$ or the normal gradient of kinetic energy at the wall. Essentially, this criterion distils down to the net squared wall shear (
$(\partial u / \partial n)^2$) over the aerofoil reaching a maximum, as shown by
\begin{equation} \left(\frac{\delta q}{\delta n}\right)_{wall} = \frac{u^2_{flow}}{2\,\delta n} = \frac{1}{2}\left|\frac{\delta u}{\delta n}\right|^2_{wall} \delta n, \end{equation}
since 
$u_{flow} = ({\delta u}/{\delta n})_{wall}\,\delta n$; the subscript ‘flow’ refers to the quantity at the first grid point from the wall in the normal direction.
Figure 11. Integral of the normal gradient of kinetic energy, 
$\partial q/ \partial n$ over the aerofoil surface, for (a) SD200, (b) SD500 and (c) NC500. This quantity shows the same trend as 
$|BEF|$.
3.6. Potential use of $BEF$ for stall control
Application of the proposed 
$BEF$ criterion in practice requires two elements: (1) flow information near the aerofoil leading edge to compute 
$BEF$, and (2) a procedure to identify when 
$\max (|BEF|)$ occurs. In this section, we propose ideas on how these could be achieved.
While measurements of pressure are possible via surface-mounted sensors, velocity information near the aerofoil leading edge has to be inferred. We argue that recent advances in the field of scientific machine learning (SciML) provide a path towards estimating the near-surface velocity field in real time. In particular, recent work based on the concept of physics-informed neural networks, for instance, Raissi, Perdikaris & Karniadakis (Reference Raissi, Perdikaris and Karniadakis2019a), has shown that the flow field information can be extracted from sparse measurements (Raissi et al. Reference Raissi, Wang, Triantafyllou and Karniadakis2019b), and can be deployed in practice (Cai et al. Reference Cai, Li, Zheng, Kong, Dao, Karniadakis and Suresh2021). We present an application of this approach to a simple problem in Appendix C.
Next, we present a way to utilise the filtered 
$|BEF|$ obtained in real time to determine the point of initiation of control efforts. Since the bulk of the 
$BEF$ contribution originates from the laminar leading edge, and we want to reduce the noisy contribution from the LSB region, we recommend using the first 1 % chord (on both sides) for 
$BEF$ calculations. 
$|BEF|_{max}(t)$ at each data point can be set equal to the maximum value of 
$|BEF|$ reached up to that time instant. The percentage change of the instantaneous 
$|BEF|$, denoted 
$|BEF|(t)$, from the maximum value 
$|BEF|_{max}(t)$, can be tracked. We call this percentage change 
$\sigma _{|BEF|}$, as given in the following equation:
\begin{equation} \sigma_{|BEF|} = 100 \times \left(\frac{|BEF|(t)}{|BEF|_{max}(t)} - 1.0\right). \end{equation}
When 
$\sigma _{|BEF|}$ falls by, say, 5 %, it is assumed that stall onset is imminent and control efforts can be deployed at the corresponding angle of attack, 
$\alpha _{control}$. The variation of 
$\sigma _{|BEF|}$ with 
$\alpha$ is shown for a couple of cases in figure 12. The variation of 
$\sigma _{LESP}$ is shown for comparison. 
$BEF$ is integrated up to 1 % chord, since its dominant contribution originates from the laminar leading edge region, as explained in § 3.4. We obtain 
$LESP$ by integrating up to 25 % chord, as defined in Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016). (The value 
$\max (LESP)$, calculated by integrating up to 1 % chord, has a smaller delay from 
$\max (|BEF|)$, close to the values in table 3.) The delay between 
$\alpha _{max(|BEF|)}$, the angle where the actual peak of 
$|BEF|$ occurs, and 
$\alpha _{control,BEF}$ will depend on the threshold chosen. A similar delay applies for any other parameter (e.g. 
$LESP$) that is tracked in this way.
Figure 12. The graphs show 
$\sigma _{BEF}$ with 
$\alpha _{control}$ when 
$\sigma ^{threshold}$ is set to 
$-5\,\%$; 
$\sigma _{LESP}$ is also shown for reference. Panel (a) NC200 (full pitch-up), (b) SD500.
To summarise, we have demonstrated that the vorticity-based 
$|BEF|$ parameter reaches its maximum within the identified stall onset regime. A major advantage of using 
$BEF$ to signal stall onset is that it is nearly independent of the integration length selected. By decomposing the contributions to 
$BEF$ from the laminar and transitional/turbulent regions in the flow, we have shown that the region very close to the leading edge is sufficient to capture the 
$BEF$ peak for identifying stall onset. We have also related the 
$|BEF|$ maximum to the instant of maximum net wall shear analytically. While the evaluation of 
$BEF$ requires higher-order derivatives, we argue that recent advances in SciML offer potential for its use in practical applications.
4. Conclusion
The current work was focused on characterising the onset of leading-edge-type dynamic stall in pitching aerofoils. In some wind turbine and helicopter applications, we would like to avoid DSV formation completely, if possible (Doligalski Reference Doligalski1994). An important factor in implementing flow control to prevent DSV formation is that there is a certain critical period before DSV formation when stall control is most effective. Therefore, the motivation to characterise dynamic stall derives from the need to deploy control efforts in a timely manner, before DSV formation.
In attempting to characterise dynamic stall onset, any proposed criterion needs to correlate with observed flow phenomena that are agreed upon as signalling DSV initiation. Keeping in mind flow control effectiveness, we defined a stall onset regime, setting its lower bound to be the first time when a spike reaching zero appears in the negative 
$C_f$ distribution of the LSB (similar to the 
$C_f$ signature criterion, Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). We set the upper bound to the time when the DSV is observed as having formed, since this marks the end of the critical period for flow control efforts. Our analysis is based on span-averaged results from wall-resolved large-eddy simulations of aerofoils undergoing a pitch-up motion at a constant rate, pivoted about the quarter-chord point at Mach number 0.1, and 
$Re_c$ between 200 000 and 500 000. Our dataset comprises two different aerofoils (SD7003 and NACA0012) at two different values of 
$Re_c$ (200 000 and 500 000).
We proposed the vorticity-based boundary enstrophy flux (
$BEF$) parameter to characterise leading edge dynamic stall onset. The parameter 
$|BEF|$ reaches its maximum within the identified bounds of the stall onset regime for all our cases and is nearly independent of the integration region used for its evaluation, which makes its use appealing. It is an aggregate parameter that encompasses inviscid effects as well as small-scale boundary layer effects. We found that the dominant contribution to 
$BEF$ is from the laminar leading edge region, which has large clockwise vorticity and favourable pressure gradient due to the pitching motion. A relatively small contribution originates from the transitional/turbulent LSB region, due to LSB-induced counter-clockwise vorticity and adverse pressure gradient. We also found – by relating analytically 
$BEF$ to the total pressure and applying simplifying assumptions – that 
$\max (|BEF|)$ corresponds to the instant of maximum net shear prevailing at the wall.
Finally, we compared 
$\max (|BEF|)$ with the most suitable standalone criterion that can be derived from the existing literature, namely, the maximum of the leading edge suction parameter (Ramesh et al. Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014). We found that 
$BEF$ and 
$LESP$ capture the instability point in different ways. The parameter 
$BEF$, reflecting the accumulation of wall shear, captures the onset of instabilities in the LSB, and reaches its maximum earlier. The pressure-based 
$LESP$ tracks the rotation of the force vector acting on the aerofoil, moving away from the chord-wise direction as the DSV coalesces. The parameter 
$|BEF|$ reaches its peak between 
$0.3^{\circ }$ and 
$0.8^{\circ }$ in rotation earlier than 
$LESP$. We therefore identified 
$BEF$ as a more conservative dynamic stall onset criterion.
Acknowledgements
We acknowledge the assistance by Drs M. Visbal and D. Garmann from the Air Force Research Laboratory in providing initial CFD solutions. We also acknowledge fruitful discussions with Professor P. Durbin on the topic.
Funding
Funding for this research is provided by the National Science Foundation (grants CBET-1554196 and 1935255). Co-author Sharma acknowledges the support provided to him by the 2019 AFOSR Summer Faculty Fellowship. Computational resources are provided by NSF XSEDE (Grant no. TG-CTS130004) and the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Additional results
A.1. NC500 timeline
The relatively gradual stall of NC500, extending over a wide region in time (purple) is shown in the timeline of figure 13. Here, 
$\max (|C_p|)$ occurs ahead of 
$\max (|BEF|)$; however, the trend of 
$BEF$ reaching a peak before 
$LESP$ is still maintained. The flow images that follow also provide an idea of the gradual development of the DSV.
Figure 13. Sequence of events within the stall onset regime for NC500. The trend of 
$\max (|BEF|)$ preceding 
$\max (LESP)$ is maintained. The large extent of the stall onset regime (purple) indicates gradual stall.
A.2. Flow images
We present vorticity contours overlaid with streamlines for a sequence of angles of attack for NC500 (figure 14) and SD500 (figure 15). The progression of the LSB as it coalesces into the DSV can be observed. Also, the contrast between the gradual stall in NC500 and abrupt stall in SD500 is evident from 
$\Delta \alpha$ between 
$C_f$ instabilities and DSV formation. There is an abrupt change in the scale of the LSB for the SD500 case. Also, due to the smaller leading edge radius of the SD7003 aerofoil, laminar separation and reattachment occurs closer to the leading edge compared to the NACA0012 aerofoil.
Figure 14. Vorticity contours overlaid with streamlines for NC500.
Figure 15. Vorticity contours overlaid with streamlines for SD500.
Appendix B. Stall onset regime identification
The lower bound of the stall onset regime is at the first time of occurrence of a spike reaching zero within the 
$C_f$ profile of the LSB region. Since these spikes in the 
$C_f$ profile are transient, a careful inspection of the profiles is necessary at each time step. For example, we show in figure 16 the intermittent appearance of these positive spikes for SD200. A positive 
$C_f$ spike (shown by red arrows) in the aft portion of the LSB is visible at angles 
$14.42^{\circ }$ and 
$14.86^{\circ }$, but not at an intermediate angle 
$14.63^{\circ }$. We consider the earliest time step at which such a spike is observed as the lower bound of the stall onset regime.
Figure 16. The graphs show 
$C_f$ on the suction surface for SD200 around onset of instabilities. Red arrows indicate the appearance of spikes reaching zero within the LSB.
The upper bound of the stall onset regime is considered to be the time when the DSV is observed to be formed from visual examination of streamlines and 
$C_p$ distributions. For example, vorticity contours overlaid with streamlines for SD200 are shown in figure 17 for different angles of attack. We observe that panels (a–c) in the figure are qualitatively different from panel (e), with the streamwise extent of the recirculation region significantly increased. Therefore, we estimate that the DSV has formed sometime between panels (d) and (e). In panel (f), the vortex has grown further.
Figure 17. Vorticity contours overlaid with streamlines for SD200 around DSV formation.
Appendix C. Flow field reconstruction from sparse measurements
Computing 
$BEF$ requires knowledge of the velocity gradients close to the surface. While a complete workflow is beyond the scope of this work, we present a simplified example to illustrate how SciML can be used to obtain the velocity field from sparse surface pressure measurements.
Building from recent work (Rao, Sun & Liu Reference Rao, Sun and Liu2020), we train a neural network to reconstruct the full 2-D flow field given sparse surface pressure measurements. We consider steady flow past a 2-D cylinder at Reynolds number 
$100$. The numerically solved flow fields (using Fluent) are shown in figure 18(b,d). From these data, we take sparse pressure measurements (at 50 uniformly distributed points around the cylinder) and train a physics informed neural network (PINN) to reconstruct the full 2-D velocity field, which is shown in figure 18(a,c). The integrated r.m.s. error between the reconstructed flow field and the true flow field was less than 
$10^{-2}$, indicating the promise of such approaches.
Figure 18. Panels (a,c) show reconstructed 2-D velocity fields using sparse pressure measurements on the surface of the cylinder; (b,d) show ‘true’ velocity fields for flow past a cylinder at Reynolds number 100.
Training the neural network for this problem took less than an hour on a laptop. Figure 19 plots the loss function that the neural network optimises, which is a sum of the measurement error (predicted pressure against measured pressure) and the residual of the Navier–Stokes equation (how well the predicted velocity field satisfies the Navier–Stokes equation). No effort was made to optimise, fine-tune and accelerate the code. This is especially promising given that GPU acceleration can reduce training times significantly, paving the way to (near) real-time estimation of near-surface flow fields.
Figure 19. Training loss as a function of epochs for the PINN. We switch from a first-order optimiser to a second-order optimiser after 100 epochs.
Appendix D. $LESP$ definition
We have used the 
$LESP$ definition from Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2016) in our calculations, as a scaled chord-wise suction force. Although the 
$LESP$ definition was updated to be taken in the camber-wise direction at the leading edge in Narsipur et al. (Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020), this amounts to an insignificant change for our cases. The 
$LESP$ values for the cambered SD7003 aerofoil are scaled by a constant of about 0.99. An additional update in the 
$LESP$ definition was the change in dynamic pressure scaling to be based on the relative velocity of the aerofoil at half-chord. This also scales our 
$LESP$ values by a constant nearly equal to 1, with insignificant changes to the location of the 
$LESP$ peak in time, which is our primary interest. Therefore, these changes have not been applied to our presented results.
Open access
