3.1. General flow features

Figure 1 shows the lift ($C_l$), drag ($C_d$) and pitching moment ($C_m$) coefficients as functions of the effective angle of attack for freestream Mach numbers $M_\infty = 0.1$ and $0.4$. The dashed blue lines indicate the phase-averaged coefficients obtained by Visbal (Reference Visbal2011) after six cycles, for $M_\infty = 0.1$. From this figure, good agreement is observed between the present results and those from Visbal, despite the differences in the numerical methods, computational grids and phase-averaging effects. At this point, it is important to reiterate that our results display the aerodynamic loadings for a single cycle, and that they lie within the cycle-to-cycle variations shown in Visbal's work. The dashed red lines, in turn, stand for the results obtained by one cycle of the $M_\infty = 0.4$ flow with the finer grid. Differences observed are due to the sensitivity of the flow to the initial conditions, as would occur due to cycle-to-cycle variations. The phase angle $\phi$ is also shown in the plots for particular instants of the motion, and the circle and cross symbols on the $C_l$ and $C_d$ plots mark the dynamic stall onset time based on different criteria (further details are provided in § 3.3). It is observed that compressibility acts on the sense to attenuate the aerodynamic loads in the hysteresis loop while maintaining their maximum and minimum values, a trend that was also observed by Sangwan et al. (Reference Sangwan, Sengupta and Suchandra2017) for a two-dimensional simulation of a pitching aerofoil. Here, we observe that the largest changes in $C_l$ and $C_d$ occur when the effective angle of attack is reduced from $\alpha _{eff} \approx 18^{\circ }$ to $10^{\circ }$, which corresponds to the time interval when the aerofoil is ceasing its downstroke motion. At $\alpha _{eff} = 8^{\circ }$, the aerofoil is at the bottom-most position for the downstroke.

Figure 1. (a) Lift, (b) drag and (c) pitching moment coefficients, as functions of the effective angle of attack for $M_\infty = 0.1$ and $0.4$. Dashed blue lines are phase-averaged results from Visbal (Reference Visbal2011), and dashed red lines are results for $M_\infty = 0.4$ with a finer grid. Circle and cross symbols represent onset time based on LESP and ${\rm \Delta} z$ criteria, respectively (see § 3.3).

These results, especially the accentuated decay of the coefficients during the post-stall period with the increasing Mach number, can be understood better from figures 2 and 3. The former exhibits the spanwise-averaged pressure coefficient ($C_p$) contours at different flow instants for the two Mach numbers investigated (figures 2a,b) and its distribution over the aerofoil suction side (figure 2c). The formation of the DSV is also shown in grey shading using a flow entropy measure defined as $(p/p_0)/(\rho /\rho _0)^{\gamma } - 1$, where $p_0 = 1/\gamma$ and $\rho _0 = 1$ are reference values for pressure and density, respectively. This entropy measure gives a better visual representation of the separated flow since it indicates the regions where entropy is changing due to viscous effects. Figure 3 presents the full history of pressure ($C_p$) and skin friction ($C_f$) coefficients along the suction side of the aerofoil.

Figure 2. Spanwise-averaged pressure coefficients at different flow instants for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$; and (c) $C_p$ over the aerofoil suction side. The evolution of the DSV can be seen in a movie (supplementary movies are available at https://doi.org/10.1017/jfm.2022.165).

Figure 3. Comparison of (a,b) pressure and (c,d) skin friction coefficients for (a,c) $M_{\infty } = 0.1$ and (b,d) $M_{\infty } = 0.4$. Auxiliary lines and markers are included to facilitate the comparison between the two flows and to represent specific events as described in the text.

From figure 2, two observations can be drawn immediately. The first is that the suction peak is reduced near the leading edge for the higher Mach number flow, as demonstrated in figure 2(c) for $t=4.5$. This condition is accompanied by a weaker, more diffuse, pressure core of the DSV. The second observation concerns the gestation period and residence time for which the DSV remains over the aerofoil. Heuristic reasoning to explain the weakening of the DSV strength with compressibility is provided by Chandrasekhara & Carr (Reference Chandrasekhara and Carr1990). They argue that for higher Mach numbers, the pressure gradient is reduced near the leading edge due to the smaller curvature of the streamlines from the earlier flow separation. At this condition, the net vorticity introduced is lower, which leads to a weaker vortex formation. With respect to the second observation, figure 2 shows that the DSV residence time is smaller for $M_{\infty } = 0.4$, and that it is more spread at $t=6.5$, spanning a larger region over the aerofoil suction side. On the other hand, the DSV is relatively concentrated at the mid-chord location for the lower Mach number case. At $t=7.5$, the DSV is already leaving the aerofoil surface for the higher Mach number flow, followed by a trailing edge vortex.

The same conclusions can be drawn from the map of $C_p$ over the aerofoil suction side, depicted in figures 3(a,b). In this $x/c$–$t$ diagram, the dark blue contours represent the signature of the DSV over the aerofoil suction side followed by the formation of the trailing edge vortex. We observe that the $C_p$ contours are lighter for $M_{\infty } = 0.4$ (see marker  4 in the figures), indicating the presence of a weaker pressure core. Moreover, with closer inspection, it is possible to measure a delay of $17^{\circ }$ in the phase angle between the formation of the trailing edge vortex for the two Mach number flows, as indicated by the horizontal black lines on marker 5. The solid and dashed horizontal lines show the initial formation of the trailing edge vortex for the $M_{\infty } = 0.1$ and $0.4$ cases, respectively. While the trailing edge vortex starts to form at $\phi = 213^{\circ }$ for $M_{\infty } = 0.4$, its appearance occurs at $\phi = 230^{\circ }$ for $M_{\infty } = 0.1$.

The $C_f$ plots in figures 3(c,d) present further information about the flow history over the suction side of the aerofoil. At about $90^{\circ }$ phase angle, an oscillatory pattern appears close to the aerofoil trailing edge as shown by marker 1. These oscillations arise due to initial shedding of von Kármán vortices from the aerofoil trailing edge, which exhibit a higher frequency at the lower Mach number flow. This vortex shedding can be visualized better in a movie (see supplementary movies). In the following moments, high-frequency $C_f$ variations can be observed spanning almost the entire chord (denoted by marker 2) in the range $120^{\circ } <\phi < 130^{\circ }$ for both flow conditions, but they are more pronounced for the $M_{\infty } = 0.1$ case. These fluctuations are due to Kelvin–Helmholtz instabilities forming during the downstroke motion. Figure 4 shows this primary instability stage of the flow at different instants using entropy contours. As is evident, initially the instabilities arise closer to the trailing edge, but they develop rapidly along the entire suction side. They appear to form at an earlier time for the lower Mach number and exhibit a more organized and denser behaviour where flow structures display a higher wavenumber.

Figure 4. Entropy contours showing the development of Kelvin–Helmholtz instabilities (here called the primary instability stage) along the aerofoil suction side for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$.

For both cases investigated, the formation of the dynamic stall vortex appears just after the Kelvin–Helmholtz instabilities approach the leading edge, at phase angle $\phi \approx 135^{\circ }$, shown by marker 3 in figure 3. In the same figure, the travelling signature of the DSV is shown by the negative $C_f$ values shown in blue contours along the aerofoil chord for intermediate phase angles $\phi$. Auxiliary black lines are included to facilitate the comparison of the mean velocity at which the DSV is being transported for each compressible regime. The lines represent an approximate path of the DSV when moving, considering the centre of the blue region. The solid line corresponds to the DSV signature of the $M_{\infty } = 0.1$ flow, while the dashed line refers to the $M_{\infty } = 0.4$ case. They are repeated in both subplots to facilitate the comparison of their slopes, and show that the DSV is being advected faster in the higher Mach number flow. The weaker signature of friction coefficient for $M_{\infty } = 0.4$ is also noticeable compared to that for $M_{\infty } = 0.1$.

The skin friction maps in figure 3 also display the presence of a flow reversal at earlier times that progressively advances over the aerofoil suction side from the trailing edge towards the leading edge as $\phi$ increases. This region appears for $\phi \leq 120^{\circ }$ in light blue contours, below markers 1, 2 and 3. The time (here also visualized as a function of the plunge cycle angle $\phi$) in which the flow reversal initiates does not appear to change under the distinct compressible regimes investigated. The main differences in $C_f$ between the two flows start after the beginning of the primary instabilities (marker 2), which are more sparse for the higher Mach number flow. As mentioned before, the transport of the DSV over the suction side is faster for $M_{\infty } = 0.4$, and its signature is more diffuse. Moreover, the point where the development of the DSV occurs changes from approximately $x/c = 0.1$ for $M_{\infty } = 0.1$ to $x/c = 0.2$ for $M_{\infty } = 0.4$ (see different positions of marker 3). This behaviour can be understood better through an inspection of the boundary layer forming over the aerofoil suction side.

3.2. Onset of dynamic stall

Figure 5 shows the spanwise-averaged $z$-vorticity contours in a short time window that marks the moment of formation of the DSV. As can be seen from the figure, ${\rm \Delta} t = 0.3$ based on the freestream velocity, which corresponds to nearly 2.4 % of the plunging cycle period. Unlike for steady separation, in dynamic stall the outer flow continues to follow the aerofoil contour despite the presence of flow reversal. As a consequence, a local shear layer forms between the displaced leading edge boundary layer and the reversed fluid layer that is, at a later stage, subjected to inflectional instabilities. This leads to the generation and growth of coherent structures that can be visualized in figure 5 and in a movie (see supplementary movies). In the figure, auxiliary lines are plotted to indicate regions where $u$- and $v$-velocity components change direction, from which we identify the presence of flow reversal and the shear layer over the leading edge.

Figure 5. Spanwise-averaged $z$-vorticity field near the leading edge ($x/c \le 0.25$) revealing the sudden boundary layer separation for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$. Dashed and solid lines represent the boundaries for which $u=0$ and $v=0$, respectively, indicating flow regions where velocity components change direction.

In previous work by van Dommelen & Shen (Reference van Dommelen and Shen1980), it is speculated that the onset of the unsteady separation phenomenon is an inviscid process, independent of the Reynolds number. This characteristic is attributed to the fact that the initial flow reversal near the surface, which is triggered by a strong adverse pressure gradient, is later governed by inertial effects on the zero-vorticity line. Along this line, formed as the fluid approaches the separation region, the convective terms dominate the viscous effects, justifying the inviscid assumption. Under the influence of the increasing adverse pressure gradient, the local reversed flow begins to accelerate rapidly upstream near the leading edge region, generating an intense shear. Because of the presence of the solid surface, fluid is propelled away from the wall. Meanwhile, the vorticity conservation yields an outward distortion of the zero-vorticity line, destabilizing the local vorticity distribution and resulting in the formation of a large vortical structure. The development of this vortex and its induced secondary separation leads to its ejection from the surface.

The process described above is visualized in figure 6 and in a movie (see supplementary movies), where the local Mach number is plotted along with the zero-vorticity lines during the dynamic stall onset. In this figure, a straight line measuring $0.06 c$ in length is plotted normal to the aerofoil surface, from which tangential and normal velocity components are computed to describe the evolution of the local velocity field. The placement of these lines is chosen in order to highlight similar behaviour of the velocity field at $t=4.4$, but since separation occurs at different chord positions depending on the freestream Mach number, they are placed at different distances from the aerofoil leading edge for each flow. The values $x/c = 0.11$ and $0.15$ are chosen for $M_{\infty }=0.1$ and $M_{\infty }=0.4$, respectively. We stress that these values are meant not to indicate the precise location of the separation, but only to help to illustrate the van Dommelen and Shen model. With that clarification, figure 7 shows the profiles of tangential ($u_t$) and normal ($u_n$) velocity components in blue and red, respectively. The corresponding vectors are also displayed in the figure. The tail of each vector is placed at the respective normal distance where its velocity components are shown.

Figure 6. Local Mach number contours near the leading edge ($x/c \le 0.25$) for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$. White dashed lines represent the zero-vorticity lines, and solid red lines represent regions of sonic flow.

Figure 7. Tangential ($u_t$, blue) and normal ($u_n$, red) velocity profiles for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$.

In figure 6, the flow continues to accelerate up to a point further downstream of the leading edge where it encounters the reversed flow accelerating along the inner layer in the opposite direction. This acceleration is also evident in the increment of the tangential velocity between the first and second columns in figure 7. In this region, a mutual interaction is apparent between the outer flow and the flow inside the zero-vorticity line to form a low-pressure core (not shown). At $M_{\infty }=0.4$, the strong acceleration of the reversed flow is sufficient to achieve sonic speeds along the zero-vorticity lines (red lines in figure 6 represent regions of $M = 1$). The flow, however, is decelerated almost isentropically, and no shock waves are observed despite the sonic flow regions. The strong shear accompanied by the low pressure on the aft portion of the separation region results in the flow reversal being ejected away from the surface, as shown by the velocity vectors in the third and fourth columns of figure 7. The interaction of the ejected flow with the outer flow creates a vortical structure that intensifies the instabilities in the region. The dynamics observed in figures 5 and 6 appear to agree well with the overall concept of the van Dommelen and Shen model, which provides the description of the flow up to the early stages of separation.

As observed from the previous figures, the onset of the DSV appears further downstream from the leading edge with an increase in Mach number, as was also reported by Chandrasekhara & Carr (Reference Chandrasekhara and Carr1990) and Benton & Visbal (Reference Benton and Visbal2020) for pitching aerofoils. In their work, compressibility caused the leading edge separation to occur at an earlier time and, consequently, at a lower effective angle of attack. Therefore, the separation for higher Mach number flows would experience a less pronounced curvature of the streamlines and a weaker pressure gradient. However, this earlier separation is not observed in our simulations at $Re = 6 \times 10^{4}$, suggesting that the pressure gradient reduction may not be a direct consequence of the lower effective angle of attack. It is important to mention that the observations from the previous authors are for higher Reynolds numbers than that considered in the present work.

For the present periodic plunging motion, the leading edge separation seems to be related to the primary instabilities developing at the trailing edge region and, later, spanning the entire aerofoil suction side as shown in figure 4. These instabilities approach the leading edge and interact with the secondary instability presented in figures 5 and 6. Notice that the term ‘secondary instability’ does not refer to the jargon commonly used in classical stability analysis. This term is used here only to differentiate the leading edge instability from the Kelvin–Helmholtz one. In particular, we notice from figure 4 that the primary instabilities reach further upstream for the $M_{\infty }=0.1$ case. On the other hand, for $M_{\infty }=0.4$, a wider favourable pressure gradient region along the leading edge causes the formation of secondary instabilities to be shifted downstream over the aerofoil suction side and slightly delayed in time (see the first two columns of figure 5). This is illustrated in figure 8, which shows the instantaneous covariant derivative of pressure along the aerofoil. The variable $\xi = \xi (x,y)$ is the curvilinear coordinate of the O-mesh topology, which goes along the aerofoil profile, and $\mathrm {d}\xi / | \mathrm {d}\xi |$ is introduced to remove the influence of the mesh orientation of the covariant derivative. The basis of the curvilinear system is also normalized to yield a scale in physical units. As discussed in the literature (Li et al. Reference Li, Komperda, Ghiasi, Peyvan and Mashayek2019), compressibility has a stabilizing effect that delays the laminar to turbulence transition in free shear flows. Here, the mechanism appears to be similar, although the plunging motion together with the surface presence lead to a more complex flow.

Figure 8. Pressure gradient contours (blue and red) in the streamwise $\xi$ direction and contours of entropy (transparent shading) for (a) $M_{\infty } = 0.1$,and (b) $M_{\infty } = 0.4$, at $t=4.25$.

Given the importance of the pressure gradient on the dynamic stall onset, we display the pressure gradient history along the $\xi$ direction over the aerofoil suction side in figure 9. The saturation level is kept high to call attention to the features marked in the figure. As can be confirmed from the movie (see supplementary movies), von Kármán vortices are shed from the aerofoil trailing edge before the primary instability stage. These vortices induce flow fluctuations that are scattered by the trailing edge generating acoustic waves that propagate upstream. These waves are perceived as the white ribbons in marker 3. The traces of the acoustic waves start at $\phi \approx 60^{\circ }$ and can be observed up to $\phi \approx 120^{\circ }$. The black arrows are tangent to their traces, and it is possible to conclude that for the $M_{\infty } = 0.4$ flow, the acoustic waves propagate upstream at a lower speed than for $M_{\infty } = 0.1$ due to the Doppler effect. Moreover, the larger thicknesses of the ribbons reflect the fact that the frequency at which the acoustic waves are generated is lower for $M_{\infty } = 0.4$. The arrows also show the path taken by a particular wave that coincides with the inception of the Kelvin–Helmholtz instabilities for each flow (the latter being represented by marker 2). The black circles in both subplots are positioned in the same spatio-temporal coordinates, emphasizing that the mechanism that creates the primary instability originates from the trailing edge at the same time for the two compressible regimes. However, since the upstream propagation speed of the acoustic wave is lower for the higher Mach number, the perturbation takes longer to reach $x/c=0.65$, where the Kelvin–Helmholtz instability starts. This chord location appears to be the same for both flows and is represented by the vertical dashed line in each subplot. In the region of marker 2, we observe the Kelvin–Helmholtz instabilities that, analogously to what happens with the acoustic waves, display a lower frequency for $M_{\infty } = 0.4$.

Figure 9. History of pressure gradient in the $\xi$ direction over the aerofoil suction side for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$. Panel (c) shows the pressure gradients computed along the horizontal green line, which corresponds to $t=4.25$.

Shih, Lourenco & Krothapalli (Reference Shih, Lourenco and Krothapalli1995) pointed out that although a flow reversal extends from the trailing edge over the aerofoil suction side, there was insufficient time for this flow to reach the leading edge region. Therefore, they postulated that the flow reversal near the leading edge and the eventual initiation of the unsteady separation were essentially local phenomena. The trailing edge, in turn, was considered responsible for influencing the separation indirectly, through the increase of global circulation around the aerofoil by shedding counter-rotating vorticity into the wake. In our compressible simulations, the importance of the trailing edge lies in what appears to be the acoustic triggering of Kelvin–Helmholtz instabilities, but it is unclear whether the trailing edge is of secondary importance to the onset of the DSV. The Kelvin–Helmholtz instabilities occur due to shear, but our results suggest that the trailing edge acoustics might act as an initial disturbance that triggers the formation of such instabilities. Moreover, at marker 1 in figure 9, we see that the onset of the DSV occurs at nearly the same time for the two regimes, as do the acoustic waves that possibly trigger the primary instability. Since these latter propagate at different speeds, this means that the acoustic disturbance does not influence directly the dynamic stall onset. On the other hand, separation occurs only when the primary instability reaches the leading edge region, and the trailing edge seems to be important in the generation of this primary instability. Other studies (Benton & Visbal Reference Benton and Visbal2020) have shown that the turbulent separation from the trailing edge can play an important role in the dynamic stall onset depending on the flow configuration. However, such trailing edge separation is not present in our results.

In order to better assess the differences in pressure gradient between the two flows, we plot the values of the $\xi$ pressure gradient component over the horizontal green line displayed in figure 9, which corresponds to $t=4.25$. The time instant is the same as that shown in figure 8, which displayed a wider region of favourable pressure gradient along the leading edge for $M_{\infty } = 0.4$. This more extensive region of favourable pressure gradient delays the formation of secondary instabilities in the higher Mach number case. From the figure, we can see that the adverse pressure gradient region (positive values) is shifted downstream from the leading edge at higher Mach number. We also notice that the magnitudes of the pressure gradient fluctuations are higher for the $M_{\infty } = 0.1$ case. In addition, high frequency oscillations and a strong peak appear close to the leading edge for $M_{\infty } = 0.1$. These are related to the beginning of the destabilization of the local vorticity distribution described by the van Dommelen and Shen model. Finally, this plot provides another interesting observation about the wavelengths of the Kelvin–Helmholtz instabilities: for the higher Mach number case, longer wavelengths are observed near the trailing edge region (say $x/c > 0.6$), while shorter ones can be seen around the mid-chord. For the lower Mach number flow, on the other hand, the wavelengths appear to be more uniform. Since our results are presented in terms of spanwise-averaged quantities, this fact possibly indicates a higher three-dimensionality of the flow structures close to the trailing edge for $M_{\infty } = 0.4$.

In a recent work by Benton & Visbal (Reference Benton and Visbal2020), simulations of a static NACA 0012 at $8^{\circ }$ angle of attack and $Re = 10^{6}$ showed the presence of a wider spectrum of three-dimensional (oblique) modes in the transition process at the leading edge as compressibility was increased. Based on their observations, we investigate the presence of three-dimensional disturbances in the flow for the current plunging aerofoil. Results are shown in figure 10, where iso-surfaces of $q$-criterion $Q=1/2( |\varOmega |^{2} - |S |^{2} )$, coloured by $x$-momentum, are plotted over the aerofoil suction side. The figure is presented for $t=4.4$, which corresponds to a moment before the ejection of the first leading edge vortex (see third column of figure 6). Results show that the higher Mach number flow is more three-dimensional along the entire aerofoil suction side. Given the quicker breakup of spanwise coherence of Kelvin–Helmholtz instabilities for $M_{\infty } = 0.4$, it could be that non-normality plays an important role in this case. Even the leading edge vortices related to the secondary instability mechanism display spanwise modes at this higher Mach number. On the other hand, the rolls of the $M_{\infty } = 0.1$ flow are mostly two-dimensional, except close to the trailing edge where transition sets in.

Figure 10. Iso-surfaces of $q$-criterion coloured by $x$-momentum for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$, computed at $t=4.4$.

3.3. Analysis of empirical criteria for dynamic stall onset

As discussed by Ramesh et al. (Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards2014), the LESP $A_0$ can be used to indicate the onset of the DSV, which occurs when the maximum value of this parameter is achieved. According to Deparday & Mulleners (Reference Deparday and Mulleners2019), the LESP can be computed directly from numerical or experimental data as

(3.2)\begin{equation} A_{0} = \text{sgn}[\cos(\lambda - \alpha)] \sqrt{ \frac{2}{{\rm \pi} } |\boldsymbol{S}_{LE} | |{\cos(\lambda - \alpha)}|} . \end{equation}

In (3.2), $\alpha$ is the effective angle of attack, and $\lambda$ is the angle formed between the incoming flow and the leading edge suction vector $\boldsymbol {S}_{LE}$. This vector is determined by integrating the pressure signals along the first 10 % of the aerofoil chord. Hence $\boldsymbol {S}_{LE}$ gives the net force from the surface pressure integration. Here, the sign function is included to yield a continuous result when $\cos (\lambda - \alpha )$ goes from positive to negative, and vice versa. A schematic representation of the terms used to evaluate $A_{0}$ according to Deparday & Mulleners (Reference Deparday and Mulleners2019) is presented in figure 11(a). In this figure, the 10 % frontal portion of the aerofoil chord is coloured blue, and the chord line is shown by a dashed line.

Figure 11. Schematic of the relevant parameters used in the calculation of (a) the LESP, and (b) an entropy measure showing the evaluation of the chord-normal shear layer height.

The temporal evolution of the LESP for the periodic aerofoil motion is shown in figure 12(a) for the two Mach numbers investigated. Relevant values of the phase angle $\phi$ are marked in the figure, and the green and purple dots are the LESP thresholds, plotted along with their respective coordinates $(t, A_{0})$. The green dots represent the first local maxima found in $A_0$, while the purple dots show the global maximum. According to the previous references, the global maximum should indicate the instant of the DSV onset. However, we find here that the first local maximum provides a better estimate for the onset. Deparday & Mulleners (Reference Deparday and Mulleners2019) show that after the onset of the DSV, the LESP is susceptible to fluctuations from the separated flow. Moreover, the secondary instabilities develop downstream 10 % of the chord in the present flows. These characteristics could impact on the calculation of the LESP, since they may lead to new maximum values. Hence the green dots will be used to provide the time instant of dynamic stall onset.

Figure 12. (a) LESP, and (b) chord normal shear layer height, for the aerofoil in periodic motion at $M_{\infty } = 0.1$ and $0.4$.

Comparing the blue and red lines from figure 12, which correspond to the values of $A_{0}$ for Mach numbers $M_{\infty } = 0.1$ and $M_{\infty } = 0.4$, respectively, we notice that the maximum suction is reduced as compressibility increases. This is a consequence of the weaker low-pressure core formed close to the leading edge during the onset of dynamic stall, as shown in figure 2, which influences directly the magnitude of the suction vector. Despite this observation, the dynamic stall onset occurs almost at the same time for both Mach numbers.

The applicability and reliability of the LESP hinge upon a prior knowledge of the critical thresholds of the leading edge suction as a function of the parameters of the aerofoil motion. However, the dependency of this suction on the kinematics and flow compressibility impairs its usage as a universal dynamic stall onset indicator for several realistic flow applications. With the aim of overcoming this limitation, the shear layer height criterion proposed by Deparday & Mulleners (Reference Deparday and Mulleners2019) was shown to be a more precise, motion-independent, indicator, although not easily accessible outside of a laboratory environment. In this previous reference, the shear layer height was evaluated through an averaged location of clockwise rotating shear layer vortices identified using an Eulerian vortex criterion from the individual snapshots of the velocity field. Here, as we know the values of all conserved variables from the numerical simulation, it is possible to estimate the chord normal shear layer height, ${\rm \Delta} z$, in a simpler manner. We compute contours of an entropy measure and select a small threshold ($0.05$ in this case) as a limit for the viscous region. A representation of the procedure is shown in figure 11(b), where the blue region indicates the portion of the fluid within the limit where entropy values change above the $0.05$ threshold. Then uniformly distributed chord-normal lines are drawn, and the distances from the aerofoil surface to the entropy contour boundary are measured. These lines are represented in figure 11(b) in yellow, and their distances are non-dimensionalized by the aerofoil chord. Here, 100 lines are used in the analysis, but only a few appear in the figure for better visualization. If a given chord-normal line intersects the limiting entropy level more than once – for instance, in the rightmost yellow line in the figure – then a separate distance is computed for each intersection, and all such distances are subsequently averaged out to obtain ${\rm \Delta} z$.

The temporal evolution of ${\rm \Delta} z$ for the two Mach number flows is shown in figure 12(b). From this figure, we see that during the lifetime of the dynamic stall vortex, the shear layer height exhibits a high-frequency oscillatory behaviour. This is a consequence of locally non-convex topological spaces formed when the flow is highly separated, causing the normal lines to have multiple intersections. This noise, however, does not degrade the capability of the method in finding the critical shear layer height since it is obtained by the intersection of two linear fits. According to Deparday & Mulleners (Reference Deparday and Mulleners2019), these fits are comprised of the primary instability stage displayed by the dark grey line in figure 12, and the vortex formation stage represented by the first slope crossing this line. To draw the linear fits, a time window $3 \le t \le 4$ is used for the dark grey line, while for the vortex formation stage, the windows $4.5 \le t \le 6$ and $4.5 \le t \le 5$ are used for $M_{\infty } = 0.1$ and 0.4, respectively. For our simulations, the instant of formation of the DSV is marked by the intersection of the green with the grey dashed lines shown in figure 12(b).

The threshold value based on the shear layer height ${\rm \Delta} z \approx 0.02$ remains the same for both Mach numbers. According to Deparday & Mulleners (Reference Deparday and Mulleners2019), this critical value of ${\rm \Delta} z$ is invariant with respect to the kinematics. For the present periodical motion, in the absence of shock waves, we find that it is also invariant to the freestream Mach number. Further analysis should be conducted to verify if shock-induced separation would alter this threshold. Figure 12 indicates that at the time instant shown in the abscissa of the first local maximum value of the LESP, a reasonable match is obtained with that computed by the critical shear layer height. However, when the global maximum value of the LESP is used, it indicates a more advanced stage characterized by the roll-up of the shear layer into a large-scale DSV as the onset point. Despite the small time lapse between results from the two different methods, both answers are reasonable and consistent with the skin friction maps of figure 3.

At this point, it is important to mention that to the present date there is no consistent definition of ‘stall onset’ in the literature. For example, McCroskey (Reference McCroskey1981) defines it as the regime where maximum lift is produced; Mulleners & Raffel (Reference Mulleners and Raffel2012), on the other hand, consider the point when the DSV detaches from the aerofoil surface. Here, we do not infer a precise moment to the stall onset, but instead consider a range of possible values starting from the time that the secondary instabilities initiate near the leading edge to that when the DSV detaches. Some studies (see Deparday & Mulleners Reference Deparday and Mulleners2019; Narsipur et al. Reference Narsipur, Hosangadi, Gopalarathnam and Edwards2020) reported that the LESP reaches its critical value well after the onset of instabilities near the leading edge. From the perspective of an effective control strategy, it is important to predict the stall onset before the DSV formation (Chandrasekhara Reference Chandrasekhara2007). Overall, our results for LESP and ${\rm \Delta} z$ criteria agree reasonably well with the features seen in figures 5 and 6 (where the secondary instabilities are shown), in the sense that the critical values lie inside or close to this time interval. With increasing Mach number, however, both criteria indicate that the onset occurs earlier in time. This is consistent with the observations from Chandrasekhara & Carr (Reference Chandrasekhara and Carr1990), who reported that increasing Mach number makes the DSV form at a lower effective angle of attack and further downstream from the leading edge. However, our simulations may not share this trend depending on how one defines the stall onset. By inspection of either figure 5 or figure 6, it is possible to see that the instabilities appear sooner for the Mach 0.1 case, but the detachment of the leading edge vortex comes after that for Mach 0.4. Hence it is possible to conclude that both criteria are more closely related to the definition of stall onset as being the time when the DSV detaches from the aerofoil surface, which could be detrimental potentially to the effectiveness of dynamic stall control.

3.4. Modal analysis

Modal decomposition is now performed using the flow snapshots to identify the relevant flow dynamics related to stall onset in terms of low-rank features. As noted earlier, several different types of such decompositions exist, each with its own benefits and shortcomings. Typically, these approaches assemble snapshots in time and extract modes and corresponding eigenvalues, whose interpretation depends on the underlying theory. The benefits of applying modal decomposition in general, and DMD in particular, to examine dynamic stall may be found in Dunne et al. (Reference Dunne, Schmid and McKeon2016), Mohan et al. (Reference Mohan, Gaitonde and Visbal2016) and Ramos (Reference Ramos2019). For DMD, the modes and corresponding eigenvalues obtained in the manner of §§ 2.2 and 2.3, and Appendix A, represent a decomposition into coherent structures, each associated with a magnitude, growth or decay rate and characterized by a specific frequency.

The transient nature of the current problem, where some phenomena such as stall onset are relatively rapid and dominant in a short time range of the imposed period, motivates the need for special care in the interpretation of the results. Indeed, as shown in this section, application of the original DMD as in the cited references is problematic for the current dataset; this motivates the adoption of the relatively more recent multi-resolution variant (mrDMD), whose properties are more suited to datasets where the condition of a zero-mean second-order stochastic process is not appropriate.

In the standard DMD algorithm, the resulting modes can be ranked by either mode frequency or amplitude. Here, we chose to rank the modes in decreasing order by their amplitudes $b_j$ from (2.10), with $j$ being the index of the mode. A collection of $1024$ snapshots of pressure coefficient is used to build the matrices $\boldsymbol {X}$ and $\boldsymbol {X}'$ (details of the algorithm are provided in Appendix A) within a time period ${\rm \pi} /k$, corresponding to a full cycle of the aerofoil motion, where $k \triangleq {\rm \pi}f c / U_{\infty } = 0.25$ is the reduced frequency. Analogously to the observation of Dunne et al. (Reference Dunne, Schmid and McKeon2016) and Mohan et al. (Reference Mohan, Gaitonde and Visbal2016), when ranked by their amplitudes, the first modes exhibit frequencies that may be associated with the harmonics of the imposed motion. Here, we leave the discussion of these harmonics to Appendix B because, despite giving a global perspective about the effects of increasing compressibility on dynamic stall, the harmonics by themselves do not elucidate the underlying physics associated with the dynamic stall onset. To fill this gap, we search for modes that emphasize the contribution to the dynamics near the leading edge.

Some selected modes are presented in figure 13 together with their respective frequencies. The modes are arranged so that similar events can be compared for the two Mach numbers. Here, we stress that the results obtained from the classical DMD approach provide no guarantee that similar dynamical mechanisms are being compared. Indeed, there is no a priori assumption that similar mechanisms should exist at the two Mach numbers. Thus, by similar events, we mean that the spatial distributions of the modes are comparable. Therefore, the rationale behind the selection of the DMD modes is the presence of a spatial distribution that could be associated with the DSV onset. After obtaining the DMD modes, we selected those whose spatial distributions are related to the primary instability stage or to the separation process. However, since there are several modes that fit these requirements, we select the one with similar spatial distribution and frequency, but largest initial amplitude. We arrive at the four frequencies for each Mach number in figure 13 by simply applying the above procedure starting from the lowest frequency. Then the pairing frequencies between the two Mach number flows appear naturally in the course of the process. Note that the growth rates of the modes are not shown because they are not representative of a linear stability analysis in this case, and saturation yields eigenvalues close to the unit circle.

Figure 13. DMD modes for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$, for the periodic motion.

The standard DMD algorithm essentially combines data from all snapshots; for example, as shown in figure 13, all modes containing coherent structure information near the leading edge also include the influence of the flow from post-stall stages. Coloured regions away from the aerofoil surface in all subplots are a remnant indication of the passage of the DSV. They arise from the fact that the extracted modes exist throughout the entire temporal window analysed, which impairs the capacity of the algorithm to highlight events that may be dominant only in a finite time window of the cycle. The signature of such events is effectively spread out over the cycle and complicates the interpretation of the DMD modes. This issue is also related to the inefficacy of SVD based approaches to extract low-rank dynamics from data with translational and rotational invariance, as discussed in the work of Kutz et al. (Reference Kutz, Fu and Brunton2016), where the mrDMD algorithm is proposed.

Despite these limitations, some useful observations can still be drawn from these modes. For instance, the intense structures close to the aerofoil surface, especially around the leading edge region, are shifted slightly downstream for the $M_{\infty } = 0.4$ case. The results also suggest that the frequencies of the structures associated with dynamic stall onset between the two Mach number flows are similar and range from $St \sim 1$ to $St \sim 40$. In figure 13, only a few modes are shown, so it does not cover all frequencies mentioned in that range. Moreover, these values are not precise as they are derived from a nonlinear collection of snapshots, and provide only a perception of the order of magnitude of the frequencies at play during the dynamic stall onset. The persistence of the flow structures and their corresponding frequencies across the compressible regimes investigated indicate a possible correlation of mode shapes from a parametric modal decomposition perspective, as reported by Coleman et al. (Reference Coleman, Thomas, Gordeyev and Corke2019).

The limitations of the standard DMD are now lifted by applying the mrDMD algorithm proposed by Kutz et al. (Reference Kutz, Fu and Brunton2016). Using the same procedure as before and the binning technique of § 2.3, some of the most interesting modes are presented in figure 14. In the multi-resolution approach, the modes are no longer ranked by their amplitude so that, instead of showing the index of a given mode, it is more suitable to discuss the resolution level and the bin where it resides. For clarity, the levels and bins are converted into their equivalent time periods ${\rm \Delta} t$ and the modes are judiciously organized in the figure so that the same period of time is assimilated between the two different Mach number flows.

Figure 14. DMD modes obtained from the mrDMD algorithm for (a) $M_{\infty } = 0.1$, and (b) $M_{\infty } = 0.4$, for the periodic motion.

Comparing figures 13 and 14, it is evident that the multi-resolution variant has a better capability to deal with the transient nature of dynamic stall and is able to extract clean and physically interpretable modes. Starting from the leftmost column of figure 14, we observe the structures that represent the Kelvin–Helmholtz instability formed along the aerofoil suction side. Note that the time period encompassed by this mode coincides with the time window when the instability develops, as shown in figure 4 and also the pressure gradient plot of figure 9. For the $M_{\infty } = 0.4$ case, acoustic waves are also seen propagating from the trailing edge due to a scattering mechanism, but these are not visible for the $M_{\infty } = 0.1$ case since the acoustic waves have a lower amplitude for this flow. The spatial distribution of the modes is also different. While a quasi-continuous train of Kelvin–Helmholtz disturbances appears as a wavepacket at $M_{\infty } = 0.1$, the structures are more sparsely distributed for $M_{\infty } = 0.4$, being more concentrated at the mid-chord. This is consistent with the observation of high frequency oscillations in the skin friction and pressure gradient maps of figures 3 and 9, respectively. It is also interesting to note that the frequencies between the two compressible cases are much alike, and this holds for all subplots.

The modes depicted in the subsequent columns of figure 14 are associated with the growth of the shear layer instabilities close to the leading edge and the separation process that leads to the formation of the DSV. In the second column, the modes represent the advection of vortical structures formed in the leading edge shear layer instability, which coalesce afterwards to form the DSV. Finally, the last two columns refer to the same time period and mark the separation process described earlier, in figures 5 and 6, where nonlinear mechanisms govern the evolution of the flow. As indicated in figure 14, the frequencies increase substantially at this stage, independently of the Mach number. While the structures observed in the last two columns are more organized and compact for $M_{\infty } = 0.1$, those appearing in the $M_{\infty } = 0.4$ flow are composed of a wider range of spatial scales. This observation may reflect the higher three-dimensionality of the $M_{\infty } = 0.4$ flow, as shown in figure 10. As discussed before, the structures are also formed downstream for the $M_{\infty } = 0.4$ case and extend over a larger portion of the aerofoil chord. For these DMD modes, one can also see that acoustic waves are generated by the hydrodynamic fluctuations near the leading edge. These waves have higher wavenumbers for the $M_{\infty } = 0.4$ flow and may contribute to the earlier breakdown into three-dimensional structures.