2.1. Problem definition

We consider a spanwise homogeneous NACA 0012 airfoil with an AoA of $\alpha _0=15^\circ$, undergoing sinusoidally plunging (or heaving) motion in a uniform incoming flow. The position $h$ of the leading edge moves vertically as

(2.1)\begin{equation} h(t) = \frac{A}{2}\cos(2{\rm \pi} t/T), \end{equation}

where $A$ is the peak-to-peak amplitude and $T$ is the plunging period. The plunging amplitude $A$ is fixed at 0.5 chord length in this work. We non-dimensionalise the problem based on the chord length, $c$, and incoming flow velocity, $U_{\infty }$ and so introduce the dimensionless plunging amplitude, $A_c = A/c = 0.5$, define the chord-based Strouhal number as $St_c = c/(U_\infty T)$ and define the Reynolds number as $Re = U_\infty c /\nu$, where $\nu$ is the kinematic viscosity of the fluid. The effective AoA of the airfoil is therefore

(2.2)\begin{equation} \alpha_e(t) = \alpha_0 +\tan^{-1}\left({\rm \pi} A_c St_c \sin(2{\rm \pi} t/T)\right), \end{equation}

which varies from $\alpha _0-\tan ^{-1}({\rm \pi} A_c St_c)$ to $\alpha _0+\tan ^{-1}({\rm \pi} A_c St_c)$. The amplitude-based Strouhal number $St_A = A / (U_\infty T) = St_c A_c$ is also used. The force coefficients are non-dimensionalised using $0.5\rho U_\infty ^2 c b$, where $\rho$ is fluid density and $b=5c$ is the spanwise length of the airfoil.

2.2. Experimental set-up

The force is measured experimentally in the water tunnel facility of the University of Bath. The test section of the water tunnel is $381\times 508\times 1530$ mm$^3$ and the flow velocity range is 0–0.5 m s$^{-1}$ with a turbulence intensity of less than 0.5 %. The wing is vertically mounted in the test section via the translation stage, which generates a linear motion perpendicular to the freestream flow (see figure 1). The set-up is identical to that used in our previous investigations (Son et al. Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022). Two end plates are used at the root and the tip of the airfoil to approximate the 2-D cases. The chord length of the wing is $c = 62.7$ mm and the aspect ratio is ${\rm AR} = 5$. The wing is subjected to a sinusoidal plunging motion, as in (2.1). Reynolds number, based on the chord length, is $Re=10\,000$ for the experiments. The wing has a NACA 0012 profile and is manufactured from PA-2200 polyamide using selective laser sintering with a polished smooth surface. A T800 carbon fibre bar with a cross-section of $25\times 5$ mm$^2$ is inserted inside the wing to provide high spanwise stiffness to the wing. A Zaber LSQ150B-T3 translation stage powered by a stepping motor provides a sinusoidal plunging motion with an accuracy of 2 % of the peak-to-peak amplitude. Force measurements are performed with an ATI Mini 40 force/torque sensor. The data are collected for 50 periods of plunging motion and the data acquisition frequency is set for each case to have 2000 data points in a plunging period. For the stationary cases, the data are collected at 1000 Hz for 60 seconds. For the mean lift and drag, the data were time-averaged whereas for the time-dependant lift and drag, a band-stop filter was applied to the data to eliminate the vibrations at the natural frequencies of the wings. Then a moving-average filter was applied for 100 points and the data were phase-averaged for 50 plunging cycles.

Figure 1. Experimental set-up.

2.3. Numerical method

Numerical simulations are conducted at Reynolds numbers from 400 to $10\,000$. An iLES technique, which converges to the DNS for sufficient grid resolution, is adopted. At Reynolds numbers lower than or equal to 800, the numerical results have the DNS resolution. For higher Reynolds numbers, simulations should be considered as iLES where the high-wavenumber modes are damped by the spectral vanishing viscosity (Moura et al. Reference Moura, Aman, Peiró and Sherwin2020).

The governing equations are the incompressible Navier–Stokes (N-S) equations and the continuity equation. To avoid the need for a moving grid, the relative velocity is solved in the body frame of reference. They are expressed as

(2.3a)$$\begin{gather} \partial_t \boldsymbol{u}+ \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} =-\rho^{-1}\boldsymbol{\nabla} p + \nu\nabla^2\boldsymbol{u} - \ddot{h}\boldsymbol{e}_y, \end{gather}$$

