12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

THREE-DIMENSIONAL TRANSITION OF A LOW REYNOLDS NUMBER
FLOW PAST A PLUNGING NACA 0012 AIRFOIL AT POST-STALL ANGLEOF ATTACK

An-Kang Gao
Department of Aeronautics
Imperial College London
London SW7 2AZ, UK
ankang.gao@imperial.ac.uk

Spencer J. Sherwin
Department of Aeronautics
Imperial College London
London SW7 2AZ, UK
s.sherwin@imperial.ac.uk

Chris D. Cantwell
Department of Aeronautics
Imperial College London
London SW7 2AZ, UK
c.cantwell@imperial.ac.uk

ABSTRACT

tex (Wang, 2000) and that the critical frequency of the drag-tothrust transition follows that of the von Kármán, to reverse vonKarman, vortex street transition (Godoy-Diana et al., 2008).
The lift enhancement of a post-stall airfoil by the plunging motion was studied by Cleaver et al. (2011) both experimentallyand numerically, however, the explanation of this lift enhancement through vortex dynamics has not been fully explored andso this is one of the focus of this paper.

The two-dimensional to three-dimensional transition of aflow past a plunging NACA 0012 airfoil at a Reynolds number of Re = 400, based on the chord length c, and an angleof attack of 15 degrees was investigated using global linear
stability analysis and spanwise-homogeneous direct numericalsimulation (DNS). The peak-to-peak plunging amplitude was
fixed at A/c = 0.5 and the Strouhal number was varied from
Stc = 0.10 to Stc = 1.00. This parameter regime encompasses
flow phenomena of leading-edge vortex (LEV) dominated flow
(0.10 ≤ Stc ≤ 0.19), almost vanishing LEV-trailing-edge vortex (TEV) interaction (0.22 ≤ Stc < 0.5), strong previous cycle LEV-TEV interaction (0.49 ≤ Stc ≤ 0.95) and aperiodic
flow (Stc ≥ 0.99). For the periodic baseflow, Floquet stability analysis was conducted. Below a Strouhal number of 0.5,
the Floquet multiplier is smaller than the static airfoil whichindicates the plunging motion stabilises the two-dimensional
baseflow. For higher frequencies, a period-doubling mode appears, which has a peak Floquet multiplier around a spanwise
wavelength of 2c. This unstable mode also dominates in threedimensional direct numerical simulations (DNS). Finally, ashort-wave mode becomes unstable at Stc = 0.95, which generates more small-scale vorticies in the DNS result.

An obvious starting point for numerical simulations is apurely two-dimensional (2-D) flow. Visbal (2009) undertook
large-eddy simulations of a plunging low-aspect-ratio SD7003
airfoil at Reynolds numbers in the range 10 000 to 40 000 andshowed that the three-dimensional (3-D) force is almost thesame as the 2-D force at a reduced frequency of Stc = 1.25
and a normalised amplitude of A/c = 0.1, where c is thechord length. Zurman-Nasution et al. (2020) conducted direct numerical simulations of a flapping NACA 0016 airfoil ata Reynolds number of 5300 and found the 2-D force only to
be valid for a range of Strouhal numbers around Stc = 0.3 forpure plunging motions.
The transition boundary and mechanism of 2-D to 3-Dflow are clearly crucial for assessing the validity of 2-D simulations. Moriche et al. (2016) studied the global instabilityof a combined plunging and pitching NACA 0012 airfoil ata Reynolds number of 1000 and an amplitude of 2c. They
found a long-wave subharmonic unstable mode, which has apeak Floquet multiplier at a wavelength of 4.078c and resembles mode A of the cylinder flow. Sun et al. (2018) studied
the instability of a plunging NACA 0015 airfoil at zero angleof attack, a Reynolds number of 1700, amplitudes from 0 to
0.3c and reduced frequencies from 4 to 8.4. They found threeunstable modes: mode A, a quasi-periodic mode, and mode B,
but the quasi-periodic mode is not dominant in the transition
region.

