Computational Aerodynamics and

mustardarchaeologistMechanics

Feb 22, 2014 (3 years and 3 months ago)

217 views

Computational Aerodynamics and
Aeroacoustics for Wind Turbines
Computational Aerodynamics and Aeroacoustics for
Wind Turbines
Wen Zhong Shen
Fluid Mechanics
Department of Mechanical Engineering
TECHNICAL UNIVERSITY OF DENMARK



























Published in Denmark by
Technical University of Denmark

Copyright © Wen Zhong Shen
All Rights reserved

Section of Fluid Mechanics
Department of Mechanical Engineering
Technical University of Denmark
Nils Koppels Allé, Building 403, DK-2800 Kgs. Lyngby, Denmark
Phone + 45 4525 4300, Telefax + 45 4588 4325
E-mail: info.fm@mek.gtu.dk
WWW: http://www.mek.dtu.dk/

Publication Reference Data
Shen, Wen Zhong
Computational Aerodynamics and Aeroacoustics for Wind Turbines
Doctor Thesis
Technical University of Denmark, Section of Fluid Mechanics
Printed: October 2009
ISBN 978-87-89502-87-8
Keywords: wind turbine; rotor aerodynamics; wakes; aeroacoustics; Navier-Stokes equations; vorticity-
velocity formulation; SIMPLEC; computational fluid dynamics (CFD); computational aero
acoustics (CAA); flow-acoustic splitting technique; large eddy simulation (LES); actuator disc;
actuator line; actuator surface; tip loss correction

Computational Aerodynamics and Aeroacoustics i

Preface

This thesis summarizes my work on computational aerodynamics and aero-acoustics
for wind turbines. Nineteen papers have been selected and listed on page v. Seventeen
of them were published, or are in press, in international journals and the other two are
conference papers.
The work presented here is divided into two parts: an aerodynamic part
(Computational Fluid Dynamics) and an aero-acoustic part (Computational Aero
Acoustics) for wind turbines. The main objective of the research was to develop new
computational tools or approaching techniques for analysing wind turbine flows. A
few papers deal with applications of Blade Element Momentum (BEM) theory to
wind turbines.
The research has been carried out mainly at the Fluid Mechanics Section,
Department of Mechanical Engineering, Technical University of Denmark in the
period from 1996 to 2008. The friendly and cooperative research atmosphere in the
section provided a very good environment. Most of the work presented in the thesis
has been sponsored by the European Commission, the Danish Research Council or the
Danish Energy Agency.
During the preparation of this thesis, my main collaborator Professor Jens Nørkær
Sørensen gave me much encouraging support. Here I wish to give him my sincere
acknowledgements for leading me into the field of wind energy and constantly
providing inspiring advice and support. Moreover, I appreciate the contributions of
the following colleagues or students from DTU or abroad (in alphabetical order):
Christian Bak, Martin O. L. Hansen, Ta Phuoc Loc, Jess A. Michelsen, Robert
Mikkelsen, Xavier Munduate, Niels N. Sørensen, Jianhui Zhang and Weijun Zhu.
Furthermore, I would like to express my acknowledgements to the committee
members (Prof Steen Krenk, Prof Pierre Sagaut and Prof Laszlo Fuchs) for their
constructive suggestions. Secretary Ruth S. Vestergaard and Vivi Jessen are also
acknowledged for helping with printing the thesis. Many other researchers and PhD
students working in the Fluid Mechanics Section are also thanked for all kinds of help
and discussion. Finally, special thanks go to my wife Xu and two daughters (Julia and
Anna) for enjoyable family life.
ii Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics iii

Abstract


To analyse the aerodynamic performance of wind turbine rotors, the main tool in use
today is the 1D-Blade Element Momentum (BEM) technique combined with 2D
airfoil data. Because of its simplicity, the BEM technique is employed by industry
when designing new wind turbine blades. However, in order to obtain more detailed
information of the flow structures and to determine more accurately loads and power
yield of wind turbines or cluster of wind turbines, it is required to resort to more
sophisticated techniques, such as Computational Fluid Dynamics (CFD). As computer
resources keeps on improving year by year (about ten times every five years from
statistics over the last twenty years), CFD has now become a popular tool for studying
the aerodynamics of wind turbines.
The present thesis consists of 19 selected papers dealing with the development and
use of CFD methods for studying the aerodynamics and aero-acoustics of wind
turbines. The papers are written in the period from 1997 to 2008 and numbered
according to the list in page v. The work consists of two parts: an aerodynamic part
based on Computational Fluid Dynamics and an aero-acoustic part based on
Computational Aero Acoustics for wind turbines. The main objective of the research
was to develop new computational tools and techniques for analysing flows about
wind turbines. A few papers deal with applications of Blade Element Momentum
(BEM) theory to wind turbines.
In most cases the incompressible Navier-Stokes equations in primitive variables
(velocity-pressure formulation) are employed as the basic governing equations.
However, since fluid mechanical problems essentially are governed by vortex
dynamics, it is sometimes advantageous to use the concept of vorticity (defined as the
curl of velocity).

In vorticity form the Navier-Stokes equations may be formulated in
different ways, using a vorticity-stream function formulation, a vorticity-velocity
formulation or a vorticity-potential-stream function formulation.
In [1] - [3] two different vorticity formulations were developed for 2D and 3D
wind turbine flows. In [4] and [5] numerical techniques for avoiding pressure
oscillations were developed when solving the velocity-pressure coupling system in the
in-house EllipSys2D/3D code, which originally was developed in a cooperation
between DTU (Michelsen, 1992) and Risø (Sørensen, 1995). In [6] – [8] different
actuator disc techniques combined with CFD are presented. This includes actuator
disc, actuator line and actuator surface techniques, which were developed to simulate
flows past one or more wind turbines. In [9] and [10] a tip loss correction method that
improves the conventional models was developed for use in combination with BEM
or actuator/Navier-Stokes computations. A simple and efficient technique for
determining the angle of attack for flow past a wind turbine rotor was developed in
[11], and in [12] tunnel wall corrections for wind tunnels with closed or open test
sections were developed.
The second part of the thesis deals with Computational Aero-Acoustics (CAA).
With the spread of wind turbines near urban areas, there is an increasing need for
accurate predictions of aerodynamically generated noise. Indeed, noise has become
one of the most important issues for further development of wind power, and the
iv Wen Zhong Shen

ability of controlling and minimising noise emission may be advantageous when
competing on the world energy market.
To predict generation and propagation of aerodynamic noise, it is required to
solve the compressible Navier-Stokes equations. As the scales of the flow and the
acoustic waves are quite different (about 1/M, M=Mach number=U
￿
/c), it is difficult
to resolve them together at the same time. Hardin and Pope proposed a non-linear
two-step (viscous incompressible flow and inviscid acoustic perturbation) splitting
procedure for computational aero-acoustics that is suitable for both generation and
propagation. The advantage of the splitting approach, as compared to the acoustic
analogies, is that the source strength is obtained directly and that it accounts for sound
radiation as well as scattering.
In [13] and [14] an inconsistency in the original formulation of Hardin and Pope
1994 was analysed and a consistent formulation was proposed and applied to laminar
flows. An aero-acoustic formulation for turbulent flows was in [15] developed for
Large Eddy Simulation (LES), Unsteady Reynolds Averaged Navier-Stokes
Simulation (URANS) and Detached Eddy Simulation (DES). In [16] a collocated grid
/ finite volume method for aero-acoustic computations was developed and
implemented in the EllipSys2D/3D code. In [17] and [18] three dimensional flow-
acoustic computations were carried out. Finally, the aero-acoustic formulation using
high order Finite Difference schemes (Dispersion Relation Preserving (DRP) /
Optimized Compact schemes) was developed in [19] and implemented in the
EllipSys2D/3D code.
Computational Aerodynamics and Aeroacoustics v
List of thesis papers
[1] W. Z. Shen and J. N. Sørensen. Quasi-3D Navier-Stokes model for a rotating airfoil.
Journal of Computational Physics, 150: 518-548, 1999.
[2] W. Z. Shen and T. P. Loc. Numerical methods for unsteady 3D Navier-Stokes
equations in velocity and vorticity form. Computers & Fluids, 26: 193-216, 1997.
[3] M. O. L. Hansen, J. N. Sørensen and W. Z. Shen. Vorticity-velocity formulation of
the 3D Navier-Stokes equations in cylindrical coordinates. International Journal for
Numerical Methods in Fluids, 41:29-45, 2003.
[4] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. Improved Rhie-Chow interpolation
for unsteady flow computations. AIAA Journal, 39: 2406-2409, 2001.
[5] W. Z. Shen, J. A. Michelsen, N. N. Sørensen and J. N. Sørensen. An improved
SIMPLEC method on collocated grid for steady and unsteady flow computations.
Numerical Heat Transfer, Part. B, 43: 221-239, 2003.
[6] J. N. Sørensen, W. Z. Shen and X. Munduate. Analysis of wake states by a full-field
actuator disc model. Wind Energy, 1: 73-88, 1998.
[7] J. N. Sørensen and W. Z. Shen. Numerical modeling of wind turbine wakes. Journal
of Fluid Engineering, 124: 393-399. 2002.
[8] W. Z. Shen, J. N. Sørensen and J. H. Zhang. The actuator surface model: A new
Navier-Stokes based model for rotor computations. Journal of Solar Energy
Engineering, 131: 011002 (9 pages), 2009.
[9] W. Z. Shen, R. Mikkelsen, J. N. Sørensen and Ch. Bak. Tip loss corrections for wind
turbine computations. Wind Energy, 8: 457-475, 2005.
[10] W. Z. Shen, J. N. Sørensen and R. Mikkelsen. Tip loss correction for
actuator/Navier-Stokes computations. Journal of Solar Energy Engineering, 127:
209-213, 2005.
[11] W. Z. Shen, M. O. L. Hansen and J. N. Sørensen. Determination of the angle of
attack for rotor blades. Wind Energy, 12: 91-98,