(2.3b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

where the $-\ddot {h}\boldsymbol {e}_y$ term captures the effect of the plunging motion. On the airfoil surface, a no-slip boundary condition $\boldsymbol {u}=0$ is imposed. In the far field, the flow velocity is prescribed to be uniform in space,

(2.4)\begin{equation} \boldsymbol{u} = U_{\infty}\boldsymbol{e}_x - \dot{h}\boldsymbol{e}_y \end{equation}

and accounts for the sinusoidal motion of the frame of reference.

The momentum equation (2.3a) and the continuity equation (2.3b) are solved using a high-order velocity-correction scheme (Karniadakis, Israeli & Orszag Reference Karniadakis, Israeli and Orszag1991). A second-order implicit–explicit scheme is applied for time integration, in which the nonlinear advection term is treated explicitly and the viscous term implicitly. At least $20$ plunging periods are simulated to ensure that the flow is fully developed. The last 10 periods are used for analysis.

For 2-D simulations, the computational domain is spatially discretised with an elemental decomposition using the spectral/hp element method, implemented in the open-source code Nektar++ (Cantwell et al. Reference Cantwell2015; Moxey et al. Reference Moxey2020). The size of the computational domain is $[-40c, 60c]\times [-40c, 40c]$. The domain is tessellated by 9921 quadrilaterals and 64 triangles, where the height of the layer of elements adjacent to the airfoil surface is $0.005c$. The airfoil trailing edge is rounded with a radius of $0.00145c$. The 2-D computational domain and spectral element mesh are shown in figure 2. Within each element, the solution is represented in terms of a tensor product of modified Jacobi polynomials of order $P$, with Gauss–Lobatto–Legendre points for performing numerical integration. The grid convergence study in Appendix A.1 indicates that using $P=4$ is sufficient for the range of flow configurations in this study. A reduced mesh with 3020 elements and an expansion order of $P=5$ used by Son et al. (Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022), where the grid convergence has been tested, is also applied here for $Re=10\,000$ simulations.

Figure 2. Two-dimensional computational domain and spectral/hp element mesh used for the numerical simulations. (a) Global view showing the full domain. (b) Close-up view of mesh around the airfoil.

For 3-D simulations, the 2-D computational mesh described previously is extended to a 3-D computational discretisation using a Fourier expansion with 128 planes in the $z$-coordinate. Spectral dealiasing is applied to reduce the nonlinear error. The highest wavenumber that can be captured is therefore 63. The grid convergence in the $z$-direction has been assessed in Son et al. (Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022) and is also presented in § A.3.

The boundary conditions are enforced numerically as follows. On the wall, a no-slip condition is used for the velocity, while a high-order Neumann condition is used for the pressure (Orszag, Israeli & Deville Reference Orszag, Israeli and Deville1986). On the upstream boundary a Dirichlet condition, imposing (2.4), is used for velocity, while a zero-Dirichlet condition is used for pressure. The downstream boundary uses a high-order outflow boundary condition (Dong, Karniadakis & Chryssostomidis Reference Dong, Karniadakis and Chryssostomidis2014). The cross-stream boundaries are periodic. A spectral vanishing viscosity (Moura et al. Reference Moura, Aman, Peiró and Sherwin2020) is applied with a coefficient of 0.1 to stabilise the simulation.

In the visualisation of flow fields at $Re\ge 1000$ where the numerical results have the LES resolution, a Gaussian filter is applied to remove small-scale vortices that are affected by numerical errors and only show the more physically realistic large-scale vortices. Based on the performance of the spatial filter in the work of Son et al. (Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022), the standard deviation of the Gaussian filter is taken as $0.04c$, which damps the wave amplitude by half and more at a wavelength smaller than $0.2c$.

2.4. Linear stability analysis

In linearising the N-S equations (2.3), the flow is considered as an infinitesimal perturbation $\tilde {\boldsymbol {u}}$ around a 2-D time-periodic base flow $\boldsymbol {U}$, with period $T$,

(2.5a,b)\begin{equation} \boldsymbol{u} = \boldsymbol{U} + \tilde{\boldsymbol{u}}, \quad p = P + \tilde{p}. \end{equation}

Since the base flow $\boldsymbol {U}$ satisfies the N-S equations (2.3), eliminating these terms and perturbation terms of higher than linear order leaves the linearised N-S equations,

(2.6a)$$\begin{gather} \partial_t \tilde{\boldsymbol{u}} + \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\boldsymbol{u}} + \tilde{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U} =-\rho^{-1}\boldsymbol{\nabla} \tilde{p} + \nu\nabla^2\tilde{\boldsymbol{u}}, \end{gather}$$

(2.6b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}} = 0. \end{gather}$$

describing the evolution of the perturbation. For the linearised problem, the condition $\tilde {\boldsymbol {u}}=0$ is imposed on all boundaries.

The base flows for the required flow configurations are computed by solving (2.3) using DNS. For each configuration, the solution over one time period is stored using 128 snapshots. During linear stability analysis, a fourth-order Lagrange interpolation is used to generate the instantaneous base flow $\boldsymbol {U}$ at the required time through the cycle.

Based on Floquet theory, solutions to (2.6) can be expressed as

(2.7)\begin{equation} \tilde{\boldsymbol{u}}(x, y, z, t_0 + n T) = \mu^n\tilde{\boldsymbol{u}}(x, y, z, t_0). \end{equation}

for integer $n$, where $\mu$ are the Floquet multipliers. For $|\mu |>1$, the corresponding mode grows exponentially and is therefore unstable.

The disturbance may be 3-D and due to the symmetry of the problem a further simplification is possible, whereby the solution can be represented using only a half Fourier mode in the $z$ direction,

(2.8a)$$\begin{gather} \tilde{\boldsymbol{u}}(x,y,z,t) = \left(\hat{u} \cos(\beta z), \hat{v} \cos(\beta z), \hat{w} \sin(\beta z)\right), \end{gather}$$

(2.8b)$$\begin{gather}\tilde{p}(t, x, y, z) = \hat{p}\cos(\beta z), \end{gather}$$

where $\hat {\boldsymbol {u}} = (\hat {u}, \hat {v}, \hat {w})$, and $\hat {p}$ are functions of $(x, y, t)$ only and $\lambda =2{\rm \pi} \beta ^{-1}$ is the wavelength in the $z$ direction.