Introduction
For insects flying at low Reynolds numbers, the flapping
wing is a common form of motion. By flapping its wings, theinsect can achieve a much higher lift than the classical static
wing (Ellington et al., 1996) and a greater maneuverability(Alexander, 1986). This motivates the study of this mechanism and to develop smaller bio-inspired aerial vehicles.
As the simplest form of flapping, the pure plunging motion has been studied extensively. Thrust generation of aplunging airfoil at zero angle of attack (AoA) was studied by
Triantafyllou et al. (1991); Wang (2000); Godoy-Diana et al.
(2008) and Andersen et al. (2017). It was found that this gainin thrust is caused by the suction force of the leading-edge vor-

This study aims to explore the 2-D to 3-D transition mechanism of a very low Reynolds number (Re = 400) flow past a
post-stall NACA 0012 airfoil. The plunging amplitude is fixed
at 0.5c which is between the amplitude used by Moriche et al.

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Here, µ is the Floquet multiplier and |µ| > 1 indicates an unstable mode. β = 2π/λ is the wavenumber and λ is the wavelength. A modified Arnoldi algorithm with an iteration tolerance of 10−5 is used to find the most unstable eigenvalue andeigenmode (Barkley et al., 2008). The periodic baseflow is
captured by 128 snapshots over one time period and interpolated using 4th-order Lagrange interpolation.

(2016) and Sun et al. (2018). The plunging frequency is increased gradually from below the natural shedding frequency
to chaotic flow, to cover the complete range of flow regimes.
Direct numerical simulations are also conducted. The link between linear stability analysis and direct numerical simulations
is further discussed.

Numerical methods
The airfoil is treated as a rigid body with a fixed angle ofattack (AoA) of 15 degrees. To avoid the need for a moving
grid, the absolute velocity u is solved in the body frame ofreference. The governing equations in dimensionless form are∇ · u = 0,

Results
The time-averaged lift and drag coefficients at Re = 400,
A/c = 0.5, are shown in Figure 1. For comparison, experimental force coefficients of Cleaver et al. (2011) at Re = 10000,
A/c = 0.4, and implicit large-eddy simulation results of Son
et al. (2022) at Re = 10000, A/c = 0.5, are also presented.
For Stc < 0.5, both the lift and the drag forces increase
with plunging frequency. Around Stc = 0.32 (k = 1), the lowReynolds number lift agrees with high-Reynolds number results and the lift is insensitive to the plunging amplitude. The
low-Reynolds number drag is larger than the high-Reynolds
number results due to stronger friction.
A jump in the forces exists around Stc = 0.5 with a hysteresis from Stc = 0.49 to 0.5. The lift on the left side is 1.2
times higher than the lift on the right side of the jump. On theleft branch, the 3-D flow is the same as the 2-D flow and therefore no 3-D data is given there. On the right branch, the 3-Dflow deviates from the 2-D flow. This indicates the 2-D flow
becomes unstable in the low-lift right hand branch.
For Stc > 0.5, the lift from the numerical simulations issmaller than the experimental result of Cleaver et al. (2011).
However, the 3-D lift is higher than the 2-D lift. The dragdecreases with plunging frequency and transitions to thrust
around Stc = 0.76 (k = 2.4). The 3-D low-Reynolds number
drag agrees with the high-Reynolds number result.

(1)
−1 2

∂t u + (u − UB ) · ∇u = −∇p + Re

∇ u.

(2)

In this study, Re = 400 is based on the uniform incoming flow
U∞ = 1 and the chord length c = 1. The airfoil moves transversely under the motion
UB = −πA f sin(2π f t)ey ,(3)

where f is the plunging frequency and A = 0.5c is the peakto-peak amplitude.
Two dimensionless frequencies can be defined. These are
the reduced frequency, which is usually used to measure theflow unsteadiness, given byk = π f c/U∞ ,(4)and the Strouhal number,
Stc = k/π = f c/U∞ .

(5)

