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, DK2800 Kgs. Lyngby, Denmark
Phone + 45 4525 4300, Telefax + 45 4588 4325
Email: 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 9788789502878
Keywords: wind turbine; rotor aerodynamics; wakes; aeroacoustics; NavierStokes equations; vorticity
velocity formulation; SIMPLEC; computational fluid dynamics (CFD); computational aero
acoustics (CAA); flowacoustic 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 aeroacoustics
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 aeroacoustic 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 1DBlade 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 aeroacoustics 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 aeroacoustic 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 NavierStokes equations in primitive variables
(velocitypressure 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 NavierStokes equations may be formulated in
different ways, using a vorticitystream function formulation, a vorticityvelocity
formulation or a vorticitypotentialstream 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 velocitypressure coupling system in the
inhouse 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/NavierStokes 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 AeroAcoustics (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 NavierStokes 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 nonlinear
twostep (viscous incompressible flow and inviscid acoustic perturbation) splitting
procedure for computational aeroacoustics 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 aeroacoustic formulation for turbulent flows was in [15] developed for
Large Eddy Simulation (LES), Unsteady Reynolds Averaged NavierStokes
Simulation (URANS) and Detached Eddy Simulation (DES). In [16] a collocated grid
/ finite volume method for aeroacoustic computations was developed and
implemented in the EllipSys2D/3D code. In [17] and [18] three dimensional flow
acoustic computations were carried out. Finally, the aeroacoustic 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. Quasi3D NavierStokes model for a rotating airfoil.
Journal of Computational Physics, 150: 518548, 1999.
[2] W. Z. Shen and T. P. Loc. Numerical methods for unsteady 3D NavierStokes
equations in velocity and vorticity form. Computers & Fluids, 26: 193216, 1997.
[3] M. O. L. Hansen, J. N. Sørensen and W. Z. Shen. Vorticityvelocity formulation of
the 3D NavierStokes equations in cylindrical coordinates. International Journal for
Numerical Methods in Fluids, 41:2945, 2003.
[4] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. Improved RhieChow interpolation
for unsteady flow computations. AIAA Journal, 39: 24062409, 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: 221239, 2003.
[6] J. N. Sørensen, W. Z. Shen and X. Munduate. Analysis of wake states by a fullfield
actuator disc model. Wind Energy, 1: 7388, 1998.
[7] J. N. Sørensen and W. Z. Shen. Numerical modeling of wind turbine wakes. Journal
of Fluid Engineering, 124: 393399. 2002.
[8] W. Z. Shen, J. N. Sørensen and J. H. Zhang. The actuator surface model: A new
NavierStokes 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: 457475, 2005.
[10] W. Z. Shen, J. N. Sørensen and R. Mikkelsen. Tip loss correction for
actuator/NavierStokes computations. Journal of Solar Energy Engineering, 127:
209213, 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: 9198,
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: 18901894, 2006.
[13] W. Z. Shen and J. N. Sørensen. Comment on the aeroacoustic formulation of
Hardin and Pope. AIAA Journal, 37: 141143, 1999.
[14] W. Z. Shen and J. N. Sørensen. Aeroacoustic modelling of lowspeed flows.
Theoretical and Computational Fluid Dynamics, 13: 271289, 1999.
[15] W. Z. Shen and J. N. Sørensen. Aeroacoustic modeling of turbulent airfoil flows.
AIAA Journal, 39: 10571064, 2001.
[16] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. A collocated grid finite volume
method for aeroacoustic computations of lowspeed flows. Journal of Computational
Physics, 196: 348366, 2004.
[17] W. Z. Shen, J. A. Michelsen and J. N. Sørensen. Aeroacoustic computations of wind
turbines. AIAA paper, No. 20020043, Reno, Nevada, USA, 2002.
[18] W. Z. Shen, W. J. Zhu and Sørensen. Aeroacoustic computations for turbulent
airfoil flow. AIAA Journal, 47: 15181527, 2009.
[19] W. J. Zhu, W. Z. Shen and J. N. Sørensen. Computational AeroAcoustic using
highorder finite differences schemes. Journal of Physics: Conference Series 75, doi:
10.1088/17426596/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 [13].......................................................................................5
A quasi3D NavierStokes model for a rotating airfoil [1]....................................5
Vorticityvelocity formulations [23]....................................................................8
3. Velocitypressure coupling [45]...............................................................................9
4. Actuator models [68]..............................................................................................13
The actuator disc / NavierStokes model [6].......................................................13
The actuator line / NavierStokes model [7]........................................................15
The actuator surface/NavierStokes model [8]....................................................15
5. Tip loss corrections [910].......................................................................................19
6. Miscellaneous aerodynamic problems [1112]........................................................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 [1315]...............................................29
Improvement on an existing AeroAcoustic formulation [1314].......................29
Aeroacoustic formulation for turbulent flow [15]..............................................30
8. Secondorder finite volume/aeroacoustic code with applications [1618].............33
9. High order finite difference aeroacoustic 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, DaubeLoc 1978 and Tannehill et al. 1997), Finite Volume methods (for
example, Patankar 1980, FerzigerPeric 1996 and AkeselvollMoin 1996), Finite
Element methods (for example, CiarletLions 1996, DhattTouzot 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 spectralelement method of KarniadakisSherwin
2005 and the finite differencevortex method of Shen 1993, Guermond et al. 1993,
ShenLoc 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 NavierStokes equations. Since the relative velocity (U
) on wind
turbine blades is small when compared to the speed of sound (c340 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 NavierStokes equations obtained from the
compressible NavierStokes 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 (velocitypressure 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 vorticitystream function,
vorticityvelocity and vorticitypotentialstream function formulations.
In [13] two vorticity formulations were developed for 2D and 3D wind turbine
flows. In [45], some numerical techniques for avoiding pressure oscillations were
developed when solving the velocitypressure coupling system used in the inhouse
EllipSys2D/3D code developed at DTU by Michelsen 1992 and at RISØ by Sørensen
1995. In [68] the actuator techniques including actuator disc, actuator line and
actuator surface techniques were developed to simulate flows past one or more wind
turbines. In [910] a tip loss correction method that improves the classical correction
models was developed for tip losses for BEM or actuator/NavierStokes 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 AeroAcoustics (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
semiempirical model based on a series of wind tunnel measurements (for example,
Brooks et al. 1989, Lowson 1993, FuglsangMadsen 1996 and Zhu et al. 2005). A
semiempirical model can give an overall 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 WilliamsHawkings 1968, Lyrintzis 1994 and FarassatSucci
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 coworkers (Bogey, Bailly
and Juvé 2003, Bogey and Bailly 2005, 2006) studied the sound generation from a
turbulent jet by solving the compressible NavierStokes 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 nonlinear twostep (viscous incompressible flow and inviscid acoustic
perturbation) splitting procedure (SP) for computational aeroacoustics 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 coworkers (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 coworkers (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 sourceterm 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 [1314] 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 aeroacoustic formulation for turbulence flows was developed in [15] for Large
Eddy Simulations (LES), Unsteady Reynolds Averaged NavierStokes Simulations
(URANS) and Detached Eddy Simulations (DES). In [16] a collocated grid / finite
volume method for aeroacoustic computations was developed and implemented in
the inhouse EllipSys2D/3D code. In [1718] three dimensional flowacoustic
computations were carried out. Finally, the aeroacoustic 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 [13]
The concept of vorticity has been used to solve the incompressible NavierStokes
equations since the sixties of the last century. It was first formulated through the use
of the stream function to solve twodimensional (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, vorticityvelocity formulations were studied very
intensively by e. g. Daube 1992, Gatski 1991, StellaGuj 1989, GujStella 1993,
Hansen 1994, Quartapelle 1993, Speziale 1987 and in [23].
In order to study the rotating effects on a rotating wind turbine blade, a quasi3D
model in a rotating frame of reference was developed in [1] and solved using the
formulation of the vorticitystream function. In the paper [2], a general 3D vorticity
velocity formulation was developed and in [3] a 3D vorticityvelocity formulation
based on the CauchyRiemann equations was developed in cylindrical coordinates.
The developed vorticityvelocity based code was used to simulate flows past a wind
turbine in [67], with the turbine modelled by Actuator Disc and Actuator Line
techniques.
A quasi3D NavierStokes 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 quasi3D NavierStokes model was developed in [1]. The
hypotheses used in the model are based on Sears 1950 and FogartySears 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 quasi3D NavierStokes equations are obtained
and formulated into a vortictystream functionspanwise 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 semiimplicit, 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 OverRelaxation (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 secondorder 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 oneequation model of BaldwinBarth 1990 is used. From
the early European project VISCWIND (see Sørensen 1999, Sørensen and Nygreen
2001), different turbulence models including the BaldwinLomax 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 SpalartAllmaras 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.
Vorticityvelocity formulations [23]
In [2] several vorticityvelocity formulations were discussed and followed with a
suggested formulation. In the paper, it was proved that the new formulation is
equivalent to the NavierStokes equations in velocitypressure 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 velocitypressure 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 vorticityvelocity formulation based on solving the CauchyRiemann
equations for the velocity was proposed. The formulation was used to solve the 3D
NavierStokes 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 nonaxisymmetric 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. Velocitypressure coupling [45]
The EllipSys2D/3D code developed as a collaborative work between DTU and Risø,
Michelsen 1992 and Sørensen 1995, is based on the velocitypressure 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 velocitypressure
formulation. Therefore, a velocitypressure formulation needs also to be considered.
The EllipSys code uses cellcentred finitevolume 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 predictorcorrector method, Rhie
Chow 1983 proposed a procedure using a momentumbased 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 velocitypressure 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 RhieChow interpolation with
t=0.001, c)
original RhieChow interpolation with
t=0.005, d) original RhieChow 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 RhieChow 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 timesteps. 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 DoormaalRaithby 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 timelike 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 timestep 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 [1619].
12 Wen Zhong Shen
Computational Aerodynamics and Aeroacoustics 13
4. Actuator models [68]
The concept of an actuator disc used together with a NavierStokes solver was
introduced early by SørensenMyken 1992 and SørensenKock 1995. Later a similar
NavierStokes based actuator disc model was also proposed by Madsen 1996. The
actuator disc/NavierStokes model was used to study various wake states in [6]. An
actuator line/NavierStokes model was proposed to study wind turbine wakes in [7].
An actuator surface/NavierStokes 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 casedependent.
The actuator disc / NavierStokes model [6]
In [6] the actuator disc/NavierStokes model was implemented in the finite difference
code based on vorticityvelocity 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 onedimensional 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 / NavierStokes 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 /NavierStokes model
was introduced in [7] for studying threedimensional flow field about wind turbines.
The model solves the threedimensional NavierStokes 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/NavierStokes 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 twodimensional flow
past a Vertical Axis Wind Turbines (VAWT). The method, a socalled 2D actuator
surface (AS) technique, consists of a twodimensional NavierStokes 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 predefined
functions are curvefitted 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
twobladed 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 LeishmanBeddoes
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 tipspeed 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 threedimensional NavierStokes 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 tipspeed
ratio
5.2
, computed by the actuator surface technique, from [8].
18 Wen Zhong Shen
Computational Aerodynamics and Aeroacoustics 19
5. Tip loss corrections [910]
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 onedimensional 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 WilsonLissaman 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
WilsonLissaman 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 tipspeedratio (<5.4). For high tipspeedratio (>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/NavierStokes 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/NavierStokes model, the threedimensional
NavierStokes 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/NavierStokes and actuator line/NavierStokes 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 threedimensional NavierStokes 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 [1112]
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 aeroelastic 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 crosssection 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 NavierStokes
computations on flows past the Tellus 95 kW wind turbine by HansenJohansen 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 JohansenSø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, MikkelsenSørensen 2002 derived
a simple model based on a onedimensional 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.0a
and
3b
. 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 [1315]
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, aeroacoustic formulations of the splitting
technique proposed originally by HardinPope 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 aeroacoustic formulation were given and
some flow cases were considered. In [15], a general form of the aeroacoustic
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 AeroAcoustic formulation [1314]
For the original formulation given by HardinPope 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 aeroacoustic 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 aeroacoustic code
uses the semiimplicit AdamsBashforth/Crank Nicolson scheme at time level
1/2n t
for time discretization and second order central difference scheme in
space, except the convectivedeformation terms that are discretized by the second
order upwind QUICK scheme. The TamWebb 1993 acoustic radiation conditions are
applied at the farfield boundary and convective boundary conditions are used at the
outlet boundary. The aeroacoustic code was used together with the vorticitystream
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.
Aeroacoustic formulation for turbulent flow [15]
In order to deal with turbulent flows, the aeroacoustic 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 timeaveraging. The final aeroacoustic 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
, 20deg 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 BaldwinBarth
32 Wen Zhong Shen
oneequation 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. Secondorder finite volume/aeroacoustic code with
applications [1618]
In order to study aeroacoustics in more general cases, the aeroacoustic formulation
was implemented in the NavierStokes code EllipSys2D/3D. Since EllipSys2D/3D
uses a finite volume discretization on a collocated grid in multiblocks, the same
methodology was used in [16] for discretizing the acoustic equations. The features
for the acoustic code are as follows:
1.
Semiimplicit CrankNicolson 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 RhieChow 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 5level multigrid 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 NavierStokes
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, airfoilgust computations were carried out.
Short and long vortical gust waves were considered and comparisons to the inviscid
results of LockardMorris 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 aeroacoustic 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 0deg incidence, from [16].
In order to simulate noise generated from a wind turbine, the aeroacoustic
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 0deg 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 aeroacoustic
model was coupled to a SubGridScale 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 SheldahlKlimas
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 SheldahlKlimas
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 aeroacoustic 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 aeroacoustic
tool to treat such problems. The choice here is highorder finite difference schemes
such as the Dispersion Relation Preserving (DRP) schemes introduced by TamWebb
1993 and the compact schemes introduced by Lele 1992. A highorder formulation of
the aeroacoustic code based on the splitting aeroacoustic technique was developed
in [19]. The code uses finite difference discretizations and is implemented in the
EllipSys2D/3D code.
In [19], the aeroacoustic equations shown in [1315] were rewritten in a
conservative matrix form
Q E F G
S
t x y z
The equation is solved using a 4
th
order RungeKutta 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 farfield. 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 BogeyBailly 2004 are needed to damp out oscillations. For
more details about the highorder aeroacoustic code, the reader is referred to Zhu
2007.
Figure 27: Seven point stencils used for interior, wall and farfiled points, from [19].
In order to determine the performance of the highorder aeroacoustic 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 AeroAcoustics (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 1DBlade 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 quasi3D 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 / NavierStokes (AD/NS)
[6], Actuator Line / NavierStokes (AL/NS) [7] and Actuator Surface / NavierStokes
(AS/NS) [8] models can be employed. These models are much faster to use than full
CFD computations.
Concerning aeroacoustics 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
semiempirical code developed by Zhu et al. 2005 can be used. Since this code is
semiempirical, 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 AeroAcoustics (CAA). The
developed CAA model is based on the splitting technique developed in [1315] and
can be used for predicting noise from airfoils and wind turbines. Both secondorder
finite volume methods [16] and highorder 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
NavierStokes equations in confined axisymmetric geometries,
Journal of
Computational Physics
125
, 454463.
Arakawa A.: 1966, Computational design for longterm numerical integration of the
equations of fluid motion: twodimensional incompressible flow,
Journal of
Computational Physics
1
, 119143.
Baldwin B. S. and Barth T. J.: 1990,
A oneequation turbulence transport model for
high Reynolds number wallbounded flows
, NASA TM 202847.
Baldwin B. S. and Lomax H.: 1978,
Thin layer approximation and algebraic model
for separated turbulent flow
, AIAA paper 78257.
Béchara W. Bailly C., Lafon P. and Candel S.: 1994, Stochastic approach to noise
modeling for free turbulent flows,
AIAA Journal
32
, 455463.
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
, 10001007.
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
, 13441358.
Bogey C., Bailly C.and Juvé D.: 2002, Computational of flow noise using source
terms in linearized Euler’s equations,
AIAA Journal
40
, 235243.
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
, 273297.
Brooks T. F., Pope D. S. and Marcolini M. A.: 1989,
Airfoil selfnoise 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
, SpringerVerlag, New York.
Chorin A. J.: 1973, Numerical study of slightly viscous flow,
Journal of Fluid
Mechanics
57
, 785796.
Ciarlet P. G. and Lions J. L.: 1996,
Handbook of numerical analysis2: Finiteelement
methods
, NorthHolland, Amsterdam.
Colonius T.: 1997,
Lectures on computational aeroacoustics
, Lectures Series 199707,
von Karman Institute for Fluid Dynamics, RhodeSaintGenese, 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
, 651678.
Daube O.: 1992, Resolution of the 2D NavierStokes equations in velocityvorticity
form by means of an influence matrix technique,
Journal of Computational
Physics
103
, 402414.
De Vries O.: 1979,
Fluid dynamic aspects of wind energy conversion
, AGARD
Report AG243.
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
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment