WIND ENERGY
Wind Energ.
2011; 14:799–819
Published online 4 February 2011 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/we.458
SPECIAL ISSUE PAPER
Review of computational fluid dynamics for wind
turbine wake aerodynamics
B. Sanderse
1,2
,S.P.vanderPijl
2
and B. Koren
2,3
1
Energy Research Centre of the Netherlands (ECN), Petten, The Netherlands
2
Centrum Wiskunde & Informatica (CWI), Amsterdam, The Netherlands
3
Mathematical Institute, Leiden University, Leiden, The Netherlands
ABSTRACT
This article reviews the state-of-the-art numerical calculation of wind turbine wake aerodynamics. Different computational
fluid dynamics techniques for modeling the rotor and the wake are discussed. Regarding rotor modeling, recent advances
in the generalized actuator approach and the direct model are discussed, as far as it attributes to the wake description. For
the wake, the focus is on the different turbulence models that are employed to study wake effects on downstream turbines.
Copyright © 2011 John Wiley & Sons, Ltd.
KEYWORDS
wind energy; wake aerodynamics; CFD; turbulence modeling; rotor modeling
Correspondence
B. Sanderse, Energy Research Centre of the Netherlands (ECN), Petten, The Netherlands.
E-mail: sanderse@ecn.nl
Received 30 November 2009; Revised 20 August 2010; Accepted 4 January 2011
1. INTRODUCTION
During the last decades, wind turbines have been installed in large wind farms. The grouping of turbines in farms introduces
two major issues: reduced power production, because of wake velocity deficits, and increased dynamic loads on the blades,
because of higher turbulence levels. Depending on the layout and the wind conditions of a wind farm, the power loss of a
downstream turbine can easily reach 40% in full-wake conditions. When averaged over different wind directions, losses of
approximately 8% are observed for onshore farms and 12% for offshore farms (see e.g. Barthelmie et al.
1,2
).
When studying power losses and blade loading, wind turbine wakes are typically divided into near and far wakes.
3
The near wake is the region from the turbine to approximately one or two rotor diameters downstream, where the turbine
geometry directly affects the flow, leading to the presence of distinct tip vortices. Tip and root vortices lead to sharp
gradients in the velocity and to peaks in the turbulence intensity. For very high tip speed ratios, the tip vortices form a
continuous vorticity sheet: a shear layer. The turbine extracts momentum and energy from the flow, causing a pressure
jump and consequently an axial pressure gradient, an expansion of the wake and a decrease of the axial velocity. In the
far wake, the actual rotor shape is only felt indirectly, by means of the reduced axial velocity and the increased turbulence
intensity. Turbulence is the dominating physical process in the far wake, and three sources can be identified: atmospheric
turbulence (from surface roughness and t hermal effects), mechanical turbulence (from the blades and the tower) and wake
turbulence (from tip vortex breakdown). Turbulence acts as an efficient mixer, leading to the recovery of the velocity deficit
and a decrease in the overall turbulence intensity. Far downstream, the velocity deficit becomes approximately Gaussian,
axisymmetric and self-similar. Wake meandering, the large-scale movement of the entire wake, might further reduce the
velocity deficit, although it can considerably increase fatigue and extreme loads on a downwind turbine. It is believed to be
driven by the large-scale turbulent structures in the atmosphere.
4–6
The distinction between near and far wakes is also apparent when classifying the existing numerical models for wind
turbine wake aerodynamics (see Table I). The first and simplest approach is an analytical method that exploits the self-
similar nature of the far wake to obtain expressions for the velocity deficit and the turbulence intensity. The second,
Blade Element Momentum (BEM) theory, uses a global momentum balance together with two-dimensional (2D) blade
Copyright © 2011 John Wiley & Sons, Ltd.
799
Review of CFD for wind turbine wake aerodynamics B. Sanderse, S. P. v. d. Pijl and B. Koren
Ta b l e I . Classification of models.
Method Blade model Wake model
Kinematic Thrust coefficient Self-similar solutions
BEM Actuator disk + blade element Quasi one-dimensional momentum theory
Vortex lattice, vortex particle Lifting line/surface + blade element Free/fixed vorticity sheet, particles
Panels Surface mesh Free/fixed vorticity sheet
Generalized actuator Actuator disk/line/surface Volume mesh, Euler/RANS/LES
a
Direct Volume mesh Volume mesh, Euler/RANS/LES
a
RANS, Reynolds-averaged Navier–Stokes; LES, large eddy simulation.
elements to calculate aerodynamic blade characteristics. The vortex-lattice and vortex-particle methods assume invis-
cid, incompressible flow and describe it with vorticity concentrated in sheets or particles. Panel methods similarly
describe an inviscid flow field, but the blade geometry is taken into account more accurately, and viscous effects
can be included with a boundary-layer code; the wake follows as in vortex-wake methods. These four methods
have been extensively discussed in previous reviews, such as Vermeer et al.,
3
Crespo et al.,
7
Snel
8,9
and Hansen
et al.
10
The last two methods, the generalized actuator disk method and the direct method, are relatively new and
are commonly called computational fluid dynamics (CFD) methods. In this s urvey, we will review these methods,
discuss their ability to predict wind turbine wakes and give an outlook into possible future developments.
The paper is organized as follows. First, we will consider the Navier–Stokes equations and discuss their use to predict
turbulent flows (Section 2). Section 3 then discusses rotor modeling, Section 4 discusses wake modeling and Section 5
deals with the question how to verify and validate wind turbine wake CFD codes.
2. GOVERNING EQUATIONS
2.1. The incompressible Navier–Stokes equations
It is reasonable to assume that the flow field in wind turbine wakes is incompressible, since the velocities upstream and
downstream of a turbine placed in the atmosphere are typically in the range of 5–25 m s
1
. Only when calculating the
aerodynamics at blade tips that compressibility effects may be important. Since in most calculations of wind turbine wakes
the rotor is not modeled directly (which will be discussed later), the incompressible Navier–Stokes equations are a suitable
model to describe the aerodynamics of wind turbine wakes:
ru D 0;
@u
@t
C .u r/u D
1
rp C r
2
u
(1)
supplemented with initial and boundary conditions, which will be discussed in Section 2.3. In the case of a non-neutral
atmosphere, the Boussinesq approximation is typically employed to account for buoyancy effects, and an extra equation
for the temperature has to be solved. The effect of the rotation of the Earth, given by the Coriolis term, is typically neglected
in many wake studies but can have an effect when computations involve large wind turbines and wind farms (see e.g.
Porté-Agel et al.
11
).
Although this set of equations provides a complete model for the description of turbulent flows, it is not easily solved.
The difficulty associated with turbulent flows is the presence of the non-linear convective term, which creates a wide range
of time and length scales.
12
For example, in the atmospheric boundary layer (ABL), the largest turbulent scales are of
the order of 1 km, whereas the smallest scales are of the order of 1 mm.
13
Inside the blade boundary layers, the scales
are even smaller. The range of scales depends on the Reynolds number (Re), the dimensionless parameter that indicates
the ratio of convective forces to viscous forces in the flow. Large values of the Reynolds number, encountered in blade
and wake calculations, lead to a large range of scales, making computer simulations extremely expensive. Resolving
all scales in the flow, so-called direct numerical simulation (DNS), is therefore not feasible. Turbulence models need to
be constructed, modeling the effect of the unresolved small scales based on the behavior of the large scales. However,
even with the cost reduction provided by a turbulence model, one cannot resolve both the boundary layers on the turbine
blades and the turbulent structures in the wake. This necessitates a simplified representation of the wind turbine
in case of wake calculations and a simplified representation of the wake in case of blade calculations.
800
Wind Energ.
2011; 14:799–819 © 2011 John Wiley & Sons, Ltd.
DOI: 10.1002/we
B. Sanderse, S. P. v. d. Pijl and B. Koren Review of CFD for wind turbine wake aerodynamics
2.2. Turbulence modeling
A l arge number of turbulence models have been constructed in the last decennia (see e.g. Wilcox,
14
Sagaut
15
and Geurts
16
).
This section will discuss the two most important methodologies in turbulence modeling for wind turbine wakes, RANS and
LES, their applicability and their limitations.
2.2.1. RANS
Reynolds-averaged Navier–Stokes methods aim for a statistical description of the flow. Flow quantities such as velocity
and pressure are split in an average and a fluctuation, the so-called Reynolds decomposition:
u.x;t/D
u.x/ C u
0
.x;t/ (2)
The averaging procedure, ensemble averaging, is such that
u.x/ D u.x/ and u
0
.x;t/ D 0. The Reynolds decomposition
[equation (2)] is substituted into the Navier–Stokes equations, which are then averaged, resulting in
12
@u
@t
C .
u r/u D
1
r
p C r
2
u r.u
0
u
0
/ (3)
The term
u
0
u
0
is called the Reynolds stress tensor, which appears as a consequence of the non-linearity of the convective
term, and represents the averaged momentum transfer because of turbulent fluctuations. The Reynolds stresses can be inter-
preted as turbulent diffusive forces. In wind turbine wakes, they are much larger than the molecular diffusive forces r
2
u,
except near solid boundaries. In order to close the system of equations, a model is needed to express the Reynolds stresses
in terms of mean flow quantities.
A widely adopted approach of modeling the Reynolds stresses exploits the Boussinesq hypothesis
17
(not to be confused
with the Boussinesq approximation mentioned earlier). Based on an analogy with laminar flow, it states that the Reynolds
stress tensor can be related to the mean velocity gradients via a turbulent ‘eddy’ viscosity
T
:
u
0
u
0
D
T
r
u C .ru/
T
(4)
so that the RANS equation (3) becomes
@
u
@t
C .
u r/u D
1
rNp Cr
. C
T
/.ru C .ru/
T
/
(5)
This approach of modeling the effect of turbulence as an added viscosity is widely used for turbulent flow simulations. It is
very useful as an engineering method, because the computational time is only weakly dependent on the Reynolds number.
However, the validity of the Boussinesq hypothesis is limited. In contrast to ,
T
is not a property of the fluid but rather
a property of the type of flow in question. Since eddies are fundamentally different from molecules, there is no sound
physical basis for equation (4),
14
and DNS calculations have indeed not shown a clear correlation between u
0
u
0
and ru.
18
The Boussinesq hypothesis is therefore inadequate in many situations, for example, for flows with sudden changes in mean
strain rate (e.g. the shear layer of the wake), anisotropic flows (e.g. the atmosphere) and three-dimensional (3D) flows.
12,14
Many different methods have been suggested to calculate
T
, typically called zero-equation (algebraic closure, such as
mixing length), one-equation and two-equation models (see e.g. Wilcox
14
and the references therein). The k model is
an example of a two-equation model often encountered in wind energy wake applications; the k ! model [with shear
stress transport (SST) limiter] is more convenient near blade surfaces. In the k model, two additional partial differen-
tial equations are introduced, one for the turbulent kinetic energy k and one for the turbulent diffusion . They contain a
number of constants that have been determined by applying the model to some very general flow situations (e.g. isotropic
turbulence decay or flow over a flat plate).
A fundamentally different approach is the Reynolds stress model (RSM),
19
also called differential second-moment clo-
sure model, which does not rely on the Boussinesq hypothesis. In the RSM, all the components of the Reynolds stress tensor
are modeled, which makes it suitable for anisotropic flows. However, it leads to six additional partial differential equations
(PDEs), making the approach expensive. Moreover, these PDEs contain terms that have to be modeled again, and often,
closure relations resembling the Boussinesq hypothesis are still employed. Lastly, the disappearance of the (stabilizing)
eddy viscosity term can lead to numerical problems.
2.2.2. LES
In recent years, LES is receiving more attention in the wind energy wake community, because of its ability to han-
dle unsteady, anisotropic turbulent flows dominated by large-scale structures and turbulent mixing. This is a significant
Wind Energ.
2011; 14:799–819 © 2011 John Wiley & Sons, Ltd.
DOI: 10.1002/we
801
Review of CFD for wind turbine wake aerodynamics B. Sanderse, S. P. v. d. Pijl and B. Koren
advantage over RANS methods, but the drawback is that the computational requirements for LES are much higher than
for RANS. In LES, the large eddies of the flow are calculated, whereas the eddies smaller than the grid are modeled with
a subgrid-scale model. This is based on the assumption that the smallest eddies in the flow have a more or less universal
character that does not depend on the flow geometry. Mathematically, this scale separation is carried out by spatially fil-
tering the velocity field, splitting it in a resolved (also called large scale, simulated or filtered) velocity and an unresolved
(small scale) part. In general, this filtering operation is defined as a convolution integral:
eu.x;t/D
Z
u.;t/G.x ;/d (6)
where G.x ;/ is the convolution kernel, depending on the filter width . The subgrid velocity is then defined as the
difference between the flow velocity and the filtered velocity:
u
0
.x;t/D u.x;t/ eu.x;t/ (7)
This decomposition resembles that of Reynolds averaging but with the difference that in general
e
eu ¤Qu and
e
u
0
¤ 0.
Applying the filtering operation to the Navier–Stokes equations (and assuming certain properties of the filter) leads to the
following equation:
@eu
@t
C .eu r/eu D
1
rep C r
2
eu r.fuu eueu/ (8)
As in the RANS equation (3) a new term appears, the subgrid-scale (SGS) stresses. These stresses represent the effect of the
small (unresolved) scales on the large scales. A widely used model to calculate these stresses is the Smagorinsky model,
20
which again employs the Boussinesq hypothesis:
SGS
D fuu eueu D
SGS
reu C .reu /
T
(9)
Possible ways to calculate the subgrid-scale eddy viscosity
SGS
are to use an analogy of the mixing-length formulation
and to use one-equation or two-equation models involving kinetic energy and turbulent dissipation. Because of the use of
the Boussinesq hypothesis, similar limitations as in RANS are encountered. A large number of other subgrid-scale models
have therefore been proposed, for example, dynamic models, regularization models and variational multiscale models (see
e.g. Sagaut
15
and Geurts
16
).
In contrast to RANS, where the computational cost is only weakly dependent on Re, the computational cost of LES
scales roughly with Re
2
. Near solid boundaries, where boundary layers are present, LES is extremely expensive because
it requires refinement in three directions, whereas RANS only requires refinement in the direction normal to the wall. A
possibility is to employ a hybrid approach: RANS to resolve the attached boundary layers and LES outside the wall region,
so-called detached eddy simulation (DES).
21,22
Since equations (4) and (9) both have a similar form, this switch between
RANS and LES can be made by changing the eddy viscosity based on a wall-distance function.
2.3. Boundary conditions
One of the ongoing challenges in CFD simulations of wind turbine wakes, especially when comparing with experimental
data, is the prescription of inflow conditions that mimic all relevant characteristics of the atmosphere, such as the sheared
velocity profile, the anisotropy of the turbulence and the instationary nature of the inflow. In early CFD simulations, uni-
form, laminar inflow profiles were used, but LES simulations showed later that both the presence of the shear inflow profile
(the ABL)
23
as well as the turbulence in the incoming flow
24
have a pronounced effect on the flow field behind the rotor.
For R ANS simulations, Monin–Obukhov similarity theory (see e.g. Panofsky and Dutton
25
) can be used to prescribe
profiles for the velocity components and turbulence quantities, which are independent of time. For LES, unsteady inflow
data are necessary, and basically, two different types of methods exist: synthesized inlet methods and precursor simulation
methods.
26
The synthesis technique was employed in the models of Mann,
27,28
Veer s
29
and Winkelaar.
30
They are based
on the construction of spectral tensors, which model the frequency content of the wind field, and have the advantage that
certain parameters like turbulent length or time scales can be directly specified. Mann’s model combines rapid distortion
theory (based on a linearization of the Navier–Stokes equations) and the von Kármán spectral tensor to generate a 3D
incompressible atmospheric turbulence field, which is homogeneous, stationary, Gaussian and anisotropic but does not
include the effect of the ground boundary. Troldborg used this model to create a turbulent inflow profile.
24,31
This profile
is introduced close to the rotor (one rotor radius upstream) to prevent the decay of turbulent fluctuations before they reach
the rotor. In order to ensure a divergence-free velocity field, the turbulent fluctuations are introduced via body forces in
the momentum equation, following Mikkelsen et al.
32
In the second category, a separate precursor simulation is made as
802
Wind Energ.
2011; 14:799–819 © 2011 John Wiley & Sons, Ltd.
DOI: 10.1002/we
B. Sanderse, S. P. v. d. Pijl and B. Koren Review of CFD for wind turbine wake aerodynamics
an inflow model for a successor simulation. The advantages of this method are that the ground boundary is taken into
account and that the resulting turbulent velocity field is a solution to the non-linear Navier–Stokes equations. However, it
cannot be easily manipulated to obtain certain desired turbulence characteristics, and it is computationally more expensive.
In the work of Bechmann,
33
Meyers et al.
34
and Stovall et al.,
35
such a method was used. In Stovall et al.,
35
roughness
blocks were introduced to generate the turbulence. When a desired turbulence level is achieved, the blocks are removed,
and the turbines are introduced in the domain. The model of Bechmann can handle the anisotropy of the atmosphere, except
near the surface, where it switches to a RANS solver. As was stated earlier by Vermeer,
3
Bechmann observed again that
‘imposing realistic and even experimentally obtained inflow conditions was found extremely difficult’.
At other boundaries of the computational domain, the prescription of boundary conditions is not straightforward either.
The most physical model for the presence of the ground is a body-fitted mesh with a no-slip boundary condition, but such a
condition does not allow the prescription of a ground roughness and would require very fine grids near the surface. There-
fore, a common approach in both RANS and LES models is to use wall functions to prescribe wall friction. Porté-Agel
et al. prescribed the surface shear stress and heat flux by adopting a time-dependent variant of Monin–Obukhov similarity
theory.
11
An alternative to no-slip conditions are slip conditions,
23,36
so that a laminar shear profile can be imposed with-
out depending on the building-up of a boundary layer associated with the Reynolds number of the atmosphere. A turbulent
inflow is difficult to prescribe, because the turbulent viscosity of the ABL is much (10
6
times) larger than the molecular vis-
cosity, making the effective Reynolds number of the blades much too low. At the upper boundary, symmetry conditions are
often prescribed, but care should be taken when the height of the ABL is of the same order as the turbine height, something
which is not unusual for modern-size turbines, especially at night time. In streamwise direction, periodicity conditions can
be used to simulate infinite wind farms or to generate a turbulence field with the precursor method discussed before, so that
an ABL can form without the need for a very large domain. Periodicity conditions are also convenient from a numerical
point of view because they allow the application of fast Poisson solvers like the fast Fourier transform.
3. ROTOR MODELING
To solve the RANS equation (3) or LES equation (8) in the near and far wakes of a wind turbine, a representation of the
blades is necessary. Basically, two approaches exist: the generalized actuator disk approach, in which the blades are repre-
sented by a body force (Section 3.1), and the direct approach, in which the presence of the blades is taken into account by
discretizing the actual blades on a computational mesh (Section 3.2).
3.1. Generalized actuator disk modeling
In many near and far wake calculations, the rotor is represented by an actuator disk or an actuator line. Such a representation
circumvents the explicit calculation of the blade boundary layers, reducing computational cost and easing mesh generation.
The actuator disk exerts a force on the flow, acting as a momentum sink. This force is explicitly added to the momentum
equation (1):
Z
˝
@u
@t
d˝ C
Z
@˝
uu n dS D
Z
@˝
1
p n dS C
Z
@˝
ru n dS C
Z
A\˝
f dA (10)
which are written in weak form, because the force leads to a discontinuity in pressure. Apart from this momentum sink, one
should also introduce sources of turbulence corresponding to the mechanical turbulence generated by the blades. Currently,
three different approaches for prescribing the force term f exist: the actuator disk, the actuator line and the actuator surface
models (see Figure 1).
AL
AD
AS
Figure 1. Illustration of the actuator disk (AD), line (AL) and surface (AS) concept.
Wind Energ.
2011; 14:799–819 © 2011 John Wiley & Sons, Ltd.
DOI: 10.1002/we
803