2-D baseflow
The wake pattern of the 2-D baseflow is shown in Figure 2. Based on the streamwise advection length of the LEV
and the vortex interactions, the Strouhal number range can beclassified into four regions.
Region 1 corresponds to 0.10 ≤ Stc ≤ 0.19 (0.3 ≤ k ≤0.6). In this Stc range, the streamwise advection length of theLEV during a half period is longer than the chord length c. One
strong trailing-edge vortex (TEV) sheds when the LEV passes
the trailing edge and several small TEVs also shed between
two successive LEVs.
Region 2 corresponds to 0.22 ≤ Stc < 0.5 (0.7 ≤ k <1.57). In this range, the streamwise advection length of theLEV during a half period is shorter than c but the advection length during one period is longer than c. The interaction between the LEV and the TEV is very weak. For
0.22 ≤ Stc ≤ 0.38 (0.7 ≤ k ≤ 1.2), the TEV is positioned nearly
equidistant between the current LEV and the previous LEV(pLEV). Interaction between the LEV and the pLEV is also
weak. Two single vortices are shed in each period and no vortex pair is formed. Vortices in the wake align along a straightline. In this Stc range, the flow is insensitive to the Reynolds
number and the plunging amplitude, and the 2-D flow is stable, as can be seen by a comparison between the present vorticity fields at Re = 400 and the PIV results of Son et al.
(2022) at Re = 10000 in Figure 3. The stable 2-D flow aroundk = 1 is also confirmed by Zurman-Nasution et al. (2020). For
0.41 ≤ Stc < 0.5 (1.3 ≤ k < 1.57), the pLEV interacts with the
TEV and forms a vortex pair in the wake. In this vortex pair,

A second-order velocity correction scheme (Karniadakis
et al., 1991) is used to transform equations (1) and (2) into
a set of Helmholtz equations (one for each space dimension)and one pressure Poisson equation. The spectral/hp element
method, which is implemented in the open-source code, Nektar++ (Cantwell et al., 2015), is used to solve equations (1) and(2).
The computational domain spans −40 ≤ x ≤ 60, −40 ≤y ≤ 40 and contains 9985 quadrilateral and 64 triangular elements, with 5 polynomial modes in each dimension. The
height of the first layer elements on the body surface is enforced to be 0.005c. For the spanwise-homogeneous 3-D simulation, a Fourier expansion with de-aliasing is applied. The
spanwise length is 10c and 64 Fourier planes are used (32complex modes). A spectral vanishing viscosity (SVV) witha DG-mimic kernel (Moura et al., 2020) is also applied in thepolynomial expansion to improve numerical stability. Basedon the parameters in Son et al. (2022), the above resolution is
sufficient.
The linearised equation (2) with infinitesimal 3-D disturbance is also solved in order to explore the linear stability ofthe periodic baseflow. According to Floquet theory, an eigenmode can be assumed to have the following formu(x, y, z,t0 + nT ) = µ n û(x, y,t0 ) cos(β z),(6)v(x, y, z,t0 + nT ) = µ n v̂(x, y,t0 ) cos(β z),(7)nw(x, y, z,t0 + nT ) = µ ŵ(x, y,t0 ) sin(β z).

(8)212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Stc

k

λ /c

0.50

1.57

2 ± 0.06-1.0569

0.64

2.00

2 ± 0.07-1.6519

0.80

2.50

1.8 ± 0.3

-1.7205

0.95

3.00

1.8 ± 0.3

-2.0912

A spanwise probe line was placed at x = 0.16, y = 0.16 asshown by the white dot in figure 5. The period-doubling phenomenon is clearly shown by the time evolution of the spanwise velocity w on this probe line in Figure 6. The alternative change of the sign of w along the spanwise direction indicates the bending wave on the LEV. At Stc = 0.64 (k = 2) andStc = 0.80 (k = 2), the flow evolution pattern is regular and periodic. At Stc = 0.95 (k = 3), the long-wave period-doubling
pattern is still the most significant spanwise pattern but there
are also many aperiodic small scales.
To measure the energy distribution along the spanwise
wavenumber, the time-averaged power spectral density of w
on the probe line is shown in Figure 7. At Stc = 0.64 (k = 2)
and Stc = 0.80 (k = 2), the spanwise disturbance mainly distributed on the most energetic wavelength λ = 2c (β = π) and
its higher harmonics. This wavelength 2c agrees with the linear stability result. At Stc = 0.95 (k = 3), the most energetic
wavelength is λ = 1.7c but the power spectral density showsa continuous distribution across all wavenumbers. This is inagreement with linear stability results that the mode is unstable for a wide range of wavelengths at Stc = 0.95 (k = 3).

µTable 1. The most unstable wavelength and the corresponding Floquet multiplier at Re = 400, A/c = 0.5.

