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 (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 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.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 [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

## Comments 0

Log in to post a comment