J. Fluid Mech. (2008), vol. 602, pp. 175–207.
c
2008 Cambridge University Press
doi:10.1017/S0022112008000864 Printed in the United Kingdom
175
Direct numerical simulations of forced
and unforced separation bubbles on an
airfoil at incidence
L. E. JONES, R. D. SANDBERG AND N. D. SANDHAM
Aerodynamics and Flight Mechanics Research Group, School of Engineering Sciences, University of
Southampton, Southampton, SO17 1BJ, UK
(Received 21 August 2007 and in revised form 23 January 2008)
Direct numerical simulations (DNS) of laminar separation bubbles on a NACA-0012
airfoil at Re
c
=5× 10
4
and incidence 5
◦
are presented. Initially volume forcing is
introduced in order to promote transition to turbulence. After obtaining sufficient
data from this forced case, the explicitly added disturbances are removed and the
simulation run further. With no forcing the turbulence is observed to self-sustain, with
increased turbulence intensity in the reattachment region. A comparison of the forced
and unforced cases shows that the forcing improves the aerodynamic performance
whilst requiring little energy input. Classical linear stability analysis is performed
upon the time-averaged flow field; however no absolute instability is observed that
could explain the presence of self-sustaining turbulence. Finally, a series of simplified
DNS are presented that illustrate a three-dimensional absolute instability of the two-
dimensional vortex shedding that occurs naturally. Three-dimensional perturbations
are amplified in the braid region of developing vortices, and subsequently convected
upstream by local regions of reverse flow, within which the upstream velocity
magnitude greatly exceeds that of the time-average. The perturbations are convected
into the braid region of the next developing vortex, where they are amplified further,
hence the cycle repeats with increasing amplitude. The fact that this transition process
is independent of upstream disturbances has implications for modelling separation
bubbles.
1. Introduction
Under an adverse pressure gradient a boundary layer may separate, leading to
reverse (upstream) fluid flow. Within the separated region disturbances are strongly
amplified, typically leading to transition to turbulence. The resultant turbulent flow
enhances mixing and momentum transfer in the wall-normal direction, and causes
the boundary layer to reattach. This system of laminar separation, transition and
turbulent reattachment is referred to as a laminar (or transitional) separation bubble
(LSB), and is typically associated with flows at low to moderate Reynolds numbers.
When present on an airfoil, laminar separation bubbles have a marked effect upon
aerodynamic performance. Drag forces are typically increased, and the presence of
a separation bubble may determine stall behaviour (Gault 1957). The phenomenon
of bubble bursting, where a small increase in incidence leads to a sudden increase in
bubble length, causes a dramatic loss in aerodynamic performance and hence is an
important consideration in low-Reynolds-number airfoil design.
176 L. E. Jones, R. D. Sandberg and N. D. Sandham
Using the results of Gaster (1967), Horton (1969) was the first to describe the
time-averaged structure of a laminar separation bubble, and proposed an empirical
model for predicting bubble behaviour. Despite refinements such as the use of the
e
n
transition prediction method, modelling of low-Reynolds-number effects and the
dependence on background turbulence levels, present day models do not adequately
predict bubble bursting or unsteady behaviour.
More recently, advances in understanding of laminar separation bubbles have been
made via numerical methods. The first numerical simulations of separation bubbles
were limited either to two-dimensional analysis (Pauley, Moin & Reynolds 1990), or
only studied primary/linear instability and did not resolve transition (Pauley 1994;
Rist 1994). Being less computationally expensive, linear stability analysis could be
performed before fully resolved direct numerical simulations (DNS) were possible.
Hammond & Redekopp (1998) performed local analysis of separated boundary layer
profiles in order to determine whether absolute instability could be observed, as it
had been for separated shear layers (Huerre & Monkewitz 1985) and bluff-body
wakes (Hannemann & Oertel 1989). For certain profiles local absolute instability
was observed, depending on both the maximum reverse flow and the height of
the reverse flow region. Similar criteria were investigated and confirmed later by
Rist & Maucher (2002). Hammond & Redekopp (1998) found that for profiles at
Re
δ
∗
=10
3
, a minimum reverse flow velocity of 20 % was required to observe local
absolute instability. Theofilis (2000) performed global linear stability analysis of a
laminar separation bubble, specifying a spanwise wavenumber, β, and computing the
resultant two-dimensional disturbance eigenvectors. A temporally growing stationary
mode was observed, associated with unsteadiness at the reattachment point but not
affecting the separation point. Growth rates were found to be significantly lower than
those associated with amplification of convective instabilities within the shear layer;
however Theofilis (2003) suggests that the existence of this mode may potentially be
relevant to the phenomenon of vortex shedding from separation bubbles as observed
by Pauley et al. (1990).
The first numerical simulations to fully resolve transition to turbulence within a
laminar separation bubble were conducted by Alam & Sandham (2000) and Spalart
& Strelets (2000). Alam & Sandham (2000) performed DNS of a laminar separation
bubble on a flat surface, induced by upper boundary transpiration. Performing linear
stability analysis on analytic velocity profiles similar to those observed in the DNS,
Alam & Sandham found that reverse flow greater than 15 % would be required
in order to sustain absolute instability, compared to an observed reverse flow of
only 4–8 %. As a result, it was concluded that the transition process was driven
by convective instability. Spalart & Strelets (2000) conducted DNS of a laminar
separation bubble for the purpose of assessing turbulence models. No unsteadiness
was introduced and inflow disturbances were less than 0.1 % of the free-stream
velocity; however transition to turbulence was still observed. As a result the study
stated that entry-region disturbances (referring to Tollmien–Schlichting, or TS, type
waves) could be discarded as the mechanism behind transition; however the study
also stated that the magnitude of reverse flow present was unlikely to be sufficient
to sustain absolute instability. In the study of Spalart & Strelets three-dimensionality
was observed to occur rapidly, with no clear regions of primary or secondary
instability, whereas Alam & Sandham observed -vortex-induced breakdown. Hence
the first two fully resolved DNS of laminar separation bubbles apparently observed
different instability mechanisms leading to transition, and different transitional
behaviour.
Forced and unforced separation bubbles on an airfoil at incidence 177
Marxen et al. (2003) performed a combined DNS and experimental study of an LSB
on a flat plate. Periodic two-dimensional disturbances were introduced upstream of
separation, and three-dimensionality was introduced via a spanwise array of spacers.
The separated shear layer was observed to roll up to form vortices which subsequently
broke down to turbulence. The same configuration was studied further by Lang, Rist
& Wagner (2004) and again by Marxen, Rist & Wagner (2004) in order to quantify the
respective roles of two-dimensional and three-dimensional disturbances. Marxen et al.
concluded that transition was driven by convective amplification of a two-dimensional
TS wave, which also determined the length of the bubble, and that the dominant
mechanism behind transition is an absolute secondary instability in a manner first
proposed by Maucher, Rist & Wagner (1997), and investigated further by Maucher,
Rist & Wagner (1998).
It is clear that numerical simulations of separation bubbles can differ markedly in
behaviour. Unlike wake and shear layer flows, the presence of regions of absolute
instability within laminar separation bubbles has not yet been confirmed by linear
stability analysis of DNS or experimental velocity profiles, and stability characteristics
of separation bubbles not well defined in all cases. With continued advances in
computing power, it is now possible to perform DNS of laminar separation bubbles
on full airfoil configurations. This contrasts with previous numerical studies, which
have been limited to bubbles on flat plates or other simplified geometries in order to
reduce the computational cost. The advantage of studying full airfoil configurations
is that the bubble can interact more strongly with the potential flow, in particular via
the Kutta condition at the trailing edge. The bubble will be closer in nature to those
observed under flight conditions, and the influence of the bubble behaviour on the
aerodynamic performance of the airfoil can be observed directly.
The purpose of the current study is to investigate the dependence of bubble
behaviour on the presence of boundary layer disturbances, and to investigate the
role of instability mechanisms in separation bubble transition. First, data from both
forced and unforced three-dimensional simulations of a laminar separation bubble
on a NACA-0012 airfoil will be compared. Classical linear stability analysis will then
be performed upon the time-averaged flow fields obtained from the DNS, in order
to determine whether any local absolute instability is present. Finally, a series of
computationally inexpensive simulations will be presented, intended to explain the
self-sustained transition to turbulence observed in the first part of the study.
2. Direct numerical simulations
2.1. Simulation geometry
The chosen airfoil geometry is a NACA-0012, extended to include a sharp trailing
edge and rescaled to unit chord length. The coordinate system for the curvilinear
C-type grids used in all simulations is given in figure 1. Grids are equidistantly spaced
in the z-direction for three-dimensional simulations. The two parameters governing
domain size are the wake length W , and the domain radius R. The airfoil chord is
used as the reference length scale and the coordinate system is defined such that the
trailing edge is located at (x, y)=(1, 0).
2.2. Governing equations
All simulations were run at a Reynolds number based on airfoil chord of Re
c
=5×10
4
,
and Mach number M =0.4 unless otherwise stated. A compressible-flow formulation
was chosen so that sound waves originating at the trailing edge could also be studied
178 L. E. Jones, R. D. Sandberg and N. D. Sandham
R
R
W
ξ
η
Figure 1. Topology of the computational domain.
(Sandberg, Sandham & Joseph 2007). The DNS code directly solves the unsteady,
compressible Navier–Stokes equations, written in curvilinear form as
∂ Q
∂t
+
∂ E
∂ξ
+
∂ F
∂η
+
∂ G
∂z
=
∂ R
∂ξ
+
∂ S
∂η
+
∂ T
∂z
. (2.1)
The conservative vector Q, inviscid flux vectors E, F and G, and the viscous vector
terms R, S and T are defined as
Q =
⎛
⎜
⎜
⎜
⎜
⎝
ρ
ρu
ρv
ρw
E
t
⎞
⎟
⎟
⎟
⎟
⎠
, E =
⎛
⎜
⎜
⎜
⎜
⎜
⎝
ρU
ρuU + pξ
x
ρvU + pξ
y
ρwU
(E
t
+ p)U
⎞
⎟
⎟
⎟
⎟
⎟
⎠
, (2.2)
F =
⎛
⎜
⎜
⎜
⎜
⎝
ρV
ρuV + pη
x
ρvV + pη
y
ρwV
(E
t
+ p)V
⎞
⎟
⎟
⎟
⎟
⎠
, G =
⎛
⎜
⎜
⎜
⎜
⎝
ρw
ρuw
ρvw
ρww + p
(E
t
+ p)w
⎞
⎟
⎟
⎟
⎟
⎠
, (2.3)
R =
⎛
⎜
⎜
⎜
⎜
⎜
⎝
0
τ
xx
ξ
x
+ τ
xy
ξ
y
τ
yx
ξ
x
+ τ
yy
ξ
y
τ
zx
ξ
x
+ τ
zy
ξ
y
Q
x
ξ
x
+ Q
y
ξ
y
⎞
⎟
⎟
⎟
⎟
⎟
⎠
, S =
⎛
⎜
⎜
⎜
⎜
⎜
⎝
0
τ
xx
η
x
+ τ
xy
η
y
τ
yx
η
x
+ τ
yy
η
y
τ
zx
η
x
+ τ
zy
η
y
Q
x
η
x
+ Q
y
η
y
⎞
⎟
⎟
⎟
⎟
⎟
⎠
, T =
⎛
⎜
⎜
⎜
⎜
⎜
⎝
0
τ
xz
τ
yz
τ
zz
Q
z
⎞
⎟
⎟
⎟
⎟
⎟
⎠
, (2.4)
where ρ is the fluid density, u, v and w are velocity components in the Cartesian x, y
and z directions, p is the pressure, and E
t
is the total energy per unit volume, defined
as
E
t
= ρe +
1
2
ρ(uu + vv + ww), (2.5)
where
e =
T
γ (γ − 1)M
2
. (2.6)
Forced and unforced separation bubbles on an airfoil at incidence 179
Metric terms are defined as
ξ
x
=
y
η
J
,η
x
= −
y
ξ
J
,ξ
y
= −
x
η
J
,η
y
=
x
ξ
J
, (2.7)
noting that terms ξ
z
and η
z
are both equal to zero for computational grids with no
spanwise variation, as used in the current study, and the Jacobian J is defined as
J = x
ξ
y
η
− x
η
y
ξ
. (2.8)
The contravariant velocities U and V are defined as
U = ξ
x
u + ξ
y
v, V = η
y
v + η
x
u, (2.9)
and the stress terms τ
ij
as
τ
ij
=
µ
Re
∂u
i
∂x
j
+
∂u
j
∂x
i
−
2
3
µ
Re
∂u
k
∂x
k
δ
ij
. (2.10)
The terms Q
i
comprise the conduction and work terms of the energy equation,
Q
i
= −q
i
+ u
j
τ
ij
, (2.11)
where
q
i
=
−µ
(γ − 1)M
2
RePr
∂T
∂x
i
. (2.12)
Viscosity is calculated using Sutherland’s law,
µ = T
3/2
1+C
T + C
,C=0.368
˙
6, (2.13)
and finally, the perfect gas law relates p, ρ and T
p =
ρT
γM
2
. (2.14)
The primitive variables ρ,u,v and T have been non-dimensionalized by the free-
stream conditions and dimensionless parameters Re, Pr and M are defined using
free-stream (reference) flow properties. The ratio of specific heats is specified as
γ =1.4 and the Prandtl number as Pr =0.72.
2.3. Numerical method
Fourth-order-accurate central differences utilizing a five-point stencil are used for
spatial discretization. Fourth-order accuracy is extended to the domain boundaries by
use of a Carpenter boundary scheme (Carpenter, Nordstr
¨
om & Gottlieb 1999). No
artificial viscosity or filtering is used. Instead, stability is enhanced by entropy splitting
of the inviscid flux terms in combination with a Laplacian formulation of the viscous
terms (Sandham, Li & Yee 2002). The explicit fourth-order-accurate Runge–Kutta
scheme is used for time-stepping, and all cases were run with a constant time step of
t =1× 10
−4
. Appropriate boundary conditions must be applied to avoid unphysical
reflections. At the free-stream (η
+
) boundary, where the only disturbances likely to
reach the boundary will be in the form of linear waves, an integral characteristic
boundary condition is applied (Sandhu & Sandham 1994). At the downstream exit
boundary (ξ
±
), which will be subject to the passage of nonlinear fluid structures,
a zonal characteristic boundary condition (Sandberg & Sandham 2006) is applied
for increased effectiveness. At the airfoil surface an adiabatic, no-slip condition is
applied.