the pLEV is stronger than the TEV. Therefore, their induced
motion leads to this vortex pair moving upwards.
Region 3 corresponds to 0.49 ≤ Stc ≤ 0.95 (1.53 ≤ k ≤3.00). In this range, the 2-D flow is still periodic but the advection length of the LEV during one period is shorter than thechord length c. The TEV has enough time to grow as the pLEV
starts interacting with it. Therefore, the TEV is stronger in thevortex pair, forming a downward motion pLEV-TEV vortex
pair in the wake. The induced velocity of this vortex pair also
reduces the strength of the current LEV. As a result, the lift,
which is correlated with the LEV, is also reduced.
Region 4 corresponds to Stc ≥ 0.99 (k ≥ 3.1). In thisrange, the 2-D baseflow becomes aperiodic.

Conclusions
The two-dimensional (2-D) to three-dimensional (3-D)
transition of a plunging airfoil at a post-stall angle of attackof 15 degrees, a Reynolds number of 400 and a plunging amplitude of 0.5c was explored using linear stability analysis andspanwise-homogeneous direct numerical simulations.
For Stc < 0.5, the lift grows with plunging frequency andthe 2-D flow is stable. The lift drops to less than half ofthe highest value at around Stc = 0.5, with a hysteresis from
Stc = 0.49 to Stc = 0.5. For St ≥ 0.49, the 2-D flow is unstable.
A period-doubling mode with a wavelength of approximately
twice the chord length first becomes unstable. This mode isalso the dominating spanwise flow structure in the fully developed 3-D flows.

Linear stability result
The linear stability results are shown in Figure 4. To test
the convergence of the numerical result, a refined mesh with 7
polynomial modes and a baseflow with 256 snapshots are also
used and are shown as red circles. The relative difference iswithin 1%.
For Stc < 0.5 (k < 1.57), no absolute unstable mode is observed. For Stc ≥ 0.49 (k ≥ 1.53), a periodic-doubling mode(µ is a negative real number) first becomes unstable. The
most unstable wavelength and the corresponding Floquet multiplier are shown in table 1. For all four frequencies explored,
the most amplified mode has a wavelength of about twice
the chord length. At Stc = 0.5 (k = 1.57) and 0.64 (k = 2),
this long-wave mode has a narrow bandwidth. At Stc = 0.80
(k = 2.5), two short-wave stable modes with wavelength ofabout 0.7c and 0.3c appear. At Stc = 0.95 (k = 3), modes in awide wavelength range λ ≥ 0.25c are unstable.

Acknowledgements
The authors acknowledge the Engineering andPhysical Sciences Research Council (EPSRC) Grant
No.
EP/S029389/1 and UK Turbulence Consortium
EP/R029326/1. This work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).
The authors would also like to acknowledge HPC support
from Imperial College Research Computing Service (DOI:
10.14469/hpc/2232).

3-D result
Since the characteristic wavelength of the most amplified
mode is around 2c, the spanwise length of the 3-D simulationis set as 10c, which is 5 times the wavelength. The flow is firstsimulated for 20 time periods using 32 planes in the spanwise
direction. Then 64 planes are used to continue the simulation
for another 20 time periods.
Evolution of the spanwise vorticity at Stc = 0.64(k = 2),
Re = 400, A/c = 0.5, is shown in Figure 5. Both the 2-D flow
and the spanwise-averaged flow are shown. At t/T = 0.75
when the LEV has just formed, the LEV is very strong andthere is no significant displacement along the vortex column.
At t/T = 1, the LEV in the 2-D flow has an almost circularshape, but the LEV in the 3-D flow are stretched in a 45-degree
direction. The 3-D effect makes the spanwise averaged LEV
appear to diffuse faster (through 3D deformations) and reduces
the adverse effect of the pLEV on the lift. Therefore the 3-D
lift is higher than the 2-D lift.