The task of computing the Floquet multiplier that has the maximum absolute value can be formulated as an eigenvalue problem and solved using a Krylov-based iterative Arnoldi method (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008). While the disturbance field is necessarily simulated using the entire computational domain, only the near-body region is used for the calculation of the Floquet multipliers. The maximum dimension of the Krylov subspace is 16 and the iterative tolerance is set to $10^{-4}$. Appendix A.2 demonstrates that the eigenvalue results are independent of the chosen domain restriction.

2.5. Verification and validation

The time-averaged lift and drag coefficients, non-dimensionalised by $0.5 \rho U_\infty ^2 c b$, from experiments and numerical simulations at $Re=10\,000$, are listed in table 1. The lift and drag coefficients from numerical simulations roughly agree with those from the experimental measurements, although a relative difference of about 10 % is observed. The discrepancy in the drag coefficient is larger when the drag is close to zero or transits to thrust for $St_c\ge 0.64$, which may be caused by the large transient drag at high plunging frequencies.

Table 1. Time-averaged force coefficients of the plunging airfoil from numerical simulation and experimental measurement at $Re=10\,000$, $A/c=0.5$ and $\text {AoA} = 15^\circ$.

The time-dependent force coefficients at $St_c=0.32$ and $St_c=0.64$ are also compared in figure 3, where lines represent experimental results and symbols are numerical results. Good agreement of the force is found when the position of the airfoil is around the equilibrium point, however, notable discrepancies exist when the airfoil is at the lowest ($t/T=0.5$) and highest ($t/T=0$) positions. At $St_c=0.32$, $t/T=0.5$, the numerical simulation exhibits small oscillations of the force, but the experimental data is very smooth. At $St_c=0.64$, $t/T=0$ and $t/T=0.5$, the magnitude of lift from the numerical simulation is lower than that from the experimental measurement by 20 %. However, because the signs of these discrepancies are different at the highest ($t/T=0$) and at the lowest ($t/T=0.5$) position, the relative difference of the time-averaged force is 10 %.

Figure 3. Time-dependent (a) lift and (b) drag coefficients at $Re=10\,000$, $A/c=0.5$ and $\text {AoA} = 15^\circ$. Phase averaging is applied.

In addition, the spanwise vorticity field is compared between the numerical simulation and experimental data from Son et al. (Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022), see figure 4. Both phase averaging and spanwise averaging are applied. In experiments, because the laser is placed above the airfoil, the flow field below the airfoil is not illuminated and therefore the corresponding vorticity field is not available. Generally, the vortex patterns and locations of main vortices, such as the LEV, the secondary vortex caused by the LEV–wall interaction and the TEV, agree between the numerical simulation and experimental results. In experiments, however, the shear layer is much thicker and the size of the vortex identified by the contour plot is also larger.

Figure 4. Phase-averaged and spanwise-averaged vorticity field from numerical simulation and experimental measurements at $Re=10\,000$, $A/c=0.5$ and $\text {AoA} = 15^\circ$: (a,e) $St_c=0.16$, $t/T=0$; (b,f) $St_c=0.16$, $t/T=0.25$; (c,g) $St_c=0.32$, $t/T=0$; (d,h) $St_c=0.32$, $t/T=0.5$. The position of the airfoil is shown by a black dot on a black line. Panels (a–d) are numerical results. Panels (e–h) are experimental results from Son et al. (Reference Son, Gao, Gursul, Cantwell, Wang and Sherwin2022).

4.1. Floquet results

Linear stability analysis is conducted using the periodic 2-D base flows to identify the onset of the three-dimensionality and the transitional nature of the flow. A convergence study of the base flow and the sensitivity of the Floquet multipliers is reported in Appendix A and validation against previously published results for the NACA airfoils at AoA of $20^\circ$ at $Re=500$ are provided in Appendix B. Since the 2-D base flow deviates significantly from the spanwise-averaged 3-D flow for $St_c>0.70$, the linear stability analysis mainly considers cases where $0.50^+\leq St_c\leq 0.70$, which corresponds the optimal amplitude-based Strouhal number $St_A$ of 0.25 to 0.35.

The absolute value of the Floquet multiplier $|\mu |$ as a function of the spanwise wavenumber $\beta =2{\rm \pi} \lambda ^{-1}$ at $Re=225$, $Re=400$ and $Re=800$ is shown in figure 16, where $\beta$ is scaled by $c^{-1}$. The most amplified dimensionless wavelength $\lambda c^{-1}$ and the Floquet multiplier are also provided in table 2. For the static airfoil shown by the green line in figure 16(a), only a short-wavelength (high-wavenumber) period-doubling mode is observed to be unstable, which is consistent with findings of Deng et al. (Reference Deng, Sun and Shao2017) at $\text {AoA}=15^\circ$. The most amplified wave length is $0.4c$ at $Re=800$. For plunging cases in regime IV, a long-wavelength (low-wavenumber) period-doubling mode is the only unstable mode. The growth rate of this period doubling mode increases with the plunging Strouhal number $St_c$ and the Reynolds number as can be seen from figures 16(b) and 16(c) where we plot the Floquet multipliers at $St_c = 0.50^+$ and $St_c = 0.64$. The most amplified wavelength, which is around $2c$, increases with the Reynolds number. At slightly higher $St_c$, a short wave mode also appears but this mode is stable under the current Reynolds number.

Figure 16. Floquet multiplier vs dimensionless wavenumber $\beta c$ from linear stability analysis. (a) Varying with $St_c$ at $Re=800$. From high to low, $St_c = 0.0, 0.70, 0.64$ and $0.50^+$. (b) Varying with Reynolds number at $St_c=0.50^+$. From high to low, $Re=800, 400$. (c) Varying with Reynolds number at $St_c=0.64$. From high to low, $Re=800,400,225$.

