Three-dimensional transition and force characteristics of low-Reynolds-number flows past a plunging airfoil

Abstract

The three-dimensional (3-D) transition of the leading-edge vortex (LEV) and the force characteristics of the plunging airfoil are investigated in the chord-based Strouhal number $St_c$ range of 0.10 to 1.0 by means of experimental measurements, numerical simulations and linear stability analysis in order to understand the spanwise instabilities and the effects on the force. We find that the interaction pattern of the LEV, the LEV from a previous cycle (pLEV) and the trailing-edge vortex (TEV) is the primary mechanism that affects the 3-D transition and associated force characteristics. For $St_c \leq 0.16$, the 3-D transition is dominated by the LEV–TEV interaction. For $0.16 < St_c \leq 0.44$, the TEV lies in the middle of the LEV and the pLEV and therefore vortex interaction between them is relatively weak; as a result, the LEV remains two-dimensional up to a relatively high Reynolds number of $Re = 4000$ at $St_c = 0.32$. For $0.44 < St_{c} \leq 0.54$, and at relatively low Reynolds numbers, the pLEV and the TEV tend to form a clockwise vortex pair, which is beneficial for the high lift and stability of the LEV. For $0.49 \leq St_c$, the pLEV and TEV tend to form an anticlockwise vortex pair, which is detrimental to the lift and flow stability. In the last $St_c$ range, vortex interaction involving the LEV, the TEV and the pLEV results in an unstable period-doubling mode which has a wavelength of about two chord-lengths and the 3-D transition enhances the lift.

1. Introduction

A geometrically two-dimensional (2-D) airfoil undergoing plunging motion, or experiencing an unsteady incoming flow, is a classical problem in unsteady aerodynamics. It is encountered not only in biological locomotion (Shyy et al. 2010) and micro-aerial vehicle control (Gursul, Cleaver & Wang 2014; Greenblatt & Williams 2022), but also with traditional large-aspect-ratio wings in a gust (Jones, Cetiner & Smith 2022) as well as helicopter rotor blades (McCroskey 1982). These flows are characterised by flow separation at the leading and the trailing edges. The concentrated vortices, which are formed by the rolling up of separated boundary layers and modulated by the motion of the airfoil, play crucial roles in generating unsteady force and momentum. A thorough understanding of the evolution characteristics of these vortices is important to uncover the underlying mechanisms producing unsteady forces and thereby help improve the aerodynamic performance of a wing inboard of the tip region.

The vortex structures in the wake have been observed to be closely correlated with the unsteady force (Zhang 2017). For example, a symmetric airfoil undergoing plunging can generate wake vortices, with increasing flapping frequency, that are in the form of a Bérnad–von Kármán vortex street, a reverse Bérnad–von Kármán vortex street, a deflected vortex street or aperiodic patterns (Lewin & Haj-Hariri 2003). The drag-to-thrust transition is found to be correlated with the transition of Bérnad–von Kármán to reverse Bérnad–von Kármán vortex street, although the thrust transition takes place at a slightly higher frequency (Godoy-Diana, Aider & Wesfreid 2008; Andersen et al. 2017). Analogously, the deflected vortex wake indicates a non-zero mean lift (Cleaver, Wang & Gursul 2012; Deng et al. 2016). However, the wake pattern itself does not always provide a reliable indicator of the force. Indeed, Floryan, Van Buren & Smits (2020) found that in some situations the change in the wake patterns results in no change of the force, especially when the wake contains self-interacting vortex patterns.

Heuristically, the effect of a vortex structure decreases when it moves away from the body. Chang (1992) expressed the force exerted on a body by the fluid as the dot product of a potential flow and the Lamb vector. The potential flow, which is first introduced by Quartapelle & Napolitano (1983), represents the importance of a fluid element at a specified spatial location. The velocity of this auxiliary potential flow decays as $O(r^{-2})$ in 2-D space or $O(r^{-3})$ in three-dimensional (3-D) space when the distance to the body $r$ goes to infinity. Lee et al. (2012) applied Chang's method to the impulsively started finite plate. They found that flow structures within three chord lengths of the plate contribute 99 % lift. Using the finite-domain vorticity moment theory, Li & Lu (2012) found the vorticity moment of body-connecting vortical structures make the largest contribution to the force. Gao et al. (2019) subsequently studied the pressure Poisson equation of the incompressible flow and observed that the source term of the pressure Poisson equation, within three chord lengths, contributes almost all of the pressure force of a 2-D airfoil. It can therefore be concluded that the near-body flow structures are at the heart of force generation. Therefore, in the following analysis, we focus our study on resolving and analysing the near-body flows.

The near-body vortex system around the plunging airfoil includes leading-edge vortex (LEV), trailing-edge vortex (TEV) and the interaction of these vortices with the viscous wall that can lead to additional secondary vortices (Lewin & Haj-Hariri 2003; Wu, Ma & Zhou 2006). The primary LEV is found to be the key to high lift (Ellington et al. 1996; Eldredge & Jones 2019) and thrust (Wang 2000). The LEV usually moves downstream after its formation. Under high-frequency and low-amplitude plunging motion, however, the LEV can move upstream and circumscribe the leading edge to the opposite side of the airfoil (Lewin & Haj-Hariri 2003; Gao et al. 2019). Further, the interaction between the LEV and the TEV can also greatly affect the thrust and propulsion efficiency (Lewin & Haj-Hariri 2003; Choi, Colonius & Williams 2015).

A plunging airfoil at a large angle of attack (AoA) generates strong LEVs and contains rich vortex evolution phenomena and therefore is adopted in this study. The plunging motion also enhances the lift of the post-stall airfoil and therefore provides promising control methods (Gursul et al. 2014). The salient dimensionless parameters that describe this flow are the AoA, the Reynolds number $Re$, the peak-to-peak amplitude $A_c = A/c$ and the chord-length-based Strouhal number $St_c =c / (U_\infty T)$. We note that the amplitude-based Strouhal number $St_A = A / (U_\infty T)$ is also often adopted instead of $St_c$.