REFERENCES
Alexander, D. E. 1986 Wind tunnel studies of turns by flying
dragonflies. Journal of Experimental Biology 59 (3), 81.
Andersen, A, Bohr, Tomas, Schnipper, Teis & Walther,Jens Honore 2017 Wake structure and thrust generation of a
flapping foil in two-dimensional flow. Journal of fluid mechanics 812.
Barkley, D., Blackburn, H. M. & Sherwin, S. J. 2008 Direct optimal growth analysis for timesteppers. International
Journal for Numerical Methods in Fluids 57 (9), 1435–
1458.
Cantwell, Chris D, Moxey, David, Comerford, Andrew, Bolis, Alessandro, Rocco, Gabriele, Mengaldo, Gianmarco,
De Grazia, Daniele, Yakovlev, Sergey, Lombard, J-E,
Ekelschot, Dirk et al. 2015 Nektar++: An open-source312th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
mimicking spectral vanishing viscosity for high reynolds
number flows. Journal of Computational Physics 406,
109112.
Son, Onur, Gao, A-K, Gursul, Ismet, Cantwell, Chris D, Wang,
Zhijin & Sherwin, Spencer J 2022 Leading-edge vortex dynamics on plunging airfoils and wings. Journal of Fluid Mechanics 940.
Sun, Liping, Deng, Jian & Shao, Xueming 2018 Threedimensional instabilities for the flow around a heaving foil.
Physical Review E 97 (1), 013110.
Triantafyllou, M. S., Triantafyllou, G. S. & Gopalkrishnan, R.
1991 Wake mechanics for thrust generation in oscillating
foils. Physics of Fluids A: Fluid Dynamics 3 (12), 2835–
2837.
Visbal, Miguel R 2009 High-fidelity simulation of transitionalflows past a plunging airfoil. AIAA journal 47 (11), 2685–
2697.
Wang, Z. J. 2000 Vortex shedding and frequency selection inflapping flight. Journal of Fluid Mechanics 410 (410), 323–
341.
Zurman-Nasution, AN, Ganapathisubramani, Bharathram &Weymouth, GD 2020 Influence of three-dimensionality on
propulsive flapping. Journal of Fluid Mechanics 886.

spectral/hp element framework. Computer physics communications 192, 205–219.
Cleaver, D. J., Wang, Z., Gursul, I. & Visbal, M. R. 2011 Liftenhancement by means of small-amplitude airfoil oscillations at low reynolds numbers. AIAA Journal 49 (9), 2018–
2033.
Ellington, Charles P, Van Den Berg, Coen, Willmott, Alexander P & Thomas, Adrian LR 1996 Leading-edge vortices in
insect flight. Nature 384 (6610), 626–630.
Godoy-Diana, Ramiro, Aider, Jean-Luc & Wesfreid, José Eduardo 2008 Transitions in the wake of a flapping foil. Physical Review E 77 (1), 016308.
Karniadakis, George Em, Israeli, Moshe & Orszag, Steven A
1991 High-order splitting methods for the incompressiblenavier-stokes equations. Journal of computational physics
97 (2), 414–443.
Moriche, M, Flores, O & Garcı́a-Villalba, M 2016 Threedimensional instabilities in the wake of a flapping wing at
low reynolds number. International Journal of Heat andFluid Flow 62, 44–55.
Moura, Rodrigo C, Aman, Mansoor, Peiró, Joaquim & Sherwin, Spencer J 2020 Spatial eigenanalysis of spectral/hp
continuous galerkin schemes and their stabilisation via dg-

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 1. Time-averaged lift and drag coefficient scaled by 0.5ρU∞2 c of the plunging airfoil. In present simulations, Re = 400. In Son
et al. (2022); Cleaver et al. (2011), Re = 10000.

Figure 2.

Wake patterns of the 2-D plunging airfoil at Re = 400, t/T = 0.5.

Figure 3. Upper panels: instantaneous vorticity field of the 2-D plunging airfoil at k = 1, Re = 400. Lower panels: phase-averaged
vorticity field of the plunging airfoil at k = 1, Re = 10000 from Son et al. (2022)

Figure 4.

Floquet multiplier of the 2-D plunging airfoil at Re = 400. Red circles are the refined results to show numerical convergence.

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 5. Upper panels: instantaneous vorticity fields of the 2-D plunging airfoil at k = 2, Re = 400. Lower panels: spanwise-averaged
vorticity field of the 3-D plunging airfoil at k = 2, Re = 400.

Figure 6.

Figure 7.

Evolution of the spanwise velocity on the probe line (x = 0.16, y = 0.16) for (a) k = 2, (b) k = 2.5, (c) k = 3.

Power spectrum of the spanwise velocity on the probe line (x = 0.16, y = 0.16). From low to high, k = 2, k = 2.5 and k = 3.

6