Table 2. Comparison of the most amplified wavelength with Floquet multiplier in brackets, $\lambda c^{-1}(\mu )$, for the static case and the plunging cases at $St_c = 0.50^+$ and $St_c = 0.64$.

The wavelength and period-doubling feature predicted by the linear stability is confirmed by the current and previous (Gao et al. Reference Gao, Sherwin and Cantwell2022) 3-D simulations, which have spanwise wavelengths of $1.7c$ or $2.5c$ in this parametric range. Since the spanwise length of the flow is $5c$, only spanwise wavelengths of $\lambda =5c/n$, such as $5c$, $2.5c$, $1.7c$ and $1.25c$, can be efficiently captured.

Figure 17 shows the dimensionless amplitude function of the spanwise vorticity, $\hat {\omega }_z^* = (\partial _x\hat {v}-\partial _y\hat {u}) c/U_\infty$, of the most unstable mode at $Re=800, St_c=0.50^+, \lambda =2.2c$ ($\beta = 2.9c^{-1}$). The effective AoA of the airfoil varies from $-23^\circ$ to $53^\circ$. The dimensionless spanwise vorticity of the base flow, $\varOmega _z^*=\varOmega _z c/U_\infty$, is also shown by solid lines (positive value) and dashed lines (negative value). The period-doubling feature of the instability can be seen in figure 17(a,c), which are separated by one period of the plunging motion. When the airfoil is at the lowest displacement at $t/T=0.5$, the disturbance vorticity is almost invisible in the LEV region, R0, whereas it is very strong in the pLEV region, R2, as well as the TEV region. The disturbance vorticity in each concentrated vortex has two subregions of different signs, which indicates the concentrated vortex experiences a bending displacement in this unstable mode. The bending orientations are different for the pLEV and the TEV. Even for the LEV itself, the bending orientation also changes with time. When the airfoil is at the highest displacement at $t/T=1$, disturbance vorticity appears in the LEV region, R1. Figure 18 shows the 3-D reconstruction of the base flow and the most unstable mode, where (2.8a) was used for the reconstruction of the perturbations. The kinetic energy of the disturbance flow is 9 % of that of the base flow in the region $-0.2c\leq x\leq 3c,\ -c\leq y\leq 1.4c$. The bending displacement of the counter-rotating pLEV–TEV pair, which resembles the Crow-type instability, and the period-doubling feature are clearly shown.

Figure 17. Amplitude function of the spanwise vorticity of the most unstable mode at $Re=800, St_c=0.50^+, \lambda =2.2c (\beta = 2.9c^{-1})$: (a) $t/T=0.5$, (b) $t/T=1$ and (c) $t/T = 1.5$. Black lines are the contours of the spanwise vorticity of the base flow, $\varOmega _z^*=-3, -1, -0.3, 0.3, 1, 3$ with negative values represented by dashed lines. Green rectangles show the LEV regions.

Figure 18. Three-dimensional reconstruction of the most unstable mode at $Re=800, St_c=0.50^+, \lambda =2.2c (\beta = 2.9c^{-1})$. Isosurfaces are the dimensionless $Q^*=50$ coloured by the streamwise vorticity. Inset plots in (a,c) shows the pLEV–TEV pair from a different view.

4.2. Linear stability mechanism

To further investigate the instability mechanism of the long-wavelength, period-doubling mode, the enstrophy of the disturbance flow, which is a localised quantity, is analysed. The governing equation of the enstrophy can be obtained as the inner product of disturbance vorticity $\tilde {\boldsymbol \omega }=\boldsymbol {\nabla }\times \tilde {\boldsymbol {u}}$ and the curl of (2.6a). For a 2-D base flow, it is