Experimental studies of this problem include water tunnel studies by Gursul & Ho (1992) who investigated the aerodynamic force of a static NACA 0012 airfoil submerged in an unsteady incoming flow at a Reynolds number of $50\,000$ and AoA of $20^\circ$. They found the time-averaged lift is greatly enhanced around an optimal frequency of $St_c=0.25$, where the lift can be higher than three times the steady lift. The optimal frequency is independent of the oscillating amplitude of the incoming flow but the peak lift increases with this amplitude. This high lift was attributed to a vortex pair which formed around the trailing edge. Cleaver et al. (2011) studied the lift enhancement by a plunging airfoil with a peak-to-peak amplitude of $A_c\in [0.05, 0.4]$ and a Strouhal number of $St_c\in [0, 3]$. The Reynolds number was $10\,000$ and the AoA was $15^\circ$, while the spanwise length of their airfoil was three times the chord length. Their frequency-lift curves show a double-peak structure where the first peak arose around $St_c=0.5$ for all amplitudes. The second peak has a higher lift and occurs in the amplitude-based Strouhal number interval $St_A\in [0.4, 0.5]$. The drop in lift after the first peak was caused by an adverse interaction between the LEV from a previous period and the TEV. Using the same airfoil geometry and under the same Reynolds number, Cleaver et al. (2012) observed a bifurcation of the frequency-lift curve, from a lift-increasing branch to a lift-decreasing branch, for an AoA of ${\le }10^\circ$. For an ${\rm AoA} \geq 12.5^\circ$, bifurcation was not observed. Son et al. (2022) conducted combined volumetric particle image velocimetry (PIV) measurements and numerical simulations of the vortical structures around a plunging airfoil at $Re=10\,000$, where the plunging amplitude was $A_c=0.5$, the AoA was $15^\circ$ and the spanwise length $b$ was five chord lengths. At Strouhal numbers $St_c = 0.32, 0.64, 0.95$, they found 3-D vortex structures behind the trailing edge of the airfoil and argued that the three-dimensionality was caused by a Crow-type instability of the LEV–TEV or LEV–wall interaction vortex pair.

Numerical simulations are widely used in the study of the vortex-dominated flows generated by the plunging airfoil. Many of these numerical simulations have considered 2-D flow because of their relatively low computational cost. For the plunging airfoil at a large AoA, 2-D simulations have been performed at a much lower Reynolds number, which assumes the flow is Reynolds number independent. Andro & Jacquin (2009) studied the plunging NACA 0012 airfoil at $Re = 1000$, $St_A=0.27$ and an AoA of $15^\circ$. The $St_c$ varies from 0 to 2 by adjusting the plunging period and the amplitude simultaneously. They found the first lift peak occurs at $St_c\approx 0.4$. However, due to the decreasing amplitude with $St_c$, the second lift peak is not apparent in their results. Choi et al. (2015) studied the horizontal surging and vertical plunging of a NACA 0012 airfoil at $Re \le 500$. For the plunging wing, the lift reaches a peak value around $St_c=0.5$ and then drops almost vertically to less than one-third of the maximum value. This dramatic drop in the lift was not observed in experiments where the lift loss is only about 10 %. Choi et al. (2015) found the sudden loss of lift was due to the formation of an adverse vortex pair whose induced velocity was directed upstream. The second lift peak is very weak or not present. Gao, Sherwin & Cantwell (2022) conducted 2-D simulations at $Re=400$, using the same airfoil geometry as Choi et al. (2015). They found the lift discontinuity around $St_c=0.5$ is hysteretic in the range of $0.49\leq St_c\leq 0.5$. Although 2-D simulation provides some key features of the vortex patterns, they exhibit many disagreements with experimental results, especially for $St_c$ around 0.5 and above.

Biglobal linear stability analysis is a useful tool to assess the validity of 2-D simulations and for obtaining the critical Reynolds number and the dominating mode in the 3-D transition. Meneghini et al. (2011), He et al. (2017) and Deng, Sun & Shao (2017) studied the linear stability of the 2-D periodic flow around a static airfoil at a post-stall AoA. They found two unstable modes whose spanwise wavelengths are 0.2 and 2 chord lengths, respectively. The short-wavelength mode first becomes unstable. At an AoA of $15^\circ$, the critical Reynolds number of the NACA 0015 airfoil is $Re=730$ (Deng et al. 2017). At an AoA of $20^\circ$, the critical Reynolds number was $Re=561$ (He et al. 2017). Moriche, Flores & García-Villalba (2016) studied the linear stability of the 2-D periodic flow around a flapping airfoil at $Re=1000$. They found a long-wavelength period-doubling unstable mode, which has a wavelength of $4.078c$ and resembles mode A of the cylinder flow. Similarly, Sun, Deng & Shao (2018) studied the linear stability of a plunging NACA 0015 airfoil at zero AoA and $Re=1700$. They identified a long-wavelength and a short-wavelength unstable mode whose wavelengths were around 1 and 0.2 chord lengths, respectively. Gao, Sherwin & Cantwell (2020) and Gao et al. (2022) studied the linear stability of a plunging NACA 0012 airfoil at AoA of $15^\circ$, $Re=400$, $A_c=0.5$, $St_c \in (0, 1)$. For $St_c \geq 0.5$ where the lift is low, they found a period-doubling unstable mode whose wavelength was about two chord lengths. For $St_c \leq 0.5$ where a high lift is achieved, no linearly unstable mode was obtained. The linear stability analysis partially answered the 2-D to 3-D transition mechanisms, but it was still insufficient to uncover the transition mechanisms at lower frequencies and the features of the 3-D flows.

Fully 3-D numerical simulations clearly make the fewest assumptions and therefore potentially provide the most salient information about the flow characteristics, although they also come at a notable computational expense. Visbal (2009) studied the flow past a plunging SD7003 airfoil using high-fidelity implicit large-eddy simulations (iLES) at the chord-based Strouhal number of $St_c=1.25$, $A_c=0.1$ and ${\rm AoA} = 4^\circ$. The spanwise length of the airfoil is $0.2c$. At $Re=40\,000$, an abrupt breakdown of the LEV due to the spanwise instability is observed, however, the force is almost unaffected by the flow transition. Visbal (2011) studied the 3-D transition of the plunging SD7003 airfoil with respect to the Reynolds number from $Re=1000$ to $Re=120\,000$, at $St_c=0.080$, $A_c=1.0$ and ${\rm AoA} = 8^\circ$. The spanwise length of the airfoil is $0.4c$. The 3-D transition of the flow is found at $Re=5000$. A delay in the shedding of the LEV and the TEV is observed with the increasing Reynolds number, which results in a phase shift in the force coefficients. Zurman-Nasution, Ganapathisubramani & Weymouth (2020) performed 3-D iLES of a NACA 0016 airfoil at $St_c = 1.875$, $A_c\in [0, 0.64]$, $Re=5300$ and ${\rm AoA} = 10^\circ$ with a spanwise length of $6c$. They observed the 2-D flow was valid only around the peak-to-peak amplitude-defined Strouhal number of $St_A=0.3$. Away from this Strouhal number, three-dimensionality can drastically change the wake patterns and aerodynamic forces. Son et al. (2022) performed iLES of a NACA 0012 airfoil of length $5c$ at $A_c=0.5$, $Re=10\,000$ and ${\rm AoA} = 15^\circ$, at Strouhal numbers of $St_c= 0.32$, 0.64 and 0.95. All three cases show 3-D vortex patterns. Gao et al. (2022) performed 3-D direct numerical simulations (DNS) of a NACA 0012 airfoil of spanwise length $10c$ at $A_c=0.5, Re=400$, ${\rm AoA} = 15^\circ$ and $St_c\in (0, 1)$. The flow was found to be 3-D for the low-lift regime of $St_c\ge 0.49$, but 2-D for the high-lift regime of $St_c\leq 0.5$. The dominating wavelengths in 3-D results, which are around two chord lengths, agree with the biglobal instability results. In addition to the periodically plunging motion, Moriche et al. (2021a) also studied the 3-D transition of a NACA 0012 airfoil undergoing a plunging-type manoeuvre, i.e. the plunging duration of the airfoil is one period. They found the 3-D effect is notable only when the plunging velocity is large enough so that the LEV interacts with the small vortex pair formed below it.