2009.
[12] J. N. Sørensen, W. Z. Shen and R. Mikkelsen. Wall correction model for wind
tunnels with open test section. AIAA Journal, 44: 1890-1894, 2006.
[13] W. Z. Shen and J. N. Sørensen. Comment on the aeroacoustic formulation of
Hardin and Pope. AIAA Journal, 37: 141-143, 1999.
[14] W. Z. Shen and J. N. Sørensen. Aeroacoustic modelling of low-speed flows.
Theoretical and Computational Fluid Dynamics, 13: 271-289, 1999.
[15] W. Z. Shen and J. N. Sørensen. Aeroacoustic modeling of turbulent airfoil flows.
AIAA Journal, 39: 1057-1064, 2001.
[16] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. A collocated grid finite volume
method for aeroacoustic computations of low-speed flows. Journal of Computational
Physics, 196: 348-366, 2004.
[17] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. Aeroacoustic computations of wind
turbines. AIAA paper, No. 2002-0043, Reno, Nevada, USA, 2002.
[18] W. Z. Shen, W. J. Zhu and Sørensen. Aero-acoustic computations for turbulent
airfoil flow. AIAA Journal, 47: 1518-1527, 2009.
[19] W. J. Zhu, W. Z. Shen and J. N. Sørensen. Computational Aero-Acoustic using
high-order finite differences schemes. Journal of Physics: Conference Series 75, doi:
10.1088/1742-6596/75/1/012084 (11 pages), 2007.
vi Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics vii
Contents

Preface.............................................................................................................................i
Abstract.........................................................................................................................iii
List of thesis papers.......................................................................................................v
Contents.......................................................................................................................vii
1. Introduction................................................................................................................1
2. Vorticity formulations [1-3].......................................................................................5
A quasi-3D Navier-Stokes model for a rotating airfoil [1]....................................5
Vorticity-velocity formulations [2-3]....................................................................8
3. Velocity-pressure coupling [4-5]...............................................................................9
4. Actuator models [6-8]..............................................................................................13
The actuator disc / Navier-Stokes model [6].......................................................13
The actuator line / Navier-Stokes model [7]........................................................15
The actuator surface/Navier-Stokes model [8]....................................................15
5. Tip loss corrections [9-10].......................................................................................19
6. Miscellaneous aerodynamic problems [11-12]........................................................25
Determination of AOA for a rotating wind turbine [11]......................................25
Wall correction model for wind tunnels with open test section [12]...................27
7. Development of Aeroacoustic formulations [13-15]...............................................29
Improvement on an existing Aero-Acoustic formulation [13-14].......................29
Aero-acoustic formulation for turbulent flow [15]..............................................30
8. Second-order finite volume/aero-acoustic code with applications [16-18].............33
9. High order finite difference aero-acoustic code [19]...............................................37
10. Concluding remarks...............................................................................................39
References....................................................................................................................41
Dansk Resumé.............................................................................................................47

viii Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics 1
1. Introduction

Computational Fluid Dynamics (CFD) started in the early sixties of the last century
and is now becoming a mature tool for predicting flow fields of different kinds of
fluid dynamic problems. During the past 30 years, many computational techniques
and methods have been developed: Finite Difference methods (for example, Arakawa
1966, Daube-Loc 1978 and Tannehill et al. 1997), Finite Volume methods (for
example, Patankar 1980, Ferziger-Peric 1996 and Akeselvoll-Moin 1996), Finite
Element methods (for example, Ciarlet-Lions 1996, Dhatt-Touzot 1984 and Girault-
Raviart 1986), Spectral methods (for example, Canuto et al. 1988) and Vortex
methods (for example, Chorin 1973, Leonard 1980 and Raviart 1985), and many
hybrid methods (for example, the spectral-element method of Karniadakis-Sherwin
2005 and the finite difference-vortex method of Shen 1993, Guermond et al. 1993,
Shen-Loc 1995 and 1997). Among those, Finite Difference / Finite Volume methods
are particularly useful in computational fluid dynamics because of their simplicity of
implementation and efficiency for parallel computations on multiprocessor computers
with Message Passing Interfaces (MPI).
The first part of the thesis deals with Computational Fluid Dynamics (CFD). For
the flow past a wind turbine, the governing equations are known to be the
compressible Navier-Stokes equations. Since the relative velocity (U
￿
) on wind
turbine blades is small when compared to the speed of sound (c￿340 m/s) and the
flow and acoustic sound wave scales are very different (about c/U
￿
times), it is very
difficult to solve them together within a reasonable accuracy. If only the flow is of
interest, one can solve the incompressible Navier-Stokes equations obtained from the
compressible Navier-Stokes equations by using the assumption that the change in
density can be neglected.
In order to get a flow solution, it is natural to solve the incompressible Navier-
Stokes equations in primitive variables (velocity-pressure formulation). Since many
fluid mechanical problems essentially are governed by vortex dynamics, it can be
advantageous to use the concept of vorticity (defined as the curl of velocity). There
are different formulations based on vorticity, such as the vorticity-stream function,
vorticity-velocity and vorticity-potential-stream function formulations.
In [1-3] two vorticity formulations were developed for 2D and 3D wind turbine
flows. In [4-5], some numerical techniques for avoiding pressure oscillations were
developed when solving the velocity-pressure coupling system used in the in-house
EllipSys2D/3D code developed at DTU by Michelsen 1992 and at RISØ by Sørensen
1995. In [6-8] the actuator techniques including actuator disc, actuator line and
actuator surface techniques were developed to simulate flows past one or more wind
turbines. In [9-10] a tip loss correction method that improves the classical correction
models was developed for tip losses for BEM or actuator/Navier-Stokes computations.
A simple and efficient technique for determining the angle of attack for flow past a
wind turbine was developed in [11]. In [12] a tunnel wall correction method for wind
tunnels with closed or open test sections was studied.
The second part of the thesis deals with Computational Aero-Acoustics (CAA). It
is well known that noise can be very disturbing on human daily life; hence there is an
increasing need to accurately predict aerodynamically generated noise. For wind
turbines, the noise has become an important issue for further development, and
reduction control of noise emission can have great advantages in the world energy
market.
2 Wen Zhong Shen