(4.1)\begin{equation} \partial_t \tilde{\mathcal{E}} =- \boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\mathcal{E}} - \tilde{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\varOmega_z\tilde{\omega}_z + \tilde{\boldsymbol\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U}\boldsymbol{\cdot}\tilde{\boldsymbol\omega} + \varOmega_z\partial_z\tilde{\boldsymbol{u}}\boldsymbol{\cdot}\tilde{\boldsymbol\omega} + \nu\tilde{\boldsymbol\omega}\boldsymbol{\cdot} \nabla^2\tilde{\boldsymbol\omega}. \end{equation}

Here $\tilde {\mathcal {E}}=0.5\tilde {\boldsymbol \omega }{\cdot }\tilde {\boldsymbol \omega }$ is the enstrophy of the disturbance flow and $\varOmega _z$ is the vorticity of the base flow. The physical meanings of the terms on the right-hand side are listed in table 3.

Table 3. Physics interpretation of right-hand side terms of enstrophy equation (4.1).

In the following, we calculate the magnitude of each term integrated over regions of the flow during the evolution of the instability to analyse the dominating mechanism in flow instability. To this end, it is convenient to consider two Eulerian integration regions. The first is a region around the leading edge, $\{(x,y)\,|\,x^2+y^2\le 0.36c^2\}$, which only contains the LEV during the early stage of the downward half period. The second is a near-body region, $\{(x,y)\,|\,-c\le x \le 1.4c, -2.5c\le y\le 2.5c\}$, which contains the near-body vortices, including the LEV, TEV and pLEV, and excludes vortices from earlier periods. The integral of each term in (4.1) within these two regions are shown in figure 19 where figure 19(a) shows the leading-edge region and figure 19(b) shows the near-body region. In both the leading-edge and the near-body regions, the term $-\tilde {\boldsymbol {u}}\boldsymbol{\cdot}\boldsymbol {\nabla }\varOmega _z\tilde {\omega }_z$ is the main contribution to the change in enstrophy as indicated by the growth of the red solid line and the ${\rm d}\tilde {\mathcal {E}}/{\rm d} t$ line (in black). This term represents the mechanism where the vorticity of the base flow is displaced by the inducing of the disturbance vorticity. Figure 19(b) shows the distribution of the induction term at $t/T = 0.625$ when it reaches the maximum magnitude. In the bending direction of the pLEV and TEV, the induction term $-\tilde {\boldsymbol {u}}\boldsymbol{\cdot}\boldsymbol {\nabla }\varOmega _z\tilde {\omega }_z$ is positive which further enhances the bending displacement; perpendicular to the bending direction, it is negative. We further observe that viscous diffusion is the main stabilisation effect. The advection term $-\boldsymbol {U}\boldsymbol{\cdot}\boldsymbol {\nabla }\tilde {\mathcal {E}}$ becomes notable when the vortex leaves the integration region but otherwise, all other terms are very weak.

Figure 19. Evolution of the terms in the disturbance enstrophy equation (4.1). (a) Integration over the leading-edge region $x^2+y^2\le 0.36c^2$. (b) Integration over the near-body region $-c\le x\le 1.4c,\ -2.5c\le y\le 2.5c$. The subplot shows the distribution of the induction term at $t/T=0.625$.

The eigenvector of the most unstable mode can also be activated in different zones to elucidate the key process that sustains this absolute instability. To achieve this, the disturbance velocity in a specific region is retained whilst the initial disturbance velocity outside of this region is set as zero. We then consider the enstrophy growth triggered by different parts of the initial eigenvector and compare it with the growth when using the full eigenvector.

As shown in figure 17(a), regions R0 and R2 are selected when the airfoil is at the lowest position. A further region R1, as indicated in figure 17(b), is selected when the airfoil is at the highest position. Using the disturbance flow activated in these three regions as well as in the whole domain, we plot the integral of the perturbation enstrophy in the leading-edge and near-body regions over three base flow periods in figure 20.

Figure 20. Evolution of the enstrophy with different initial conditions. (a) Integration over the leading-edge region. (b) Integration over the near-body region.

Since the masked initial conditions do not satisfy incompressibility, there is a short adjustment period in the initial evolution stage. However, after this initial adjustment, the regions R1 and R2 efficiently capture the whole enstrophy, particularly in the near-body region. In the leading-edge region, R2 captures the whole domain slightly better than region R1, however, we observe that the initial condition in region R0 leads to a reduction in the enstrophy magnitude over the first 1.5 periods.

In terms of the leading-edge region, the fast decay of disturbance triggered by R0 implies that the disturbance in the LEV does not grow by itself. Instead, pLEV is the more important region captured by R1 and R2, since the disturbance vorticity in these regions induces the LEV in region R0 to deform thereby triggering the next stage of instability.

This absolute instability mechanism can therefore be summarised as follows. Starting from $t/T=0.5$, when the airfoil is at the lowest position, the disturbance in the pLEV (R2 of figure 17a) is key since it induces a disturbance in the current LEV. During the downward motion of the first half period, the Crow-like instability of the counter-rotating pLEV–TEV pair amplifies the disturbance. The displacement of the pLEV continues to grow and the vortex centreline of the TEV also bends backwards. The new LEV is also influenced by the displacement of the pLEV–TEV vortex pair, leading to a bending motion in a direction opposite to the motion of pLEV, which results in the period doubling of the disturbed flow. The key amplification mechanism in this unstable mode is therefore the Crow instability of the pLEV–TEV vortex pair, and the self-sustaining mechanism of the disturbance is the interaction of the LEV with the pLEV–TEV vortex pair which requires the pLEV–TEV vortex pair being sufficiently close to the LEV during the upwards half period. This also explains why the base flow in regime III does not experience the same instability, since the pLEV–TEV vortex pair does not remain sufficiently close to the LEV.

5.1. Decomposition of the aerodynamic force

The force element theory (Chang Reference Chang1992; Lee et al. Reference Lee, Hsieh, Chang and Chu2012), which is based on the volume integral of the projection of the Lamb vector ($\boldsymbol \omega \times \boldsymbol {u}$) and originates from the work of Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983) and Chang (Reference Chang1992) is a powerful tool for identifying the main vortex structures that affect the aerodynamic force exerted on the body. Applications of this theory can be found in the literature such as Martín-Alcántara, Fernandez-Feria & Sanmiguel-Rojas (Reference Martín-Alcántara, Fernandez-Feria and Sanmiguel-Rojas2015), where they found the thrust is mainly generated by the growing LEV, and Moriche et al. (Reference Moriche, Sedky, Jones, Flores and García-Villalba2021b), where they found the main LEV contributes most to the force during the manoeuvre of the airfoil. Martín-Alcántara & Fernandez-Feria (Reference Martín-Alcántara and Fernandez-Feria2019) compared the force element theory with theories based on the integral of the Lamb vector and the momentum, and they found the force element theory always achieves better accuracy in predicting the total force.