The transition to 3-D flow has a notable effect on the aerodynamic performance of the plunging airfoil at post-stall AoAs. Clarifying the transition boundary and the aerodynamic force characteristics is important for both flow physics and flow controls such as lift enhancement or alleviation. However, the transition map in the Strouhal number range $St_c\in (0,1)$ is still rather unexplored and a significant discrepancy exists around $St_c=0.50$, between the numerical simulations at Reynolds numbers of order $10^2$ and the experimental results at $Re=10\,000$. In this paper, we take a complementary approach to study the 3-D transition and the force characteristics of the flow around the plunging airfoil at a post-stall AoA of $15^\circ$. The main objectives are to identify the spanwise instabilities, their origin, evolution with Reynolds number and their effect on the aerodynamic force in a chord-based Reynolds number from $Re=400$ to $Re=10\,000$. We consider data from experimental measurements, 2-D and 3-D simulations, to identify flow regimes over the Strouhal number range of $St_c\in [0.10, 1.0)$ and explore the 3-D flow features in each regime. Subsequently, we explore the instability mechanisms using biglobal stability analysis and 3-D numerical simulations. Finally, we discuss the vortex origin of the lift force using a novel Galilean invariant force decomposition theory. Initial findings of our study were presented in Gao et al. (2020, 2022).

The paper is organised as follows. The experimental set-up, the numerical approach and the validations are given in § 2. An overview of the flow transition features, including the time-averaged force, regimes classification, the 3-D vortex patterns and the spanwise wavelength, are presented in § 3.1. In § 3.2, the flow evolution with Reynolds number in each of the first four regimes is explored. In § 4, we examine the biglobal linear instability for the fourth regime with 2-D periodic base flows and the relevant instability mechanisms. Then, in § 5, we discuss the vortex origin of the lift force. Finally, in § 6, we summarise the main conclusions.

2. Methodology

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. 2022). 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. 2020).

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 1991). 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. 2015; Moxey et al. 2020). 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. (2022), 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. (2022) 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 1986). 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 2014). The cross-stream boundaries are periodic. A spectral vanishing viscosity (Moura et al. 2020) 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. (2022), 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 2008). 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. (2022), 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. (2022).

3. Results

3.1. Overview of the flow transition features

3.1.1. Time-averaged force

The time-averaged lift and drag coefficients, non-dimensionalised by $0.5 \rho U_\infty ^2 c b$, are shown in figure 5. The parametric space is $\alpha _0 = 15^\circ$, $0\leq St_c < 1$, $400\leq Re\leq 10\,000$ and $A_c=0.5$. In the experimental measurements, the spanwise length of the airfoil is $5c$ and end plates are attached on both ends. In the numerical simulations, the spanwise length of the airfoil is also $5c$, whereas a periodic boundary condition is applied on both ends.

Figure 5. Time-averaged (a) lift and (b) drag coefficients of the plunging airfoil at $\text {AoA} = 15^\circ$. Roman numerals show the regimes classification. Black vertical lines indicate the starting $St_c$ of each regime.

At $Re=10\,000$, the 3-D numerical simulation captures the varying trend of the experimental results, as indicated by the red plus symbols and the solid red line in figure 5. Two lift peaks, previously observed by Cleaver et al. (2011), are captured. However, unlike the results of Cleaver et al. (2011), where the second peak has a much higher lift for an airfoil with a spanwise length of $3c$ and plunging amplitude of $A_c\leq 0.4$, the two peak values in the current study are almost the same. The first peak occurs at $St_c=0.48$ and is followed by a rapid drop in the lift. For $St_c<0.48$, the lift force continuously grows with increasing plunging frequency/Strouhal number. The second peak is flatter in nature and maintains a high lift value in the $St_c$ range of 0.73 to 0.83.

The 3-D lift at a lower Reynolds number of $Re=800$ is indicated by the black solid circles in figure 5(a), which also has two lift peaks. The first peak occurs at $St_c=0.53^-$ and achieves a 16 % higher lift than the $Re=10\,000$ cases, however, the lift coefficient for $St_c>0.50^+$ is slightly lower than the $Re=10\,000$ cases, resulting in a sharper drop of lift around $St_c=0.5$. The flow around $St_c=0.5$ is found to be hysteretic since increasing $St_c$ from a lower value or decreasing $St_c$ from a higher value allows the flow to maintain the high or low lift value of the previous state, meaning that both states can be stably achieved at $St_c=0.5$. To distinguish cases in this hysteric region, superscripts ‘$-$’ and ‘$+$’ are used to denote cases from the left (high lift) and right (low lift), respectively. Hysteresis is not observed at $Re=10\,000$.

Two-dimensional simulations were also performed at $Re = 400$ and $Re=800$ and the force is shown by open symbols in figure 5. In this case, there is only one lift peak, which is around $St_c=0.5$. Similar to the 3-D results at $Re=800$, the lift drop around $St_c=0.5$ is hysteretic. The hysteresis process happens in a Strouhal number interval of $St_c\in [0.48, 0.50]$ for $Re=400$, and in the interval $St_c\in [0.49, 0.54]$ for $Re=800$. At $Re=800$ and for $0.10\leq St_c\leq 0.54^-$, the lift force of the 2-D simulation is identical to that of the 3-D simulation, which indicates the flow is 2-D within this parametric range.

The validity of 2-D simulations can be assessed by the comparison of the time-averaged lift. At low plunging frequency, $St_c \leq 0.16$, the 2-D lift coefficient at $Re=400$ and $Re=800$ is slightly smaller than the 3-D lift coefficient at $Re=10\,000$. However, 2-D simulations at $Re=10\,000$, as indicated by empty circles in figure 5(a), significantly over-predict the lift. In the Strouhal number range of $St_c\in (0.16, 0.44]$, the lift coefficients of the 2-D simulations and the 3-D simulations at both $Re=800$ and $Re=10\,000$ are almost the same. This phenomenon seems similar to the finding of Zurman-Nasution et al. (2020) who reported that 2-D flow can exist in a medium Strouhal number range of $St_A\approx 0.3$ for a NACA 0016 airfoil at an AoA of $10^\circ$ and $Re=5300$. Nevertheless, as we shall demonstrate, the flow at $Re=10\,000$ is 3-D. In the Strouhal number range of $St_c\in (0.44, 0.54^-]$, the transition characteristic from high lift to low lift and the transition Strouhal number both depend on the Reynolds number. As already mentioned, hysteresis exists in this Strouhal number range for $Re=800$ but not for the case when $Re = 10\,000$. For $St_c> 0.49^+$, the 2-D lift and 3-D lift are different. The 3-D lifts at $Re=800$ and $Re=10\,000$ have the same increasing or decreasing trends, but their magnitudes differ.

The time-averaged drag coefficients are shown in figure 5(b). Two different behaviours exist. For $St_c<0.5^-$ which corresponds to an increasing lift regime, the drag coefficient at $Re=10\,000$ is at first almost independent of the plunging frequency for $St_c\le 0.22$, before slightly increasing. This varying trend is in agreement with the results of Cleaver et al. (2011). Drag coefficients at $Re= 400$ and $Re=800$ keep increasing with $St_c$. For $St_c>0.50^+$, the drag decreases and a drag-to-thrust transition occurs around $St_c=0.76$. Hysteresis also exists on the drag curve around $St_c=0.50$ at $Re=400$ and $Re=800$.

Based on the varying trends of the 3-D time-averaged lift at $Re=800$ and $Re=10\,000$, five regimes can be identified in the frequency range of $St_c\in [0.10, 1.00)$. They are regime I $[0.10, 0.16]$, regime II $(0.16, 0.44]$, regime III $(0.44, 0.54^-]$, regime IV $[0.49^+, 0.80]$ and regime V $(0.80, 1)$. The regimes are labelled by Roman numerals in figure 5(a), with the starting $St_c$ indicated by black vertical lines. The force varies monotonically in each regime, except for $St_c$ around 0.50, which corresponds to regime III and the start of regime IV, where a decrease of lift occurs. Regime I is separated because the 2-D lift differs from the 3-D lift at $Re=10\,000$ in this Strouhal number range. The classification criterion for regime III and regime IV is the near-field vortex pattern which is discussed in the following and also in Appendix D. The upper boundary of regime III and the lower boundary of regime IV vary with the Reynolds number.

3.1.2. Vortex patterns

Due to the large AoA of the airfoil, the effective AoA during the upward motion is much smaller than during the downward motion, which results in a much weaker vortex shedding from the lower surface. Therefore, the LEV in the following discussions refers to the one shed from the upper surface if not otherwise specified.

The spanwise-averaged vortex patterns at $Re=800$ and $Re=10\,000$ are shown in figure 6, where phase averaging is also applied to reduce statistical errors of aperiodic flows. The time snapshot at $t/T=0.5$, when the airfoil just finishes the downward motion and is at the lowest displacement, is selected, because the LEV and the TEV are strong at this moment. Driven by the distance the LEV advects during the downward motion, four distinct vortex-shedding patterns can be identified in the Strouhal number range of $0.10 \leq St_c < 1$.

Figure 6. Typical near-field vortex patterns at the lowest displacement ($t/T=0.5$), coloured by the spanwise vorticity. Phase-averaged and spanwise-averaged 3-D simulation results are shown.

In regime I where $0.10\le St_c\leq 0.16$, the LEV advects past the trailing edge when the airfoil is at its lowest displacement ($t/T=0.5$). The separated shear layer that feeds the LEV is cut off by the TEV which is also formed during the downward half period. The remaining shear layer accumulates before the TEV and forms another concentrated vortex, V2.

In regime II where $0.16 < St \le 0.44$, the LEV has not reached the trailing edge when the airfoil is at the lowest displacement ($t/T=0.5$) but the LEV from the previous period (pLEV) has passed the trailing edge when the airfoil is at the highest displacement. When the primary TEV forms at $t/T=0.5$, the separation between the primary TEV and current LEV as well as the previous LEV is so large that there is little vortex interaction among them. A secondary vortex is generated between the wall and the LEV because of the LEV–wall interaction.

In regime III where $0.44< St\leq 0.54^-$, the LEV has not reached the trailing edge when the airfoil is at the lowest displacement ($t/T=0.5$) and the pLEV has only just passed the trailing edge when the airfoil is at the highest displacement ($t/T=0$). When the primary TEV forms at $t/T=0.5$, the pLEV is still very close to the TEV and pairs with it. The fully formed TEV has a weaker circulation than the pLEV. As a result, the induced motion of the pLEV–TEV pair moves the vortex pair upwards and away from the airfoil, resulting in a separate vortex system above the airfoil. This vortex pattern is observed at $Re=800$ but is not found at $Re=10\,000$, which indicates that this regime only exists under a certain Reynolds number.

In regime IV where $0.49^+\leq St_c\leq 0.80$, the LEV has also not reached the trailing edge when the airfoil is at the lowest displacement ($t/T=0.5$) and the pLEV has not passed the trailing edge either when the airfoil is at the highest displacement. The TEV forms in the presence of the pLEV and again pairs with the pLEV. The main difference of this vortex pair to the previous case is that the TEV has a stronger circulation than the pLEV, which induces a velocity against the incoming flow and a downward motion on the vortex pair. This vortex pattern exists at both $Re=800$ and $Re=10\,000$. This different behaviour of the pLEV–TEV pair, which is caused by the relative strength of the pLEV and the TEV, is the classification criterion for regime III and IV.

In regime V where $0.80< St_c<1$, the advection length of the LEV is even shorter. The vortex-shedding pattern in this regime is similar to that in regime IV, but the LEV is stronger due to the larger plunging velocity. This regime has a very chaotic vorticity field and the input power is high whereas the lift force starts decreasing, therefore, this regime is not discussed in detail in this work.

The 2-D vortex pattern bears a similar feature near the airfoil in regimes I, II and III, but due to the lack of a 3-D effect, the concentrated vortices are preserved for a longer distance in the wake. In regimes IV and V, the 2-D vortex patterns are drastically different from the 3-D results. Since the 2-D vortex patterns have been studied extensively in previous work such as Lewin & Haj-Hariri (2003), Choi et al. (2015) and Deng, Sun & Shao (2015), the relevant results are placed in Appendix C.

The three-dimensionality of the flow and the streamwise vorticity are also analysed. The 3-D vortical structures, identified by the $Q$-criterion and coloured by the streamwise vorticity, are shown in figures 7 and 8, where a dimensionless value of $Q^*=Qc^2U_\infty ^{-2}=10$ is chosen.

Figure 7. Instantaneous 3-D vortical structures ($Q^*=10$) coloured by the streamwise vorticity when the airfoil is at the lowest displacement: (a) $Re=800, St_c=0.10$; (b) $Re=800, St_c=0.32$; (c) $Re=10\,000, St_c=0.10$; (d) $Re=10\,000, St_c=0.32$.

Figure 8. Instantaneous 3-D vortical structures ($Q^*=10$) coloured by the streamwise vorticity when the airfoil is at the lowest displacement: (a,e) $Re=800, St_c=0.50^+$; (b,f) $Re=800, St_c=0.73$; (c,g) $Re=10\,000, St_c=0.50^+$; (d,h) $Re=10\,000, St_c=0.73$. Two successive periods are shown with panels (a–d) taken at $t/T=0.5$ and panels (e–h) taken at $t/T=1.5$.

In figure 7, the instantaneous 3-D vortical structures at $St_c=0.10$ and $St_c=0.32$ are shown, which correspond to regime I and regime II, respectively. At $Re=800$ (figures 7(a) and 7(b)), both Strouhal numbers do not show a 3-D effect. At $Re=10\,000$ and $St_c=0.10$ (figure 7c), although the LEV can be observed in the spanwise-averaged vorticity field, see figure 6, it breaks down into small-scale vortices and is almost invisible by the isosurface of $Q^*=10$. The vortex V2 also breaks down into small-scale vortices. Meanwhile, the TEV remains almost straight. At $Re=10\,000$ and $St_c=0.32$ (figure 7d), the pLEV and LEV are both almost straight. The TEV, which is relatively weak, has a larger displacement. The streamwise vortices that formed from the LEV–wall interaction and are shown by the blue and red spots on the LEV, reveal a short wave pattern.

The 3-D vortical structures at $St_c=0.50^+$ and 0.73 are shown in figure 8. Two successive periods are presented to show the period-doubling feature of the flow. At $Re=800$ and $St_c=0.50^+$, a spanwise periodic 3-D vortex pattern is observed. Its spanwise wavelength is $2.5c$ and its period is twice the plunging period. The long-wavelength 3-D displacement pattern of the pLEV–TEV pair resembles the Crow instability of an unequal-strength counter-rotating vortex pair. Further downstream, the crest of the pLEV from a previous period forms an upward-oriented hairpin vortex, which is caused by the stretching of the stronger previous TEV. At $Re=10\,000$ and $St_c=0.50^+$, the spanwise wavelength is $2.5c$ and the period-doubling feature can also be observed from the pLEV. However, the flow here contains many irregular small-scale vortices. The current LEV–wall interaction also induces a short-wavelength 3-D pattern. At $St_c=0.73$, the flow field has much richer vortical structures, nevertheless, the long-wavelength spanwise pattern and period-doubling feature can still be observed from the displacement of the current LEV. In the middle of the span, the LEV bends to the right at $t/T=0.5$, see figure 8(d), whereas it bends to the left at $t/T=1.5$, see figure 8(h).

3.1.3. Spanwise wavelength

At $Re=800$ and $St_c\le 0.73$, the spanwise-varying pattern of the vortical structures is periodic or almost periodic and therefore the wavelength can be measured directly from the flow field. At higher frequencies or at $Re=10\,000$, however, the flow structures are too chaotic to identify the spanwise wavelength. To evaluate the characteristic spanwise wavelength in these circumstances, the energy distribution of the spanwise velocity $w$ in the spanwise wavenumber $\beta$ space is calculated, see figure 9(a). For a 2-D flow, the spanwise velocity $w$ is zero. Therefore, the energy of $w$ can be used to represent the 3-D energy of the flow and its peak energy corresponds to the characteristic spanwise wavelength of the flow. Since the LEV and TEV are the main concern in this work, a near-body fluid domain $-c\le x\le 2.5c$, $-3c\le y\le 3c$ and one phase point at $t/T=0.5$ is selected for the spatial integration. Phase averaging is applied on the dimensionless wavenumber–energy ($\beta c - E(w^*)$) curve to reduce the statistical error. It should be noted that because the spanwise length of the airfoil is $5c$, only wavelengths of $5c/n$, such as $5c$, $2.5c$, $1.7c$, $1.25c,\ldots\,$, can be captured.

Figure 9. (a) Energy distribution of the spanwise velocity vs the spanwise wavenumber at $Re=10\,000, t/T=0.5$, the number indicates the $St_c$ of each line. (b) The most energetic wavelength $\lambda /c$. The background colour indicates the regime classification.

With increasing $St_c$, the 3-D energy of the flow decreases in regime I and increases in regime IV. The 3-D energy is low in regime II and it reaches the lowest level around $St_c=0.32$ where the pLEV–TEV and LEV–TEV interactions are weak. In regime V, more 3-D energy is distributed to small-scale structures, indicating a more chaotic flow field.

In regime I, the 3-D energy does not have a clear peak in the wavenumber space and the most energetic wavelength, which decreases with the plunging frequency, is $2.5c$ at $St_c=0.10$ and $1.25c$ at $St_c=0.16$. Inside regime II, the 3-D energy distribution has a flat peak. The most energetic wavelength, which is $0.83c$ for $0.32\leq St_c\leq 0.38$ at $Re=10\,000$, is much shorter. A second dominating wavelength of $2.5c$ also exists, and its energy is lower than the most energetic wavelength by only 4 % at $St_c=0.32$. In regimes IV and V, the flow is dominated by a long-wavelength mode which has a sharp peak and a wavelength of $2.5c$ or $1.7c$. With increasing $St_c$, the wavelength in regimes IV and V is also likely to become shorter. A summary of the most energetic wavelength $\lambda =2{\rm \pi} \beta ^{-1}$ is shown in figure 9(b), where background colours show the regimes classification. Since regime III does not exist at $Re = 10\,000$ and the LEV in regime III is still 2-D at $Re=800$, no data is presented for regime III in figure 9(b) and therefore its range is marked by the dashed lines.

3.2. Flow evolution with Reynolds number

To further clarify the 3-D transition features of the flow in the first four regimes, one Strouhal number from each regime is selected and the flow evolution with the Reynolds number is analysed. It should be noted that, once the flow becomes 3-D, the aerodynamic force and the flow field from 2-D simulations are obviously no longer physically realisable.

3.2.1. Regime I, $St_c = 0.10$

At a Strouhal number of $St_c=0.10$ in regime I, the effective AoA of the airfoil is always positive and it varies from $6^\circ$ to $24^\circ$. The 2-D to 3-D transition of the LEV occurs between $Re = 1000$ and $Re=2000$, as shown by figure 10. An increase in lift is observed after the flow transitions to a 3-D state, however, both before and after the flow transition, the lift force is almost independent of the Reynolds number. The variation of drag is non-monotonic. Before the 3-D flow transition, the drag decreases slowly with increasing Reynolds number. In the 2-D-to-3-D transition process, an increase in drag is observed. To examine the 3-D effect on the force, the lift and drag coefficients of the 2-D simulations are also indicated. The 2-D simulation predicts higher lift and drag than the 3-D simulation. The lift force of the 2-D simulation also keeps increasing instead of reaching a constant value. At $Re = 2000$, the relative errors of the lift and drag in 2-D simulations are $16\,\%$ and $10\,\%$, respectively. At $Re = 10\,000$, the relative errors of the lift and drag increase to $47\,\%$ and $32\,\%$, respectively.

Figure 10. Flow evolution with Reynolds number at $St_c=0.10$ in regime I. (a) Time-averaged lift and drag coefficients. The black vertical line shows the onset Reynolds number of the 3-D flow. (b) 3-D vortical structures ($Q^*=10$) coloured by $\omega _x^*$ at $Re=1000$. (c) 3-D vortical structures coloured by $\omega _x^*$ ($Q^*=10$) at $Re=2000$. The spanwise vorticity on the plane $z/c=0$ is also shown.

The 3-D vortical structures identified by the $Q$-criterion $(Q^*=Qc^2U_\infty ^{-2}=10)$ and coloured by the streamwise vorticity are shown in figure 10(b,c). At $Re=1000$, the flow is 2-D and only straight vortex columns are observed. At $Re=2000$, the LEV, the TEV and the concentrated vortex upstream of the TEV, all become 3-D. The energy distribution of the spanwise velocity on the spanwise wavenumber space is computed in the near-body domain $-c\le x\le 2.5c$, $-3c\le y\le 3c$. At $Re=2000$, results show that the 3-D energy distributes across wavelengths from $0.71c$ to $2.5c$, without a strong dominating wavelength being identifiable. At a higher Reynolds number of $Re=6000$ and $Re=10\,000$, the flow becomes more chaotic and the dominating wavelength is $1.7c$.

To uncover the flow physics behind the variation in force, the dimensionless spanwise-averaged vorticity fields ($\omega _z^*$) at $St_c=0.10$ are shown in figure 11, where both 2-D and the spanwise-averaged 3-D simulation results are provided. We first consider the 2-D simulations. At $Re=2000$, the LEV rolls up earlier than the LEV at $Re=1000$, because of the stronger fluid inertia against the viscosity. The change of LEV formation time at $Re=2000$ results in a stronger LEV being shed and higher lift. At higher Reynolds numbers of $Re=6000$ and $10\,000$, the free shear layer rolls up into smaller vortices above the airfoil, which enhances the lift in a similar mechanism as the LEV. In 3-D simulations and at $Re=2000$, the concentrated vortex before the TEV breaks down and annihilates with the TEV, resulting in a weaker TEV and a lower lift than the 2-D results. The free shear layer above the airfoil does not roll up into a concentrated vortex but instead disperses in a broad region. The 3-D simulation results at higher Reynolds numbers of $Re = 6000$ and $Re=10\,000$ do not change significantly compared with $Re=2000$, which explains why the time-averaged lift force is almost constant for $Re\ge 2000$ in 3-D simulations.

Figure 11. Comparison of the phase-averaged and spanwise-averaged $\omega _z^*$ field in 2-D and 3-D simulations: (a,e) $Re=1000$; (b,f) $Re=2000$; (c,g) $Re=6000$; (d,h) $Re=10\,000$. The airfoil is at the lowest displacement. Here $St_c=0.10$. Panels (a–d) are 2-D simulations; panels (e–h) are 3-D simulations. The ‘2-D’ label in panels (b–d) highlights these are 2-D simulations of 3-D flows and therefore are not physically realisable.

3.2.2. Regime II, $St_c = 0.32$

At a Strouhal number of $St_c=0.32$ in regime II, the effective AoA of the airfoil varies from $-12^\circ$ to $42^\circ$. The 2-D-to-3-D transition of the LEV at $t/T=0.5$ occurs between $Re = 4000$ and $Re=5000$ with dominating wavelengths of $1.7c$ and $0.21c$, as shown in figure 12. Also shown in this figure are the 3-D vortical structures identified by the $Q$-criterion and coloured by the streamwise vorticity, where $Q^*=10$ is once again applied. The smaller LEV (sLEV) indicated by the black arrow in figure 12(c) has the most significant 3-D displacement whereas the primary LEV is almost straight. After the flow transitions to 3-D, the time-averaged force, which is mainly influenced by the primary LEV, is still very close to that of the 2-D simulation, with a relative difference smaller than $10\,\%$. Although the LEV becomes 3-D from $Re=5000$, the 3-D transition of the TEV occurs between $Re=3000$ and $Re=4000$. Flow transition in the far wake is observed at $Re=2000$.

Figure 12. Flow evolution with Reynolds number at $St_c=0.32$ in regime II. (a) Time-averaged lift and drag coefficients. The black vertical line shows the onset Reynolds number of the 3-D LEV. The larger empty symbols have a refined mesh. (b) 3-D vortical structures ($Q^*=10$) coloured by $\omega _x^*$ at $Re=4000$. (c) 3-D vortical structures coloured by $\omega _x^*$ ($Q^*=10$) at $Re=5000$. The spanwise vorticity on the plane $z/c=0$ is also shown.