In order to predict noise generated from wind turbines, an easy way is to employ a
semi-empirical model based on a series of wind tunnel measurements (for example,
Brooks et al. 1989, Lowson 1993, Fuglsang-Madsen 1996 and Zhu et al. 2005). A
semi-empirical model can give an over-all noise level at a specific observer position.
Since the model is experimentally based, there are many questions that may be raised
regarding accuracy and applicability. Another kind of prediction methods are based
on the Lighthill acoustic analogy resulting from linearized wave equations combined
with sound sources, such as monopoles, dipoles and quadrupoles (for example,
Lighthill 1952, Ffowcs Williams-Hawkings 1968, Lyrintzis 1994 and Farassat-Succi
1983 and Farassat and Brentner 1988).
In order to accurately predict noise, it is natural to solve the compressible Navier-
Stokes equations because it contains laws for the acoustic generation and propagation.
For subsonic turbulent flows, the cost of a Direct Numerical Simulation (DNS) for
both flow and acoustics is proportional to
3 4
Re/M
(according to Colonius 1997),
where
Re/UD
￿
￿￿
is the Reynolds number and /M U c
￿
￿
is the Mach number.
Hence, DNS is limited to low Reynolds number flows. In order to consider sound
radiation from high Reynolds number flows, Bogey and co-workers (Bogey, Bailly
and Juvé 2003, Bogey and Bailly 2005, 2006) studied the sound generation from a
turbulent jet by solving the compressible Navier-Stokes equations using Large Eddy
Simulation. As the scales of flow and acoustic waves are quite different (about 1/M),
it is rather difficult to resolve them together at the same time. Hardin and Pope 1994
proposed a non-linear two-step (viscous incompressible flow and inviscid acoustic
perturbation) splitting procedure (SP) for computational aero-acoustics that is suitable
for both generation and propagation. The advantage of the splitting approach is that
the source strength is obtained directly and that it accounts for sound radiation as well
as scattering. Alternatively, Béchara et al. 1994 developed an approach using
linearized Euler’s equations (LEE) to study the sound generation from shear mean
flows. Subsequently, discussions on the different source terms in the LEE approach
were carried out by Bogey et al. 2002. The third approach, referred as the nonlinear
disturbance equations (NLDE), was developed by Morris and his co-workers (Morris
et al. 1997, 1998). The NLDE equations were obtained from the compressible Navier-
Stokes equations by subtracting the mean compressible flow equations from the
instantaneous equations. The technique was combined with a zonal RANS/LES
approach for noise source predictions by Terracol and co-workers (Labourasse and
Sagaut 2002, Terracol et al. 2005 and Terracol 2006). Another approach, denoted the
acoustic perturbation equations (APE), was developed by Ewert and Schröder 2003,
where the concepts of enthalpy and the source-term filtering are used to obtain
acoustic sources from a combination of acoustic and vorticity sources. Subsequently,
the APE approach was used to study the noise generation from a turbulent jet by
Gröschel et al. 2008.
In [13-14] an inconsistency of the original formulation of Hardin and Pope 1994
was analysed and a consistent formulation was proposed and applied to laminar flows.
An aero-acoustic formulation for turbulence flows was developed in [15] for Large
Eddy Simulations (LES), Unsteady Reynolds Averaged Navier-Stokes Simulations
(URANS) and Detached Eddy Simulations (DES). In [16] a collocated grid / finite
volume method for aero-acoustic computations was developed and implemented in
the in-house EllipSys2D/3D code. In [17-18] three dimensional flow-acoustic
computations were carried out. Finally, the aero-acoustic formulation using high order
Finite Difference schemes (Dispersion Relation Preserving (DRP) / Optimized
Computational Aerodynamics and Aeroacoustics 3
Compact schemes) was developed in [19] and implemented in the EllipSys2D/3D
code.
4 Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics 5
2. Vorticity formulations [1-3]

The concept of vorticity has been used to solve the incompressible Navier-Stokes
equations since the sixties of the last century. It was first formulated through the use
of the stream function to solve two-dimensional (2D) problems. The advantage of
such a formulation is that the continuity equation is satisfied completely, see Daube-
Loc 1987. For three dimensional (3D) flows, the stream function is a vector and has
three components. Since it was difficult to find the correct boundary conditions at
walls for the stream function, vorticity-velocity formulations were studied very
intensively by e. g. Daube 1992, Gatski 1991, Stella-Guj 1989, Guj-Stella 1993,
Hansen 1994, Quartapelle 1993, Speziale 1987 and in [2-3].
In order to study the rotating effects on a rotating wind turbine blade, a quasi-3D
model in a rotating frame of reference was developed in [1] and solved using the
formulation of the vorticity-stream function. In the paper [2], a general 3D vorticity-
velocity formulation was developed and in [3] a 3D vorticity-velocity formulation
based on the Cauchy-Riemann equations was developed in cylindrical coordinates.
The developed vorticity-velocity based code was used to simulate flows past a wind
turbine in [6-7], with the turbine modelled by Actuator Disc and Actuator Line
techniques.
A quasi-3D Navier-Stokes model for a rotating airfoil [1]
For an airfoil of a wind turbine blade in rotation it is important to model centrifugal
and Coriolis forces, since these forces may have a significant influence on the
aerodynamics, as compared to a two dimensional airfoil. The Coriolis force in the
chordwise direction acts as a favourable pressure gradient that tends to delay
boundary layer separation and increase the lift force coefficient. In order to take into
account rotating effects, a quasi-3D Navier-Stokes model was developed in [1]. The
hypotheses used in the model are based on Sears 1950 and Fogarty-Sears 1950
equations for the inviscid velocity distributions along a rotating blade. The hypotheses
are summarized as
,
,
0.
u u
z z
v v
z z
w
z
￿
￿
￿
￿
￿
￿
￿
￿
￿

In order to satisfy the mathematical identity
￿ ￿
0curl grad f ￿￿
, the total pressure
must satisfy the following property, which then becomes an assumption for the static
pressure
￿ ￿
2
2 2 2
2
2 2 2
/2 ( )/2/2
/2 ( )/2/2
p u r r
p u r r
z z
￿ ￿ ￿ ￿ ￿￿
￿ ￿ ￿ ￿￿
￿
￿
￿
￿ ￿ ￿
￿ ￿


.
Based on these hypotheses, the quasi-3D Navier-Stokes equations are obtained
and formulated into a vorticty-stream function-spanwise velocity form. The
formulation has been expressed in general curvilinear coordinates and solved on an
orthogonal staggered mesh generated from a conformal mapping. The vorticity
equations are discretized in time using a semi-implicit, combined Adams-
Bashforth/Crank Nicolson scheme that has second order accuracy and second order
6 Wen Zhong Shen

central differences in space, except that the convective terms are discretized by the
upwind QUICK scheme. The resulting matrix is solved iteratively using the
Successive Line Over-Relaxation (SLOR) method. The stream function and spanwise
velocity equations are solved by first applying Fourier transformation in the periodic
tangential direction and discretizing in the frequency domain with second-order finite
difference discretizations in the radial direction.

Figure 1: Streamlines for flow around a NACA 63
2
-415 airfoil at incidence of 20
o
and
Reynolds number of Re=1.5x10
6
, (a) 2D, (b) k=6, (c) k=4, from [1].

For turbulent flow, the one-equation model of Baldwin-Barth 1990 is used. From
the early European project VISCWIND (see Sørensen 1999, Sørensen and Nygreen
2001), different turbulence models including the Baldwin-Lomax 1978, Spalart-
Allmaras 1994 and k-￿ Baseline and Shear Stress Transport models of Menter 1994
Computational Aerodynamics and Aeroacoustics 7
were compared for turbulent flow past an airfoil. The outcome is that the Baldwin-
Barth model works as well as the Spalart-Allmaras and k-￿ SST models. For more
details about the performances of different turbulence models, the reader is referred to
Sørensen 1998. As an illustration, streamlines for flow past a NACA 63
2
-415 airfoil at
an angle of attack of 20
o
and a Reynolds number of 1.5 x10
6
are plotted in Figure 1.
From the figure, it is seen that the separation bubble becomes smaller when the
position of the airfoil section is closer to the rotor centre (measured by the parameter
k=r/chord where r is the distance to the rotor centre). The normal force coefficient is
increased significantly and a small increase in the tangential force coefficient is also
seen.


Figure 2: Cn and Ct versus angle of attack for a NACA 63
2
-415 airfoil at
6
Re 1.5 10
￿ ￿
,
from [1].

In Figure 2, the normal and tangential force coefficients are plotted. From the
figure, it is seen that both force coefficients increase when the airfoil section is closer
to the rotor centre. As a conclusion, the model is capable of determining the
8 Wen Zhong Shen

qualitative behaviour for airfoils subject to rotation. By using the model it is possible
to derive airfoil data for use in engineering predictive codes.
Vorticity-velocity formulations [2-3]
In [2] several vorticity-velocity formulations were discussed and followed with a
suggested formulation. In the paper, it was proved that the new formulation is
equivalent to the Navier-Stokes equations in velocity-pressure variables. The outcome
of the formulation is that: (1) the vorticity transport equation written in full curl form
can control the transport of vorticity with a divergence free vorticity field:
￿ ￿
u
t
￿
￿ ￿ ￿
￿
￿￿￿ ￿ ￿ ￿ ￿￿￿￿
￿
￿
￿ ￿￿
;
(2) the velocity field can be found by using two components of the velocity Poission
equations and the continuity equation in order to reduce computing cost. The vorticity
code based on the proposed formulation on a staggered mesh had been validated
against the velocity-pressure based ONERA Pegase code by Sagaut 1995 for the cases
of the flow past a cube and the experimental data in Schlichting 1968 for the case of
flow past a sphere.
In [3] a vorticity-velocity formulation based on solving the Cauchy-Riemann
equations for the velocity was proposed. The formulation was used to solve the 3D
Navier-Stokes equations in cylindrical coordinates. Special treatments are needed at
the polar singularity axis. The code was validated against experimental data for the
cases of axisymmetric and non-axisymmetric flows in a pipe and in a closed cylinder.
As an illustration, streamlines for the flow past a spinning sphere are plotted in
Figure 3. From the figure, it is seen that the treatment of polar singularity is smooth
and consistent.


Figure 3: Flow past a rotating sphere at Re=500 and
/
U a
￿
￿￿
, from [2].
Computational Aerodynamics and Aeroacoustics 9
3. Velocity-pressure coupling [4-5]