One inconvenience of the force element theory is that the Lamb vector depends on the choice of frame of reference, making the contour plot of the volume force element, i.e. the projection of Lamb's vector, as well as the force decomposition result, not Galilean invariant (Gao & Wu Reference Gao and Wu2019). This problem becomes more serious in moving-body problems where no special frame of reference exists. To overcome this inconvenience, a modified force-element theory (Gao et al. Reference Gao, Zou, Shi and Wu2019) that is based on the weighted integral of the second invariant of the velocity gradient tensor is adopted in this work.

In the modified force element theory, the force can be expressed as

(5.1)\begin{equation} \boldsymbol{F}\boldsymbol{\cdot}\boldsymbol{e}_j = \int_{{Vol}_\infty}{2\rho Q \phi_j \,{\rm d} V} + \oint_{\partial B}{\rho\phi_j\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{a} \,{\rm d} S} + \oint_{\partial B}{-\rho\nu\phi_j\boldsymbol{n}\boldsymbol{\cdot}\nabla^2\boldsymbol{u} \,{\rm d} S} + \boldsymbol{F}_{f}\boldsymbol{\cdot}\boldsymbol{e}_j, \end{equation}

where $\boldsymbol {e}_j$ is the unit vector in the $j$th direction, $\partial B$ is the body surface and ${Vol}_\infty$ is the whole flow domain. Furthermore, $Q=-0.5\boldsymbol {\nabla }\boldsymbol {u}:\boldsymbol {\nabla }\boldsymbol {u}$ is the second invariant of the velocity gradient tensor, $\boldsymbol {a}$ is the acceleration of the Lagrangian points on the body surface in the inertial frame of reference, $\boldsymbol {n}$ is the unit inward-pointing normal vector on the body surface and $\boldsymbol {F}_{f}$ is the viscous friction force.

The function $\phi _j$ introduced in the work of Quartapelle & Napolitano (Reference Quartapelle and Napolitano1983) is an auxiliary harmonic function which satisfies $\nabla ^2\phi _j = 0$ in the fluid, $\boldsymbol {\nabla }\phi _j\boldsymbol{\cdot}\boldsymbol {n} = -\boldsymbol {e}_j\boldsymbol{\cdot}\boldsymbol {n}$ on the body surface, and $\phi _j=0$ in the far field at infinity. Function $\phi _j$ represents the sensitivity of the aerodynamic force on the spatial location and it only depends on the body geometry in a flow of constant density. The derivation and application of (5.1) can be found in the work of Gao et al. (Reference Gao, Zou, Shi and Wu2019), Gao & Wu (Reference Gao and Wu2019) and Menon & Mittal (Reference Menon and Mittal2021a,Reference Menon and Mittalb). Here, we briefly explain the physical meaning of each term on the right-hand side of (5.1).

Since the pressure field in incompressible flow satisfies $\nabla ^2 p = 2\rho Q$, the first term on the right-hand side of (5.1), which is denoted by $\boldsymbol {F}_{Q}$, represents the force contributed by the source of the pressure Poisson equation,

(5.2)\begin{equation} \boldsymbol{F}_{Q}\boldsymbol{\cdot}\boldsymbol{e}_j = \int_{{Vol}_\infty}{2\rho Q \phi_j \,{\rm d} V}, \end{equation}

with a positive $Q$ generating an attractive pressure force and a negative value generating a repulsive pressure force on the body. The corresponding force coefficient, scaled by $0.5\rho U_\infty ^2$, is denoted by $(C_{D, Q}, C_{L,Q})$ and the integrand $2\phi _j Q$ was named the volume force element by Chang (Reference Chang1992). Since $Q$ is widely used in the identification of vortical structures, this term was named vortex-induced force by Menon & Mittal (Reference Menon and Mittal2021a). For 3-D simulation results, the spanwise-averaging, which is denoted by an overline, is applied to the $Q$ field

(5.3)\begin{equation} \bar{Q} = \int{Q(x,y,z) \,{\rm d} z}, \end{equation}

which, as well as $2\rho \bar {Q}\phi _j$, can be shown as a contour plot on the $x$–$y$ plane.

The second term on the right-hand side of (5.1), which is denoted by $\boldsymbol {F}_a$, represents the force contributed by the body acceleration. The corresponding force coefficient scaled by $0.5\rho U_\infty ^2$ is denoted by $(C_{D, a}, C_{L,a})$. For the translational motion of a rigid body where the acceleration $\boldsymbol {a}$ does not depend on the spatial location, the body-acceleration-related force $\boldsymbol {F}_a$ equals the added-mass force of the potential flow, which is

(5.4)\begin{equation} \boldsymbol{F}_{a} = \boldsymbol{e}_j\oint_{\partial B}{\rho\phi_j\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{a} \,{\rm d} S} =-a_i M_{ij} \boldsymbol{e}_j, \end{equation}

with

(5.5)\begin{equation} M_{ij} =-\rho\oint_{\partial B}{\phi_j n_i\,{\rm d} S} \end{equation}

being the added-mass matrix. For the NACA 0012 airfoil with a chord length of $c=1$ and at ${\rm AoA}=15^\circ$, $M_{ij}$ is calculated as

(5.6)\begin{equation} M = \rho\left(\begin{array}{@{}cc@{}} 0.062066 & 0.192912 \\ 0.192912 & 0.730332 \end{array}\right). \end{equation}

In a periodic flow, the added-mass force has a zero-mean value and it has no contribution to the time-averaged force.

The third term on the right-hand side of (5.1), which is denoted by $\boldsymbol {F}_{p, vis}$, represents the viscous effect in the Neumann-type pressure boundary condition on the airfoil surface $\partial B$,

(5.7)\begin{equation} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla} p =-\rho\boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{a} + \rho\nu\boldsymbol{n}\boldsymbol{\cdot}\nabla^2\boldsymbol{u}. \end{equation}

The corresponding force coefficient, scaled by $0.5\rho U_\infty ^2$, is denoted by $(C_{D,p,vis}, C_{L,p,vis})$. The fourth term $\boldsymbol {F}_{f}$ is the friction force

(5.8)\begin{equation} \boldsymbol{F}_{f} = \oint_{\partial B}{\rho\nu(\boldsymbol{n}\times\boldsymbol\omega) \,{\rm d} S}. \end{equation}

The force coefficient of the friction force, scaled by $0.5\rho U_\infty ^2$, is denoted by $(C_{D,\,f}, C_{L,\,f})$. The third and the fourth terms are small in the current Reynolds number range and the total force is dominated by the vortex-induced force $\boldsymbol {F}_Q$ and the added-mass force $\boldsymbol {F}_a$.

It should be noted that although the vortex-induced force introduced in (5.2) can be used to evaluate the impact of vortex structures on the force, it can not reveal how the vortex structures are formed, which is discussed in detail in § 3. Since (5.1) is derived from the exact Navier–Stokes equations, for the DNS result at $Re=800$, the force computed using (5.1) equals the wall-stress integral result with high accuracy. For the wall-resolved iLES results, however, a large error may be experienced in the volume-integral term $\boldsymbol {F}_Q$. Therefore, in the following analysis, the $\boldsymbol {F}_Q$ is only computed at $Re=800$.

5.2. Quantitative analysis of the lift force

In the four terms on the right-hand side of (5.1), $\boldsymbol {F}_a$ is determined by the instantaneous wall acceleration; the other three terms are determined by the instantaneous velocity field and its spatial gradients. As a result, vortex structures identified from the instantaneous velocity field can only reflect the origin of the force component $\boldsymbol {F} - \boldsymbol {F}_a$, which is mainly due to the vortex-induced force $\boldsymbol {F}_Q$. For the lift components, $\boldsymbol {e}_y$, which is the unit vector of the $y$-axis, and $\phi _y$ are used in (5.1).

The force decomposition based on (5.1) is shown in figure 21. As shown in figure 21(a), the viscous components $C_{f}+C_{L,p,vis}$ are very small at $Re=800$, and they become even smaller at $Re=10\,000$. Therefore, the lift component $C_L - C_{L, a} = C_{L, Q} + C_{f}+C_{L,p,vis}$ is almost entirely due to the vortex-induced lift $C_{L, Q}$. Therefore, the $C_L-C_{L,a}$, which is shown in figure 21(b), can be used to represent the vortex-induced lift.

Figure 21. Time-dependent lift coefficients at $Re=800$ and $10\,000$: (a) viscous effect, $C_{L,\,f} + C_{L, p, vis}$; (b) $C_L - C_{L,a} = C_{L, Q} + C_{L,\,f} + C_{L, p, vis}$; (c) added-mass force, $C_{L,a}$.

The vortex-induced lift reaches the maximum value around $t/T=0.25$ when the airfoil has the maximum downward velocity. At $St_c=0.73$, a negative lift peak is observed around $t/T=0.75$ when the airfoil has the maximum upward velocity. At lower Strouhal numbers, this negative lift peak is, however, not found due to the missing LEV from the lower surface. In terms of the plunging Strouhal number, cases with a higher plunging frequency have a larger vortex-induced lift. Regarding the Reynolds number, cases at $Re=10\,000$ have larger vortex-induced lift than cases at $Re=800$, except at $St_c=0.32$, where both cases have the same time-averaged lift. These phenomena are in agreement with the time-averaged lift shown in figure 5(a).

The added-mass force is calculated using (5.4) from potential-flow theory and, as a result, it does not depend on the Reynolds number, see figure 21(c). The added-mass force varies as a cosine function that achieves the maximum lift at $t/T=0.0$ and the minimum lift at $t/T=0.5$, which correspond to the highest displacement and the lowest displacement, respectively. Although the vortex-induced lift and the added-mass force have the same magnitude, they do not cancel each other out in the total lift because of the phase difference of ${\rm \pi} /2$. The maximum total lift $C_L$ occurs between the positive peak of the vortex-induced lift and that of the added-mass force.

The phase-averaged maximum and minimum values of $C_L-C_{L,a}$ and $C_{L,a}$ are shown in figure 22. The solid black lines show the added-mass force, which is symmetric about the $y=0$ line. The symbols show the maximum and minimum values of $C_L-C_{L,a}$, which is biased to the positive value. In flow regimes I, II and III, the maximum vortex-induced lift grows with the $St_c$ while the minimum $C_{L}-C_{L,a}$ is almost a constant. A drop of vortex-induced lift is observed around $St_c=0.5$ when the flow switches from regime III to regime IV. In regimes IV and V, the maximum vortex-induced lift keeps increasing in a similar trend with the added-mass force. The minimum vortex-induced lift starts decreasing with the $St_c$ in regime IV and it decreases faster in regime V. Cases at $Re=10\,000$ have a slightly higher vortex-induced lift in terms of the absolute value.