A sudden increase of lift is observed at $Re=4000$ and $Re=8000$. Further simulation on a refined grid shows this change of lift is grid-independent. This change of lift is related to a change in the self-interaction of the vortex system around the primary LEV. At $Re=4000$ and $t/T=0.5$, the sLEV whose vortex centre is located at $(0.12c, 0.18c)$, as indicated by the black arrow in figure 12(b), is still in close proximity to the aerofoil. At $Re=5000$ and $t/T=0.5$, however, this sLEV which is now centred at $(0.23c, 0.30c)$, has already lifted up through the induced motion of the primary LEV and it even merges with the primary LEV at $Re=10\,000$. The advection speeds of the primary LEV are also different. At $t/T=0.75$, the vortex centre of the LEV is $(1.02c, 0.01c)$ at $Re=4000$ and $(1.22c, 0.04c)$ at $Re=5000$. The interaction among the sLEV, the secondary vortex and the primary LEV at $Re=4000$ and $Re=8000$ delays the advection of the primary LEV causing it to remain for a longer duration above the airfoil and therefore produce an increased time-averaged lift.

3.2.3. Regime III, $St_c = 0.50^-$

At a Strouhal number of $St_c=0.5$, the effective AoA of the airfoil varies from $-23^\circ$ to $53^\circ$. Simulations were conducted from $Re=800$ to $Re=10\, 000$ using the flow field at $St_c=0.50^-, Re=800$ as an initial condition. The time-averaged force is shown in figure 13, where the 3-D vortical structures identified by the $Q$-criterion and coloured by the streamwise vorticity are also shown. The flow is 2-D at $Re=2000$. The lift and drag coefficients decrease with increasing $Re$, for $Re\le 6000$, and the lift and drag coefficients of the 3-D simulation remain consistent with those of the 2-D simulation. For $Re>6000$, 3-D simulations have a notably lower lift and are clearly different from the 2-D simulation. The vortical structures also show strong three-dimensionality.

Figure 13. Flow evolution with Reynolds number at $St_c=0.50$. The flow field at $Re=800, St_c=0.50^-$ is used as the initial condition of the simulation. (a) Time-averaged lift and drag coefficients. The black vertical line shows the onset Reynolds number of a 3-D LEV. (b) 3-D vortical structures ($Q^*=10$) coloured by $\omega _x^*$ at $Re=2000$. (c) 3-D vortical structures coloured by $\omega _x^*$ ($Q^*=10$) at $Re=7000$. The spanwise vorticity on the plane $z/c=0$ is also shown.

The 3-D vortical structures identified by the $Q$-criterion and coloured by the spanwise vorticity are shown in figure 14 plotted for an isocontour of $Q^*=10$. The corresponding spanwise-averaged vorticity field is also shown. The 2-D-to-3-D transition of the pLEV–TEV pair occurs between $Re = 2000$ and $3000$. The pLEV has the most significant 3-D displacement with dominant wavelengths of $1.3c$ and $0.36c$. At $Re=5000$, strong 3-D flow is observed in the pLEV–TEV pair, especially in the pLEV, but the spanwise-averaged vortex pattern still bears a similar characteristic to the $St_c=0.50^-, Re=800$ case. At $Re=6000$, however, the self-induced motion of the pLEV–TEV becomes anticlockwise, which is a feature similar to the $St_c=0.50^+, Re=800$ case. The dominating wavelength also becomes $2.5c$. The flow at $Re=6000$ is a marginal case, where the LEV is still 2-D at $t/T=0.5$. For Reynolds numbers $Re\ge 7000$, the anticlockwise moving pLEV–TEV vortex pair is observed. The strong blockage effect caused by this vortex pair weakens the LEV and results in a lower lift and a 3-D LEV. Therefore, regime III only exists below a critical Reynolds number, which is $Re\leq 5000$ at $St_c=0.50^-$.

Figure 14. (a–d) Instantaneous vortical structures ($Q^*=10$) coloured by $\omega _z^*$ at different Reynolds numbers and (e–h) the corresponding phase-averaged and spanwise-averaged $\omega _z^*$ fields: (a,e) $Re=3000$; (b,f) $Re=5000$; (c,g) $Re=6000$; (d,h) $Re=7000$.

3.2.4. Regime IV, $St_c =0.73$

At a Strouhal number of $St_c=0.73$, the effective AoA of the airfoil varies from $-34^\circ$ to $64^\circ$. The time-averaged force for Reynolds numbers from $Re=400$ to $Re=10\,000$ is shown in figure 15(a). The lift coefficient increases, and the drag coefficient decreases, gradually with the Reynolds number. For $Re\geq 5000$, the lift coefficient becomes almost constant. In the 2-D simulation, strong oscillations exists in the time-averaged force curve. The 2-D forces do not agree with the 3-D results either.

Figure 15. Flow evolution with Reynolds number at $St_c=0.73$ in regime IV. (a) Time-averaged lift and drag coefficients. (b) 3-D vortical structures ($Q^*=10$) coloured by $\omega _x^*$ at $Re=800$. (c) 3-D vortical structures coloured by $\omega _x^*$ ($Q^*=10$) at $Re=2000$. The spanwise vorticity on the plane $z/c=0$ is also shown.

The 3-D vortical structures identified by the $Q$-criterion and coloured by the streamwise vorticity are shown in figure 15(b,c) for a value of $Q^*=10$. Even at a very low Reynolds number of $Re=800$, the flow shows strong three-dimensionality. The flow is very chaotic at higher Reynolds numbers. At $Re=800$, $Re=2000$ and $Re=10\,000$, the dominating wavelength is $2.5c$. At $Re=6000$, the dominating wavelength is $1.7c$. The flow transition mechanism and the vortex dynamics explanation of the force are presented in the next two sections.

4. Flow transition mechanisms

In regime IV, the three-dimensionality of the flow is dramatically different from the previous regimes, with a much lower transition Reynolds number. The 3-D effect also has a significant impact on the vortex evolution and the aerodynamic force. In flying or swimming, the optimal Strouhal number, which is between 0.25 to 0.35 in the amplitude-based Strouhal number $St_A$ (Triantafyllou, Triantafyllou & Grosenbaugh 1993), lies in regime IV which has an amplitude-based Strouhal number of $0.25\le St_A\le 0.40$ ($0.49\le St_c\le 0.80$). Therefore, the instability mechanism in regime IV deserves further investigation.

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. (2017) 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. 2022) 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. Vortex origin of the lift force

5.1. Decomposition of the aerodynamic force