The EllipSys2D/3D code developed as a collaborative work between DTU and Risø,
Michelsen 1992 and Sørensen 1995, is based on the velocity-pressure formulation. As
a vorticity formulation can satisfy the continuity equation or mass conservation
completely or to a high degree by using the stream function or by discretizing directly
the continuity equation, it is preferable to use it for solving two dimensional flows.
For three dimensional flows, the number of governing equations for a vorticity
formulation is increased to 6 (3 vorticity and 3 velocity equations) or 7 (plus the
continuity equation) and it is computationally more expensive than a velocity-pressure
formulation. Therefore, a velocity-pressure formulation needs also to be considered.
The EllipSys code uses cell-centred finite-volume discretizations and the multiblock
strategy. It is suitable for parallel computations with Message Passing Interface (MPI).
In order to couple the velocity and the pressure in a predictor-corrector method, Rhie-
Chow 1983 proposed a procedure using a momentum-based interpolation for the cell
face mass fluxes in the continuity equation so that the mass fluxes are driven by
1
￿

pressure difference across the faces. Hence, the velocity-pressure decoupling cannot
occur.
Figure 4: Pressure contours for flow past a circular cylinder at Re=40: a) new Rhie-
Chow interpolation with
￿
t=0.001, b) original Rhie-Chow interpolation with
￿
t=0.001, c)
original Rhie-Chow interpolation with
￿
t=0.005, d) original Rhie-Chow interpolation
with
￿
t=0.01, from [4].
10 Wen Zhong Shen


Unfortunately, oscillations may appear when a small time step is used in unsteady
computations. In order to solve the problem, an improved Rhie-Chow interpolation
scheme was proposed in [4]. The improved interpolation scheme calculates the cell
face fluxes using the information obtained from the previous iteration. As an
illustration, pressure contours are plotted in Figure 4 for flow past a circular cylinder
with different time-steps. From the figure, it is seen that oscillations died out.
In order to couple the velocity and pressure more consistently, a SIMPLEC
(Consistent SIMPLE) scheme on collocated grid was developed in [5]. It was shown
in the paper that the original SIMPLEC scheme by Van Doormaal-Raithby 1984
cannot be used directly and is inconsistent for computations on collocated grids. A
remedy was proposed by adding a small term in the central terms on both sides when
computing cell face fluxes:
* * * * *
* * * * *
( ) ( )
( ) ( )
m m m m
P P P P P i i P P P
EWNS
m m m m
P P P P P i i P P P
EWNS
A u A u u Au A A u y p y p
A v A v v Av A A v x p x p
￿
￿ ￿ ￿
￿ ￿ ￿ ￿
￿ ￿
￿ ￿
￿ ￿ ￿ ￿ ￿ ￿ ￿
￿ ￿ ￿ ￿ ￿ ￿ ￿
￿
￿

where
*
(1 )/
P
u P u
A A￿ ￿￿ ￿ and ￿ is a parameter quantifying the amount added in the
interpolation. The remedy term
m
P
P
A u￿￿ can be interpreted in two meanings: (1) a
term related to time-like derivative
* *
P
P
A
u
; or (2) a term related to spatial
derivative
*
P
P
A u
. In the former case, the remedy term can be expressed as
* m
P P
A
u￿￿
￿

where
￿ ￿
/1
u u
￿
￿￿ ￿￿ ￿
￿
obtained from the definition of
*
P
A
.
For unsteady computations, the remedy can be introduced in two different ways:
m
t P
Au
￿
￿
￿
or
m
P
P
A
u￿￿
. Both ways were tested in [5]. Remark that the revised scheme
becomes the original SIMPLEC when 1
￿ ￿
￿
. It was also shown in the paper that the
solution obtained with the new scheme is time-step independent when a steady
solution is resulted from the unsteady solver.

Figure 5: Plot of number of iterations using the revised SIMPLEC for flow past a NACA
0012 at Re=3x10
6
and
￿
=4
o
, from [5].

Computational Aerodynamics and Aeroacoustics 11
As a basic test case, the flow past a circular cylinder at Re=40 was used for testing
both steady and unsteady versions of the revised SIMPLEC scheme. Results show
that the proposed scheme is consistent and can damp out oscillations successfully. It
has also been applied to turbulent flows past a NACA 0012 at a Reynolds number of
3x10
6
and an incidence of 4
o
. In Figure 5 the number of iterations is plotted. From the
figures, it can be seen that the revised SIMPLEC scheme retains the property of the
original scheme.
The revised SIMPLE and SIMPLEC schemes are implemented in the
EllipSys2D/3D code. The code will be used later for flow [8] [12] and acoustic
computations [16-19].

12 Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics 13
4. Actuator models [6-8]

The concept of an actuator disc used together with a Navier-Stokes solver was
introduced early by Sørensen-Myken 1992 and Sørensen-Kock 1995. Later a similar
Navier-Stokes based actuator disc model was also proposed by Madsen 1996. The
actuator disc/Navier-Stokes model was used to study various wake states in [6]. An
actuator line/Navier-Stokes model was proposed to study wind turbine wakes in [7].
An actuator surface/Navier-Stokes model was proposed to solve 2D flow past a
Vertical Axis Wind Turbine (VAWT) in [8] and 3D flow past a Horizontal Axis Wind
Turbine (HAWT) by Shen et al. 2007. The idea of introducing such actuator models is
to simplify the process of solving the boundary layers of the rotor and focus more on
the rotor generated wakes. One obvious limitation of the models is that the detailed
near wall boundary layer information cannot be obtained. Another limitation is that
the blade loading is not found directly by solving the flow near walls but calculated
with predetermined airfoil data that may be case-dependent.
The actuator disc / Navier-Stokes model [6]
In [6] the actuator disc/Navier-Stokes model was implemented in the finite difference
code based on vorticity-velocity formulation. Constant loading with different thrust
coefficients was considered. The wake states comprise wind turbine brake state,
turbulent wake state, vortex ring state, hover state and propeller state. It was shown
that the turbulent wake and vortex ring states are unstable regimes for a rotor with a
constant loading and that these states settle to a steady state after a complicated
transient phase. The results were compared to classic one-dimensional momentum
theory and experiments.
In Figure 6, the thrust coefficient is plotted against the axial interference factor.
From the figure, it is seen that the computed results follow the classic momentum
theory for thrust coefficient
Ct
<0.9 and for higher
Ct
they follow the Glauert
empirical relation (based on the data of helicopter experience), see Eggleston-
Stoddard 1987. In Figure 7, the dimensionless induced velocity
//2
v v T A
￿
￿

where
￿
,
T
and
A
are the air density, thrust and area covered by the wind turbine rotor
is plotted against the dimensionless undisturbed velocity
//2 4/
o o t
V V T A C￿￿ ￿
.
From the figure, it is seen that the results are in good agreements with the momentum
theory.

14 Wen Zhong Shen


Figure 6: Thrust coefficient as a function of axial interference factor, from [6].


Figure 7: Dimensionless induced velocity versus dimensionless undisturbed velocity,
from [6].
Computational Aerodynamics and Aeroacoustics 15

Figure 8: Computed vorticity field showing the formation of the wake structure at a
wind speed of 10 m/s, from [7].
The actuator line / Navier-Stokes model [7]
As an actuator disc model considers an averaged wake, this may have problems in
describing correctly a turbine wake and wake interactions between turbines in a wind
farm. In order to overcome these shortcomings, an actuator line /Navier-Stokes model
was introduced in [7] for studying three-dimensional flow field about wind turbines.
The model solves the three-dimensional Navier-Stokes equations with a loading
distributed along lines representing the blade force. The technique gives detailed
information about wind turbine wakes. The flow past a 500 kW Nordtank wind
turbine equipped with three LM19.1 blades was studied. In Figure 8, the absolute
vorticity is plotted for a 500 kW Nordtank at a wind speed of 10 m/s. From the figure,
the wake structure behind the wind turbine is clearly seen.
The actuator surface/Navier-Stokes model [8]
In order to capture the structure near rotor blade correctly, an actuator surface/Navier-
Stokes model has been introduced by Zhang 2004 and in [8] for two-dimensional flow
past a Vertical Axis Wind Turbines (VAWT). The method, a so-called 2D actuator
surface (AS) technique, consists of a two-dimensional Navier-Stokes solver in which
the pressure distribution is represented by body forces that are distributed along the
chord of the airfoils. The distribution of body force is determined from a set of pre-
defined functions that depend on angle of attack and airfoil shape. The pre-defined
functions are curve-fitted using pressure distributions obtained from a viscous-
inviscid interactive code (XFOIL). In Figure 9, streamlines and pressure are plotted
for flow past a NACA 0015 airfoil at a Reynolds number of 10
6
and an angle of attack
of 10
o
. From the figure, it is seen that the actuator surface model predicts the flow
field about an airfoil quite well.
16 Wen Zhong Shen


Figure 9: Streamline and pressure plot for flow past a NACA 0015 airfoil at Reynolds
number of 10
6
and angle of attack of 10
o
with (a) RANS; (b) Actuator Surface model,
from [8].