Figure 22. The phase-averaged maximum and minimum lift coefficients at $Re=800$. The solid black line shows the added-mass force $C_{L,a}$. Symbols show $C_L-C_{L,a}$.

To further identify the vortex structures that are responsible for the vortex-induced lift $C_{L,Q}$, the integration of the volume lift element is computed in the truncated fluid domain $-4c\leq x \leq x_w$, $-4c\leq y\leq 4c$, at $t/T=0.25$ and $t/T=0.75$ when the vortex-induced lift has the maximum or the minimum value. The added-mass force equals zero at these two phases. Dependence of this volume integration result, denoted by $C_{L, Q}(x_w)$, on the downstream location $x_w$ of the integration domain is shown in figures 23 and 24. For $x_w\leq -0.5c$, the vortex-induced lift $C_{L, Q}(x_w)$ is almost zero. For $x_w\ge 2.5c$, the total lift $C_{L, Q}(x_w)+C_{L,p,vis}+C_{L,\,f}$ hardly changes and equals the wall-stress integral result, which indicates that vortex structures in the range of $x<-0.5c$ or $x>2.5c$ have a negligible contribution to the lift force.

Figure 23. Dependence of the lift decomposition results on the integration domain $-4c\leq x \leq x_w, -4c\leq y\leq 4c$ at $Re=800 \textrm{ and } t/T=0.25$: (a) $St_c=0.10$; (b) $St_c=0.32$; (c) $St_c=0.50^-$. Here $C_{L,a}=0$. The top row is $\omega _z$; the middle row is $2\rho \phi _y Q$; the bottom row is the force coefficients vs $x_w$.

Figure 24. Dependence of the lift decomposition results on the integration domain $-4c\leq x \leq x_w, -4c\leq y\leq 4c$ at $Re=800$: (a) $St_c=0.50^+, t/T=0.25$; (b) $St_c=0.73, t/T=0.25$; (c) $St_c=0.73, t/T=0.75$. Here $C_{L,a}=0$. The top row is $\overline {\omega _z}$; the middle row is $2\rho \phi _y \bar {Q}$; the bottom row is the force coefficients vs $x_w$.

The distribution of the spanwise vorticity $\omega _z$, scaled by $U_\infty c^{-1}$, and the volume lift element $2\rho \phi _y Q$, or $2\rho \phi _y \bar {Q}$ when scaled by $\rho U_\infty ^2 c^{-1}$, are also shown in figures 23 and 24, where their $x$-axes are aligned with the $x_w$-axis in the bottom $x_w-C_L$ plot. In the centre of a concentrated vortex where rotation is stronger than the strain rate, $Q$ is positive and it contributes a suction vortex-induced force. At the edges of the airfoil and the surroundings of the vortices where the strain rate is stronger, $Q$ is negative and it contributes a repulsive vortex-induced force. In the boundary layer and straight free shear layer, $Q$ is very small and has little contribution to the force.

We first consider three cases from regimes I to III in figure 23 where $t/T$ is 0.25 and the vortex-induced lift reaches the maximum value. All three cases have stable 2-D flows. The strong strain-rate field from the lower surface induces a positive lift when the fluid passes the sharp leading edge. This lift is, however, cancelled out by the strain field on the left side of the LEV from the upper surface because the $\phi _y$ has different signs on the upper and on the lower surfaces. The LEV suction is the main cause of the vortex-induced lift and it grows with the $St_c$. The strain-rate-dominated region on the right side of the LEV reduces the lift. At $St_c=0.10$ in regime I, the LEV-induced TEV is very close to the airfoil and, therefore, significantly contributes to the lift. At $St_c=0.32$ in regime II, the TEV is at a downstream location of $x=1.5c$ and almost has no contribution to the lift. At $St_c=0.50^-$ in regime III, the pLEV–TEV pair has a small positive contribution to the lift.

Two cases at $St_c=0.50^+$ and $St_c=0.73$ in regime IV are shown in figure 24(a,b), at $t/T=0.25$. Both cases are 3-D and therefore the spanwise-averaged vorticity and $Q$ are used. At $St_c=0.50^+$, $t/T=0.25$ in regime IV, the suction peak of the LEV is reduced by 50 % compared with the case at $St_c=0.50^-$ in regime III. On the left side of the TEV, the strong strain-rate field, which stretches the pLEV making it become 3-D, also contributes a strong negative lift. The TEV at the moment of $t/T=0.25$ is still very close to the airfoil and contributes a relatively strong lift. The net contribution of the pLEV–TEV pair is positive. The chaotic wake has little contribution to the lift. Since more lift is generated around the trailing edge, the centre of pressure also moves to the trailing edge by $0.05c$, compared with the case at $St_c=0.50^-$. A similar situation is observed at $St_c=0.73, t/T=0.25$. It can be concluded that both the LEV and the TEV are the main cause of the lift in regime IV.

The case at $St_c=0.73, t/T=0.75$ in regime IV is shown in figure 24(c) to explore the cause of the negative lift peak at $t/T=0.75$. At this moment, the airfoil is moving upward and a small LEV is formed on the lower surface. The negative volume lift element in the centre of this LEV has a magnitude of $O(500)$ while its surrounding positive value has an order of $O(50)$. As a result, this LEV from the lower surface contributes a strong negative lift via vortex suction. The LEV on the upper surface also has a strong positive lift contribution, which is however largely cancelled out by the strain-rate field on both sides of the upper LEV. The net contribution of the upper LEV is not strong at this time instant.