The force element theory (Chang 1992; Lee et al. 2012), 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 (1983) and Chang (1992) 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 (2015), where they found the thrust is mainly generated by the growing LEV, and Moriche et al. (2021b), where they found the main LEV contributes most to the force during the manoeuvre of the airfoil. Martín-Alcántara & Fernandez-Feria (2019) 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 2019). 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. 2019) 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 (1983) 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. (2019), Gao & Wu (2019) and Menon & Mittal (2021a,b). 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 (1992). Since $Q$ is widely used in the identification of vortical structures, this term was named vortex-induced force by Menon & Mittal (2021a). 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.

6. Conclusions

This work conducts a parametric study of the plunging airfoil using experimental measurements, spanwise homogeneous DNS and iLES. The airfoil has a profile of a NACA 0012 airfoil at a post-stall AoA of $15^\circ$, a spanwise length of 5 chord lengths and plunges at a peak-to-peak amplitude of 0.5 chord length. We focus on the lift characteristics, 3-D transition properties and the relevant flow mechanisms in the Reynolds number range of $Re=400$ to $Re=10\,000$ and chord-based Strouhal number range of $St_c=0.10$ to $St_c=1.0$. The results have a DNS resolution for $Re\leq 800$ and an iLES resolution for $Re\ge 1000$. The numerical results are carefully validated and good agreement of the lift force between the numerical simulations and the experimental measurement is achieved.

Our primary finding is that the interaction pattern of the LEV, LEV from the previous cycle (pLEV) and the TEV, which is determined by the advection length of the LEV in one-half period or one period, is the main factor that affects the 3-D transition of the LEV and the force features. Based on the vortex interaction pattern, five regimes of different lift and transition properties are identified as $St_c$ is varied.

1. (1) Regime I, $0.10 \leq St_c \leq 0.16$. The lift increases with the $St_c$ and the LEV has been advected past the trailing edge when the airfoil is at the lowest displacement. The transition of the LEV is dominated by the LEV–TEV interaction. For $St_c=0.10$, the 3-D transition of the flow occurs at $Re=2000$ and the dominant wavelength is $1.7c$ at $Re=10\,000$. Once the LEV becomes 3-D, the lift first increases with, and then becomes independent of, the Reynolds number. The 2-D simulation over-predicts the lift and drag due to the strong concentrated vortices formed by the free shear layer on the upper surface.

2. (2) Regime II, $0.16 < St_c \leq 0.44$. The lift increases with the $St_c$ and the TEV is in the middle of the LEV and the pLEV. Vortex interactions between the TEV and the LEV as well as the pLEV are very weak due to the large distance between two nearby vortices, and therefore the transition of the LEV is dominated by the LEV–wall interaction. For $St_c=0.32$, the 3-D transition of the LEV occurs at $Re=5000$ and the dominant wavelength is $0.83c$ at $Re=10\,000$. Compared with other regimes, the time-averaged force of a 3-D simulation, even after the flow transits to 3-D, is still very close to that of the 2-D simulation, with a relative difference smaller than $10\,\%$. With varying Reynolds numbers, the LEV and the secondary vortex it induces may form self-organised vortex patterns, which results in oscillations of the time-averaged force with respect to the $Re$.

3. (3) Regime III corresponds to $0.44 < St_c \leq 0.54^-$ and only exists below a critical Reynolds number. A very high lift, which is $25\,\%$ higher than the lift peak of regime II, can be achieved in this regime. The pLEV pairs with the TEV, forming a clockwise moving vortex pair, which moves upward and leaves an isolated LEV on the upper surface of the airfoil. Although the pLEV–TEV pair can be 3-D, the LEV is always 2-D and the 2-D simulation has approximately the same time-averaged force as the 3-D simulation. Above a critical Reynolds number of $Re=5000$ for $St_c=0.50^-$, the pLEV–TEV pair breaks down and an anticlockwise pLEV-LEV pair is formed instead, which is the feature of regime IV. The transition of the vortex pattern between regime III and regime IV is caused by the different relative strengths of the vortices in the pLEV–TEV pair and it is hysteretic. In order to distinguish the two regimes, superscripts ‘$-$’ and ‘$+$’ are used to denote regime III and regime IV, respectively.

4. (4) Regime IV, $0.49^+ \leq St_c \leq 0.80$. The lift first drops to a very low value and then gradually recovers to the peak value of regime II. At the beginning of this regime, the pLEV and the TEV form an anticlockwise moving vortex pair, which induces a velocity directing upstream that blocks the incoming flow and, as a result, reduces the strength of the LEV. For $St_c=0.50^+$, the flow becomes 3-D from $Re=400$ with a dominant wavelength of 2 chord lengths and the 3-D flow has a period equal to twice the plunging period. Biglobal linear stability analysis reveals that the Crow-type instability between the pLEV and the TEV is the main disturbance-amplification mechanism and that the interaction between the LEV and the pLEV-LEV pair is the reason for the doubled period in the disturbed flow. The wavelength, which is about two chord lengths, predicted by the biglobal linear stability agrees with 3-D simulation results. At higher $St_c$, the pLEV breaks down and the blockage effect of the pLEV–TEV pair is relieved. As a result, the time-averaged lift increases to the peak value of regime II.

5. (5) Regime V, $0.80 < St_c \leq 1.00$, the lift decreases with $St_c$ and part of the vorticity in the LEV circumvents the leading edge, moving to the lower surface during the upward half period. Due to the loss of lift and the highly chaotic characteristics of the flow, this regime is not discussed in depth in this work.

A quantitative analysis of the vortex origin of the lift force is also performed using a modified force-element theory, which is based on the weighted integral of the second invariant of the velocity gradient tensor and has a Galilean invariant volume force element. It is found that the time-dependent lift is dominated by the added-mass force and the vortex-induced lift. These two lift components have the same order of magnitude but the phase of the added-mass force is ahead of the vortex-induced lift by ${\rm \pi} /2$. Meanwhile, the time-averaged lift is mainly contributed by the vortex-induced force. In regimes I, II and III, the LEV is the main cause of lift. In regime IV, both the LEV and the TEV have significant contributions to the lift force.

This work clarifies the discrepancy between 2-D numerical studies at the Reynolds number of order $10^2$ and experimental studies at a Reynolds number of $Re=10\,000$. The vortex pattern transition between regime III and regime IV is clarified and the instability mechanism in regime IV is fully uncovered. This work provides physical insights into the validity of 2-D simulations and the Reynolds number dependence of the flow past a plunging airfoil.