In [8], the actuator surface technique was also applied
to
compute the flow past a
two-bladed vertical axis wind turbine equipped with NACA 0012 airfoils. In order to
compare with experiments, the normal and tangential force coefficients based on wind
speed and airfoil chord have been computed. Figure 10 displays the tangential force
coefficient acting on the turbine blade as a function of the azimuth position. From the
figure, it is seen that the actuator surface method with the Leishman-Beddoes
dynamic stall model is in fairly good agreement with the experimental data while the
vortex method for the Darrieus Turbine with dynamical stall model (VDART2) of
Stickland et al. 1981 in
Paraschivoiu 2002
fails to reproduce the hysteresis of the
tangential force. A snapshot of the particle tracers at a dimensionless time of 34 is
shown in Figure 11. From the figure two main series of vortices are seen: one series
on the upper side of the wake where the vortices rotate clockwise and another series
on the lower side where the vortices rotate counter clockwise. Between these main
vortices a series of secondary vortices is observed.

Computational Aerodynamics and Aeroacoustics 17


Figure 10: Tangential force coefficient for flow past a VAWT at a tip-speed ratio of 2.5,
from [8].

The technique has been further developed for flow past a Horizontal Axis Wind
Turbine (HAWT) in Shen et al., 2007. The method, denoted as the actuator surface
technique, consists of a three-dimensional Navier-Stokes solver and a body force
distributed along the blade surfaces of a wind turbine. The force at each airfoil section
is obtained from tabulated lift and drag coefficients, which depend on the local angle
of attack and resulting velocity. The local force is further distributed by a set of pre-
determined functions that depends on angle of attack and airfoil shape. The actuator
surface technique was applied to compute the flow past a Nordtank 500kW turbine.
Detailed comparisons with the previously developed Actuator Line model show that
the new model determines the flow structures in the region close to the rotor blades
more correctly.

Figure 11: Particle tracers for flow past a Vertical Axis Wind Turbine at a tip-speed-
ratio
5.2￿
￿
, computed by the actuator surface technique, from [8].

18 Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics 19
5. Tip loss corrections [9-10]

To take into account the difference between the physics of an actuator disc with
infinitely many blades and an actual wind turbine or propeller with a finite number of
blades, Prandtl 1927 introduced the concept of tip loss and showed that the circulation
of a real rotor tends to zero exponentially when approaching the blade tip. Glauert
1963 used the Blade Element Momentum (BEM) theory as a computational tool to
predict aerodynamic loading and power (see also Hansen 2000). The theory is based
on one-dimensional momentum theory in which forces are distributed continuously in
the azimuth direction, corresponding to an infinite number of blades with no tip loss.
In order to make BEM computations more realistic, Glauert 1963 corrected the
induced velocity in the momentum equations by exploiting that the ratio between the
averaged induced velocity and the induced velocity at the blade position tends to zero
at the tip by the expression developed by Prandtl. A refined tip loss model was later
introduced by Wilson-Lissaman 1974 who suggested that the mass flow through the
rotor disc should be corrected in the same manner as the induced velocity in the wake.
This, however, leads to a formulation in which the orthogonality of the induced
velocity to the relative velocity at the blade element is not satisfied. In order to satisfy
the orthogonality condition, De Vries 1979 refined further the tip correction of
Wilson-Lissaman by correcting the mass flux in the tangential momentum equation in
the same way as in the axial momentum equation.

Figure 12: Comparison of normal forces computed by BEM code with the present tip
correction and with Glauert’s tip correction for flows past the NREL experimental rotor
at different wind speeds: (a) 7 m/s (
￿
=5.4); (b) 10 m/s (
￿
=3.9); (c) 13 m/s (
￿
=2.9); (d) 15
m/s (
￿
=2.5), from [9].

Theoretical analyses and comparisons with measurements were carried out in [9]
and showed that the existing tip loss corrections are inconsistent and fail to predict
20 Wen Zhong Shen

correctly the physical behaviour in the proximity of the tip. In order to overcome the
difficulty, a new tip correction was proposed in [9] where 2D force coefficients
should include further the tip loss effects at the blade tip (corrected 2D force
coefficients). This relation between uncorrected and corrected force coefficients was
proposed as
n
r
n
CFC
1
￿
,
t
r
t
CFC
1
￿
,
where the function F1 is
￿
￿
￿
￿
￿
￿ ￿
￿￿
￿
)
sin2
)(
exp(cos
2
1
1
R
R
rRB
gF
￿￿
,
￿ ￿
1.0)21(125.0exp ￿￿￿￿
￿
Bg
.

The F1 function was adjusted using the experimental data for flows past the NREL
rotor and the Swedish WG 500 rotor. A possible better fit can be obtained in the
future when more experimental data or CFD rotor computations for flow past rotors
with different rotor tips and rotor blades. A validation or comparison of the F1
function with the CFD computations in Hansen and Johansen 2004 for flows past the
Tellus 95 kW turbine shows good agreement (not shown here).
In order to show the performance, the tip loss correction models were used to
compute flows past the NREL experimental rotor and the Swedish WG 500 rotor. In
Figure 12, the new tip correction model is seen to predict correctly the forces on the
blade at wind speeds of small tip-speed-ratio (<5.4). For high tip-speed-ratio (>5.4),
the Swedish WG 500 rotor was used to test different correction models. The normal
force coefficients with the present tip loss correction and Glauert’s model were
calculated and compared to experimental data in Figure 13. From the figure, it is seen
that the present model correctly predicts the force changes near the tip.

Computational Aerodynamics and Aeroacoustics 21

Figure 13: Comparison of normal force coefficients computed by BEM code with the
new tip correction and Glauert’s tip correction models for flows past the Swedish WG
500 rotor at different tip speed ratios, (a) 5.44; (b) 8.24; (c) 14, from [9].

22 Wen Zhong Shen


Figure 14: Comparison of normal forces computed by the AD/NS and AL/NS models,
BEM with the new and classical tip corrections for flows past the NREL experimental
rotor at wind speeds of 10 and 15 m/s corresponding to tip speed ratios TSR of 3.79 and
2.53, respectively, from [10].

Since the actuator disc/Navier-Stokes models use the actuator disc principle, tip
loss corrections are necessary for such model (see Mikkelsen et al. 2001 and
Mikkelsen 2003). In the actuator line/Navier-Stokes model, the three-dimensional
Navier-Stokes equations are solved by introducing the loading from each of the blades
in the computations. From this point of view, it does not need any correction as used
in actuator disc models. However, 2D airfoil data are used everywhere on the blade
and the flow in the region near the blade tip does not look like 2D flow. In [10], the
actuator disc/Navier-Stokes and actuator line/Navier-Stokes models with tip loss
corrections were discussed. To illustrate the performance, the same test cases were
considered. In Figure 14, the normal forces for flows past the NREL experimental
rotor at wind speeds of 10 m/s and 15 m/s are compared. From the figure, it is seen
that the AD/NS works almost on the whole blade as the BEM with the new tip loss
correction. In the region near the tip, all three models work alright. It should be noted
that the AL/NS solves the three-dimensional Navier-Stokes equations with the loading
distributed along blade lines. Hence, it predicts the forces more correctly. In Figure 15,
the Swedish WG 500 rotor is considered again. From the figure, it is seen that the
AD/NS and AL/NS can work correctly with the new tip loss function,
F
1
.
Computational Aerodynamics and Aeroacoustics 23

Figure 15: Comparison of normal force coefficients computed by the AD/NS and AL/NS
models and BEM code with the new and classical tip corrections for flows past the
Swedish WG 500 rotor at tip speed ratio TSR of 5.44 and 8.24, from [10].

24 Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics 25
6. Miscellaneous aerodynamic problems [11-12]

In this section, two miscellaneous aerodynamic problems were treated. In [11], a
technique that determines Angle of attack (AOA) for rotating wind turbines was
proposed. In [12], a wall correction model for wall interference on wind turbine rotors
in wind tunnels with open test section was proposed.
Determination of AOA for a rotating wind turbine [11]
For a 2D airfoil the angle of attack (AOA) is defined as the geometrical angle
between the flow direction and the chord. The concept of angle of attack is widely
used in aero-elastic engineering models (i.e. FLEX by Øye 1996 and HAWC by
Petersen 1996) as an input to tabulated airfoil data that normally are established
through a combination of wind tunnel tests and corrections for the effect of Coriolis
and centrifugal forces in a rotating boundary layer. For a rotating blade the flow
passing by a blade section is bended due to the rotation of the rotor, and the local flow
field is influenced by the bound circulation on the blade. As a further complication,
3D effects from tip and root vortices render a precise definition of the AOA difficult.
In [11] two simple methods are proposed to compute correctly the angle of attack for
wind turbines in standard operations and in general flow conditions.

Method 1
consists of 6 steps:
Step 1: Determining initial flow angles using the measured velocity at every
monitor point
￿ ￿
Ni,1￿

￿ ￿
.
,/tan
2
,
2
,
0
,
,,
10
iizirel
iizi
VVV
VV
￿
￿
￿
￿￿
￿
￿

where
￿
￿
iiz
VV
,,
,
￿
are the velocity components in the axial and azimuth
directions, respectively and N is the number of cross sections.
Step 2: Estimating lift and drag forces using the previous angles of attack and the
measured local blade forces
￿
￿
iiz
FF
,,
,
￿
in the axial and azimuth directions,
￿ ￿
Ni,1￿

,,
,,
cos sin,
sin cos,
n n n
i z i i i i
n n n
i z i i i i
L F F
D F F
￿
￿
￿
￿
￿
￿
￿ ￿
￿ ￿

where n is the number of iteration.
Step 3: Computing associated circulations from the estimated lift forces as
n
irel
n
i
n
i
VL
,
/￿￿
.
Step 4: Computing the induced velocity created by the bound vortices using
circulations
￿ ￿
￿ ￿
￿ ￿
3
0
1
1
,,( )/
4
b
n
R
n n n n n
in r z j j j j
j
u x u u u y x y x y dr
￿
￿
￿
￿ ￿ ￿ ￿ ￿ ￿
￿
￿
￿
￿ ￿ ￿ ￿ ￿ ￿

where R and
b
n are the rotor radius and number of blades, respectively.
Step 5: Computing the new flow angle and the new velocity
￿
￿
n
i
n
iziz
uVuV
,,,
,
￿￿
￿￿

by subtracting the induced velocity generated by bound vortex
26 Wen Zhong Shen

￿
￿
￿
￿
￿ ￿ ￿ ￿
.
,tan
2
,,
2
,,
1
,
,,,,
11
n
ii
n
iziz
n
irel
n
ii
n
iziz
n
i
uVuVV
uVuV
￿￿
￿￿
￿
￿￿￿￿
￿￿￿
￿
￿￿

Step 6: If not converged, go to Step 2.
When the convergence is reached, the angle of attack and force coefficients can be
determined:
iii
￿
￿
￿
￿￿,
iirel
i
id
cV
D
C
2
,
,
2
￿
￿
,
iirel
i
il
cV
L
C
2
,
,
2
￿
￿
.
where
i
￿
is the sum of pitch and twist angles and
i
c is the chord.

Method 2
consists of the following 4 steps:
Step 1: Determine the edge velocity on the whole blade surface from the measured
pressure coefficient
￿ ￿
2
2
p
p
p
C
U￿
￿
￿
￿
￿ using the steady Bernoulli equation
(assumed to be true) at the edge of the boundary layer
2
2
1 1
2 2
p
u p U
￿
￿
￿ ￿
￿ ￿ ￿
￿

￿ ￿
2
2
1
p
p
U U
U
￿
￿
￿
￿
￿
￿
￿ ￿

Step 2: Determine the local bound circulation on the cross section as the velocity
jump over the boundary layer with the sign based on the edge velocity as
shown in Figure 2,
( ) ( 0)
r r
wall
x
U e U e
￿
￿ ￿￿ ￿ ￿ ￿
￿ ￿ ￿ ￿

Step 3: Compute the self induction at a monitor point in the cross-section from the
local bound circulation
￿ ￿ ￿ ￿
￿ ￿ ￿ ￿
3
1
1
,,
4
b
n
i
in r z
i
surface
y x y
u x u u u d dr
x y
￿
￿
￿
￿
￿
￿ ￿
￿ ￿
￿
￿
￿￿
￿ ￿ ￿ ￿
￿ ￿
￿ ￿

where
￿ ￿
,r￿
are the coordinates in the tangential and radial directions,
respectively, of the local coordinate system,
y
￿
,on the blade surface.
Step 4: Compute the flow angle in a section from the velocity measured at the
monitor point in which the self induction from the bound vortices is
subtracted
￿ ￿ ￿ ￿
￿ ￿ ￿ ￿
1
2 2
tan
z z
rel z z
V u V u
V V u V u
￿ ￿
￿ ￿
￿
￿
￿ ￿ ￿
￿ ￿ ￿ ￿

The Angle of Attack and force coefficients are then determined as
￿
￿
￿
￿￿
,
2
2
d
rel
D
C
V c￿
￿,
2
2
l
rel
L
C
V c￿
￿.
Computational Aerodynamics and Aeroacoustics 27


Figure 16: Drag and lift coefficients at r/R=65.7%, from [11].

The techniques were used to determine the AOA using the data of Navier-Stokes
computations on flows past the Tellus 95 kW wind turbine by Hansen-Johansen 2004.
As an illustration, the lift and drag coefficients at r/R=65.7% were plotted in Figure
16. From the figure, it is seen that our approach works as well as the Averaged
Technique (AT) by Johansen-Sørensen 2004 based on the averaged velocity fields at
upstream and downstream. It should be remarked that for unsteady flow the AT may
have difficulties to find the area to make averaging whereas the proposed technique
does not have such problem.
It is worth noting that in the second technique no iteration is required and the
monitor point can be chosen to be closer to the cross section. On the other hand, the
difficulty of using the second technique is to find the separation point where the local
circulation changes sign. The problem may be solved from the distribution of skin
friction when CFD solutions are used.
Wall correction model for wind tunnels with open test section [12]
It is known that experimental data measured in a wind tunnel are influenced greatly
by tunnel walls. This is also true for wind turbine measurements carried out in a wind
tunnel. For a wind turbine in a closed wind tunnel, Mikkelsen-Sørensen 2002 derived
a simple model based on a one-dimensional momentum approach to study the
blockage effect. For wind tunnels with open test section, a momentum change (loss or
increase) from the boundary between the tunnel and open test section was used in [12]
when considering the axial momentum conservation
￿
￿￿￿
CS
x
AdVVT
￿
￿
￿
~
.
Based on the results of an AD/NS model for flows past an actuator disc with
constant loading in a wind tunnel with open test section, the momentum loss can be
curve fitted with a simple function
b
T
o
Ca
SV
T
￿￿
￿
2
2/1
~
￿

where
03.0￿a
and
3￿b
. Using this expression, the velocity at the disc plane can be
calculated with a simple formula
28 Wen Zhong Shen

￿ ￿
￿ ￿
￿￿￿￿￿￿
￿￿￿￿
￿￿￿
￿
￿
223
~
2
u

where
1
1
￿
￿￿
b
T
Ca￿. Since this expression contains
T
C, the solution of the system of
equations can be found iteratively.

Figure 17: Ratio of corrected and uncorrected stream velocities as function of thrust
coefficient and ratio of tunnel to rotor radius in a solid wall wind tunnel (left) and a
wind tunnel with open test section, from [12].

In order to illustrate the model, the BEM with tunnel corrections are compared to
AD/NS computations in the same tunnel conditions in Figure 17. From the figure, it is
seen that the 1D BEM model can predict correctly the influence of wall blockage in
the MEXICO wind tunnel.

Computational Aerodynamics and Aeroacoustics 29
7. Development of Aeroacoustic formulations [13-15]

As aerodynamic problems are important for the technological development of wind
turbines, noise problems are more critical for further development of wind power in
the society. The noise problems come directly from the public acceptance for
installing wind turbines as noise has negative impacts on human daily life. Because of
the scale differences in flow and acoustics at low Mach number, classic compressible
solution techniques are difficult and computational expensive to employ for low-
speed flows. As a consequence, it is preferable to solve them separately using a
splitting technique. In this section, aero-acoustic formulations of the splitting
technique proposed originally by Hardin-Pope 1994 are discussed. In [13], a comment
on the original formulation was given and followed with an improved formulation. In
[14], more discussions about the improved aero-acoustic formulation were given and
some flow cases were considered. In [15], a general form of the aero-acoustic
formulation for both laminar and turbulent flows was given where a turbulence model
was introduced to cope with turbulent flow simulations.
Improvement on an existing Aero-Acoustic formulation [13-14]
For the original formulation given by Hardin-Pope 1994, several issues were
discussed:
1.

The variable for hydrodynamic density correction
1
￿
is not needed in the aero-
acoustic formulation and can be integrated together with the fluctuation
￿
￿.
2.

The original formulation does not contain any source term: i.e.
(
1
0,0,0
i
u p
￿
￿￿ ￿ ￿￿ ￿ ￿ ￿ ) is the trivial solution.
3.

It is not consistent when the flow becomes isentropic.

Based on these points, a new formulation was proposed in [13]. The compressible
variables are decomposed as
0
i i i
u U u
p
P p
￿
￿ ￿
￿
￿ ￿
￿
￿ ￿
￿￿ ￿

where
￿ ￿
0
,,
i
U P
￿
are the incompressible variables and
￿ ￿
,,
i
u p
￿
￿ ￿ ￿
are the acoustic
variables. The aero-acoustic formulation is
￿ ￿
￿ ￿
0
2
0
i
i
i
i j j ij i j
j j
f
t x
f
f
U u p U u
t x x
p P
c
t t t
￿
￿ ￿
￿
￿ ￿￿
￿ ￿
￿ ￿
￿ ￿ ￿
￿ ￿￿ ￿ ￿
￿ ￿ ￿ ￿ ￿
￿ ￿
￿ ￿ ￿
￿ ￿￿ ￿ ￿
￿ ￿ ￿
￿ ￿ ￿

where
i i i
f
u U
￿
￿
￿ ￿
￿ ￿ and
￿ ￿ ￿ ￿
2
0
/c P p
￿
￿ ￿
￿ ￿
￿ ￿ ￿
. From the acoustic equations
above, the terms
￿ ￿
0 i j
j
U u
x
￿
￿
￿
￿
￿
and
P
t
￿
￿
￿
are considered to be the source terms. The
first source term is related to the vortex movements while the second is related to the
changes in incompressible pressure.
30 Wen Zhong Shen

The formulation was implemented in [14] using a finite difference discretization
on an orthogonal mesh generated with a conformal mapping. The aero-acoustic code
uses the semi-implicit Adams-Bashforth/Crank Nicolson scheme at time level
￿ ￿
1/2n t￿￿
for time discretization and second order central difference scheme in
space, except the convective-deformation terms that are discretized by the second-
order upwind QUICK scheme. The Tam-Webb 1993 acoustic radiation conditions are
applied at the far-field boundary and convective boundary conditions are used at the
outlet boundary. The aero-acoustic code was used together with the vorticity-stream
function based code.

Figure 18: Contours of fluctuating density for an impulsively started cylinder at Re=200
and M=0.2 at t= 1, 2, 4, 6, 8 and 10, from [14].

Flows past a circular cylinder both for steady and unsteady cases were considered.
As an illustration, the fluctuating density for an impulsively started cylinder is plotted
in Figure 18. From the figure, it is seen that a strong starting acoustic wave propagates
out smoothly from the cylinder.

Aero-acoustic formulation for turbulent flow [15]
In order to deal with turbulent flows, the aero-acoustic formulation was further
developed. Since Direct Numerical Simulations (DNS) are too expensive for today’s
computers, computations using turbulence models for Reynolds Averaged Navier-
Stokes Simulations (RANS) or Large Eddy Simulations (LES) are often made for
industrial flows. Based on the fact that the flow induced acoustic noise has a very
Computational Aerodynamics and Aeroacoustics 31
limited influence back on the flow at low Mach number, the following assumptions
are made
1.

The turbulent fluctuation of an acoustic signal is negligible.
2.

There is no back interaction from acoustics to flow.
3.

The acoustic length scale is much bigger than the turbulent length scale.
These assumptions can avoid a further turbulence modelling for the acoustic signals
and the time average or filtering of an acoustic variable is equal to be itself. Thus, the
basic decomposition of the compressible solution can be written:
* *
* *
* *
0 0
,
,
.
i i i i i i
u U u U U u
p
P p P P p
￿ ￿ ￿ ￿ ￿
￿
￿ ￿ ￿ ￿ ￿
￿
￿ ￿ ￿ ￿ ￿
￿ ￿ ￿ ￿

where the bar denotes time-averaging. The final aero-acoustic formulation for
turbulent flows read
￿ ￿
*
* * * * *
0
* *
2
0
2
3
i
i
j
i i
i j j i j ij t
j j j i
f
t x
U
f U
f U u U u p k
t x x x x
p P
c
t t t
￿
￿ ￿ ￿ ￿￿
￿
￿￿
￿ ￿
￿ ￿
￿ ￿￿ ￿
￿
￿ ￿￿ ￿ ￿ ￿
￿ ￿
￿ ￿ ￿ ￿ ￿ ￿ ￿
￿ ￿￿ ￿
￿ ￿
￿ ￿
￿ ￿
￿ ￿ ￿ ￿ ￿
￿ ￿
￿ ￿ ￿ ￿
￿ ￿
￿ ￿
￿ ￿ ￿
￿ ￿ ￿
￿ ￿ ￿

where
￿ ￿
* * *
0i i i
f
u U￿ ￿ ￿￿ ￿ ￿ and
￿ ￿
￿ ￿
2 * *
0
/c P p
￿
￿ ￿￿ ￿ ￿.
As compared to the laminar acoustic formulation, the additional terms,
￿ ￿
*
2/3/
i
k x￿￿ ￿ and
*
/(//
j i j j i
x
U x U x￿
￿ ￿
￿ ￿ ￿ ￿ ￿￿ ￿
￿ ￿
, appear in the acoustic velocity
equations.

Figure 19: Instantaneous plot of fluctuating pressure for turbulent flow past a NACA
0015 airfoil at Re=1.5x10
6
, 20-deg incidence, M=0.2 and t=420, from [15].

The formulation was implemented using the similar technique as shown in the
previous section for laminar flows. The code was tested for turbulent flow past a
NACA 0015 airfoil at Re=1.5x10
6
, 20
o
incidence and M=0.2 using the Baldwin-Barth
32 Wen Zhong Shen

one-equation model. In Figure 19, the instantaneous fluctuating pressure at t=420 is
plotted. From the figure, it is seen that the noise is dominated by the Strouhal
frequency (tonal noise). The directivity pattern is plotted in Figure 20. From the figure,
it is seen that the noise propagates preferable in the normal direction of the airfoil.

Figure 20: Directivity pattern of NACA 0015 airfoil noise radiation at Re=1.5x10
6
, 20-
deg incidence, M=0.2 and r=12, from [15].

Computational Aerodynamics and Aeroacoustics 33
8. Second-order finite volume/aero-acoustic code with
applications [16-18]

In order to study aero-acoustics in more general cases, the aero-acoustic formulation
was implemented in the Navier-Stokes code EllipSys2D/3D. Since EllipSys2D/3D
uses a finite volume discretization on a collocated grid in multi-blocks, the same
methodology was used in [16] for discretizing the acoustic equations. The features
for the acoustic code are as follows:
1.

Semi-implicit Crank-Nicolson scheme in time (less dissipative than the second
order backward differentiation used in EllipSys2D/3D).
2.

It consists of predictor and corrector steps.
3.

Improved Rhie-Chow interpolation.
4.

The pressure correction equation for the acoustic pressure perturbation was
derived and used for correcting the acoustic pressure.
5.

The pressure correction equation is solved by a 5-level multi-grid solver.
In [17], it was found that adding the viscous term in the acoustic velocity equations
can stabilize the acoustic solution. This problem is directly a consequence of the
incompressible/acoustic splitting technique. The main reason is that an unsteady
incompressible solution does not provide correctly a density fluctuation which is a
part of an unsteady compressible flow solution. Thus the incompressible solution
(pressure and velocity) is not the exact solution of the compressible Navier-Stokes
equations even for flows at small Mach numbers. This inconsistency can destabilize
completely the acoustic equations and make the code divergent.
The flow past a circular cylinder and a NACA 0015 airfoil were used to test the
code. It was found that the solution is grid independent when the viscous term is
included. In the second part of the paper, airfoil-gust computations were carried out.
Short and long vortical gust waves were considered and comparisons to the inviscid
results of Lockard-Morris 1998 showed very good agreement.

Figure 21: Instantaneous plot of fluctuating pressure for flow past a NACA 0015 at
Re=300, M=0.2 and t=240, from [16].
34 Wen Zhong Shen

To illustrate the performance of the aero-acoustic formulation, the fluctuating
pressure is plotted for flow past a NACA 0015 airfoil in Figure 21. The fluctuating
velocity and pressure for a short vortical gust interacting with a NACA 0012 airfoil
are plotted in Figures 22 and 23, respectively.


Figure 22: Instantaneous
v v￿
plot for a NACA 0012 airfoil at k=7.85, M=0.5, Re=5000
and 0-deg incidence, from [16].


In order to simulate noise generated from a wind turbine, the aero-acoustic
equations in cylindrical coordinates was also implemented in the code. In [17], the
wind turbine was modelled by an actuator line technique that distributes loading along
blade lines. The computations were carried out at a wind speed of 15 m/s. Without
perturbation, the incompressible solution is almost steady and the acoustic waves
attenuate. In order to analyse how the acoustic waves are generated and propagated,
the rotor is oscillating in the axial direction at a frequency of twice the rotating
frequency and amplitude of 10% of the inflow wind speed.


Figure 23: Instantaneous
p
p￿
plot for a NACA 0012 airfoil at k=7.85, M=0.5, Re=5000
and 0-deg incidence, from [16].

Computational Aerodynamics and Aeroacoustics 35

Figure 24: Directivity pattern of an oscillating Nordtank 500 kW turbine equipped with
three LM19.1 blades at a wind speed of 15 m/s, Mach number of 0.1 and a distance of 6
radii, from [17].

The directivity pattern is plotted in Figure 24. From the figure, it is seen that the
noise propagates preferably in the downstream direction of about 70
o
around the rotor
axis. It can be also seen that it is less noisy in the area close to the rotor plane.
In order to study the noise generated from turbulent flows, the aero-acoustic
model was coupled to a Sub-Grid-Scale turbulence model for Large Eddy Simulations
[18]. The mixed model by Loc 1994 was used in the flow computations.


Figure 25: Comparison of lift and drag coefficients to the data from Sheldahl-Klimas
1981 for flow past a NACA 0015 at Re=160 000, from [18].

36 Wen Zhong Shen


Figure 26: Noise spectrum for flow past a NACA 0015 airfoil at a Reynolds number of
160 000, Mach number of 0.2, angle of attack of 8
o
and a distance of 12 chords from the
airfoil centre, from [18].

In order to check the performance of the turbulence model, the force coefficients
are plotted in Figure 25 and compared to the experimental data of Sheldahl-Klimas
1981. Good agreements were found before stall. There is a discrepancy in the stall
region that may be attributed to the double stall phenomenon. In Figure 26, the
acoustic noise spectrum was compared to the experimental data from Brook et al.
1989. Good agreements were found. A parametrical study of the noise pattern for
flows at angles of attack between 4
o
and 12
o
shows that the noise level is small for
angles of attack below 8
o
, increases sharply from 8
o
to 10
o
and reaches a maximum at
12
o
.
Computational Aerodynamics and Aeroacoustics 37
9. High order finite difference aero-acoustic code [19]

Since the difference between the smallest and biggest scales is quite large for
turbulent flows, it is required to develop a more efficient computational aero-acoustic
tool to treat such problems. The choice here is high-order finite difference schemes
such as the Dispersion Relation Preserving (DRP) schemes introduced by Tam-Webb
1993 and the compact schemes introduced by Lele 1992. A high-order formulation of
the aero-acoustic code based on the splitting aero-acoustic technique was developed
in [19]. The code uses finite difference discretizations and is implemented in the
EllipSys2D/3D code.
In [19], the aero-acoustic equations shown in [13-15] were re-written in a
conservative matrix form
Q E F G
S
t x y z
￿ ￿ ￿ ￿
￿ ￿ ￿ ￿
￿ ￿ ￿ ￿

The equation is solved using a 4
th
order Runge-Kutta method in time and 4
th
, 6
th
and
8
th
order DRP or optimized compact schemes in space. The boundary conditions are
the inviscid boundary condition at wall boundaries and radiation boundary condition
for outgoing acoustic waves at the far-field. Introducing a ghost cell inside the wall,
the boundary equations are discretized as follows:
￿

When explicit DRP schemes are used, the 4
th
order backward DRP schemes
are applied at the points near or at the boundaries, see Figure 27.
￿

When compact schemes (for example, the 6
th
optimized compact scheme) are
used, 4
th
order compact boundary scheme is used at the first layer and 5
th
order
at the second and third layers.
Since high order schemes are usually unstable for disturbances, high order filter
schemes proposed by Bogey-Bailly 2004 are needed to damp out oscillations. For
more details about the high-order aero-acoustic code, the reader is referred to Zhu
2007.

Figure 27: Seven point stencils used for interior, wall and far-filed points, from [19].

In order to determine the performance of the high-order aero-acoustic code, the
flow past a NACA 0012 airfoil at Re=100 000,
M
=0.2 and ￿=5
o
is considered. In
38 Wen Zhong Shen

Figure 28, the acoustic pressure is plotted. From the figure, it is seen that the acoustic
waves are propagating smoothly out from the airfoil.

Figure 28: 3D view of the acoustic pressure for flow past a NACA 0012 airfoil at Re=100
000, M=0.2 and
￿
=5
o
, from [19].
Computational Aerodynamics and Aeroacoustics
39

10. Concluding remarks

In this dissertation I have presented some of my main contributions within the field of
wind turbine aerodynamics, based of about thirteen years of research. Most of the
work was focussed on development and use of numerical techniques, i. e.
Computational Fluid Dynamics (CFD) and Computational Aero-Acoustics (CAA),
applied on wind turbine rotors. A part of work concerned improvements of classical
rotor aerodynamics. An example of this is the improvement of the 1D-Blade Element
Momentum (BEM) technique for predicting wind turbine performances, where a new
tip loss correction model using 2D airfoil data near the blade tips was developed in [9].
To improve the quality of the BEM technique, airfoil data need to be corrected more
precisely with respect to the rotational effects caused by the centrifugal and Coriolis
forces. The quasi-3D technique developed in [1] was developed for such purpose. If
detailed numerical and experimental data are available, airfoil data can be extracted
directly by using the technique of determination of the angle of attack developed in
[11]. In order to predict the performance of wind turbines located in wind farms, the
developed actuator techniques including the Actuator Disc / Navier-Stokes (AD/NS)
[6], Actuator Line / Navier-Stokes (AL/NS) [7] and Actuator Surface / Navier-Stokes
(AS/NS) [8] models can be employed. These models are much faster to use than full
CFD computations.
Concerning aero-acoustics of wind turbines, the noise issue is an important point
for further development of wind power in the society. To design low noise blades, the
semi-empirical code developed by Zhu et al. 2005 can be used. Since this code is
semi-empirical, there are many issues related to the applicability of the model for
general wind turbine flows that need to be checked. In order to precisely predict flow
induced noise, the natural choice is Computational Aero-Acoustics (CAA). The
developed CAA model is based on the splitting technique developed in [13-15] and
can be used for predicting noise from airfoils and wind turbines. Both second-order
finite volume methods [16] and high-order DRP/optimized compact finite difference
methods [19] have been developed and used for predictions. For turbulent flows,
aeroacoustic computations combing with flow solutions obtained from LES (Large
Eddy Simulations) or DES (Detached Eddy Simulations) are in [18] shown to be a
reliable way for noise predictions.
40 Wen Zhong Shen

Computational Aerodynamics and Aeroacoustics
41
References
Akselvoll K. and Moin P.: 1996, An efficient method for temporal integration of the
Navier-Stokes equations in confined axisymmetric geometries,
Journal of
Computational Physics

125
, 454-463.
Arakawa A.: 1966, Computational design for long-term numerical integration of the
equations of fluid motion: two-dimensional incompressible flow,
Journal of
Computational Physics

1
, 119-143.
Baldwin B. S. and Barth T. J.: 1990,
A one-equation turbulence transport model for
high Reynolds number wall-bounded flows
, NASA TM 202847.
Baldwin B. S. and Lomax H.: 1978,
Thin layer approximation and algebraic model
for separated turbulent flow
, AIAA paper 78-257.
Béchara W. Bailly C., Lafon P. and Candel S.: 1994, Stochastic approach to noise
modeling for free turbulent flows,
AIAA Journal

32
, 455-463.
Bogey C. and Bailly C.: 2004, A family of low dispersive and low dissipative explicit
schemes for noise computations,
Journal of Computational Physics

194
, 194-
214.
Bogey C. and Bailly C.: 2005, Effects of inflow conditions and forcing ob subsonic
jet flows and noise,
AIAA Journal

43
, 1000-1007.
Bogey C. and Bailly C.: 2006, Computation of a high Reynolds number jet and its
radiated noise using large eddy simulation based on explicit filtering,
Computers
& Fluids

35
, 1344-1358.
Bogey C., Bailly C.and Juvé D.: 2002, Computational of flow noise using source
terms in linearized Euler’s equations,
AIAA Journal

40
, 235-243.
Bogey C., Bailly C.and Juvé D.: 2003, Noise investigation of a high subsonic,
moderate Reynolds number jet using a compressible LES,
Theoretical and
Computational Fluid Dynamics

16
, 273-297.
Brooks T. F., Pope D. S. and Marcolini M. A.: 1989,
Airfoil self-noise and prediction
,
NASA Reference Publication
1218
, National Aeronautics and Space
Administration, USA.
Canuto C., Hussaini M. Y., Quarteroni A. and Zang T. A.: 1988,
Spectral Methods in
Fluid Dynamics
, Springer-Verlag, New York.
Chorin A. J.: 1973, Numerical study of slightly viscous flow,
Journal of Fluid
Mechanics

57
, 785-796.
Ciarlet P. G. and Lions J. L.: 1996,
Handbook of numerical analysis-2: Finite-element
methods
, North-Holland, Amsterdam.
Colonius T.: 1997,
Lectures on computational aeroacoustics
, Lectures Series 1997-07,
von Karman Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium.
Daube O. and Loc T. P.: 1978, Etude numérique d’écoulements instationnaires de
fluide visqueux incompressible autour de corps profilés par une méthode
combinée d’ordre O(h
4
) et O(h
2
),
Journal de Mécanique

17
, 651-678.
Daube O.: 1992, Resolution of the 2D Navier-Stokes equations in velocity-vorticity
form by means of an influence matrix technique,
Journal of Computational
Physics

103
, 402-414.
De Vries O.: 1979,
Fluid dynamic aspects of wind energy conversion
, AGARD
Report AG-243.
Dhatt G. and Touzot G.: 1984,
The finite element method displayed
, John Wiley &
Son.
Eggleston D. M. and Stoddard F. S.: 1987,
Wind turbine Engineering Design
, Van
Nostrand Reinhold, New York.
42 Wen Zhong Shen