Johnson Ben Final E-Thesis (Master Copy)

Computational Fluid Dynamics (CFD) Modelling
of Renewable Energy Turbine Wake Interactions
by
Benjamin Michael Carver Johnson
A thesis submitted in partial fulfilment for the requirements for the degree of
Doctor of Philosophy at the University of Central Lancashire
May 2015
STUDENT DECLARATION FORM
Concurrent registration for two or more academic awards
I declare that while registered as a candidate for the research degree, I have not been a registered
candidate or enrolled student for another award of the University or other academic or professional
institution
___________________________________________________________________________________
Material submitted for another award
I declare that no material contained in the thesis has been used in any other submission for an
academic award and is solely my own work
___________________________________________________________________________________
Signature of Candidate _______________________________________________________
Type of Award PhD
School School of CEPS

Contents
Abstract 17
Acknowledgements 19
1 Introduction 31
1.1 Background/motivation . . . . . . . . . . . . . . . . . . . . . . . . . 32
1.2 Aim and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
1.3 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4 Wind energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
1.4.1 Wind-farm . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
1.4.2 Wind turbine wake . . . . . . . . . . . . . . . . . . . . . . . 38
1.5 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . 38
1.5.1 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
1.5.2 Reynolds number . . . . . . . . . . . . . . . . . . . . . . . 41
1.5.3 Reynolds-Averaged Navier-Stokes . . . . . . . . . . . . . . 42
1.5.4 Large Eddy Simulations . . . . . . . . . . . . . . . . . . . . 43
2 Theory 45
2.1 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . 45
2.1.1 Reynolds-Averaged Navier-Stokes . . . . . . . . . . . . . . 46
3
2.1.3 Detached Eddy Simulation . . . . . . . . . . . . . . . . . . 61
2.1.4 Direct Numerical Simulations . . . . . . . . . . . . . . . . . 62
2.1.5 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . 63
2.2 Rotor models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.2.1 Porous disc . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.2.2 Actuator disc method . . . . . . . . . . . . . . . . . . . . . 70
2.2.3 Actuator line method . . . . . . . . . . . . . . . . . . . . . . 76
3 Literature review 80
3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.1.1 One-dimensional momentum theory . . . . . . . . . . . . . . 80
3.1.2 Betz limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.1.3 Ideal Horizontal Axis Wind Turbine (HAWT) . . . . . . . . 83
3.2 Previous experimental work . . . . . . . . . . . . . . . . . . . . . . 85
3.2.1 NREL Unsteady Aerodynamics Experiment (UAE) project . 85
3.2.2 Model EXperiments In COntrolled conditions (MEXICO) project 87
3.2.3 Mexnext project . . . . . . . . . . . . . . . . . . . . . . . . 88
3.3 Previous numerical work . . . . . . . . . . . . . . . . . . . . . . . . 88
3.3.1 Blade models . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.3.2 Vortex wake methods . . . . . . . . . . . . . . . . . . . . . . 89
3.3.3 Wake models and computer programs . . . . . . . . . . . . . 90
3.4 Previous computational work . . . . . . . . . . . . . . . . . . . . . . 91
3.4.1 Software/codes . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.4.2 Reynolds-Averaged Navier-Stokes . . . . . . . . . . . . . . . 93
3.4.3 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3.4.5 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . 96
4
3.5 Rotor modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
3.5.1 Actuator disc method . . . . . . . . . . . . . . . . . . . . . 98
3.5.2 Actuator line method . . . . . . . . . . . . . . . . . . . . . . 99
3.5.3 Actuator surface . . . . . . . . . . . . . . . . . . . . . . . . 100
3.5.4 Full Rotor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4 Methods 103
4.1 One-dimensional momentum theory study . . . . . . . . . . . . . . . 103
4.2 Porous disc study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.2.1 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.2.3 Resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.2.4 Code comparison . . . . . . . . . . . . . . . . . . . . . . . 110
4.2.5 Cell type study . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.3 Actuator disc study . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.3.1 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.4 Two turbine study . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.4.1 Current work . . . . . . . . . . . . . . . . . . . . . . . . . . 118
4.4.2 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.4.4 Repeatability study . . . . . . . . . . . . . . . . . . . . . . 122
4.5 Rotor diameter study . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5 Results 126
5.1 One-dimensional momentum theory results . . . . . . . . . . . . . . 126
5.2 Porous disc Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.2.1 Influence of the boundary conditions . . . . . . . . . . . . . . 128
5
5.2.2 Wake predictions . . . . . . . . . . . . . . . . . . . . . . . . 128
5.2.3 Cell type study results . . . . . . . . . . . . . . . . . . . . . 132
5.2.4 Code comparison results . . . . . . . . . . . . . . . . . . . . 135
5.3 Actuator disc study results . . . . . . . . . . . . . . . . . . . . . . . 144
5.4 Two turbine study results . . . . . . . . . . . . . . . . . . . . . . . . 149
5.4.1 Repeatability study results . . . . . . . . . . . . . . . . . . . 150
5.5 Rotor diameter study results . . . . . . . . . . . . . . . . . . . . . . 151
5.5.1 First turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
5.5.2 60m spacing wake . . . . . . . . . . . . . . . . . . . . . . . 154
5.5.3 60m spacing second turbine velocity . . . . . . . . . . . . . 155
5.5.4 60m spacing second turbine power . . . . . . . . . . . . . . 158
5.5.5 30m spacing wake . . . . . . . . . . . . . . . . . . . . . . . 162
5.5.6 30m spacing second turbine velocity . . . . . . . . . . . . . 164
5.5.7 30m spacing second turbine power . . . . . . . . . . . . . . 165
5.5.8 Power comparison between spacings . . . . . . . . . . . . . 169
6 Discussion 171
6.1 Computational techniques . . . . . . . . . . . . . . . . . . . . . . . 171
6.1.1 RANS and LES . . . . . . . . . . . . . . . . . . . . . . . . . 172
6.1.2 Cell type . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
6.1.4 ANSYS-CFX and ANSYS-Fluent . . . . . . . . . . . . . . . 175
6.2 Rotor models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
6.2.1 Near wake discrepancy . . . . . . . . . . . . . . . . . . . . . 176
6.2.2 Actuator disc . . . . . . . . . . . . . . . . . . . . . . . . . . 178
6.2.3 Actuator line . . . . . . . . . . . . . . . . . . . . . . . . . . 179
6.2.4 Ground affect and shear . . . . . . . . . . . . . . . . . . . . 181
6
6.3 Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
6.3.1 Reynolds number . . . . . . . . . . . . . . . . . . . . . . . 183
6.4 Rotor diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
6.4.1 Wake . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
6.4.2 Power . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
6.4.3 Power Coefficient . . . . . . . . . . . . . . . . . . . . . . . . 189
6.4.4 Overall performance . . . . . . . . . . . . . . . . . . . . . . 194
6.4.5 Wind-farm . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
7 Conclusions & Recommendations 196

7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198
7.1.1 Porous disk in water . . . . . . . . . . . . . . . . . . . . . . 198
7.1.2 Influence of the boundary conditions . . . . . . . . . . . . . . 199
7.1.3 Actuator disc in air . . . . . . . . . . . . . . . . . . . . . . . 199
7.1.4 Two turbine interactions . . . . . . . . . . . . . . . . . . . . 200
7.1.5 Influence of rotor diameter . . . . . . . . . . . . . . . . . . . 201
7.2 Recommended work . . . . . . . . . . . . . . . . . . . . . . . . . . 202
References 204
A Reynolds number 226
7
List of Tables
4.1 Mesh densities using linear and quadratic cells . . . . . . . . . . . . 111

4.2 Various mesh densities and tip refinement. . . . . . . . . . . . . . . . 121
4.3 Computational time for various mesh densities . . . . . . . . . . . . 122
5.1 Thrust coefficient results to 3 s.f. . . . . . . . . . . . . . . . . . . . 130

5.2 Mesh densities using linear and quadratic cell types . . . . . . . . . . 133
5.3 Mesh variation in the repeatability study . . . . . . . . . . . . . . . . 150
5.4 Velocity at the first turbine to 3 s.f. . . . . . . . . . . . . . . . . . . 154
5.5 Approximate power extracted by the first turbine to 3 s.f. . . . . . . . 154
5.6 Power difference of the first turbines compared with the baseline to 3
s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
5.7 Average velocity at the second turbine to 3 s.f. . . . . . . . . . . . . 159
5.8 Velocity difference at the second turbine compared to the baseline to 3
s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
5.9 Difference in power extracted by the second turbines compared to the
baseline to 3 s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
5.10 Comparison in the combined power extracted by the 7.5m and 12.5m
configurations with the baseline configuration to 3 s.f. . . . . . . . . 161
5.11 Power deficit of the second turbine to 3 s.f. . . . . . . . . . . . . . . 162
5.12 Average velocity at the second turbine with a spacing of 30m to 3 s.f. 165
8
5.13 Velocity difference at the second turbine compared to that behind the
baseline to 3 s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
5.14 Difference in power of the second turbines compared to 10m turbine
to 3 s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
configurations with the baseline configuration to 3 s.f. . . . . . . . . 169
5.16 Power deficit of second turbine compared to the first turbine to 3 s.f. . 169
5.17 Comparison of the power extracted by the turbines with a spacing of
60m and 30m to 3 s.f. . . . . . . . . . . . . . . . . . . . . . . . . . 170
6.1 New CT values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174

6.2 Average velocity difference between this study and Réthoré et al. . . 179
6.3 Power of 9m configuration . . . . . . . . . . . . . . . . . . . . . . . 189
6.4 New power coefficient values to 3 s.f. . . . . . . . . . . . . . . . . . 191
6.5 Approximate power extract by the first turbine with a variable CP to 3
s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
6.6 Power difference of the first turbines compared with the baseline using
a variable CP to 3 s.f. . . . . . . . . . . . . . . . . . . . . . . . . . . 191
6.7 CP values calculated for the second turbine referred to as CP2 . . . . . 192
6.8 Approximate power extract by the second turbine for different CP val-
ues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192
6.9 Difference in power of the second turbines compared to 10m turbine . 193
6.10 Power deficit of the second turbine . . . . . . . . . . . . . . . . . . . 193
6.11 Combined power of the first and second turbines for two variable CP
values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
9
configurations with the power extracted by the 10m configuration for
two variable CP values . . . . . . . . . . . . . . . . . . . . . . . . . 194
10
List of Figures
1.1 Types of wind turbine . . . . . . . . . . . . . . . . . . . . . . . . . . 36

1.2 Illustration of the resolved and modelled parts of the flow . . . . . . 40
2.1 Difference in grid resolution using a wall function . . . . . . . . . . 57

2.2 Control volumes of (a) the cell-centred scheme, (b) the cell-vertex
scheme (overlapping) and (c) the cell-vertex scheme (dual control vol-
umes) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.3 Geometry of two methods used to define disc location a) algebraically
b) geometrically . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
2.4 Aerofoil cross-section showing force vectors . . . . . . . . . . . . . 71
3.1 1D momentum theory . . . . . . . . . . . . . . . . . . . . . . . . . . 82

3.2 Aerofoil forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.1 1D momentum theory . . . . . . . . . . . . . . . . . . . . . . . . . 104

4.2 Normalized velocity profiles showing different mesh densities (a) along
the centre line and (b) 7 diameters downstream of the disc . . . . . . 107
4.3 Normalized velocity at the inlet of this study (solid line) and experi-
mental data (o) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.4 Difference in cell type . . . . . . . . . . . . . . . . . . . . . . . . . 111
11
4.5 Representation of the blade used in the MEXICO project where the red
zones represent the transition zones between aerofoils . . . . . . . . 112
4.6 Mexico study wind tunnel . . . . . . . . . . . . . . . . . . . . . . . 113
4.7 Mesh used in the actuator disc study (a) cross-section through the do-
main and (b) detailed view . . . . . . . . . . . . . . . . . . . . . . . 114
4.8 Centre line velocity for the 10m/s inlet with the 15 × 106 cell mesh
(dashed line), 30 × 106 cell mesh (solid line), numerical study (o) and
experimental data (.) . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.9 Two turbine domain with 3 diameter spacing on top and 6 diameter
spacing on the bottom . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.10 Mesh refinement zones . . . . . . . . . . . . . . . . . . . . . . . . . 120
4.11 Mesh comparison (a) 70m behind the first turbine and (b) detail view
of the velocity deficit . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.12 Velocity difference between the 22 million and 37 million cell mesh
densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.1 Pressure and velocity profiles along the centre line of the rotor given
by (a) the 1D momentum theory and (b) the actuator disc method . . 127
5.2 Velocity profile with no disc of the channel flow (solid line) and the
duct flow (dashed line) at 24m, 31m and 40m from the inlet . . . . . 128
5.3 Velocity along the centre line showing the channel (solid line), duct
(dashed line), numerical (x) and experimental data (o) for the CT values
of a) 0.61, b) 0.86 and c) 0.97 for the experimental study. The inlet 1
simulations are shown black and inlet 2 simulations are in red . . . . . 130
12
5.4 Turbulence intensity along the centre line showing the channel (solid
line), duct (dashed line), numerical (x) and experimental data (o) for
the CT values of a) 0.61, b) 0.86 and c) 0.97 for the experimental study.
The inlet 1 simulations are shown black and inlet 2 simulations are in
red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.5 Normalized velocity of the channel (solid line), duct (dashed line) and
experimental data (o) at different CT values of a) 0.61, 0.62, b) 0.86,
0.91 and c) 0.97, 0.99 for this study and the experimental study respec-
tively. The inlet 1 simulations are shown black and inlet 2 simulations
are in red. Each figure represents the distance downstream of the disc
corresponding to 8R, 14R, 22R, 30R and 40R respectively . . . . . . 133
5.6 Turbulence intensity at 8R, 14R, 22R, 30R and 40R downstream of the
disc, with the channel (solid line), duct (dashed line) and experimental
data (o) at different CT values of a) 0.61, 0.62, b) 0.86, 0.91 and c)
0.97, 0.99 for this study and the experimental study respectively. The
inlet 1 simulations are shown black and inlet 2 simulations are in red . 134
5.7 Velocity profile 7 diameters behind the disc of the channel simulations
with the quadratic cells (solid line) and linear cells (dashed line) . . . 136
5.8 Velocity profile 7 diameters behind the disc of the duct simulations
with the quadratic cells (solid line) and linear cells (dashed line) . . . 137
5.9 Normalized turbulence intensity behind the disc . . . . . . . . . . . . 138
5.10 Normalized velocity profile through the centre of the domain with
ANSYS-CFX (red), ANSYS-Fluent (blue) and experimental data (o)
for the channel (solid line) and duct (dashed line) simulations . . . . 138
x
5.11 Normalized velocity profile behind the disc at 2R = 4 with ANSYS-
CFX (red), ANSYS-Fluent (blue) and experimental data (o) for the
channel (solid line) and duct (dashed line) simulations . . . . . . . . . 140
13
5.12 Normalized velocity profile behind the disc with ANSYS-CFX (red),
ANSYS-Fluent (blue) and experimental data (o) for the channel (solid
line) and duct (dashed line) simulations . . . . . . . . . . . . . . . . 141
5.13 Turbulence intensity through the centre of the domain with ANSYS-
CFX (red), ANSYS-Fluent (blue) and experimental data (o) for the
channel (solid line) and duct (dashed line) simulations . . . . . . . . 142
x
5.14 Turbulent intensity behind the disc at 2R = 4 with ANSYS-CFX (red),
ANSYS-Fluent (blue) and experimental data (o) for the channel (solid
line) and duct (dashed line) simulations . . . . . . . . . . . . . . . . 142
5.15 Turbulence intensity behind the disc with ANSYS-CFX (red), ANSYS-
Fluent (blue) and experimental data (o) for the channel (solid line) and
duct (dashed line) simulations . . . . . . . . . . . . . . . . . . . . . 143
5.16 Velocity profile through the centre of the domain at 1.84797m or z/R =
0.82132 for this study (solid line), experimental data (.) and previous
numerical study (o) for three inlet velocities . . . . . . . . . . . . . . 145
5.17 Velocity variation towards the inlet for this study (solid line), experi-
mental data (.) and previous numerical study (o) . . . . . . . . . . . 146
5.18 Normalized velocity profile at x/2R = −0.4 with this study (blue) and
previous numerical study (red o) . . . . . . . . . . . . . . . . . . . . 147
5.19 Normalized velocity profiles of this study (blue) and previous numeri-
cal study (red o) at various locations . . . . . . . . . . . . . . . . . . 148
5.20 Normalized velocity profiles at 50m, 60m and 70m downstream re-
spectively of the first turbine with this study (solid line), previous nu-
merical data at 5% (o) and 10% (x) ambient turbulence and experimen-
tal data (*) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5.21 Velocity difference along yz plane between two simulations used in the
repeatability study . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
14
5.22 Detailed view of the velocity difference . . . . . . . . . . . . . . . . 151
5.23 Velocity profiles at the first turbine for each wind speed of 7m/s, 10m/s,
15m/s and 20m/s respectively and the different colours represent the
different turbine configurations; 7.5m turbine (blue), 10m turbine (green)
and 12.5m turbine (red) . . . . . . . . . . . . . . . . . . . . . . . . 153
5.24 Velocity at 70m behind the first turbine for each wind speed of 7m/s
10m/s, 15m/s and 20m/s respectively and the different colours rep-
resent the different turbine configurations; 7.5m turbine (blue), 10m
turbine (green) and 12.5m turbine (red) . . . . . . . . . . . . . . . . 156
5.25 Velocity profile through the centre of the domain offset by 4m with
velocity behind the 7.5m (blue), 10m (green) and 12.5m (red) turbines 157
5.27 Velocity profiles at the second turbine for each wind speed of 7m/s
10m/s, 15m/s and 20m/s respectively and the different colours rep-
resent the different turbine configurations; 7.5m turbine (blue), 10m
turbine (green) and 12.5m turbine (red) . . . . . . . . . . . . . . . . 158
5.28 Power extracted by the second turbine (7.5m turbine (blue), 10m tur-
bine (green) and 12.5m turbine (red)) at different wind speeds . . . . 160
5.29 Combined power extracted by the first and second turbines at different
wind speeds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
5.30 Velocity profiles 70m behind the first turbine for each wind speed of
7m/s 10m/s, 15m/s and 20m/s respectively with the different colours
representing the different turbine configurations; 7.5m turbine (blue),
10m turbine (green) and 12.5m turbine (red) . . . . . . . . . . . . . 163
15
5.33 Velocity profiles at the second turbine for each wind speed of 7m/s
10m/s, 15m/s and 20m/s respectively and the different colours repre-
sent the different turbine behind the 7.5m turbine (blue), 10m turbine
(green) and 12.5m turbine (red) . . . . . . . . . . . . . . . . . . . . 166
5.34 Power extracted by the second turbine (7.5m turbine (blue), 10m tur-
bine (green) and 12.5m turbine (red)) at different wind speeds . . . . 167
5.35 Combination of the power extracted by the first and second turbines at
different wind speeds . . . . . . . . . . . . . . . . . . . . . . . . . . 168
6.1 Force distribution of (a) actuator disc and (b) actuator line methods . 180
6.2 Isosurface of the vorticity (a) actuator disc and (b) actuator line meth-
ods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
6.3 Log relationship between the Reynolds number and rotor diameter . . 184
6.4 Velocity contours on the YZ plane for the (a) 7.5m configuration, (b)
10m configuration and the (c) 12.5m configuration . . . . . . . . . . 186
6.5 Streamlines viewed along (a) the yz−plane and (b) the xy−plane . . . 187
6.6 Power of each turbine at different wind speeds. Blue is the 7.5m tur-
bine, green is the 10m turbine, red is the 12.5m turbine and light blue
is the 9m turbine. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
16
Abstract
This thesis presents Computational Fluid Dynamics (CFD) simulations of renewable

turbines akin to those used for wind, hydro, and tidal applications. In line with sim-
ilar studies in the literature and to significantly reduce computational expense, the
models developed took the form of actuator discs with the solution of incompressible
Reynolds-Averaged Navier-Stokes equations with the k-ω SST turbulence models.
Simulations were initially conducted of a single turbine in water and air and then
two axially aligned turbines to study the flow field interactions. These models were
compared with previous theoretical, experimental and numerical data evident in the
literature. Generally, good agreement was found between these models and other anal-
ogous data sources in terms of velocity profiles in the far wake. The actuator disc
method was underpinned using the transient actuator line method, which showed good
agreement from a quantitative and qualitative viewpoint. However, it required signifi-
cant additional computational time when compared to the actuator disc method.
Each of the models were developed and solved using complimentary commercially
available CFD codes, ANSYS-CFX and ANSYS-Fluent. For this type of study, a crit-
ical evaluation of these codes was in all probability performed for the first time, where
it is shown that for the studies investigated in this thesis ANSYS-CFX performed bet-
ter than ANSYS-Fluent with respect to the computational effort (i.e. time and lines of
code).
17
The novel contributions to knowledge of the work described in this thesis is in
regard to the aerodynamic interactions between two turbines with different rotor diam-
eters and its effect on performance. It was observed that the rotor diameter of the first
turbine had a significant influence on the downstream turbines performance. Varying
the power coefficient exaggerates the differences between the first turbines but predicts
similar power extraction from the second turbines. Similar predictions with respect to
power extraction are made if the power coefficient of the second turbine is defined us-
ing the free stream velocity or the local free stream velocity even though there is less
power available for the second turbine to extract.
The work carried out as part of this thesis has been presented at two conferences
[1][2] and published within peer-reviewed journals [3][4].
18
Acknowledgements
I would like to thank my fiancée for supporting me in every way possible and moving
around the country for me to take this opportunity. My family for their continued
support, particularly my granddad for proof reading and providing feedback on my
multiple drafts.
I would like to thank my supervisors, Drs Justin Whitty and Jonathan Francis for
guiding and shaping me into the researcher I am today. I am also grateful for the other
researcher at the University who make my time far more enjoyable.
Lastly, I would like to acknowledge the Samuel Lindow Foundation who provided
me with a bursary and the opportunity to carry out this research.
19
Nomenclature
Latin
a Axial induction factor [−]
a0 Angular or tangential induction factor [−]
a1 Constant (fraction of the rate of turbulence production) [−]
A Boundary layer model constant [−]
Ar Rotor swept area [m2 ]
b Distance between a grid point and source point [m]
B Number of blades [−]
c Chord length [m]
cp Specific heat capacity [Nm/K]
C Closure coefficient [−]
C1,2 Prescribed matrix [−]
CD Drag coefficient [−]
Ci jk Turbulent transport tensor [−]
20
21
CL Lift coefficient [−]
Cp Power coefficient [−]
CS Smagorinsky constant [−]
CT Thrust coefficient [−]
CDkw Cross diffusion term [−]
D Drag force [N]
DAD Actuator Disc [N/m3 ]
DADUT Actuator Disc with Uniform Thrust definition term [N/m2 ]
DADV T Actuator Disc with Variable Thrust definition term [N/m]
DADL Actuator Disc with Lift definition term [N/m]
DADLD Actuator Disc with Lift and Drag definition term [N/m]
DAL Actuator Line definition term [N/m]
Dn Disc definition term [N/m]
d Rotor thickness/depth [m]
eL Lift unit vector [−]
eD Drag unit vector [−]
f Base term used in the definition of the actuator methods [N/m2 ]
F1 ,F2 Blending function [−]
FP Prandtl’s correction factor [−]

22
FG Glauert’s correction factor [−]
FS Shen et al. correction factor [−]
G Filter function [−]
H(x) Step function [−]
i, j Represents components in the i, jth direction [−]
k Specific turbulence kinetic energy [m2 /s2 ]
k Thermal diffusivity [m2 /s]
K Resistance coefficient [−]
l Characteristic length [m]
L Lift force [N]
Ln Actuator line area [m2 ]
LS Mixing length for the SGS [m]
ṁ Mass flow-rate [m/s]
Ma Mach number [−]
n Integer number [−]
N The blade number [−]
p Pressure [N/m2 ]
P Power [Nm/s]
PW F Wind-farm power [Nm/s]

23
PW T Wind turbine power [Nm/s]
r Radius (variable) [m]
r, θ, z Cylindrical coordinates [m, rad, m]
R Radius [m]
ℜ Function defining the volume of the rotor [−]
RR Function defining the rotor area [−]
Rr Variable rotor radius [m]
RT Function defining the rotor thickness [−]
S Mean strain tensor [N/m2 ]
S Invariant measure of the strain rate [1/s]
Si Source term in ith direction [various]
S¯i j Rate-of-strain tensor [1/s]
Sr,θ,x Source terms in cylindrical coordinates [various]
t Time [s]
T Thrust [N]
ϒ Temperature [K]
u, v, w Velocity components in the x, y, z directions [m/s]
ū, v̄, w̄ Mean velocity components in the x, y, z directions [m/s]
ūi Mean velocity component [m/s]

24
u0 , v0 , w0 Velocity fluctuation components in the x, y, z directions [m/s]
Ū Mean velocity [m/s]
U Velocity [m/s]
U0 Free stream Velocity [m/s]
Ui Velocity in i = x, y or z direction [m/s]
Ur Velocity at the rotor [m/s]
Ut Velocity at the turbine [m/s]
U∗ Friction velocity [m/s]
Ūi Large scale filtered component [−]
Uinlet Inlet velocity [m/s]
Urel Relative velocity [m/s]
V Volume of the computational cell [−]
→
−
V Velocity vector [m/s]
Vr,θ,x Velocity in cylindrical coordinates [m/s]
Vz Axial velocity [m/s]
Vθ Tangential velocity [m/s]
x, y, z Cartesian coordinate system [−]
xi, j Cartesian coordinate system [−]
(X,Y, Z)R Local coordinate system of the rotor [m]

25
yd Distance to the closest wall [m]
ys Distance to the nearest surface [m]
yw Water depth [m]
Greek
αω , βω ω ω ω
1 , β2 , σ1 , σ2 Closure coefficients of k-ω turbulence model [−]
αε , βε1 , σε1 , σε2 ,Cε Closure coefficients of k-ε turbulence model [−]
αsst , βsst sst sst sst

1 , β2 , σ1 , σ2 Closure coefficients of k-ω SST turbulence model [−]
α Local angle of attack [◦ ]
βi Initial blade position [◦ ]
βn Blade position [◦ ]
δi j Kronecker delta [−]
ς Constant of smoothing function [−]
ε Specific dissipation [m2 /s3 ]
εi j Dissipation tensor [m2 /s3 ]
ξ Thermal conductivity [W /mK]
γ Local pitch angle [◦ ]
κ von Karman constant [−]
` Turbulence length scale [m]

26
`mix Mixing length [m]
λ Tip speed ratio [−]
µ Dynamic viscosity [Ns/m2 ]
µT Eddy viscosity [m2 /s]
µsgs Eddy viscosity of the SGS [m2 /s]
ης Smoothing function [−]
η Kolmogorov length scale [m]
Ω Rotational speed [rad/s]
ω Specific dissipation rate [1/s]
ρ Density [kg/m3 ]
−u0i u0j Reynolds stresses [N/m2 ]
τi j Reynolds stress tensor [N/m2 ]
τsgs Subgrid scale stress term [N/m2 ]
ν Kinematic viscosity [m2 /s]
ψ(x,t) Flow variable [−]
ψ̄(x,t) Large scale filtered component [−]
ψ0 (x,t) Small scale sub-filtered component [−]
φsst Variables used in the k-ω SST model [−]
φω Variables used in the k-ω model (αω , ...) [−]

27
φε Variables used in the transformed k-ε model (αωε , ...) [−]
φ Angle between Urel and rotor plane [◦ ]
ϕ Blade twist [◦ ]
Θ Fluid domain [−]
∇2 Laplacian operator [−]
∆ = (∆1 ∆2 ∆3 )1/3 Filter width [m]
∆i Filter width in the ith direction [m]
ζ Viscous dissipation term [m2 /s]
∆p Change in pressure over the disc [N/m2 ]
∆t Time step [s]
∆x Spatial step [m]
Πi j Pressure-strain correlation tensor [m2 /s3 ]
σk Closure coefficient [−]
γ̇ Shear rate [1/s]
Abbreviations
1D One Dimension
2D Two Dimensions
3D Three Dimensions
28
ABL Atmospheric Boundary Layer
AD Actuator Disc
ADL Actuator Disc with Lift
ADLD Actuator Disc with Lift and Drag
ADM Actuator Disc Method
ADUT Actuator Disc with Uniform Thrust
ADV T Actuator Disc with Variable Thrust
ADV Acoustic Doppler Velocimeter
AL Actuator Line
ALM Actuator Line Method
APG Adverse Pressure Gradient
AS Actuator Surface
BE Blade Element
BEM Blade Element Momentum
BSL Baseline
CEL CFX Expression Language
CFD Computational Fluid Dynamics
CFL Courant–Friedrichs–Lewy
DES Detached Eddy Simulation

29
DNS Direct Numerical Simulation
DNW German Dutch Wind tunnel Organization
EU European Union
FDM Finite Difference Method
FEM Finite Element Method
FV M Finite Volume Method
GE General Electric
HAW T Horizontal Axis Wind Turbine
HPC High Performance Computer
IEA International Energy Agency
KE Kinetic energy
LES Large Eddy Simulation
LLF Large Scale Low Speed Facility
MCT Marine Current Turbines
MEXICO Model EXperiments In COntrolled COnditions
NDA Non Disclosure Agreements
NREL National Renewable Energy Laboratory
PIV Particle Image Velocimetry
RANS Reynolds-Averaged Navier-Stokes

30
RSM Reynolds Stress Model
RMS Root-Mean-Square
SST Shear Stress Transport
VAW T Vertical Axis Wind Turbine
VTM Vorticity Transport Model
UDF User Defined Function
UK United Kingdom
Chapter 1
Introduction
The world is currently moving towards renewable energy sources. This growth is
driven by policy and renewable energy targets such as the European Union (EU) Re-
newables directive [5] which aims to have 20% of the EU energy needs met by re-
newable energy by 2020. Two of these sources, wind and tidal, employ turbines to
extract energy from the environment. Tidal or marine energy is very promising as it
provides predictable energy and the UK has an estimated 50% of Europe’s tidal energy
resource1 . However, it is currently in its infancy and the UK only has 9MW installed
capacity. The wind energy industry is more mature, and continues to grow despite the
economic climate. In the year 2000 there was 17.4GW of installed capacity worldwide.
This grew by an average of 27% each year until 238GW was installed by 2011 [6]. As
capacity increases the percentage of growth is likely to fall. Between 2011 and 2012
growth slowed to 19% with total installed capacity reaching 283GW [6][7]. However,
approximately 13GW was installed during this time in both the US and China (the
two largest renewable energy markets) and almost 45GW worldwide [6]. Due to the
increased maturity and knowledge available, the majority of the work described in this
thesis will be focused on wind energy.
1 http://www.renewableuk.com/en/renewable-energy/wave-and-tidal
31
CHAPTER 1. INTRODUCTION 32
The UK wind industry, which is leading the offshore market, grew by 29% in 2012
to 8.4GW [6] and now contributes as much as 12% 2 to UK energy requirements.
This growth and the trend in wind-farms increasing in size, to the point where they are
sometimes constructed in stages, means manufacturers and developers are often unable
to create wind-farms using a single wind turbine specification. Many existing wind-
farms are also under consideration to be re-powered and/or extended. This can result in
wind-farms being made up of multiple designs and/or supplied by a number of different
companies, meaning different wind turbines with varying sizes and specifications are
installed at a single site. An example of this is the Roscoe wind-farm in Texas, USA
which uses three types of wind turbines (Mitsubishi 1MW, Siemens 2.3MW and GE
1.5MW) ranging from 350ft to 415ft in height. Such situations are seldom by design,
but are rapidly becoming more common as developers strive to meet renewable energy
policy targets by extending current wind-farms and building larger ones.
1.1 Background/motivation
A wind-farm using different wind turbines is not necessarily a negative proposition and
may, in fact, improve overall power output. Chowdhury et al. [8] discussed the factors
affecting the maximum power output of a wind-farm and included the wind turbine
rotor diameter as a variable in their optimization algorithm. The study [8] focused on
a single wind-farm case, which originally consisted of a 3 × 3 wind turbine array. To
allow thousands of configurations to be investigated, a simple velocity deficit model
[9] was used with the optimization algorithm. The algorithm found a 30% increase in
energy yield using identical wind turbines and a 43% increase using rotor diameter as
a variable. This suggests a 10% increase in energy yield of the wind-farm is possible
using wind turbines with different sizes. This provides promising insight into possible
2 http://www.gridwatch.templar.co.uk/
new optimization approaches, which have not yet been investigated. However, the
study [8] was limited to a single case using a single wind speed, hub height and blade
geometry while varying the rotor diameter and location. It is necessary to conduct a
more detailed study into how the rotor diameter influences the energy yield of a wind-
farm. To the author’s knowledge, there are no computational studies on the interactions
of different sized turbines within the literature.
1.2 Aim and objectives
This thesis aims to investigate the affect rotor diameter size has on turbine interactions,
through the use of Computational Fluid Dynamics (CFD) modelling and aerodynamic
analysis.
These will be demonstrated via achievement of the following objectives:
1. to identify current turbine models used within the literature and demonstrate
their implementation within commercially available CFD software;
2. to produce salient benchmark regimes of turbine model(s) for verification of

modelling procedures with analogous studies within the literature;
3. to predict the aerodynamic flow field interactions between two similar axially
aligned turbines;
4. to analyse the relative sensitivity of the aforementioned interactions to rotor size;
5. to critically evaluate the overall performance of two non-identical (different rotor

sizes) turbines;
6. to demonstrate if there are any benefits/limitations of using non-identical tur-

bines compared to using identical turbines;
7. to produce publishable quantitative and qualitative evidence in order to fully

describe any benefits/limitation realized.
1.3 Thesis overview
This thesis is divided into seven chapters. The opening chapter describes the motiva-
tion, key concepts and principles for the research described in this thesis. Chapter 2
discusses the theory behind the methods used in this thesis; considering the mathemat-
ics behind the main CFD methods available, focusing on the RANS and LES methods.
It concludes by detailing the rotor models employed in this thesis and their implemen-
tation within CFD codes. Chapter 3 consists of a literature review covering previous
theoretical, experimental and numerical work on turbines and their wakes, concentrat-
ing on current numerical methods. Chapter 4 describes the five CFD studies conducted
as part of this thesis, featuring simulations of a single turbine in water and air and two
axially aligned turbines. Chapter 5 details the results of the five studies, described in
Chapter 4, focusing on the interactions between two axially aligned turbines. Chapter
6 is the discussion chapter, where quantitative and qualitative evaluation of the results
obtained in this thesis is performed. Finally, Chapter 7 presents the major conclusions
of the work conducted in this thesis, highlighting the novel findings and concluding
with recommendations for further work.
1.4 Wind energy
Wind energy has always been an important resource. Windmills were built as long
as 3000 years ago and wind was very important for early travel utilizing sails [10].
For more information on wind energy’s history see the references [10][11] and their
extended reading lists.
Wind turbines work by converting the kinetic energy of the wind into usable en-
ergy, normally electrical. However, windmills convert the kinetic energy into mechan-
ical power and are commonly used to pump water or grind grain. There are two main
types of wind turbine: the Horizontal Axis Wind Turbine (HAWT), which is the most
common design today and will be the focus of this thesis; and the Vertical Axis Wind
Turbine (VAWT) such as the Darrieus and Savonius designs. Both designs are shown
in Figure 1.1. There are two variations on these types: drag-based variations which
try to catch or block the wind and use drag forces to rotate (e.g. cup anemometers)
and lift-based turbines, which use blades that work the same way as aeroplane wings,
responding to the wind by generating aerodynamic lift. Both of these forces are ex-
plained in more detail in subsection 3.1.3. The force created turns the rotor about its
axis, thus generating shaft torque which turns a generator. Turbines may or may not
have a gearbox depending on the configuration (see reference [10][11] for more de-
tail). The power extracted by a turbine is calculated from the kinetic energy in the air,
multiplied by the fraction extracted as shown in equation (1.1).
1 1 1
P = KE ×C p = ṁu2C p = (ρAr u)u2C p = ρπr2 u3C p . (1.1)
2 2 2
Where KE is the kinetic energy per sec, C p is the power coefficient, ṁ is the mass
flow-rate, u is the wind speed, ρ is the air density, Ar is the rotor swept area and r is
the rotor radius.
Equation (1.1) shows that changes in wind speed have the largest affect on power
extraction, which is why it is so important to place the wind turbine in the best location.
The power coefficient (discussed in subsections 3.1.1 and 3.1.2) varies between 0 − 0.6
and the density remains approximately equal to one, so have a smaller affect on the
power. The rotor radius is the easiest variable to change by design. It has a greater
affect on the power extracted than other variables, as the radius is squared. This is the
Figure 1.1: Types of wind turbine
main reason for the trend in wind turbine design towards bigger turbines in order to
produce more power.
1.4.1 Wind-farm
Calculating the power of a wind-farm is much more complex than a single turbine, as
a wind-farm is not just a sum of its individual parts.
PW F 6= n × PW T (1.2)
Where PW F is the wind-farm power, n is the number of turbines and PW T is the wind
turbine power [12].
Downwind turbines experience a power reduction by as much as 40% compared to
that of an isolated turbine [13][14]. When averaged over all wind directions this power
loss has been calculated at approximately 10 − 20% [15]. These power losses are due
to the wake of the upwind turbine and its interactions on the turbines downwind; this is
often referred to as the wake loss. The wake interactions also affects the fatigue loads
of the turbines. Due to these factors, wind-farm design is concerned with minimizing
wake interactions.
Wind-farm optimization
To overcome wake interactions, wind turbines need to be positioned to minimize the

affect of this phenomenon. However, this is not always practical, as it would require a
wide spaced wind-farm, creating a very low power density (power per area). The wake
of a wind turbine is extremely complex in nature and is affected by an extensive list
of variables including atmospheric turbulence and wind shear (change in wind speed
with height), as well as the wind direction. The latter can change rapidly and varies
both in time and space (location), which makes it difficult to predict and model.
There are a wide range of wind-farm optimization techniques and programs which
use various variables and assumptions. These are predominately mathematical in na-
ture, using simplified engineering models of the flow parameters with a mathematical
algorithm to determine the turbine’s location. There are a number of optimization
algorithms used including genetic algorithms, heuristic algorithms and swarm opti-
mization. For more details see reference [16].
The optimization methods have limitations, as all programs sacrifice accuracy for
speed. This enables hundreds and thousands of possible layout combinations to be
simulated and assessed. Some codes are also limited in terms of the turbine locations,
using only grid layouts or discrete points instead of continuous points in space. An-
other important limitation is the input data (which will be discussed in section 3.2).
Many optimization programs use a flow modeller to provide an estimation of the flow
field, instead of calculating it explicitly. One popular method is the Wind Atlas Anal-
ysis and Application Program, referred to as WAsP3 . WAsP uses wind mast data and
linear flow models to predict the flow field.
3 www.wasp.dk
1.4.2 Wind turbine wake
The rotor causes the wind to slow; this slower moving air is referred to as the wake.
The wake is commonly divided into two or three regions: the near wake and the far
wake, with some also considering the transition phase as an additional region [17].
This division in the wake is made by defining the near wake as the region which
is directly affected by the geometry of the rotor and the aerodynamics of the blades
(number of blades and their associated vortices). This definition is ambiguous, result-
ing in various region sizes ranging from one diameter to several, depending on the flow
conditions. Analysis in the near wake region focuses on wind turbine performance and
wake development. The far wake region is defined as beyond the near wake, where the
vortex system has broken down. Here the focus is on wake modelling and the turbine’s
influence on its surroundings.
1.5 Computational Fluid Dynamics
Fluid mechanics is the study of fluids either in motion (fluid dynamics) or at rest (fluid
statics). Computational Fluid Dynamics (CFD) is an extension of fluid mechanics; it
describes fluid motion through the use of numerical techniques using a computer. It
is concerned with the solution of the velocity field and hence, the pressure field for
given fluid properties (e.g. density, viscosity (measure of fluids resistance to flow)
and surface tension). Both liquids and gases are considered as fluids. The two main
obstacles in fluid mechanics to a workable theory are viscosity and geometry [18].
In a small number of idealized cases viscosity can be neglected, rendering the flow
fields inviscid. However, for the vast majority of cases viscosity is present and may
even vary, which increases the difficulty of constructing and using the basic equations
[18]. Moreover, it can have a destabilizing effect on the aforementioned flow field
creating a disorderly, random phenomenon known as turbulence.

Complex geometry featuring rapid changes in the geometrical bounds of the fluid,
make it difficult to find so-called ‘closed-formed analytical solutions’ of problems,
which is why most textbooks concentrate on simple geometries such as flat plates and
pipes. It is possible to solve problems featuring complex geometry and viscosity us-
ing numerical techniques on a computer which is commonly known as Computational
Fluid Dynamics (CFD).
CFD implements conservation principles, including the three fundamental conser-
vation equations of:
• Mass (continuity)
• Momentum (Newton’s second law [19])
• Energy (first law of thermodynamics [20][21])
By applying conservation principles to fluid motion it is possible to define the Navier-

Stokes (N-S) Equations, which can be given in different forms. Equation (1.3) - (1.4)
gives them in incompressible vector form. These equations describe the motion of any
Newtonian fluid.
→
−
∇· V = 0 (1.3)
→
− !
∂V → − → − − →
→ −
ρ + V · ∇ V = −∇p + µ∇2 V + S (1.4)
∂t
→
−
Where V is the velocity vector, t is time, ρ is density, p is the pressure, µ is the
→
−
viscosity, S is a source term and ∇2 is the Laplacian operator.
Incompressibility is assumed for flows when the local Mach number (the ratio of
the speed of an object to the speed of sound of the fluid) is below 0.3. For wind turbines
the flow has velocities normally around 5 − 25m/s, resulting in a Mach number of less
Figure 1.2: Illustration of the resolved and modelled parts of the flow
than 0.1, therefore the incompressible assumption is completely justifiable [22].

Although the incompressible Navier-Stokes equations implicitly describe turbulent
flow fields they are inherently difficult to solve due to the non-linear nature, specifically
→
− →
−
the convective acceleration term ( V · ∇ V ). Due to this there are a number of differ-
ent CFD methodologies, each of which numerically resolve the Navier-Stokes equa-
tions and model the resulting flows to varying degrees as shown in Figure (1.2). Four
such methods are the Reynolds-Averaged Navier-Stokes (RANS), Large Eddy Simula-
tion (LES), Detached Eddy Simulation (DES) and Direct Numerical Simulation (DNS)
which are briefly mentioned in this section with a more detailed description in Chapter
2.1. There is much research reported in the literature as well as a number of textbooks
with greater detail regarding these methods, including references [18][23][24]. The
difference between these models is their consideration of turbulence modelling.
1.5.1 Turbulence
Turbulence is characterized by random velocity variations within the flow that are gen-
erally modelled using statistical methods, based on experimental work [24]. If a fluid
is not turbulent then it is said to be laminar, and is normally characterized by steady
parallel fluid layers with no mixing. The problem with solving turbulent fluid prob-
lems is that turbulence is inherently three-dimensional, and operates at a smaller length
scale than laminar (non-turbulent) flows. The Atmospheric Boundary Layer (ABL)
can have a length scale ranging from approximately 1km to as small as 1mm [25] and
blade boundary layers can have even smaller scales. This creates a problem in trying
to resolve all the length scales. However, it is sometimes feasible and is called Direct
Numerical Simulation (DNS) which is extremely computationally expensive and is not
always viable. To overcome this a large number of turbulence models have been con-
structed for these small length scales, allowing the large scales to be resolved. There
are a large number of books and other materials available on this approach, such as
reference [24].
1.5.2 Reynolds number
The Reynolds number, Re, is a dimensionless number (has no units) that is used in fluid
mechanics to describe flow. In 1883 Osborne Reynolds [26] showed that the transition
from laminar to turbulent flow occurs at a critical value. The Reynolds number is
commonly calculated using equations (1.6)-(1.7), referred to as the ratio of inertia
force to viscous force, as shown in equation (1.5) with a proof given in Appendix A.
Laminar flows feature low Reynolds numbers. The higher the Reynolds number the
more likely a flow is to be turbulent; and the finer the small length scales (mixing)
within the flow.
Inertia f orce
Re = (1.5)
Viscous f orce
ρul
= (1.6)
µ
ul
= (1.7)
ν
Where ρ is the density, u is the speed, l is a characteristic length, µ is the dynamic

µ
viscosity and ν = ρ is the kinematic viscosity.
1.5.3 Reynolds-Averaged Navier-Stokes
The Reynolds-Averaged Navier-Stokes (RANS) method is the fastest CFD approach.

It was first described in 1895 by Osborne Reynolds [27]. It uses time averaged motion,
which introduces a new quantity known as the Reynolds stress. RANS models com-
pute this variable either through the use of an eddy-viscosity model or Reynolds Stress
model. The eddy-viscosity models vary with complexity, ranging from algebraic mod-
els to two equation models. The most commonly used in wind turbine applications
are the two equation models k-ε [28], k-ω [24] and k-ω SST [29], where k is the tur-
bulence kinetic energy per unit mass (or specific turbulence kinetic energy), ε is the
specific dissipation and ω is the specific dissipation rate and are described in section
2.1.1.
The problem with the eddy-viscosity models is that they are based on the Boussi-
nesq hypothesis [30], discussed in section 2.1.1. It is considered inadequate and invalid
for flows featuring sudden changes in mean strain rate, rotation and 3D effects [24], all
of which feature in turbine wakes. Nonetheless, it is still possible to use RANS meth-
ods without relying on the Boussinesq hypothesis by using the Reynolds Stress Method
(RSM). However, RSM requires six additional equations to be solved, making it much
more resource expensive; and similar assumptions to the Boussinesq hypothesis are
used.
1.5.4 Large Eddy Simulations
Large Eddy Simulation (LES) is able to handle unsteady anisotropic turbulent flow
much better than RANS, but requires more computational resources. LES is increas-
ingly being used due to the increase in computer power and resources that are now
available.
As the name suggests, LES resolves large eddies while the eddies smaller than the
grid size are modelled using a sub-grid scale (SGS) model. This approach is based
on the assumption that the smallest eddies have a more or less universal characteristic
that is not dependent on the flow geometry. This distinction between resolved and
modelled is done through a filtering function applied to the Navier-Stokes equations.
This introduces a new term called the SGS stress (similar to the Reynolds stress tensor
of the RANS method). This stress is used to represent the effects of the SGS modelled
eddies on the larger resolved eddies of the flow. This stress needs to be calculated and
there are a large number of models available.
Although specific models have been developed, it is possible to use turbulence
models developed for the RANS application (such as k-ε and k-ω models), which
means that problems with the Boussinesq Hypothesis are still encountered. However,
the fact that they are only used on the small scale eddies reduces the documented lim-
itations. Réthoré [31] showed that if cell size is taken small enough, then the error
made by the k-ε model and the eddy-viscosity concept are negligible. Due to this, it is
very common to use one equation models based on the Smagorinsky model [32], such
as the Smagorinsky Lilly model [33][34].
Near solid boundaries, LES becomes extremely computationally expensive. A
method used to reduce this is a hybrid approach where RANS is employed in the
near wall region and LES everywhere else, this is called a Detached Eddy Simulation
(DES).
Chapter 2
Theory
This chapter describes the theory behind the methods employed in this thesis. It cover
the main CFD methods focusing on the RANS and LES methods as well as discretiza-
tion methods. The chapter concludes with descriptions of the actuator methods (as
these are the main rotor models used within the literature), and their implementation
with CFD codes.
2.1 Computational Fluid Dynamics
This section discusses the mathematics used as part of the Computational Fluid Dy-
namics (CFD) process. For more details including the derivation of these equations see
the references [24][35][36]. This thesis considers only incompressible situations and
as such only the incompressible equations will be discussed. It should be noted that
some of the nomenclature used here is slightly different to those of the references, this
is to provide clearer differences between described models and provide consistency
within this thesis.
45
CHAPTER 2. THEORY 46
CFD uses numerical techniques and a computer to solve fluid flow problems, specif-
ically the Navier-Stokes (N-S) equations, which can be given in different forms. Equa-
tion (1.3) - (1.4) gives them in compact incompressible vector form which does not
do them justice in term of complexity, whereas equations (2.1) - (2.4) gives them in
Cartesian form. These equations describe the motion of any Newtonian fluid.
∂u ∂v ∂w
+ + =0 (2.1)
∂x ∂y ∂z
2
∂u ∂u ∂u ∂u
ρ +u +v +w = − ∂p
∂x + µ ∂ u
∂x2
+ ∂2 u
∂y2
+ ∂2 u
∂z2
+ Sx (2.2)
∂t ∂x ∂y ∂z
2
∂v ∂v ∂v ∂v
ρ +u +v +w = − ∂p
∂y + µ ∂ v
∂x2
+ ∂2 v
∂y2
+ ∂2 v
∂z2
+ Sy (2.3)
∂t ∂x ∂y ∂z
2
∂w ∂w ∂w ∂w
ρ +u +v +w = − ∂p
∂z + µ ∂ w
∂x2
+ ∂2 w
∂y2
+ ∂2 w
∂z2
+ Sz (2.4)
∂t ∂x ∂y ∂z
Where x, y, z is the Cartesian coordinate system, u, v, w is the Cartesian velocity

system, t is time, ρ is density, p is the pressure, µ is the viscosity and Sx,y,z is a source
term.
As previously mentioned the Navier-Stokes equations are inherently difficult to solve

directly due to their non-linear nature, specifically the convective acceleration terms.
There are a number of numerical methods that approximate these equations. The
most commonly used method involves using Reynolds averaging or Reynolds decom-
position which was first described in 1895 by Osborne Reynolds [27]. This splits
the instantaneous velocities into a mean velocity (ū) and a fluctuation (u0 ) such that
u = ū + u0 . In this case the mean component has to have a number of properties such as
the mean of the fluctuation is zero (ū0 = 0) and is known as a Reynolds operator. Using
Reynolds averaging on the Navier-Stokes equations and then time averaging the result
produces a time averaged version of these equations which are commonly known as
the Reynolds-Averaged Navier-Stokes (RANS) equations shown in Cartesian form in
equations (2.5)-(2.8) and in tensor form in equation (2.9).
∂ū ∂v̄ ∂w̄

+ + =0 (2.5)
∂x ∂y ∂z
2
∂ū ∂ū ∂ū ∂ū ∂ p̄ ∂ ū ∂2 ū ∂2 ū
ρ + ū + v̄ + w̄ = − +µ + +
∂t ∂x ∂y ∂z ∂x ∂x2 ∂y2 ∂z2
0 0
∂u u ∂u0 v0 ∂u0 w0
−ρ + + + Sx (2.6)
∂x ∂y ∂z
2
∂v̄ ∂v̄ ∂v̄ ∂v̄ ∂ p̄ ∂ v̄ ∂2 v̄ ∂2 v̄
ρ + ū + v̄ + w̄ = − +µ + +
∂t ∂x ∂y ∂z ∂x ∂x2 ∂y2 ∂z2
0 0
∂u v ∂v0 v0 ∂v0 w0
−ρ + + + Sy (2.7)
∂x ∂y ∂z
2
∂w̄ ∂w̄ ∂w̄ ∂w̄ ∂ p̄ ∂ w̄ ∂2 w̄ ∂2 w̄
ρ + ū + v̄ + w̄ = − +µ + 2+ 2
∂t ∂x ∂y ∂z ∂x ∂x2 ∂y ∂z
0 0
∂u w ∂v0 w0 ∂w0 w0
−ρ + + + Sz (2.8)
∂x ∂y ∂z
!
Dūi ∂ p̄ ∂u0i u0j
ρ =− + µ∆ūi − ρ + Si (2.9)
Dt ∂xi ∂x j
Where x, y, z is the Cartesian coordinate system, ū, v̄, w̄ is the mean velocity com-
ponent, u0 , v0 , w0 is the velocity fluctuation, t is time, ρ is density, p is the pressure, µ is
the viscosity, Sx,y,z is a source term and i, j represent components in the i, jth direction.
This process creates six new terms, which are commonly known as Reynolds
stresses and denoted by τi j = −ρu0i u0j . The problem is that this creates additional terms
with no additional equations, creating an incomplete system and is often referred to
as the RANS closure problem. To close the system turbulence models are used to
approximate the unknown terms.
Eddy-viscosity Turbulence models
There are many different turbulence models available, although there are two predom-
inant classes of turbulence models used with RANS which attempt to calculate or
approximate the Reynolds stresses. The first class is based on the eddy-viscosity as-
sumption and the other is the Reynolds Stress Model (RSM).
The eddy-viscosity models vary with complexity and approximate the Reynolds
stresses using the Boussinesq hypothesis also called the Boussinesq approximation
and Boussinesq eddy-viscosity assumption. The Boussinesq hypothesis, proposed by
Boussinesq [30], is a linear relation, shown in equation (2.10), between the stress tensor
(τi j ) and mean strain tensor (S) and is analogous to the linear relation between the
viscous stress tensor and viscosity [37]. It relates the Reynolds stresses to the velocity
gradients through an eddy-viscosity, in an analogy to molecular viscosity, assuming
that the Reynolds stress tensor is coincident with the mean rate tensor. The implications
of this are detailed in many textbooks and publications including Wilcox [24] which
states there is no physical basis for this assumption. In particular, Schmitt [37] analysed
the validity of the hypothesis using databases as well as DNS and LES test cases and
found that it is almost never verified. Despite this the eddy-viscosity is widely used to
model turbulence and can perform well when used correctly.
τi j = −µT S (2.10)

1 ∂ū j ∂ūi
S = + (2.11)
2 ∂xi ∂x j
Where R is the stress tensor and S is the mean strain tensor, µT is the eddy-viscosity,
ū is the mean velocity component, xi, j is the Cartesian coordinate system and i, j rep-
resent components in the i, jth direction.
There are numerous models which use the eddy-viscosity concept to different de-
grees. These models are divided into algebraic models, one equation models and two
equations models.
Algebraic models
Algebraic or zero equation turbulence models are computationally the simplest mod-
els and often compute the eddy-viscosity in terms of the mixing length. The mixing
length hypothesis was proposed by Prandtl in 1925 [38]. The hypothesis leads to defin-
ing the eddy-viscosity as:
d ū
µT = ρ`2mix (2.12)
dy
where µT is the eddy-viscosity, ρ is density and `mix is the mixing length.

The mixing length is empirical, depending on the flow and must be specified in ad-
vance. Therefore, by definition the algebraic models are incomplete turbulence models.
One equation models

One equation models calculate one turbulent transport equation, usually the turbu-
lence kinetic energy. It is technically the specific turbulent kinetic energy as it is per
unit mass, shown in equation (2.13), although it is commonly just called the turbu-
lence kinetic energy. One equation models are incomplete as they normally relate the
turbulence length scale to some flow parameter. The original one equation model was
developed by Prandtl in 1945 [39]. Prandtl calculated the specific turbulent kinetic
energy as shown in equation (2.14) and assumed that the turbulent dissipation (ε) was
proportional to k3/2 /`. The turbulence length scale still needs to be defined in advance
to close the system.
1 02
k= u +v +w
02 02 (2.13)
2

∂k ∂k ∂ūi k3/2 ∂ µT ∂k
ρ + ρū j = τi j −Cρ + µ+ (2.14)
∂t ∂x j ∂x j ` ∂x j σk ∂x j
µT = ρk1/2 ` (2.15)
Where k is the specific turbulence kinetic energy, xi, j is the Cartesian coordinate
system, ūi is the mean velocity component, t is time, ρ is density, µ is the viscosity,
µT is the eddy-viscosity, C, σk are closure coefficients, i, j represent components in the
i, jth direction and ` is the turbulence length scale.
Since Prandtl’s one equation model there have been many more models developed,
most more elaborate than his. Prandtl’s model only requires two closure coefficients
(C, σk ) whereas other have many more such as Baldwin-Barth model [40] which has
seven closure coefficients and two damping functions and the Spalart-Allmaras model
[41] which features eight closure coefficients and three damping functions.
Two equation models

Two equation models are complete, providing a closed system and are the simplest
complete system used in CFD modelling. They calculate two equations: one for the
specific turbulence kinetic energy and the other is related to the turbulence length scale.
There are many different models, the main difference being the second equation and
the variable calculated, including the turbulence length scale (`), specific dissipation
(ε) or specific dissipation rate (ω) [42][43]. The first two equation model was devel-
oped by Kolmogorov in 1942 [44] which was a k-ω model. All two equation models
feature closure coefficients which have been used to replace double and triple correla-
tion terms with algebraic expression involving known properties. If turbulence theory
was exact these coefficients would be set from first principles but it is not and the
models are developed based on dimensional analysis and observations [24]. There are
many models but the main models used today are the k-ε, k-ω and k-ω SST models,
which are discussed in more detail.
The k-ε model

The k-ε model is the most popular model and in many cases the default model. It
features equations of the specific turbulence kinetic energy (k) and the dissipation per
unit mass (ε) often just referred to as the dissipation. It was first developed by Jones &
Launder [28] and later improved by Launder and Sharma [45]. The model is so widely
used it is normally referred to as the standard k-ε model or just the k-ε model without
any references. The model defines k using equation (2.16) and ε using equation (2.17).

∂k ∂k ∂ūi ∂ µT ∂k
ρ + ρū j = τi j − ρε + µ+ ε (2.16)
∂t ∂x j ∂x j ∂x j σ1 ∂x j
2
∂ε ∂ε εε ∂ūi ε ε ∂ µT ∂ε
ρ + ρū j = α τi j −β ρ + µ+ ε (2.17)
∂t ∂x j k ∂x j k ∂x j σ2 ∂x j
system, ūi is the mean velocity component, t is time, ρ is density, µ is the viscosity, µT is
the eddy-viscosity, αε , βε , σε1 , σε2 are closure coefficients and i, j represent components
in the i, jth direction.
This model has five closure coefficients which are αε = 1.44, βε = 1.92, σε1 = 1.0,
σε2 = 1.3 and Cε = 0.09. It also has two auxiliary relations which are that ω = ε/(Cε k)
and ` = Cε k3/2 /ε. Equation 2.17 is only an approximation to the exact equation of the
turbulence dissipation [24][46] which is shown in equation (2.18). The problem is that
the exact equation is too difficult to solve as it introduces several new unknowns.
∂ε ∂ε
+ ū j = −2ν[u0i, j u0j,k ū j,k + u0k,i u0k, j ūi, j ] − 2νu0j u0i,k ūi, jk
∂t ∂x j
−2νu0i, j u0j,k u0k,i − 2ν2 u0i, j j u0i,kk (2.18)
∂ε
+ν − νu0i,k u0i,k u0j − 2νp0k u0j,k
∂x j
The k-ε model does have a number of flaws and limitations, which are well docu-
mented such as performing badly with severe pressure gradients and separation [24][29].
To this extent there have been a large number of modifications to the standard model
to improve it. The most commonly used of these are: RNG k-ε model and realizable
k-ε model. The RNG k-ε model or Renormalization Group Method k-ε model adds
a term to the turbulence dissipation equation to account for the interaction between
the dissipation and mean shear. The realizable k-ε model changes one of the closure
coefficients (Cε ) to allow it to vary based on the velocity gradients.
The k-ω model

There are a number of methods that use the specific dissipation rate (ω). The most
common k-ω model is often referred to as the standard k-ω model or the Wilcox k-
ω model [24][47]. This model defines k using equation (2.19) and ω using equation
(2.20).

∂k ∂k ∂ūi ω ∂ ω ∂k
ρ + ρū j = τi j − β1 ρkω + (µ + σ1 µT ) (2.19)
∂t ∂x j ∂x j ∂x j ∂x j

∂ω ∂ω ωω ∂ūi ω 2 ∂ ω ∂ω
ρ + ρū j = α τi j − β2 ρω + (µ + σ2 µT ) (2.20)
∂t ∂x j k ∂x j ∂x j ∂x j
system, ūi is the mean velocity component, t is time, ρ is density, µ is the viscosity, µT
is the eddy-viscosity, ω is the specific dissipation rate αω , βω ω ω ω

1 , β2 , σ1 , σ2 are closure
coefficients and i, j represent components in the i, jth direction.

This model has five closure coefficients which are αω = 5/9, βω ω
1 = 9/100, β2 =
3/40, σω ω ω
1 = 1/2 and σ2 = 1/2. It also has two auxiliary relations which are ε = β1 ωk
and ` = k1/2 /ω. It is known to perform better under adverse pressure gradients and
for flows featuring separation than the k-ε model. It is also considered more accurate
within the wall region [48]. Although it has been shown to be more sensitive to free
stream velocities at the boundary layer edge
The k-ω SST model

The k-ω shear stress transport model normally written k-ω SST or SST k-ω model is
a combination of the standard k-ω and k-ε models defined by Menter [29][49][50]. The
k-ω SST model is a variation and improvement of the baseline (BSL) model. The BSL
model varies between the k-ω and k-ε models in an almost linear manner whereas the
k-ω SST model also accounts for the turbulent shear stress. It uses the k-ω definition in
the inner part of the boundary layer and the k-ε definition outside the boundary layer.
Menter [29] created one set of equations by blending the k-ω and k-ε model equa-
tions, by first transforming the k-ε model into a k-ω formulation, shown in equation
(2.21).

∂ω ∂ω ωε ω ∂ūi ωε 2 ∂ ωε ∂ω
ρ + ρū j = α τi j − β2 ρω + (µ + σ2 µT )
1 ∂k ∂ω
+2ρσωε 2 (2.21)
ω ∂x j ∂x j
Where αωε , βωε ωε ωε ωε

1 , β2 , σ1 , σ2 are the new closure coefficients and all other vari-
ables are as previously defined. The two models are combined to create a new set of
equations for k and ω as shown in equations (2.22) and (2.23).

∂k ∂k ∂ūi sst ∂ sst
∂k
ρ + ρū j = τi j − β1 ρkω + µ + σ1 µ T (2.22)
∂t ∂x j ∂x j ∂x j ∂x j

∂ω ∂ω sst ω ∂ūi sst 2 ∂ sst
∂ω
ρ + ρū j = α τi j − β2 ρω + µ + σ2 µT
1 ∂k ∂ω
+2ρ(1 − F1 )σsst 2 (2.23)
ω ∂x j ∂x j
Where αsst , βsst sst sst sst

1 , β2 , σ1 , σ2 are the new combined closure coefficients, F1 is the
blending function and all other variables are as previously defined.

A blending function (F1 ) is used to gradually change between the two methods in
the desired regions. It is designed so that the k-ω formulation is used in the near wake
region and the k-ε formulation is used in the free shear layers. With the coefficients
also blended using equation (2.24)
φsst = F1 φω + (1 − F1 )φε (2.24)
where φsst represents the variables used in the k-ω SST model (αsst , ...), φω repre-
sents the variables used in the k-ω model (αω , ...) and φε represents the variables used
in the transformed k-ε model (αωε , ...).
The blending function is based on a set of logical arguments and utilizes both spa-
cial and fluid terms as shown in equation (2.25).
 √ ! !!4 
ωε
4ρσ2 k
k 500ν
F1 = tanh  min max ; 2 ;  (2.25)
0.09ωys ys ω CDkw y2s

ωε 1 ∂k ∂ω −10
CDkw = max 2ρσ2 , 10 (2.26)
ω ∂x j ∂x j
Where k is the specific turbulence kinetic energy, ω is the specific dissipation rate,
ys is the distance to the nearest surface, ρ is the density, ν is the kinematic viscosity,
σωε
2 is a closure coefficient and CDkw is the cross diffusion term.
The three terms in equation (2.25) correspond to arguments related to the turbulent
length scale, viscous sublayer and small free stream values. The limiter has a factor
of 10 in this formulation in line with the latest improvements [50] instead of 20 as
initially proposed [29][49].
The advantage of the k-ω SST model over the BSL model is that it accounts for the
turbulence shear stress by limiting the eddy-viscosity as shown in equation (2.27). It
is limited based on the assumption that the turbulent shear stress is proportional to the
specific turbulence kinetic energy in the logarithmic and wake regions of the turbulence
boundary layer [51].
ρa1 k
µT = (2.27)
max(a1 ω, SF2 )
 √ !!2 
k 500ν
F2 = tanh  max 2 ; 2  (2.28)
0.09ωys ys ω
Where µT is the eddy-viscosity, a1 is a constant equal to 0.31, k is the turbulence

kinetic energy, ω is the specific dissipation rate, S is the invariant measure of the strain
rate, F2 is the second blending function, ys is the distance to the nearest surface, ρ is
the density and ν is the kinematic viscosity.
The k-ω SST model is an improvement over the other two equation models for
many flow situations, although it is limited by the use of the eddy-viscosity and Boussi-
nesq hypothesis. These models assume isotropic turbulence, which is not normally
valid and can cause unrealistic results.
Near wall region
For turbulent flow conditions, the velocity profile can be split into two regions: the
viscous sublayer and the inertial layer. The location of the outer edge of the inertial
layer depends on the Reynolds number. The extent of the inertial layer grows with
increasing Reynolds number.
To fully resolve the boundary layers in the near wall region requires excessive grid
resolution. It requires the first computational cell from the wall be located within the
viscous sublayer and achieve y+ = 1, shown in equation (2.29). Wall functions are
used to significantly reduce computational expense by reducing the grid resolution and
modelling the flow using semi-empirical functions, shown in Figure 2.1. A correspond-
ing wall function grid can have a the first cell with a y+ = 100 [52]. There are many
different wall functions available, used for different turbulence models. Wall functions
will not be discussed here in detail as the topic is too vast, only a brief summary of
their implications. For more information see the references [52].
U ∗ yd
y+ ≡ (2.29)
ν
Where U ∗ is the friction velocity at the wall, yd is the distance to the nearest wall
and ν is the local kinematic viscosity of the fluid.
For the k-ω model, an analytical expression is known for the viscous sublayer and
there is an approximation for the logarithmic layer[52]. Whereas for the k-ε model
there is no analytical expression for the viscous sublayer. This means the wall function
needs to approximate the entire boundary layer. Many functions assume the first grid
point is within the fully turbulent log-law region which has a significant impact on the
result and is not always the case and is why the k-ω model is more accurate for wall
bounded simulations.
Figure 2.1: Difference in grid resolution using a wall function
Reynolds Stress Model
The Reynolds Stress Model (RSM) also known as a second moment closure model,
closes the Reynolds-Averaged Navier-Stokes equations without relying on the Boussi-
nesq hypothesis. Instead, it solves additional transport equations for all components of
the Reynolds stress tensor (τi j ). However, the derivation of the RSM requires higher
order correlations to be modelled. RSM avoids the isotropic eddy-viscosity assump-
tion allowing the full 3D nature of turbulence to be modelled. This added accuracy
requires additional computational resources to solve the additional equations.
Nearly all RSMs start from the same point, the exact differential equation of the
Reynolds stress tensor as shown in equation (2.30).

∂τi j ∂τi j ∂u¯j ∂ūi ∂ ∂τi j
+ ūk = −τik − τ jk + εi j − Πi j + ν +Ci jk (2.30)
∂t ∂xk ∂xk ∂xk ∂xk ∂xk
∂u0i ∂u0j
εi j = 2ν (2.31)
∂xk ∂xk
!
∂u0i ∂u0j
Πi j = p0 + (2.32)
∂x j ∂xi
Ci jk = pu0i u0j u0k + p0 u0i δ jk + p0 u0j δik (2.33)
Where εi j is the dissipation tensor, Πi j is the pressure-strain correlation tensor, Ci jk

is the turbulent transport tensor, δi j is the Kronecker delta and all other variables are
as previously described.
To close the equations, three tensors (dissipation tensor (εi j ), pressure-strain cor-
relation tensor (Πi j ) and the turbulent transport tensor (Ci jk )) needs to be modelled.
Rotta [53] was the first to develop a closure of the Reynolds stress equation in 1951.
There are many different approaches to modelling the tensors and some are based on
similar methods and equations used in the eddy-viscosity turbulent models. Of the
three tensors Πi j has received the most attention, this is because it has a similar mag-
nitude to the production term, so has a critical role in the flow and it is essentially
unmeasurable [24].
RSM is good for accurately predicting complex flows as it accounts for streamline
curvature, swirl, rotation and high strain rates. Although it does struggle with stability
without the eddy viscosity term and has a higher computational expense (2 − 3 times)
compared to the other turbulence models.
The Large Eddy Simulations (LES) computes the large eddies and models the smallest
eddies. The large eddies are dependent on the boundary conditions and geometry so
need to be calculated. However, the small eddies are nearly isotropic with more uni-
versal characteristics allowing them to modelled [24]. The small eddies are modelled
by subgrid scale (SGS) models.
The differentiation between the large eddies (resolved) and the small eddies (mod-
elled) of the flow field are accomplished using a filtering method. The filtering process
effectively filters out eddies smaller than the filter width. A common filtering method
decomposes a flow variable (ψ(x,t)) into the large scale filtered part (ψ̄(x,t)) and small
scale sub-filtered part (ψ0 (x,t)):
ψ(x,t) = ψ̄(x,t) + ψ0 (x,t) (2.34)
The filtered variable is defined at x0 as:
ˆ
ψ̄(x0 ) = ψ(x,t)G(∆)dx0 (2.35)
Θ
where Θ is the fluid domain, G is the filter function.

The filter function depends on the filter width ∆ = (∆1 ∆2 ∆3 )1/3 where ∆i is the filter
width in the ith direction. There are a number of different filtering functions, the most
common are the Gaussian filter (equation (2.36)), the tophat or box filter (equation
(2.37)) and cut-off filter (equation (2.38)).

6 3/2 |x0 − x|2
G = exp −6 (2.36)
π∆2 ∆2



 1 i f |x0 − x| ≤ ∆i
∆ 2
G = (2.37)


0 otherwise

3 sin ∆π (x0 − x)
G = ∏
i
(2.38)
i=1 π(x0 − x)
Whichever filtering function is used it is applied to the Navier-Stokes equations to

create filtered equations, as shown in equations (2.39) and (2.40), which are solved in
LES simulations.
∂Ūi
= 0 (2.39)
∂xi
∂Ūi ∂ 1 ∂P̄ ∂τsgs ∂2Ūi
+ (ŪiŪ j ) = − − + ν 2 + Si (2.40)
∂t ∂x j ρ ∂x ∂x j ∂x
Where Ūi is the large scale filtered component and τsgs is the subgrid scale stress
term.
Subgrid scale models
Most subgrid scale models are eddy-viscosity models using the Boussinesq hypothesis.
The most widely used models are based on the Smagorinsky model [32] which uses
an algebraic eddy-viscosity. This has been used due to its simplicity as the SGS eddy-
viscosity (µsgs ) is determined by:
q
µsgs = ρ(CS ∆) 2S¯ i j S¯i j
2
(2.41)
where S¯i j is the rate-of-strain tensor for the resolved scale, CS is the Smagorinsky
constant usually varying between 0.1 to 0.2.
Of the models based on the Smagorinsky model [32] the most commonly used is
the Smagorinsky-lilly model [33][34] which uses a modified SGS eddy-viscosity given
by:
q
µsgs = ρLS2 2S¯ i j S¯i j (2.42)
LS = min(κyd ,CsV 1/3 ) (2.43)

where LS is the mixing length for the SGS, κ is the von Karman constant, yd is the
distance to the closest wall, Cs is the Smagorinsky constant, and V is the volume of the
computational cell.
One problem with these models is that Cs is not universal. Lilly [33][34] derived a
value of 0.17 for homogeneous isotropic turbulence. However, this can cause excessive
dampening near solid boundaries and a value of around 0.1 has been found to yield the
best results for a wide range of flows [54].
LES is a transient method, requiring time discretization and definition of a time-
step. This is normally defined to maintain stability and to accurately resolve turbulence
in time. The most common condition is the Courant–Friedrichs–Lewy (CFL) condition
[55] that requires:
u∆t
≤1 (2.44)
∆x
where ∆t is the time step and ∆x is the spatial step.

For more information about LES, the SGS models and its applications see the ref-
erences [24][35][36][56][57][58].
2.1.3 Detached Eddy Simulation
Detached Eddy Simulation (DES) is a hybrid method which combines LES and RANS.
As LES becomes extremely computationally expensive within the near wall region,
one method used to reduce this expense is to use another formulation. DES uses the
LES formulation for the majority of the domain and switches to the RANS formulation
near wall boundaries. The switch between the two formulations can be accomplished
using turbulence length scales, distance variables or specified by the user within certain
zones known as zonal DES [59]. DES is only in its infancy as it was first proposed in
1997 [60], for more information on DES, its history and applications see the reference
[61].
2.1.4 Direct Numerical Simulations
Direct Numerical Simulations (DNS) solves the Navier-Stokes equations for all scales
of motion without any averaging or approximation. To capture all the turbulence scales
the small scales need to be captured as these are responsible for the turbulence dissipa-
tion and the Kolmogorov scale (η) is commonly considered as the smallest scale that
needs to be resolved, although this is not always the case [62].
This need to resolve all length scales means that the number of grid points required
are in the order of Re9/4 . Equally, the simulation must be time dependent and given that
the Kolmogorov time-scale is proportional to Re1/2 this results in an overall computa-
tional time of Re11/4 . This shows the dependence and limitation of DNS on the Re and
why only low Re flows have been simulated. For example, if Re= 106 (reasonable for
many applications such as in the aerospace industry) the grid would need to contain ap-
proximately 1013.5 (3162 billion) points, which is not possible considering the current
computational resources. In comparison, LES requires grids in the order of Re1.8 [36]
which results in LES typically requiring 5-10% the CPU time needed for analogous
DNS [24]. Despite the high computational requirements DNS is still an important tool
for CFD, particularly in regards to understanding turbulence and improving turbulence
models.
For more information on the DNS, applications and previous work see the refer-
ences [24][35][57][62].
2.1.5 Discretization
Computers do not solve the equations of motion in differential form, but in a discretized
form. The discretization process which is also referred to as meshing or grid gener-
ation, involves dividing the continuous fluid domain into discrete points (also called
nodes, or vertices) and volumes (cells or elements). There are two types of grid:
• structured grids are created with each grid point being uniquely defined by in-
dices i,j,k and are made up of hexahedral cells;
• unstructured grids are created with points that have no particular ordering. The
grid cells are commonly tetrahedral cells, although they can be made up of var-
ious cells including tetrahedral, hexahedral, prisms, pyramids and wedges such
grids are known as hybrid or mixed grids.
The advantage of a structured grid is that it is easy to relate cells and their neighbours,
which allows data to quickly and easily be calculated by adding or subtracting indices
by an integer. This is also a disadvantage when considering complex geometry that
can not be regularly divided. One approach to address this is known as the multiblock
approach which divides the geometry into smaller parts or blocks which can be dis-
cretized easier. This can create what is known as hanging nodes where nodes only
exist on one side of these block boundaries. Another approach is known as the overset
or chimera method where grids are created around boundaries independently and then
overlapped. The important part of this approach is in transferring the data between
grids within the overlapping region.
The unstructured grids are more flexible and can be generated more easily due to
the use of tetrahedral cells. However, it is common to use hexahedral and prism cells
to resolve the boundary layer. Another advantage is they can handle refinement and
coarsening with ease. The main disadvantage is that they require more complex data
structures to keep track of cells and their connections. This leads to higher memory
requirements compared to structured grids.
Once the domain has been discretized to form a grid the equations need to be
discretized. There are three main processes of discretization used; the Finite Difference
method (FDM), Finite Element Method (FEM) and Finite Volume Method (FVM).
Finite Difference Method
The Finite Difference Method (FDM) was among the first methods and was utilized by
Euler [35]. It is directly applied to the differential equations using the Taylor series ex-
pansion to discretize the flow variables. The advantage of this method is its simplicity
and ease to obtain high order approximations and hence accuracy. The disadvantage
is that it requires a structured grid and can not be applied to body fitted grids which
limits its applications. For more details on this method see the reference [35].
Finite Element Method
The Finite Element Method (FEM) was originally developed for structural analysis,
proposed in 1956 by Turner et al. [63]. It can accommodate an unstructured mesh
and the grid nodes are exclusively used to produce an arbitrarily high order of accu-
racy representation of the solution on smooth grids [35]. It works with an equivalent
integral form of the governing equations, commonly accomplished based on the weak
formulation. As this method uses an integral formulation and an unstructured grid,
it is preferable for flows featuring complex geometry and is particularly suitable for
non-Newtonian fluids. The finite element methods have a very rigorous mathemati-
cal foundation but can be shown to be mathematically equivalent to the finite volume
method for certain cases and have a higher numerical expense.
Finite Volume Method
The Finite Volume Method (FVM) was first employed in 1971 by McDonald [64] for
2D flows and directly uses the integral formulation of the Navier-stokes equations. It
discretizes the domain into control volumes on which the surface integral is calculated.
The accuracy of the discretization depends on the schemes used to define the control
volumes. As part of this methods there are two different schemes used to define the
shape and position of the control volume: the cell-centred scheme and the cell-vertex
scheme.
The cell-centred scheme has the control volumes identical to the grid cells and
the flow variables are stored at the centre of the cells (shown in Figure 2.2(a)). The
variables are calculated using the average of the fluxes computed from the values of
the cells to the left and to the right.
The cell-vertex scheme stores the variables at the grid points and can either have
the control volume as the sum of the cells sharing a grid point (overlapping control
volumes shown in Figure 2.2(b)) or the volume centred around a grid point (dual con-
trol volumes shown in Figure 2.2(c)) [35]. For the overlapping scheme the control
volumes are the cells but the variables are calculated at the vertex based on the sur-
rounding cells. The dual control volume scheme stores all variables at the node, with
the control volume constructed based on the midpoints of the cells surrounding the
node, shown in Figure 2.2(c).
Comparing the two schemes, the cell-centred scheme depend strongly on the smooth-
ness of the grid which can cause numerical errors that does not effect the vertex-centred
schemes. However, both schemes can reach second or higher order accuracy if the grid
is sufficiently smooth. The dual control volume scheme experiences difficulties at the
boundary of solid walls as the node is on the wall which leads to discretization errors
compared to the cell-centred scheme. For unstructured grids the cell-centred scheme
leads to approximately six times as many control volumes using tetrahedral and ap-
proximately three times as many when using a mixed grid compared to the cell-vertex
scheme. This suggests the cell-centred is more accurate, although it calculates the
variables from approximately half the number of points and the additional control vol-
umes result in additional computational expense, providing no clear advantage. The
cell-centred scheme requires approximately twice as much memory as it needs to store
more flow variables. Using hexahedral grids the two schemes are computationally
equivalent.
The FVM is very flexible and can be implemented on structured and unstructured
grids. It can be shown to be equivalent to both the FDM and FEM methods for certain
cases. Due to its flexibility it is the most widely used discretization process.
2.2 Rotor models
This section details the rotor models used as part of this thesis including a porous
disc and actuator methods. It begins with descriptions of each rotor model used in
this work and their underpinning theory. Followed by their implementation within
the commercially available optimized Computational Fluid Dynamics (CFD) codes
ANSYS-CFX [65] and ANSYS-Fluent [66].
2.2.1 Porous disc
A porous disc is one of the simplest models used to represent (renewable energy) tur-
bines. These commonly feature predefined momentum losses. Similar approaches can
also be described using analogous equations such as momentum sinks and drag discs.
All approaches feature the same diameter as the turbine with a thin depth/thickness.
The disc needs a thickness to avoid singularities caused by discontinuities in the pres-
sure and velocity fields.
(a) (b)
(c)
Figure 2.2: Control volumes of (a) the cell-centred scheme, (b) the cell-vertex scheme
(overlapping) and (c) the cell-vertex scheme (dual control volumes) [35]
(a) (b)
Figure 2.3: Geometry of two methods used to define disc location a) algebraically b)
geometrically
Implementation
There are two methods that can be used to define the disc location: either algebraically
as part of the solver set-up, which is described in subsection 2.2.2 or geometrically by
creating the disc as part of the geometry. Figure 2.3 shows the difference in geometry
for these two methods. If the disc is constructed geometrically then the shape and size
of the disc needs to be defined before the domain is discretized. It is then possible to
define the disc as a subdomain or cell zone for the source terms to be applied to.
ANSYS-CFX
The subdomain concept is used within ANSYS-CFX to define porous media. There
are multiple ways to impose a porous loss, for example using permeability and/or a
loss coefficient for (an)isotropic loss media. Within this thesis the loss coefficient and
directional loss model were used to model an isotropic momentum loss. In ANSYS-
CFX the directional loss model adds a momentum source term to the flow, given by:
ρ
S = −K u|u| (2.45)
2
where K is defined by the user and referred to as the resistance coefficient within
ANSYS-CFX, ρ is the density and u is the velocity. The resistance is applied as a loss
across the disc thickness and so specified by the user as Kd , where d is the thickness of
the disc. The equation features u|u| to define the source term as a vector in the direction
of the flow.
ANSYS-Fluent
ANSYS-Fluent uses cell zones analogous to the aforementioned subdomain used in

ANSYS-CFX. A cell zone can be used to define a porous media which is modelled
using an additional momentum source term, shown in equation (2.46), which is divided
into two components: a viscous loss term (linear) and inertial loss term (quadratic).
!
3 3
1
Si = − ∑ C1 µui + ∑ C2 ρui |u| (2.46)
i=1 i=1 2
1
Si = −C2 ρu|u| (2.47)
2
Where Si is the source term for the ith momentum equation, µ is the fluid viscosity,
u is the velocity and C1 and C2 are prescribed matrices. The momentum sink con-
tributes to the pressure gradient in the porous zone, creating a pressure drop that is
proportional to the fluid velocity squared in the cell [54]. For a simple homogeneous
porous media in one direction, the permeability term can be eliminated (C1 = 0) [54],
equation (2.46) becomes equation (2.47) which is analogous to equation (2.45) used
by ANSYS-CFX. Similarly, the resistance is applied as a loss across the disc thickness
C2
and so specified as d , where d is the thickness of the disc.
2.2.2 Actuator disc method
The most common way to represent the rotor is by using actuator methods, where
the blades are represented by body forces defined on either a disc, line or surface. The
Actuator Disc Method (ADM) represents a time averaged rotor using forces distributed
over a disc with the same dimensions and volume as the swept volume of the rotor.
There are a number of variations of the ADM within the literature including constantly
and variably loaded discs, as well as variations with and without rotation [67][68][69].
The forces acting on the flow are determined by the aerofoil characteristics and the
local flow velocities. The aerofoil characteristics of thrust (CT ), lift (CL ) and drag
(CD ) coefficients (discussed in subsections 3.1.1 and 3.1.3) are commonly determined
from 2D tabulated data based on the local angle of attack (α), which is defined by
equation (2.48), and the Reynolds number (§§1.5.2).The relative velocity of the flow
to the aerofoil (Urel ) is defined by equation (2.50) with the angle between the relative
velocity and the rotor plane (φ) given by equation (2.49) and shown in Figure 2.4.
α = φ−γ (2.48)

Vz
φ = tan−1 (2.49)
Ωr −Vθ
q
Urel = Vz2 + (Ωr −Vθ )2 (2.50)
Where γ is the local pitch angle, r is the radius, Vz is the axial velocity, Ω is the
rotational speed of the rotor defined by the user and Vθ is the tangential velocity.
There are a number of variations, each of which defines the force and source terms
used slightly differently. The base term of each variation is defined using equation
(2.51). It is possible to define the disc using the thrust coefficient [70][71][72] to
create a uniformly loaded disc using equation (2.52) defined as a force per unit volume.
Figure 2.4: Aerofoil cross-section showing force vectors
A non-uniformly loaded or variably loaded disc can be defined using equation (2.53)
where the variation is created by including the chord length (c) of the blades; creating a
variation with radial position which is constant over each annulus. However, it is more
common to define the disc based on the lift [73] or lift and drag coefficients [74][75]
as shown in equations (2.54)-(2.56) respectively.
1 2
f = ρU (2.51)
2 rel
1 2
DADUT = ρU CT B = fCT B (2.52)
2 rel
1 2
DADV T = ρU CT Bc = fCT Bc (2.53)
2 rel
1 2
DADL = ρU CL Bc = fCL Bc (2.54)
2 rel
1 2
DADLD = ρU Bc(CL eL +CL eD ) (2.55)
2 rel
= f Bc(CL eL +CL eD ) (2.56)
Where ρ is the density of air, Urel is the relative velocity of the rotor, CT is the thrust
coefficient, B is the number of blades, c is the chord length, CL is the lift coefficient,
CD is the drag coefficient and eL , eD are unit vectors used to define the direction of lift
and drag respectively. The subscript abbreviations stand for Actuator Disc Uniform
Thrust (ADUT), Actuator Disc Variable Thrust (ADVT), Actuator Disc Lift (ADL)
and Actuator Disc Lift and Drag (ADLD).
Whichever definition term (Dn given in equation (2.52) - (2.56)) is used, it is then
averaged over each annulus given in equation (2.57) for 3D or 2D by neglecting d in
equation (2.57). The inclusion of lift (CL ) and drag (CD ) coefficients and their unit
vectors (eL , eD ) allows for rotational effects to be included by adding tangential forces
as well as axial forces to the disc. It has been shown that the inclusion of rotation
significantly improves predictions [68][69][76]. Porté-Agel et al. [77] showed that non
rotational methods over predict turbulence intensity at the centre and under predict it
at the top of the wake. The force per volume is commonly applied to the flow in
cylindrical coordinates as a source term [73], as shown in equations (2.58)-(2.60),
within the Navier-Stokes equations as discussed in section 1.5.
Dn
DAD = (2.57)
2πrd
Sx = DAD (CL cos φ +CD sin φ) (2.58)
Sθ = DAD (CL sin φ −CD cos φ) (2.59)
Sr = 0 (2.60)
Where Dn is the disc definition term (given in equation (2.52) - (2.56)), r is the
radius, d is the depth of the rotor, CL is the lift coefficient, CD is the drag coefficient, φ
is the angle between Urel and rotor plane, defined in equation (2.49).
Tip loss correction
Use of a disc creating the assumption of an infinite number of blades neglects the fact
that towards the edge of the disc the air tends to flow around the tip of the blade sur-
face, thus reducing lift and hence power extraction at the tip. This is known as the tip
loss effect and is commonly addressed using a factor referred to as the tip-loss cor-
rection factor or tip-loss factor. There have been a number of different factors and
methods suggested within the literature [78][79][80]. It is suggested that the correc-
tion is applied to either the circulation, induced velocities, thrust coefficient or airfoil
coefficients [79].
The best known correction factor is Prandtl’s tip-loss function or correction factor
[81], given in equation (2.61), which corrects the induced velocities which affect the
loads indirectly. However, it leads to a slight over-prediction of the tip loads when
using two-dimensional data with BEM, which is due to a three-dimensional effect of
the blade geometry that is not taken into account [80]. Another popular tip correction
method is Glauert’s correction, given in equation (2.62). There are a number of semi-
empirical correction factors presented in the literature, such as Xu and Sankar [82],
Lindenburg [83] and Sant [84], which were based on the NREL phase VI experiment
[85] (§§3.2.1). The most popular correction was presented by Shen et al. [78] and
given in equation (2.63). Each tip-loss correction is applied in different ways. The
most commonly used method is to apply it to the definition equation of the disc given
in equation (2.51) or to the aerodynamics coefficients [73].

2 B R−r
FP = arccos exp − (2.61)
π 2 R sin(φ)

2 B R−r
FG = arccos exp − (2.62)
π 2 r sin(φ)

2 B R−r
FS = arccos exp −g (2.63)
π 2 R sin(φ)
g = exp[−c1 (Bλ − c2 )] + c3 (2.64)
Where FP , FG and FS are correction factors, c1,2,3 are constants with c1 , c2 fitted
from experimental data and c3 empirically such that (c1 , c2 , c3 ) = (0.125, 21, 0.1).
Implementation in ANSYS-CFX
This section describes how the Actuator Disc Method (ADM) can be implemented in
ANSYS-CFX. Similar to the porous disc there are two methods that can be used to
define the actuator disc location: either algebraically, which is described in this section
or geometrically (§§2.2.1).
It is possible to define the disc location algebraically within the software code; for
ANSYS-CFX this can be done using ANSYS-CFX Expression Language (CEL). The
main advantage of defining the disc within the code is that the rotor parameters can be
modified with ease, without the need to vary the domain or mesh. This is particularly
useful when running multiple simulations of different rotor sizes and/or different yaw
angles. One way to define an actuator disc using CEL is through the use of the built-in
step function. The step function is defined1 by equation (2.65). The function is defined
as 0.5 at x = 0 to avoid a discontinuous jump at the boundary. It should be noted that
1 There is some discussion as to the function’s value at x equals zero. The step function is commonly
defined as a discontinuous function with a value of either 1 or 0. This definition could be problematic
for CFD simulations creating a discontinuity at x = 0 and as such the value is defined as 0.5; as this is
the value obtained by the Fourier series of the step function; and is defined as such in this thesis and
within ANSYS-CFX.
the step function within ANSYS-CFX only works using dimensionless quantities.




1 f or x > 0



H(x) = 0 f or x < 0 (2.65)






0.5 i f f x = 0.
Using the step function it is possible to define the location of the disc using equa-
tions (2.66)-(2.69). These equations define a cylindrical volume with an internal value
of one and an external value of zero (also known as a cylindrical volumetric indicator
function [86]). This volume is used to define which nodes to apply the source terms
to, as described in section 3.5.1.
ℜ = RR × RT (2.66)
RR = H(R − Rr ) (2.67)
q
Rr = YR2 + XR2 (2.68)
RT = H(ZR ) × H(d − ZR ). (2.69)
Where ℜ is the function defining the volume of the rotor, RR is the function defining
the rotor area, Rr is the variable rotor radius, RT is the function defining the rotor
thickness, H(x) is the step function, R is the radius of the rotor, (X,Y, Z)R is the local
coordinate system of the rotor and d is the depth of the rotor.
The forces can be calculated within CEL or using external software [87] and im-
ported into ANSYS-CFX. If the forces are calculated within CEL then additional vari-
ables are required, specifically the tabulated variables CT , CL , CD , c and blade twist (ϕ).
It is possible to define these variables using ANSYS-CFX user functions and interpo-
lation functionality. Forces are applied in ANSYS-CFX by adding momentum source
terms. If the force is represented in cylindrical coordinates all three components are
defined as shown in equations (2.70)-(2.72):
Sz = DAD (CL cos φ +CD sin φ)ℜ (2.70)
Sθ = DAD (CL sin φ −CD cos φ)ℜ (2.71)
Sr = 0, (2.72)
where DAD is a force as described in subsection 3.5.1, CL is the lift coefficient, CD is

the drag coefficient, φ is the angle between Urel and rotor plane and ℜ is the function
defining the volume of the rotor shown in equation (2.66).
Implementation in ANSYS-Fluent
Within ANSYS-Fluent the disc can be defined in a similar manner to ANSYS-CFX.

The difference between the two is that ANSYS-Fluent does not use CEL to define ex-
pressions and functions. Instead ANSYS-Fluent uses its own User-Defined Functions
(UDFs). UDFs are written in C programming language and compiled into ANSYS-
Fluent [88]. UDFs utilize DEFINE macros and other predefined macros and functions
supplied by ANSYS-Fluent which access ANSYS-Fluent solver data and performs
other tasks [88]. It is possible to create a UDF using analogous expressions to those
detailed in section 2.2.2 to implement the ADM within ANSYS-Fluent. It should be
noted that to the author’s knowledge, there is no inbuilt step function, so it needs to be
defined as part of the UDF.
2.2.3 Actuator line method
The main problem/limitation of the actuator disc method is that it is unable to produce
tip vortices due to the way the forces are evenly distributed. The Actuator Line Method
(ALM) was created by Sørensen and Shen [89] to overcome this limitation.
The ALM works in the same way as the actuator disc except that instead of dis-
tributing the force over a stationary disc, it is defined as rotating lines, known as Ac-
tuator Lines (AL), using a smoothing function to avoid singularities. The ALM is
described by equation (2.74) where the smoothing function is commonly defined by
equation (2.75) as seen in the references [74][76][89].
1 2
DAL = ρU c(CL eL +CL eD ) (2.73)
2 rel
Fε = DAL ⊗ ης (2.74)

1 b
ης (b) = 3 3/2 exp − 2 (2.75)
ς π ς
Where ρ is the density of air, Urel is the relative velocity of the rotor, c is the chord
length, CL is the lift coefficient, CD is the drag coefficient, eL , eD are unit vectors used
to define the direction of lift and drag respectively, ης is a smoothing function, ς is a
constant that serves to adjust the strength of the function and b is the distance between
the grid point and the initial force point.
Implementation in ANSYS-CFX
The Actuator Line Method (ALM) works in the same way as the actuator disc except
that the force is distributed over rotating lines. Similarly to the disc, the lines can be
constructed either geometrically (as part of the domain) or algebraically. If the lines are
constructed geometrically either the domain or mesh must move accordingly. Although
a moving domain approach may be easier to implement, it is more computationally
expensive and has a lower accuracy [87] compared to defining the lines algebraically,
which reduces computational expense both in terms of mesh construction and solution
time.
As the ALM is a transient method, the time step needs to be considered. Troldborg
[74] showed that the movement of the blade tip during one time step should not exceed
one grid spacing. Therefore, the time step size is restricted by the tip speed ratio of
the rotor and grid spacing used at the tip. Only the tip is considered as it is the fastest
moving part of the blade.
Within ANSYS-CFX refining the disc volume into lines can be constructed using
CEL in a similar way to the ADM. The lines each need to be defined in their own
rotating coordinate system [87], shown in equations (2.76) and (2.77) with equations
(2.78) and (2.79) used to evaluate the blade position and rotational speed respectively.
c
YTransN = y − Rr sin(βn ) − sin(βn + 90) (2.76)
2
c
XTransN = x − Rr cos(βn ) − cos(βn + 90) (2.77)
2
βn = βi − β (2.78)
2π
β = Ωt (2.79)
60
Where n = {1, 2, 3} represents the blade number, Rr is the variable rotor radius, c
is the chord length (the offset distance from the leading edge being half this), βi is the
initial position of the blade,Ω is the rotational speed of the rotor and t is the time.
Using these coordinate systems it is possible to define the lines, shown in equation
(2.80) which is combined with expression (2.66) to create the volume of the blades.
This volume is used to distribute the forces over the blades in a similar way to the
smoothing function previously mentioned. DAL from equation (2.73), is averaged over
the volume of the blades, as shown in equation (2.81), to distribute the force per length.
q
Ln = H(c(r) − 2
YTransN 2
+ XTransN ) (2.80)
DAL
FAL = (2.81)
3dc(r)
Where n = {1, 2, 3} represents the blade number, c(r) is the chord length, H(x) is
the step function and d is the rotor depth.
The force per volume is defined using a momentum source in cylindrical coordi-
nates shown in equation (2.82)-(2.84) similar to the reference [89].
Sz = FAL (CL sin φ +CD cos φ)ℜ × Ln (2.82)
Sθ = FAL (CL sin φ −CD cos φ)ℜ × Ln (2.83)
Sr = 0 (2.84)
Where CL is the lift coefficient, CD is the drag coefficient and φ is the angle between
Urel and rotor plane, defined by equation (2.49) and ℜ is the function defining the
volume of the rotor shown is equation (2.66).
Implementation in ANSYS-Fluent
Within ANSYS-Fluent the lines can be defined in a similar manner to ANSYS-CFX.

ANSYS-Fluent uses User-Defined Functions (UDFs) as described in subsection 2.2.2.
It is possible to create a UDF using analogous expressions as detailed in subsection
2.2.3 to implement the ALM within ANSYS-Fluent.
Chapter 3
Literature review
This chapter consists of a literature review concentrating on work which considers

turbine wakes and their analysis. It begins with a brief review of basic turbine theory
including control volume analysis, followed by the main large scale experiments. The
next sections focus on numerical analysis ranging from semi-empirical models to full
CFD analysis and the rotor models used within them. For other reviews see references
[22, 67, 74, 90, 91, 92].
3.1 Theory
This section describes some of the theoretical work conducted into turbines and their
wakes. These theories consider idealized turbines to reduce their complexity, allowing
for limits and generalizations to be considered and calculated.
3.1.1 One-dimensional momentum theory
The first wind turbine theories were adapted from work regarding propellers of ships
and helicopters. The first people to look at wind turbines were Rankine [93] and Froude
[94] in the late 19th century. Early approaches used control volume analysis, where
80
CHAPTER 3. LITERATURE REVIEW 81
a defined volume of space (a control volume) is analysed, commonly by calculating

the flow at the boundaries. This is commonly referred to as the one-dimensional (1D)
momentum theory within the wind energy literature. It assumes a control volume made
up of stream tubes. The flow is parallel in these tubes. The wind turbine is represented
as an actuator disc, which intersects the tubes and is a discontinuous pressure/velocity
step, as shown in Figure 3.1 from reference [95]. This model relies on a number of
assumptions:
• Homogeneous, incompressible, steady fluid flow
• No frictional drag
• Static pressure far up and down stream is equal to the undisturbed ambient pres-
sure
• Uniform disc and thus, a uniform thrust and infinite number of blades
• Only axial flow and thus, a non-rotational wake
Using these assumptions and the Bernolulli equation away from the disc, it is possible
to calculate a number of variables including the power, power coefficient, thrust and
thrust coefficient given by equations (3.1)-(3.4) respectively. Power is a measure of
the energy generated by the turbine and the thrust is the axial force applied by the
turbine to the flow. For more details on this method and its derivation see references
[10, 11, 95].
1
P = ρAr u3 4a(1 − a)2 (3.1)
2
Power
CP = 1
= 4a(1 − a)2 (3.2)
2 ρAr u
3
1
T = ρAr u2 4a(1 − a) (3.3)
2
Figure 3.1: 1D momentum theory
CT = 4a(1 − a) (3.4)
u
a = 1− (3.5)
U0
Where P is the power, ρ is the density, Ar is the rotor sweep area, u is the axial
velocity, a is the axial induction factor, CP is the power coefficient, T is the thrust, CT
is the thrust coefficient and U0 is the free stream velocity.
3.1.2 Betz limit
The maximum power coefficient (C p ) can be derived by taking the derivative and set-
1
ting it equal to 0. Equation (3.6) shows the derivative after some algebra; giving a = 3
16
and a maximum C p = 27 = 0.5925, showing that the maximum efficiency of a wind
turbine is 59.3%. CT has a maximum of 1 when a = 0.5 and downstream velocity
equals zero and at maximum power CT = 89 .
dC p
= 4(1 − a)(1 − 3a) = 0 (3.6)
da
This is the theoretical maximum and is commonly known as the Betz limit after the
work done by Betz in 1920 [96], although it has been suggested that this is incorrect
[97][98][99]. Bergey [97] suggested that it should be known as the Lanchester-Betz
limit after the work done by Lanchester in 1915 [100]. Kuik [98] has further suggested
that it should be referred to as the Lanchester-Betz-Joukowsky limit after the work by
Joukowsky [101]. More recently it has been argued by Okulov & Kuik [99] that it
should be referred to as the Betz-Joukowsky limit as Lanchester did not directly accept
16 2
Froude’s work [94] and defined the limit as 27 Q where 1 < Q < 2 and suggested
Q ≈ 1.5. In practice there are many factors not considered in this limit such as wake
rotation, finite number of blades and aerodynamic drag.
3.1.3 Ideal Horizontal Axis Wind Turbine (HAWT)
The 1D momentum theory, described in section 3.1.1, was in one dimension and in-
cluded no wake rotation. This is unrealistic as all wind turbines produce a rotational
wake, which rotates in the opposite direction to the blades due to the generation of
torque. The rotation of the wake causes less energy to be extracted by the turbine and
hence, more kinetic energy remains in the wake. This added rotation to the wake is
expressed using the angular or tangential flow induction factor, a0 . The relationships
1−3a
between a and a0 can be expressed as a0 = 4a−1 . Including the rotation allows the tip
ΩR
speeds ratio to be defined as λ = u , where Ω is the rotational speed of the blade.
Actual wind turbines are not discs, but rotors made of aerofoil blades which move by
generating a pressure difference on the upper and lower surfaces, creating lift. A lift
force is defined as a force perpendicular to the oncoming wind as shown in Figure 3.2
Figure 3.2: Aerofoil forces
1 and the 2D lift coefficient is defined in equation (3.7).
L
Li f t f orce per length
Cl = = 1 l2 (3.7)
Dynamic f orce 2 ρu c
Where ρ is the density, u is the velocity, c is the chord length, L is the lift force and
l is the span.
Another important aspect of aerofoils is drag, which is defined as parallel to the
direction of the wind as shown in Figure 3.2 1 . It is caused by the viscous friction
forces on the surfaces as well as the pressure around the blade. The 2D coefficient is
defined by equation (3.8).
D
Drag f orce per length
Cd = = 1 l2 (3.8)
Dynamic f orce 2 ρu c
Where ρ is the density, u is the velocity, c is the chord length, D is the drag force and l
is the span.
1 http://www.free-online-private-pilot-ground-school.com/aerodynamics.html
3.2 Previous experimental work
This section will outline the main physical experiments that have been undertaken. A
key problem with field experiments is that it is impossible to know the entire flow field
as wind is very unsteady, continually changing in direction, magnitude and turbulence
intensity. The surrounding terrain also has an effect on any results. The alternative is to
use wind tunnel tests where the flow conditions and terrain can be defined. However,
these also experience problems, the main one being scale (as wind tunnels are far too
small for a full wind turbine). The geometry is easily scalable but the flow is not,
specifically the Reynolds number (discussed in section 6.3). This means most wind
tunnel experiments use Reynolds numbers far lower than real world conditions.
A large number of experiments have been conducted, but they have generally been
too specific in nature, concentrating on their own individual aims and objectives, taking
only limited readings and results [17]. This has a limiting effect on the ability of others
to make accurate comparisons. Although this thesis will not discuss in detail previous
experimental work there are a number of good review articles available in the literature.
The main articles are Vermeer et al. [17] and Crespo et al. [102] with a smaller section
in [74]. Vermeer et al. [17] goes into detail about various experiments, including both
near and far wake studies. Crespo et al. [102] focuses predominately on the far wake
and covers work from 1963 until 1999. Troldborg [74] covers more recent work and
concentrates on the far wake, but does have a section on the near wake as well.
There are two large scale experiments that have recently been conducted that must
be noted, as follows.
3.2.1 NREL Unsteady Aerodynamics Experiment (UAE) project
The first and largest of these was the National Renewable Energy Laboratory (NREL)
Unsteady Aerodynamics Experiment (UAE) project [103]. This project was made up
of a number of stages including an extensive 10-year field experiment phase [104].

One of the most important phases was phase VI, which consisted of placing the UAE
10m diameter wind turbine in the 24.4m×36.6m (80ft×120ft) NASA-Ames wind tun-
nel [85]. The experiment provided accurate and reliable measurements with a high
spatial and temporal resolution for approximately 1700 different wind turbine condi-
tions. This experiment was critical for a number of reasons; namely it was the first
full scale wind tunnel experiment to be conducted, providing a key source due to the
quantity and quality of data collected. This information has already been used in a
large number of comparison studies, most notable in its own blind comparison study
[105]. This comparison was used to test uncertainties in wind turbine modelling, with
all participants given the same data about geometry and inflow conditions. All partic-
ipants submitted their results before NREL compared and published the results. The
models were divided into 4 categories: performance codes, aeroelastic codes, wake
codes and CFD. This study showed a wide range of results, even for the supposedly
simple situation with no yaw or stall. Yaw is the turbine’s orientation where zero yaw
describes a turbine perpendicular to the flow. Stalled flows feature reduced lift, caused
by significant separation of the air from the blade surface. Even with no yaw or stall,
the wind turbine power was predicted ranging from 25%-175% of the measured value
[105]. This increased to 30%-275% in stall conditions [105]. Although there was no
obvious loser in the comparison (as all models showed inconsistencies) there was one
outstanding model which was a CFD model, although its application was limited to
only none yaw conditions.
3.2.2 Model EXperiments In COntrolled conditions (MEXICO)
project
The EU FP5 project Model EXperiments In COntrolled conditions commonly known

as the MEXICO project [106] used a 4.5m diameter fully instrumented wind turbine
which was setup in the 9.5×9.5m test section of the Large Scale Low Speed Facility
(LLF) facility of the German-Dutch Wind tunnel Organization (DNW). The key thing
about this experiment was the extensive flow field mapping using the stereo Particle
Image Velocimetry (PIV) technique. Its aim was to create a database of measurements
to validate and/or improve design and analysis models.
The data obtained in the MEXICO project has been used for different numeri-
cal studies featuring full rotor models [107][108][109][110][111], lifting line models
[112][113], actuator surface models [114] and actuator line models [115][116] as well
as comparative studies [117][118][119]. Carrión et al. [111] produced a clear table
showing code, solver and geometry utilized in previous numerical studies. All full
rotor CFD model simulations were conducted without inclusion of any wind tunnel
geometry. Bechmann et al. [108] and Sørensen et al. [109] modelled the whole ro-
tor and found good agreement with respect to the experimental data for axial flow.
Sørensen et al. [109] also investigated yaw flow and found large discrepancy in the
wake which was attributed to the exclusion of the nacelle and its wake in their simula-
tions. Lutz et al. [110] conducted simulations of one blade as well as the whole rotor
featuring the tower and nacelle, demonstrating their effect on the wake. Stoevesandt et
al. [107] also included the nacelle although they used a shortened version in their sim-
ulations. Grasso & Garrel [112] and Micallef et al. [113] ran lifting line simulations
on axial and yaw flow conditions, finding good agreement with the experiential data.
Micallef et al. [113] found their method was quite sensitive to input aerofoil data and
that the BEM was unsuitable for this purpose. Shen et al. [115][116] used the ALM
and LES to simulate the MEXICO rotor within the wind tunnel and in free air. Two sets
of aerofoil data as well as different turbulence models were investigated. In general,
good agreement was found although it was concluded that the original 2D aerofoil data
were not suitable and that a larger thrust than the pressure measurements predicted is
required. Réthoré et al. [120][121] compared results from actuator disc and full rotor
simulations with and without the tunnel to show the significance of the wind tunnel.
They found both the full rotor and actuator disc model could not match both the blade
loading and velocity measurements at the same time which is in agreement with Shen
et al. [115][116].
3.2.3 Mexnext project
These two large scale experiments produced a large quantity of data, far too great for
any one organization to process. Thus the International Energy Agency (IEA) Wind
Task 29 Mexnext project2 was created. Its purpose was to analyse and evaluate the
measurements and information from both the MEXICO and NASA-Ames projects.
Mexnext (Mexnext-I) ran for 3 years from June 1 2008 until the end of 2011 [117]. It
found that CFD models predicted the loads better than engineering codes, which under
predicted the loads for stalled conditions [117][118][119].
3.3 Previous numerical work
This section gives an overview of numerical methods which are commonly known as
Engineering methods. These methods have been developed from and rely on empiri-
cal data (for input and/or calibration) and are popular, as they do not require a lot of
computational resources.
2 www.mexnext.org
3.3.1 Blade models
One of the first engineering models was called the Blade Element theory or BE theory.
It analysed the forces at a section of the blade as a function of its 2D geometry. This
idea was later combined with 1D momentum theory into the strip theory, better known
as the Blade Element Momentum (BEM) theory. It is generally attributed to Glauert in
1935 [122]. The BEM theory assumes the blade is divided into small blade elements,
which are all independent of each other (no radial motion). These elements all use 2D
aerodynamic aerofoil data, calculated using local flow conditions. The momentum side
of the theory works by assuming that the pressure/momentum loss in the rotor plane is
caused by the airflow passing through the rotor plane. This theory is used to determine
the aerodynamic forces and induced velocity near the rotor in an iterative manner. It is
very well known and widely used even today due to its accuracy and speed.
3.3.2 Vortex wake methods
Vortex wake models calculate the flow field as vortices and are commonly coupled
with a lifting method to represent the blades. There are a number of lifting methods
that have been used to model a wind turbine rotor: the lifting line, the lifting surface
and the lifting panel. These methods are only briefly mentioned, for more detail see
the references [17][22][91]. Each of these methods model the blade as some arrange-
ment of distributed vortex elements. The lifting line method is the simplest model,
representing the blade as bound vortex lines along which the local vortex strength is
related to the local angle of attack. The lifting surface method represents the blades as
a surface providing more detail. A variation on this method is the lifting panel method,
which has been used extensively in the aircraft industry with only a few wind turbine
applications. The 3D nature of the lifting surface and panel methods is an advantage
over the line method, although this does come at a computational expense.
There are two types of vortex wake method: prescribed wake and free wake. The
prescribed wake approach relies on measured or empirical input to define the initial
conditions and limits. The free wake approach allows the vortex elements to con-
vect and deform freely under the influence of the velocity field. This added freedom
increases the computational expense and can lead to instability and convergence prob-
lems. The vortex wake method is good at predicting unsteady flow and yaw misalign-
ment.
3.3.3 Wake models and computer programs
There have been a large number of wake specific models and computer programs, of
which many have already been reviewed [17][22][102] and only the key ones will be
noted here. Wake models focus on the far wake where the geometry and features of
the rotor are only indirectly felt.
Wind-farms were first modelled using roughness elements and wake superposition
techniques featuring both linear [123] and squared linear [9] superposition. These
methods were very simple, using wake data from experiments and individual wake
models. The first individual wake models were the kinematic models which were
based on self similar velocity deficit profiles from experimental and theoretical work.
From these models more specific computer programs were developed. The program
EVMOD [124][125] was developed using an eddy-viscosity method, modified by an
empirical-based function to account for the lack of equilibrium. Garrad Hassan later
developed this method for wind-farm applications into EVFARM [126].
The computer program UPMWAKE was developed by Crespo et al [127][128]
using a parabolic approximation. This program does not rely on the self-similarity
assumption and allows the wake to be non-uniform, providing a more realistic predic-
tion over EVMOD. An elliptical version of this program was developed by Crespo et
al. [129] but no fundamental difference was found. UPMWAKE has been extended
to wind-farm applications in the UPMPARK program [130] and to account for the
anisotropic effects of turbulence in the UPM-ANIWAKE program [131].
The problem with all engineering models is they employ a large number of as-
sumptions and rely on empirical data for their development, calibration and input. Due
to their reliance on empirical data these models are often limited by the data they have
available to them. This input is also the main source of error in many calculations and
situations. To reduce this error it requires using as few assumptions as possible and
limiting the reliance on empirical data.
3.4 Previous computational work
The computational work detailed in this section differs from section 3.3 as it relates
to Computational Fluid Dynamics (CFD). This section outlines the types of software
codes currently being employed, followed by the methods used to solve the Navier-
Stokes equations and some of their limitations. The section ends discussing the main
methods used to model the rotor.
3.4.1 Software/codes
Due to the complex nature of the Navier-Stokes equations and the numerical methods
developed to solve them, fluid flow simulations are commonly run on High Perfor-
mance Computers (HPC) using specially developed CFD codes 3 . Some of these
codes are multipurpose codes, designed to solve the Navier-Stokes equations for a
wide range of problems. However, others are designed for specific applications such
as FLOWer [110] and Helicopter Multi-Block (HMB) [111]; or they are developed for
specific methods such as EllipSys3D which was written purely for curvilinear grids
3 http://www.cfd-online.com/Wiki/Codes
[132]. CFD codes are very difficult to write and require verification and validation
to determine their reliability and accuracy. Due to the time and effort required many
are considered in-house codes and are not available on the open market such as El-
lipSys3D. That said, there are still many codes commercially available 3 and some
open source codes available such as OpenFOAM4 . There are many advantages and
disadvantages of using both in-house and commercially available codes. One main
advantage of in-house codes is greater flexibility to change the source code compared
to commercial codes, although this is offset by a larger number of users which can
provide greater technical support.
Key codes used within the wind turbine wake sector are currently:
• EllipSys3D - an in-house code that was developed by Michelsen [133, 134] and
Sørensen [135] and is a general purpose, multi block code with a finite volume
discretization of the Navier-Stokes equations in general curvilinear coordinates.
The variables are evaluated at the cell centers to enable computations using com-
plex geometrical meshes. It is currently the most prominent code being used
within the literature to model wind turbine wakes with multiple PhD theses re-
lying on it [31][74][75][132][136].
• Ansys-CFX [65] - a commercially available general-purpose code which has

been applied to a broad range of industries and applications. It has been used to
simulate single blades [137][138] full rotors [139] and porous discs [140][141].
More recently it has been used to model wakes using actuator methods [1][2][3][87].
• Ansys-Fluent [66] - a commercially available general-purpose finite volume code

used to model turbulence, heat transfer, and reactions and has several solvers and
discretization types available. It has been used for several wind turbine studies
4 http://www.openfoam.org/
featuring single blades [142][143][144], full rotors [145] as well as actuator discs
[146][147].
• OpenFOAM - a free, open source CFD software package produced by OpenCFD

Ltd. It offers users complete freedom to customize and extend its existing func-
tionality. It has been used to simulate individual turbines such as the MEXICO
rotor [107] as well as wind farms [148].
There are other codes, but these are only being used by individual researches or insti-
tutions.
The Reynolds-Averaged Navier-Stokes (RANS) method has been the most commonly
used approach for wind turbines. This is mainly due to its speed of calculations as
previous mentioned (§§1.5.3 and §§2.1.1). There are different turbulence models used,
although the most commonly used in the literature for turbines wakes are the 2 equation
models k-ε, k-ω and k-ω SST. RSMs are rarely used due to their higher computational
cost [91]. When people first started applying CFD to turbine problems this was the
only approach available to them due to limited computer power, which today is no
longer the case.
The k-ε model
The k-ε model [28][45], described in subsection 2.1.1, is one of the most popular two
equation models for many CFD applications. The k-ε model predicts a short wake
which dissipates too quickly [149], as it over predicts wake recovery [31]. To combat
this problem a large number of modifications to the standard model have been de-
veloped. Réthoré [31] compared four different models, the El Kasmi-Masson model
[150], the Adverse Pressure Gradient (APG) model, the realizability model and a wind
turbine canopy model. It was found that although some improvement was evident,
none of them were adapted enough for wind turbine wakes [31]. Cabezón et al [151]
compared three isotropic variations of the k-ε model with the anisotropic Reynolds
Stress Method (RSM) (mentioned in subsections 2.1.1 and 2.1.1) as well as UPM-
PARK and UPMANIWAKE described in section 3.3.3. It was found that the standard
k-ε model under-estimates the wake deficit and although the modified models were
acceptable in the far wake they needed to be tuned for different cases [151].
The k-ω model
The k-ω model, detailed in Wilcox [24] and described in section 2.1.1, is an alternative
model and performs better under adverse pressure gradients and flows featuring sep-
aration. Given this the k-ω model does not outperform the k-ε model in any obvious
way as both models under predict wake effects [152].
The k-ω SST model
The k-ω shear stress transport model normally written k-ω SST or SST k-ω model
is a combination of the standard k-ω and k-ε models defined by Menter [29][49] and
discussed in section 2.1.1. It uses the k-ω definition in the inner part of the boundary
layer and gradually switches to the k-ε definition outside the boundary layer. This for-
mulation is becoming the main turbulence model used today particularly for modelling
turbine wakes [67].
3.4.3 Problems
The k-ω SST model also has problems with accuracy similar to the k-ω and k-ε models
it is based on. The problem with these models is that they are based on the Boussinesq
hypothesis [30], also called the Boussinesq eddy-viscosity assumption or the Boussi-
nesq approximation which is discussed in section 2.1.1. The Boussinesq hypothesis is
considered inadequate and invalid for flows featuring sudden changes in mean strain
rate, rotation and 3D effects [24] all of which feature in wind turbine wakes. Réthoré
[31] makes a detailed study of the k-ε model describing its limitations and associated
failing in wind turbine wakes and states that it is invalid for the majority of the wake
region up to 10 diameters downstream. Where the Boussinesq hypothesis is invalid the
turbulence models over predict the eddy-viscosity, reducing the length of the wake.
It is still possible to use RANS methods without relying on the Boussinesq Hy-
pothesis by using the Reynolds Stress Model (RSM), discussed in section 2.1.1. RSMs
model all components of the Reynolds stress tensor, making it much more resource
intensive. It is able to model the anisotropic nature of the wake more accurately then
the other RANS turbulence models within the wake [151]. However, RSMs rely on
similar assumptions and without the eddy-viscosity term for stability it can suffer from
numerical problems. Although without using the eddy viscosity RSMs are able to
model anisotropic turbulence.
The Large Eddy Simulation (LES) method (detailed in section 2.1.2) is being used
more and more for wind turbine applications due to the increase in computer power and
availability. LES is also able to handle unsteady anisotropic turbulent flow much better
than RANS, but this does come with a higher computational expense. As previous
discussed, eddy-viscosity turbulence models are used as part of LES method, which
means that the problems of the Boussinesq hypothesis are still encountered. Réthoré
[31] showed that if cell size is taken small enough then the error made by the eddy-
viscosity concept is negligible for wind turbine wakes.
As LES is computationally expensive near solid boundaries, DES of wind turbines

is promising. Johansen et al. [153] performed a DES of a single NREL phase VI blade
showing good agreement with experimental data.
3.4.5 Boundary conditions
A problem with all computational simulations is the inflow parameters. Wind turbines
experience turbulent shear inflow which varies significantly with time. Early simula-
tions used uniform, laminar inflows to reduce computational cost. More recently it
has been shown how turbulence inflow [154] and shear inflow (atmospheric boundary
layer) [155][156] affect flow behind the wind turbine.
For LES there are two ways of creating a turbulent inflow: using a synthesized
inlet and using the precursor simulations method. The synthesized inlet as detailed
in Mann [157][158] is based on creating a turbulent field by influencing the current
flow. The advantage of this is that parameters such as turbulence length scale can be
directly specified. This method was used by Troldborg [74][154] to create a turbulent
inlet. It was inserted in one rotor diameter before the rotor to prevent decay. The flow
was influenced using body forces in the momentum equation following the work by
Mikkelsen et al. [159].
The Precursor simulation method is achieved by defining a separate simulation,
which is then used as the inflow for the next simulation. The advantage of this method
is that the ground effect is taken into account. However, it cannot be easily changed
and adapted once created and it is computationally very expensive. This method has
been used in a number of papers including references [23][148][160].
When solid boundaries are not included in the simulation directly, which is nor-
mally done through an actuator method, it has been shown that the solution does not
change much when the Reynolds number is higher than 1000 [22][89][161][162]. It
was also observed that the behaviour of the wake maybe insensitive to the blades’
Reynolds number [163], which is important for many rotor models that do not con-
sider the blade geometry.
3.5 Rotor modelling
There are many ways to represent the rotor of a turbine. These range in complexity and
detail starting by representing the rotor as a simple constant momentum sink, to the full
rotor geometry including the tower and nacelle. The most obvious way to simulate a
wind turbine is to construct a model as close to real life as possible by using the full
geometry of the blade. These models are known as full rotor or direct rotor models.
The main disadvantage of this method is that it needs to resolve the blade boundary
layer which requires a fine mesh and large amount of computational power. Due to the
computational cost it is common to represent the rotor instead of simulating it directly.
The simplest approach to rotor modelling is representing the wind turbine as a
momentum sink which was the approach used in early numerical work as mentioned
in section 3.3.3. This is generally done for far wake studies where the rotor geometry
is less important, such as modelling wind-farms and other large features. Calaf et al.
[70] represented wind turbines as drag discs and used LES to study the wind turbine
array boundary layer. Prospathopoulos et al. [164] represented rotors as porous discs
on the top of two different Gaussian hills and showed how the terrain affects the wakes
and velocity recovery. They found that on a Quasi-3D hill the wake was still significant
up to 40 diameters downwind of the turbine, which is much longer than for flat terrain.
Rados et al. [152] represented wind turbines as momentum sinks when comparing two
codes using k-ε and k-ω turbulent models and found no significant difference between
the codes.
The most common way to represent the rotor is by using actuator methods, where
the blades are represented by body forces in the shape of either a disc, line or surface.
Actuator methods have been used for many years in modelling different rotors such
as helicopter rotors [165], propellers [166] and turbines [67]. This approach reduces
computational expense and mesh requirements as the blades are not represented by
solid bodies. Actuator methods take hours to complete compared to days or weeks for
full rotor simulations [167]. However, the influence of the 3D geometry of the blade is
not accounted for. There are two main types used: the actuator disc (AD) and actuator
line (AL). There is also an actuator surface (AS) method, but this is a newer method
and has not been developed as much.
3.5.1 Actuator disc method
The Actuator Disc Method (ADM) represents a time averaged rotor using forces dis-
tributed over a disc with the same dimensions and volume as the swept volume of the
rotor. There are a number of variations of the ADM within the literature including
constantly and variably loaded discs, as well as variations with and without rotation
[67][68][69].
Réthoré et al. [168] compared the AD to Conway’s analytical solution [169] and
found good agreement. Réthoré & Sørensen [170] further compared it to three more
analytical solutions showing good agreement and suggesting a method to reduce com-
putational costs. The actuator disc method has been used to model stand alone wind
turbines [75][136][161] and multiple wind turbines [136][171][172]. It has also been
used to model terrain effects, Makridis & Chick studied two wakes over a Gaussian
hill [146] as well as flow over complex terrain in a coastal region [147]. Both stud-
ies were conducted using ANSYS-Fluent and were compared with experimental data
and WAsP. They found CFD predictions were within 5% of the measurements while
WAsP was as much as 15% out. Alinot & Masson [173] used the AD to show the wake
structure in the atmospheric boundary layer using a modified k-ε turbulence model to
simulate the wind turbine operating in different types of stratification. They found that
a stable atmosphere produced the lowest power output and the unsteady atmosphere
condition produced the highest power output. Ivanell [136] used 20 AD combined
with periodic boundary conditions to simulate the Horn Revs wind-farm, showing
good agreement with measurements but the power production was found to depend
on the turbulence levels and the inflow angle [136][172].
3.5.2 Actuator line method
The main problem/limitation of the actuator disc method is that it is unable to produce
tip vortices due to the way the forces are evenly distributed. The Actuator Line Method
(ALM) was created by Sørensen & Shen [89] to overcome this limitation. The ALM
works in the same way as the actuator disc except that instead of distributing the force
over a stationary disc, it is defined as rotating lines, known as Actuator Lines (AL),
using a smoothing function to avoid singularities.
The ALM has been used for a number of studies including detailed studies by
Troldborg [74] and Mikkelsen [75]. Troldborg [74][154] used the AL combined with
LES to study how different inflow conditions affect the wake of a single wind turbine.
Troldborg [74] found that shear and ground presence affects the axisymmetric structure
of the wake development and observed minimal wake expansion vertically down. It
was also observed that with a uniform laminar inflow, the wake still becomes turbulent
especially at low velocities. Ivanell [136] used the AL method to do a stability study
on the wake structure.
The AL model has been useful in modelling turbine interactions, specifically inter-
actions of a row of three wind turbines [74][174]. Troldborg [74] showed how laminar
and turbulent inflow cases were very different. In the turbulent case the wake was un-
stable and broke up before the second turbine. However, in the laminar case this does
not happen, resulting in higher velocities at the second turbine and creating higher yaw
loads for the third turbine. Due to the transient nature of the ALM it is commonly used
with LES [74][77][154], although it has also been used with DNS [162]. Ivanell et al.
[162] used DNS on a third of the wake to provide detailed information on the wake
structure without ambient turbulence.
The ALM has been used in a number of comparison based papers to evaluate the
ADM [75][77][175]. Mikkelsen [75] compared the AD, AL and BEM methods at
simulating flow through coned rotors and yaw flow. It was found that the AD and AL
were able to model the coned and yawed rotor although the AD was found to be poor
at high yaw, with the BEM method being poor at both.
3.5.3 Actuator surface
The actuator line method was further extended into the actuator surface method by
Shen et al. [176][177]. The surface method represents a blade as a plate surface,
instead of a single line over which the forces are distributed. This requires more infor-
mation. Instead of the sectional lift and drag coefficients cl and cd used by the AD and
AL, the surface method requires knowledge of pressure and skin friction. An Actuator
surface hybrid model has been developed which replaces the blades with their mean
surfaces and a pressure jump [176][177][178]. The actuator surface model has been
used to model the NREL rotor [179][180] and MEXICO rotor [114].
3.5.4 Full Rotor
The full rotor or direct rotor method, directly models the rotor and the blades exactly.
The inclusion of the blade geometry means that compressible effects may be important
at the blades tips which could require that the compressible Navier-Stokes equations
be solved. They are very accurate, but extremely computationally expensive as the
inclusion of the blade geometry requires a fine moving/sliding [145] or overset mesh
[132][155] to capture the blade boundary layer.
There have been a number of studies using this technique. The first was done by
Sørensen & Hansen [181] and used the k-ω SST model. To reduce the computational
expense of this method some studies only simulate one blade [137][138][142][143][144][153].
One advantage of full rotor studies is they allow for the inclusion and investigation of
features such as the tower and nacelle, as well as any interactions between them and
the rotor. A number of studies have made use of this to investigate the interactions
between the rotor and the tower [132][182] and nacelle [145]. Zahle [132][182] con-
ducted a study of the NREL phase VI rotor both with and without the tower, and found
good agreement with the experimental data. He was able to show the affect of the rotor
on the tower’s shedding frequency. Zahle & Sørensen [156] have used the full rotor
and an overset grid method to investigate the effect of including ground effects and
shear inflow on a wind turbine. It was shown that including the ground effects without
a shear inflow increased the power prediction by 0.4%; adding a shear inflow increased
this further to 5.5% [156].
The main purpose of the full rotor model studies are generally towards comparison
and evaluation of simpler models. Réthoré et al. [167] compared the actuator disc and
full rotor models with the k-ω SST model in uniform and shear inflows. The study
showed the computational expense of using the different techniques, with the actuator
disc taking hours to converge whereas the full rotor model took months to converge.
It also showed that (at least for uniform inflow conditions) the difference between
the models is relatively small and does not justify the extra computational expense.
Réthoré et al. [175] went one step further and compared the AD and AL with a full
rotor simulation in the near wake region. It was found that there were small structures
created in the full rotor simulation that have a significant impact on the generation of
turbulence, which the AD and AL are unable to predict. However, this study did not
use a turbulent inflow which may have affected the overall impact.
Due to the large computational cost associated with the full rotor model, simula-
tions using this approach are commonly restricted to only one blade or at most, a single
rotor. However, there is a project to conduct a full wind farm simulated using the full
rotor model [183] although to simulate 48 turbines requires a 3840 core supercomputer
and approximately 300 million cells.
Chapter 4
Methods
This chapter details the five CFD studies carried out in this thesis to achieve the ob-
jectives set out in section 1.2. The studies increase in complexity and novelty. The
first four studies in sections 4.1−4.4 were conducted to verify (benchmark) and de-
velop the modelling methods described in section 2.2. They feature comparisons with
previous work including experimental and numerical data. Finally, section 4.5 details
novel work investigating the effect rotor diameter has on wind turbine interactions and
power generation.
4.1 One-dimensional momentum theory study
The first study featured comparison with the best known theoretical model of flow
about a turbine, namely the 1D momentum theory as discussed in section 3.1.1 and
illustrated in Figure 4.1. The theory and hence this study featured a uniform disc with-
out rotational effects. Only the porous disc and uniform actuator disc model, defined
with the thrust coefficient (CT ), were used detailed in subsection 4.2.3.
As comparison does not require any form of special domain or mesh construction,
the domain configuration defined in section 4.2 was used. The CFD simulations were
103
CHAPTER 4. METHODS 104
Figure 4.1: 1D momentum theory
conducted using the steady-state Reynolds-Averaged Navier-Stokes (RANS) equations

together with the k-ω SST turbulence model [29]. ANSYS-CFX was used to calculate
the solution using a mixed hybrid mesh, constructed with 6.2 million cells. Such a
high mesh density was used for this study as the simulation was conducted as part of
the porous disc study, described in Section 4.2.
4.2 Porous disc study
This study was a comparison with previous numerical [140][141] and experimental
[184][185] studies. The previous experiments [184][185] were conducted in a water
channel measuring 21m × 1.37m with a depth of 0.3m. Three 0.1m diameter discs of
various porosities were placed into the channel. The different porosities were used to
represent turbines with different thrust coefficients (CT ), which were measured using
a pivot arm mounted onto a load cell. The water velocities were measured at various
locations using an Acoustic Doppler Velocimeter (ADV) at a sample rate of 50Hz
and the data was averaged over 3 minutes. The previous numerical studies chosen
for comparison, are described in references [140][141]. These used ANSYS-CFX to

reproduce analogous experimental data [185]. These studies were chosen because they
were conducted using the same software as in the work presented in this thesis, hence
providing a benchmark with which to verify current modelling methods.
The numerical simulations in this study were calculated using ANSYS-Workbench,
specifically ANSYS-CFX using the steady-state solution of the Reynolds-Averaged
Navier-Stokes (RANS) equations [18], together with the k-ω SST turbulence model
[29]. This model was chosen over the k-ε model based on the literature and some
preliminary simulations which showed that the k-ω SST model performs better in flows
featuring adverse pressure gradients [24] in terms of the accuracy to predict the flow
properties. The k-ω SST model was also used in the benchmark studies [140][141].
The model domain was defined using the dimensions of the experimental channel
set-up [185]. It featured a 2m long inlet, a 3m outlet and a 0.3m deep-water column
along with a 0.1m diameter disc with a thickness of 0.001m at the centre. The flow
was assumed to be symmetrical. Allowing a symmetry plane to be created through
the centre of the disc, dividing the domain in half to create a width of just 0.685m
as opposed to the 1.37m width of the experimental channel; therefore reducing the
computational expense. As part of this study, three discs were simulated with two
different sets of boundary conditions to represent a channel and a duct. The disc was
defined with a diameter of 0.1m, a thickness of 0.001m and a uniform momentum loss
across the disc in the longitudinal (z-)direction. As part of this study three discs were
simulated with two different sets of boundary conditions to represent a channel and a
duct, each with two different inlets totaling 12 simulations.
4.2.1 Mesh
In this study a mixed hybrid mesh was constructed consisting of various mesh densities
ranging from 2.4 to 6.2 million cells. The majority of the domain was constructed out
of tetrahedral cells with a structured inflation zone of wedge cells at the boundary
of the floor and symmetry plane. The data presented corresponds to a mesh density of
approximately 6.2 million cells and 2 million nodes unless otherwise stated. Figure 4.2
shows the velocity profiles of various mesh densities along the centre line behind the
disc and at 7 diameters (7D) downstream of the disc. There was very little difference
between the predictions of the four different mesh densities showing little advantage
in refining the mesh. Figure 4.2a demonstrates a realistic velocity recovery beyond the
peak velocity drop just before z/2R = 5. Figure 4.2b shows that the main differences
between the different mesh densities are within the floor boundary layer and at the peak
velocity deficit.
The inlet velocity was defined by equation (4.1) in the same manner as the numerical
study [141] and based on the empirical data [185]:
ywU ∗
Uinlet = 2.5U ∗ ln( )+A (4.1)
ν
where Uinlet is the inlet velocity across the width of the domain, U ∗ is the friction ve-
locity, yw is the depth of the water, ν is the kinematic viscosity and A is a boundary
layer model constant. Curve fitting methods were used to define U ∗ and A. The nu-
merical paper [141] used values of U ∗ = 0.00787m/s and A = 0.197m/s which were
also used in this study. Figure 4.3 shows the inlet velocity used in this work and the
experimental study [185], normalized with a free stream velocity (U0 ) of 0.331m/s for
(a)
(b)
Figure 4.2: Normalized velocity profiles showing different mesh densities (a) along
the centre line and (b) 7 diameters downstream of the disc
the experimental study [185] and U0 = 0.33m/s in this work. The vertical height was
also normalized with the diameter of the disc (2R).
Figure 4.3: Normalized velocity at the inlet of this study (solid line) and experimental
data (o) [184]
The turbulence intensity, which is defined by equation (4.2), was described in two
different ways to define two simulations referred to in this study as inlet 1 and inlet 2.
Both inlets were set with a turbulence intensity of 5% at the inlet to produce agreement
with the experimental data [185] for y/2R > 0.5. The difference between the inlets is
that inlet 2 was also defined with a length scale of 0.3 (height of the domain). Both
these approaches are different to the reference [141] which defined the turbulent kinetic
energy and eddy dissipation.
u0
I≡ (4.2)
Ū
r
0 2
u ≡ k (4.3)
3
q
Ū ≡ Ux2 +Uy2 +Uz2 (4.4)
In equations (4.2)-(4.4) I is the turbulence intensity, u0 is the root-mean-square of

the turbulent velocity fluctuations, Ū is the mean velocity and k is the turbulent kinetic
energy and Ui = [x, y, z] is the velocity in the x, y, z directions.
The outlet was defined as a static pressure outlet with a relative pressure of 0 Pa.
The floor and far side of the domain was defined as a non-slip wall. In this study two
separate benchmark models were produced featuring different boundary conditions at
the top or roof of the domain; the first featuring an opening creating a channel and the
second featuring a non-slip wall creating a duct. These boundary condition sets were
analogues of those in the references [140][141]. Although a free-surface approach
could be considered more suitable, as the experiment was carried out in a channel
featuring water and air interactions, it was shown to only produce a 0.2% depth change
at the disc [141].
4.2.3 Resistance
The resistance was defined in the simulations using the resistance loss coefficient
in ANSYS-CFX. The work described in this study and previous numerical studies
[140][141] used identical coefficients of 1, 2 and 2.5 in separate simulations to repre-
sent the three different porous discs used in the experimental study [185]. The coeffi-
cients were derived based on the thrust coefficient measured in the experimental data
and were estimated using equation (4.5) which is a theoretical relationship between CT
and K [141].
K
CT = (4.5)
(1 + 0.25K)2
4.2.4 Code comparison
This section details a code comparison study between ANSYS-CFX and ANSYS-
Fluent. Although both codes are owned by the same parent company ANSYS they
were developed separately, and calculate the Navier-Stokes equations differently. On
first appearances the two codes differ in their implementation, with ANSYS-CFX us-
ing ANSYS-CFX Expression Language (CEL) and ANSYS-Fluent using its own User-
Defined Functions (UDFs) to customize simulations. The main difference between the
two codes is actually their method of discretization of the Navier-Stokes equations as
ANSYS-CFX uses vertex-centered scheme and ANSYS-Fluent uses a cell-centered
scheme. The difference between the two schemes is how and where the flow variables
are calculated and stored as discussed in subsection 2.1.5.
Both codes were used to simulate analogous models of a porous disc within a wa-
ter channel, specifically the disc had a resistance loss coefficient of 2 which is called
the inertial loss coefficient within ANSYS-Fluent. The predictions of both codes con-
verged to a root-mean-square residual of 1 × 10−5 and have been compared to the
experimental data [185].
4.2.5 Cell type study
This section describes a comparative study which investigated the differences between
linear and quadratic mesh cells (i.e. without and with mid-side nodes shown in Figure
4.4) [35][186]. The study was conducted using the domain and set-up described in
section 4.2 using inlet 1. The difference between linear and quadratic cells is the
number of nodes on each cell, specifically quadratic cells feature mid-side nodes as
shown in Figure 4.4. As part of the study ten different meshes were created, five with
linear cells and five with quadratic cells as shown in Table 4.1. Table 4.1 shows the
large difference in the number of nodes between the difference cell densities.
Figure 4.4: Difference in cell type
Cells Nodes Nodes

(million) linear quadratic
(million) (million)
1.36 0.466 2.64
1.57 0.513 2.96
2.42 0.694 4.22
3.44 0.959 6.08
4.38 1.38 7.81
Table 4.1: Mesh densities using linear and quadratic cells
Figure 4.5: Representation of the blade used in the MEXICO project where the red
zones represent the transition zones between aerofoils
4.3 Actuator disc study
This study reproduces the EU FP5 project ‘Model EXperiments In COntrolled COndi-
tions’ or ‘MEXICO’ project as discussed in subsection 3.2.2 and detailed in the refer-
ence [106]. The experiment was conducted using a three bladed 4.5m diameter wind
turbine set up in the 9.5×9.5m2 open test section of the Large Scale Low Speed Facility
(LLF) wind tunnel of the German Dutch Wind tunnel Organization (DNW). Data were
collected over five sections of the blades, providing loads and pressure measurements
with the flow field determined using Particle Image Velocimetry (PIV) measurements.
The blades were constructed out of three different aerofoils: DU 91-W2-250 at the root
of the blade; RISØ A1-21 at mid span; NACA 64-418 at the outer part of the blade; as
shown in Figure 4.5. For more details on the geometry see reference [106].
In this study, simulations were conducted using the actuator disc (AD) model and
results were compared to the experimental data, as well as previous numerical studies
[120][121]. The numerical simulations were calculated using ANSYS-CFX to obtain
a solution of the Reynolds-Averaged Navier-Stokes (RANS) equations [18] together
with the k-ω SST turbulence model [29].
The geometry of both the wind turbine and wind tunnel requires Non Disclosure
Agreements (NDA) and as such any figure within this document showing them has
been slightly altered and should not be considered accurate. A representation of the
wind tunnel is shown in Figure 4.6. It should be mentioned that a separate NDA is
needed for the exact geometry of the collector and nozzle, which was not obtained
Figure 4.6: Mexico study wind tunnel
for this study. Instead a representation was created using a straight square nozzle and
collector (Figure 4.6).
4.3.1 Mesh
As part of this study two mixed hybrid meshes of 15 million cells (8.8 million nodes)
and 32 million cells (11.5 million nodes) were created with tetrahedral and hexahedral
cells. The meshes consisted of a fine hexahedral refinement zone centered about the
turbine hub, with a cell size of approximately 25mm. This zone was surrounded by a
fine tetrahedral transition region containing cells of approximately 75mm and 50mm
for the coarse and fine mesh respectively, before coarser outer regions as shown in
Figure 4.7. The collector and nozzle regions were meshed using a swept hexahedral
mesh, each with 50 divisions with a bias towards the walls.
Figure 4.8 shows the velocity through the centre of the domain with the 15 million
cell mesh (dashed line), 30 million cell mesh (solid line), numerical study (o) [120]
(a)
(b)
Figure 4.7: Mesh used in the actuator disc study (a) cross-section through the domain
and (b) detailed view
Figure 4.8: Centre line velocity for the 10m/s inlet with the 15 million cell mesh
(dashed line), 30 million cell mesh (solid line), numerical study (o) [120] and ex-
perimental data (.) [106]
and experimental data (.) [106]. It shows that the 15 million cell mesh under-predicts
the velocity drop at the rotor and that the 30 million cell mesh is closer to the other
data and will be used for the rest of this study.
Three inlet velocities of 10m/s, 15m/s and 24m/s were used as part of the experiment
[106]. Preliminary simulations showed that these inlet velocities produced a higher
velocity profile than was reported in the experiment. Due to this, slightly modified
inlet velocities were used in this study, which were taken from the experimental data.
Velocities of 9.3m/s, 14.44m/s and 23.77m/s were used, as these were the average of
the first velocities given in the experimental data just after the inlet and far upstream
from the turbine. The density of the air was also defined based on the average reading
taken during the experiment, of 1.24kg/m3 .
The walls of the wind tunnel, including the nozzle and collector, were defined as
free slip walls to reduce computational expense. Preliminary simulations showed that
the different walls had little effect on the results as the wall boundary layer did not
develop significantly in the inlet.
In this study the AD method was used as described in subsection 2.2.2. The aerody-
namic coefficients were supplied for each aerofoil profile and a smoothing expression
was created in CEL to transition between the coefficients based on the radial position
along the blades. The data was obtained through a NDA and so not detailed here, see
reference [106] for more information.
4.4 Two turbine study
This study investigated the interactions between two wind turbines. The interactions
between the first two turbines are the most important as the second turbine experiences
the highest power loss [13][187][188].
There have been a number of wake interaction studies, both experimental and nu-
merical, ranging from two turbines to full wind farms including 30 turbines [189][190][191]
and 48 turbines [192] as well as infinite farms [136]. There have also been a number of
comparative studies [14][15][193]. The majority of the numerical interaction studies
use actuator methods such as the AD [160][189][194] or AL [174][195], although full
rotor models have recently been conducted [139].
The majority of studies focus on axially aligned turbines and many studies have
investigated the effect of axial or streamwise spacing on wake interactions. Choi et
al. [139] used ANSYS-CFX with the full rotor model, including nacelle and tower,
to investigate the effect of spacings ranging from 3 − 7 rotor diameters. They [139]
found that the second turbine experiences a power loss of 35-65% the power of the
first turbine for the given spacings. They suggested a five diameter spacing as a design
criteria. Yang et al. [194] and Meyer & Meneveau [160] investigated infinite aligned
wind-farm configurations using the ADM and LES. Meyer & Meneveau [160] investi-
gated the optimum spacing and included turbine cost and land cost in their calculations
and determined an approximate optimum spacing of 15 diameters.
Adaramola & Krogstad [196] and Bartl et al. [197] both conducted experimental
studies using two axially aligned 0.9m diameter turbine models at different distances.
Bartl et al [197] used spacings of 3 and 5 diameters while Adaramola & Krogstad [196]
used 3, 6 and 9 diameters. Adaramola & Krogstad [196] found a power loss of 29-45%
depending on the spacing, and that by operating the first turbine outside optimum tip
speed ratio can significantly improve the second turbines performance.
Mikkelsen et al. [174], Troldborg et al. [195] and Flecher & Brown [198] inves-
tigated the effect of streamwise (axial) and crosswind (lateral) displacement on wake
interactions, simulating full and half wake interactions. Mikkelsen et al. [174] and
Troldborg et al. [195] used the ALM in their studies. Troldborg et al. [195] investi-
gated the influence of turbulence intensity on wake structure of full and partial wake
interactions. Mikkelsen et al. [174] investigated the optimal pitch settings of the tur-
bine experiencing wind shear. Flecher & Brown [198] used the Vorticity Transport
Model (VTM) to investigate the wake structure. They found a turbine 6 diameters
apart may experience a power loss of 40-50% and that partial wake interactions cause
significant unsteadiness in the wake and fatigue reducing the life of the blades.
Multiple studies have investigated the difference between aligned and staggered
turbine layouts [189][190][199]. Wu & Porté-Agel [189] used the ADM with and
without rotation and found that a staggered layout has more homogeneous velocity
profile as well has higher velocity and lower turbulent intensity thus greater efficiency
than an aligned layout. Chamorro et al. [190] conducted an experimental study in a
wind tunnel using model turbines with a rotor diameter of 0.15m in a staggered layout.
Their analysis showed that a staggered layout has a 10% increase in power output.
Terrain effects on wake interactions have been investigated. Makridis & Chick
[146] used ANSYS-Fluent to investigated the effects of positioning two wind turbines
on top of different hills. Ozbay [199] investigated the effect hills have on wake interac-
tions using a wind tunnel and roughness elements (spires) and found that wake effects
were reduced if the turbine was placed on top of a hill with a moderate slope (12◦ ).
Tidal turbine interactions have also been investigated. Turnock et al. [200] used
ANSYS-CFX and the ADM to investigate the effect of lateral spacings using 2, 4
and 6 diameters and added a turbulence source term. They agreed with wind turbine
studies that a staggered layout configuration is beneficial. O’Doherty et al. [201]
used ANSYS-Fluent to investigate tidal contra-rotating turbine rotors spaced 1, 2 and
3 diameters apart and found the initial wake recovery is predicted to be faster than a
single turbine, although full recovery is the same.
4.4.1 Current work
Every wind-farm is individual, with different turbine size and layout as well as velocity
profiles, not to mention the terrain. The literature shows that the worst case scenario
is an axially aligned full wake interaction, and that the highest velocity deficit and
power loss is experienced by the second turbine [13][187][188]; and so that is the
configuration chosen in this thesis.
This study features two 10m NREL phase VI turbines [85] but with three blades,
axially aligned at three and six diameters apart, as they seemed to be the most widely
used spacing within the literature. Two previous studies, an experimental [197] and
a numerical [195] study, were chosen for comparison as they feature turbines spaced
three and six diameters apart. As each study uses different turbines and velocities, this
study will compare the results using normalized data in terms of free stream velocity
and turbine size. The numerical simulations were calculated using ANSYS-CFX to
Figure 4.9: Two turbine domain with 3 diameter spacing on top and 6 diameter spacing
on the bottom
obtain a solution of the Reynolds-Averaged Navier-Stokes (RANS) equations [18] to-

gether with the k-ω SST turbulence model [29]. A cylindrical domain was created with
an inlet placed three diameter (30m) in front of the first turbine and an outlet placed six
diameter (60m) behind the second turbine as shown in Figure 4.9. The ground effect
was not included as it would incur an unnecessary additional computational cost and
previous work has showed that including the ground effects without a shear inflow only
increased power production by 0.4%, and adding a shear inflow increased this further
to 5.5% [156].
4.4.2 Mesh
This particular study used a mixed hybrid mesh with unstructured tetrahedral cells
in the shear layer region at the edge of the wake and structured hexahedral cells ev-
erywhere else, as shown in Figure 4.10 which shows how the mesh was divided into
different zones of refinement. A mesh study was conducted on the six diameter domain
using five different mesh densities ranging from 2 million to 80 million cells (5 mil-
lion to 15 million nodes). The mesh densities were created by varying the refinement
within the tip region of the wake as shown in Table 4.2.
(a)
(b)
Figure 4.10: Mesh refinement zones
Nodes Cells
Tip refinement
(million) (million)
5 22 -
6 29 0.2m
11 37 0.2m
12 54 0.15m
15 80 0.125m
Table 4.2: Various mesh densities and tip refinement.
Figure 4.11(a) shows the velocity profiles of the five different meshes at 70m down-
stream of the first turbine and Figure 4.11(b) shows a more detailed view of the max-
imum velocity deficit showing the small difference between the different mesh densi-
ties. The velocities are within 5% of each other, with the greatest difference between
the coarsest two meshes of 22 million and 37 million cells. Figure 4.12 shows the
velocity difference between the 22 million and 37 million cell mesh along the yz plane
through the centre of the domain. It shows that the main difference is within the sec-
ond wake but it is only approximately 2% of the free stream velocity. The greatest
difference is at the disc location mostly due to the differences in the mesh causing a
slightly different force smoothing on the second turbine. Table 4.3 shows the differ-
ences in computational time while using 20 partitions divided between four 2.66GHz
Intel Xeon X5550 computers, which were part of UCLan’s High Performance Com-
puter (HPC) Wildcat cluster1 .
The mesh with 22 million cells and 5 million nodes was chosen for the remainder
of the study, with a six diameter spacing, as the difference between the five mesh den-
sities was small in terms of the velocity profile (within 5% as shown in Figure 4.12)
and it was determined that the other meshes densities were not worth the additional
computational expense as shown in Table 4.3. The simulations featuring a three diam-
eter spacing were conducted with a similar dimensioned domain resulting in an overall
1 www.hpcavf.uclan.ac.uk
Nodes (million) Cells (million) Wall clock time File size

5 22 < 1hour 2.4GB
6 29 ≈ 1 hour 3.1GB
11 37 1.5 hours 5GB
12 54 2 hours 6GB
15 80 3 hours 8GB
Table 4.3: Computational time for various mesh densities
shorter domain reducing the cell and nodes to 17 million cells and 3 million nodes.
For this study four wind speeds of 7m/s, 10m/s, 15m/s and 20m/s were simulated
as these represent important speeds, such as cut in and rated speeds, and there are
significant data from NREL Phase VI project [85] for these speeds. The inlet was
defined with a constant velocity and turbulence intensity of 5%. The walls of the
domain were defined as free slip walls and the outlet was defined as a static pressure
outlet with a relative pressure of zero.
The actuator disc is used in this study for both turbines as described in subsections
2.2.2 and 3.5.1, with the force exerted calculated using equations (2.70)-(2.72). The
aerodynamic coefficients used (such as lift and drag) were taken from the NREL Phase
VI project data [85]. The actuator disc was used instead of the actuator line due to the
reduced computational resources required and the number of simulations desired.
4.4.4 Repeatability study
A repeatability study was conducted as part of this study, where the mesh generation
and solution were calculated multiple times using the same parameters. For this study
the mesh containing 22 million cells and 5 million nodes was used with a wind speed
of 15m/s, as this wind speeds showed the greatest variation in the single actuator study
(a)
(b)
Figure 4.11: Mesh comparison (a) 70m behind the first turbine and (b) detail view of
the velocity deficit
Figure 4.12: Velocity difference between the 22 million and 37 million cell mesh den-
sities
(§5.3). The study involved defining the mesh and boundary conditions of four separate
simulations and clearing all generated data before recreating the mesh and calculating
each simulation individually.
4.5 Rotor diameter study
This study investigated the influence rotor diameter size has on wind turbine wakes
and the aerodynamic performance on a downwind turbine. It has been proposed by
Steinbuch et al. [202] and Johnson & Thomas [203], as well as shown experimentally
by Adaramola & Krogstad [196], that the second turbine’s performance can be signif-
icantly improved by operating the first turbine outside optimum tip speed ratio. This
would require monitoring and learning by trial and error [203] to get the most out of
the wind-farm. It would be more beneficial if the same effect could be achieved with-
out the constant need for monitoring. Chowdhury et al. [8] showed that a wind-farm
featuring different size turbines could produce 10% more power than one with identi-
cal turbines. This particular study [8] was limited to a single case using an analytical
flow model [204]. Due to this it is necessary to conduct a more detailed study into how
the rotor diameter effects the energy yield of a wind-farm.
To investigate the effect that rotor diameter has on wake interactions, six turbine
configurations were investigated each featuring two axially aligned turbines. The first
turbine featured three different diameters of 7.5m, 10m and 12.5m while the second
turbine always featured a 10m diameter rotor. These turbine sizes were chosen with
the 10m NREL phase VI turbine [85] as a starting point due to the amount of data
and previous work available. The other sizes where chosen to represent possible size
variations of full size turbines (i.e. 75m, 100m and 120m) which are currently on the
market that are rated with similar power allowing for wind-farms to be constructed
using different rotor diameters without significant loss of power.
The turbines in this study are all based on the 10m NREL phase VI turbine [85]
with three blades. All aerofoil coefficients and sizes are extrapolated from the original
10m NREL phase VI turbine [85]. This study has investigated an axially aligned layout
featuring two spacings of 30m and 60m for the same reasons as described in section
4.4. Apart from the different rotor diameters all other aspects are in line with the two
turbine study as described in section 4.4.
Chapter 5
Results
This chapter details the results obtained from the five studies detailed in Chapter 4.
5.1 One-dimensional momentum theory results
This section details the results of the comparative study carried out between CFD sim-
ulations carried out using ANSYS-CFX and the one-dimensional momentum theory,
as described in section 4.1.
The porous disc and actuator disc methods have been compared with the one-
dimensional (1D) momentum theory, as described in subsection 3.1.1 and the reference
[95], specifically the pressure and velocity profiles along the centre line of the rotor.
The 1D momentum theory, also known as the simple actuator theory, is an applica-
tion of the momentum equation applied to an idealized turbine. It uses control volume
analysis to consider an infinitely thin frictionless disc with a constant momentum sink
within an inviscid and incompressible fluid. Figure 5.1(a) shows the pressure and ve-
locity profiles given by the 1D momentum equation. Figure 5.1(b) shows the pressure
and velocity profiles produced by a constantly loaded actuator disc without rotation,
defined using the thrust coefficient (CT ). The overall profiles are in good agreement,
126
CHAPTER 5. RESULTS 127
(a) (b)
Figure 5.1: Pressure and velocity profiles along the centre line of the rotor given by (a)
the 1D momentum theory and (b) the actuator disc method
with only the magnitudes of the plot changing, depending on the CT value and the
characteristics of the disc as expected.
5.2 Porous disc Results
This section details the results of the porous disc study which was detailed in section
4.2. Simulations were carried out of a porous disc within a water channel and compared
to the results obtained by previous numerical [140][141] and experimental [184][185]
studies. The experimental data [184][185] was measured using an Acoustic Doppler
Velocimeter (ADV).
All the data presented in this study was produced using ANSYS-CFX and calcu-
lated with a root-mean-square residual of 1 × 10−5 applied to the RANS equations,
which was in line with the reference [141]. The velocities were compared to both
previous numerical and experimental results which were normalized using the free
stream velocity of the flow described between 0.5 < y/2R < 2.5 at the inlet which was
0.331m/s in the experimental study [185], 0.337m/s in the reference [141] and 0.33m/s
Figure 5.2: Velocity profile with no disc of the channel flow (solid line) and the duct
flow (dashed line) at 24m, 31m and 40m from the inlet
in the work described in this study. This range was used as there was no experimental
data beyond this range. The vertical height was also normalized using the diameter of
the disc (2R).
5.2.1 Influence of the boundary conditions
Before detailing the results of the wake study, a comparison of the effects of the bound-
ary types is needed. To do this the same domain and mesh was setup excluding the
momentum loss to observe how the velocity profiles develop without the influence of
the discs. Figure 5.2 shows the differences between the channel (solid lines) and duct
(dashed lines) velocity profiles as they develop through the domain. The figure shows
how the channel flow is almost unchanged as the inlet was defined with a channel ve-
locity profile. The duct profile changes significantly as expected with the additional
wall boundary causing a sharp decrease in velocity at the top of the domain which
forces the central velocity to increase to maintain the same mass flow rate.
5.2.2 Wake predictions
The experimental study [185] conducted experiments using discs with different poros-
ity measurements to represent different values of the thrust coefficient (CT ). CT is a non
dimensional variable used to describe a rotor’s characteristics. The greater CT value

the greater the wake expansion and turbulence levels within the wake due to increased
momentum deficit. The CT values of the porous discs in the experimental study [185]
were measured using a pivot arm attached to a load cell. The thrust coefficient can be
described numerically using equation (5.1). It requires the thrust (T ) to be estimated
which can be achieved in a number of ways. In the reference [141] the thrust coeffi-
cient was estimated from the results using equation (5.2) to define the thrust. However
in this work, as in the reference [140], the thrust was calculated using equation (5.3),
T
CT = 1 2
(5.1)
2 ρU0 Ar
ρ
T = KUr |Ur |Ar (5.2)
2
T = ∆pAr (5.3)
where T is the thrust, ρ is the density, U0 the free stream velocity, Ar is the disc
area, K is the resistance coefficient, Ur is the velocity at the rotor and ∆p is the change
in pressure over the centre of the disc.
Table 5.1 shows the CT values of the previous experimental and numerical studies
compared to the values calculated in this study. It shows that the channel simulations
conducted as part of this study have the closest CT values compared to the experimental
data. It also shows that the CT values predicted by the duct simulations have the least
agreement with the experimental data [185].
Figure 5.3 shows the velocity profiles along the centre line of the domain and shows
good agreement for all simulations in terms of the velocity characteristics. The velocity
deficit is over-predicted when compared to the numerical [141] and experimental data
[185] for inlet 1, which has a delayed velocity recovery and appears to be offset from
Study
Disc number Experimental Numerical Channel Duct
[185] [141]
1 0.61 0.65 0.60 0.64
2 0.86 0.91 0.86 0.93
3 0.97 0.98 0.93 1.00
Table 5.1: Thrust coefficient results to 3 s.f.
(a) (b) (c)
Figure 5.3: Velocity along the centre line showing the channel (solid line), duct (dashed
line), numerical (x) [141] and experimental data (o) [185] for the CT values of a) 0.61,
b) 0.86 and c) 0.97 for the experimental study [185]. The inlet 1 simulations are shown
in black and inlet 2 simulations are in red
the other data sets. Inlet 2 shows a much better prediction of the experimental data
[185] and both inlets show the duct had increased velocity recovery than the channel
simulation (discussed in §§6.1.3).
Figure 5.4 shows the turbulent intensity along the centre line of the domain and
the difference between inlet 1 and 2 (described in §§4.2.2), with inlet 1 having a lower
starting turbulence intensity and subsequent peak. Both inlets show very little change
in turbulent intensity just behind the rotor and then an almost linear increase up to
the maximum intensity. This increase is most likely due to the presence of the wake
edge shear layer, with the maximum turbulent intensity indicating the merger of the
(a) (b) (c)
Figure 5.4: Turbulence intensity along the centre line showing the channel (solid line),
duct (dashed line), numerical (x) [141] and experimental data (o) [185] for the CT
values of a) 0.61, b) 0.86 and c) 0.97 for the experimental study [185]. The inlet 1
simulations are shown in black and inlet 2 simulations are in red.
layers and subsequent end of the near wake region (discussed in §§6.2.1). The near
wake region of the flow, which is defined behind the disc up until the wake edge shear
layers meet at the centre line of the wake, creating a peak in turbulence. The near
wake region varies in distance generally from 2 to 8 diameters downstream, depending
on the disc geometry and flow conditions. Figure 5.4 shows that the model is able to
predict the intensity accurately far down stream of the disc. Although, it is unable to
predict the peak in the turbulence intensity behind the rotor both in terms of magnitude
and location (discussed in §§6.2.1).
Figures 5.5 and 5.6 show the velocity and turbulence intensity profiles at vari-
ous distances downstream of the disc and for the three different discs. The distances
downstream correspond to 8R, 14R, 22R, 30R and 40R downstream of the disc. These
locations where chosen as they coincide with the locations used in the experimental
study [185].
For all profiles in Figure 5.5 agreement was achieved1 for the majority of the profile
characteristics, such as the locations of highest and lowest velocities with the main
1 at least from a qualitative viewpoint
numerical discrepancy at the maximum velocity deficit for all simulations. While the
initial velocity drop is over-predicted at the centre, the free stream and floor boundary
layer features are predicted well. Figure 5.5 shows that the duct simulations predicted
a smaller velocity deficit at the centre of the wake and higher velocity towards the
boundaries than the channel simulations (discussed in §6.1.3).
Figure 5.5 shows how the velocity deficit of the experimental data recovers more
than the numerical data with inlet 1 simulations showing the least recovery. The ex-
perimental data seems to have an almost constant velocity by 22R downstream and
completely recovers by 30R downstream, whereas all numerical simulations still show
some velocity deficit at 40R (discussed in §6.2.1).
Figure 5.6 shows the turbulence intensity of the models in comparison with the
experimental data. There is little difference between the solid and dashed lines repre-
senting the channel and duct flows respectively; with all models predicting intensities
below that of the experiment data (discussed in §§6.2.1). The figure shows how the
experimental data peaks earlier and higher than the modelled numerical data. Beyond
approximately 22R downstream of the rotor the modelled and experimental data are
very close and almost the same by 40R .
5.2.3 Cell type study results
This section details the results of a comparative study which investigates the differ-
ences between linear and quadratic mesh cells (without and with mid-side nodes,
shown in Figure 4.4)[35][186]. The study was conducted using inlet 1. This study
compared 10 different mesh densities, shown in Table 5.2, using both ANSYS-CFX
and ANSYS-Fluent. Table 5.2 shows the quadratic mesh calculations took 48-186%
longer to solve then their linear equivalent.
Figures 5.7 and 5.8 shows the velocity profiles of the different mesh densities
(a)
(b)
(c)
Figure 5.5: Normalized velocity of the channel (solid line), duct (dashed line) and
experimental data (o) [185] at different CT values of a) 0.61, 0.62, b) 0.86, 0.91 and
c) 0.97, 0.99 for this study and the experimental study respectively [185]. The inlet 1
simulations are shown in black and inlet 2 simulations are in red. Each figure represents
the distance downstream of the disc corresponding to 8R, 14R, 22R, 30R and 40R
respectively
Nodes ×106 Time [mins] Iterations

Cells linear quadratic linear quadratic linear quadratic
×106
1.36 0.466 2.64 11 17 63 104
1.57 0.513 2.96 13.5 20 63 105
2.42 0.694 4.22 19.5 31 62 106
3.44 0.959 6.08 21 60 54 151
4.38 1.38 7.81 30 75 54 149
Table 5.2: Mesh densities using linear and quadratic cell types
(a)
(b)
(c)
Figure 5.6: Turbulence intensity at 8R, 14R, 22R, 30R and 40R downstream of the
disc, with the channel (solid line), duct (dashed line) and experimental data (o) [185]
at different CT values of a) 0.61, 0.62, b) 0.86, 0.91 and c) 0.97, 0.99 for this study and
the experimental study respectively [185]. The inlet 1 simulations are shown in black
and inlet 2 simulations are in red
7 diameters behind the disc for both the duct and channel simulations predicted by
ANSYS-Fluent. Figure 5.7(a) shows the velocity profile of all 10 mesh density. It
shows that the largest variation is at the highest velocity deficit at the centre of the
wake and that for each cell density there is little change in velocity profile. Figure
5.7(b) is a refined view of the velocity profiles showing that for each cell density the
quadratic cell version predicted a higher velocity deficit although no quadratic cell
density predicts a higher velocity deficit then a finer linear cell density. Figure 5.8
shows very similar results to the channel simulations. Figure 5.9 shows the turbulence
intensity through the centre of the domain behind the disc and that there is a variation
in the turbulence predicted by each mesh although this is less than 3%.
5.2.4 Code comparison results
This section details the results of a code comparison between ANSYS-CFX and ANSYS-
Fluent. Both codes simulated analogous models of a porous disc within a water chan-
nel, specifically the disc with a resistance of 2.
Figure 5.10 shows the normalized velocity profile through the centre of the domain
of ANSYS-CFX (red) and ANSYS-Fluent (blue) simulations for the channel (solid
line) and duct (dashed line) simulations as well as the experimental study (o) [185].
Figure 5.10 shows that ANSYS-CFX and ANSYS-Fluent predict a very similar veloc-
ity profile although both over-predict the velocity deficit compared to the experimental
data. The similarity between the two profiles is closest just behind the disc and it di-
minishes as the wake develops. ANSYS-CFX predicts the highest velocity deficit of
X
all data sets at 2R = 2 and predicts the longest wake recovery. The duct simulation
performed by ANSYS-Fluent seems to predict a velocity profile closest to the experi-
mental data [185] within the far wake.
Figures 5.11 and 5.12 show normalized velocity profiles behind the disc at various
(a)
(b)
Figure 5.7: Velocity profile 7 diameters behind the disc of the channel simulations with
the quadratic cells (solid line) and linear cells (dashed line)
(a)
(b)
Figure 5.8: Velocity profile 7 diameters behind the disc of the duct simulations with
the quadratic cells (solid line) and linear cells (dashed line)
Figure 5.9: Normalized turbulence intensity behind the disc
Figure 5.10: Normalized velocity profile through the centre of the domain with
ANSYS-CFX (red), ANSYS-Fluent (blue) and experimental data (o) [185] for the
channel (solid line) and duct (dashed line) simulations
locations. They show the results from ANSYS-CFX (red) and ANSYS-Fluent (blue)
for the channel (solid line) and duct (dashed line) simulations as well as the experimen-
tal study (o) [185]. Figure 5.11 shows the velocity profile four diameters behind the
disc and the difference between the channel and duct velocity profiles. ANSYS-CFX
predicts a larger difference in velocity within the boundary layer of the two simula-
y
tions of approximately 10% 2R = 0.5 away from the edge. Whereas, ANSYS-Fluent
y
predicts an almost identical velocity profile until 2R = 0.5 when the velocities sharply
diverges and predicts a higher velocity at the roof of the channel simulation. The chan-
nel simulation of ANSYS-CFX seems to be the closest to the experimental data [185]
at this location. Figure 5.12 shows the velocity profiles at various locations behind the
disc and that the differences between the channel and duct simulations has more of an
impact downstream of the disc with the development of the boundary layer of the duct
simulation. Figures 5.12(c) and 5.12(d) particularly show that as the boundary layer
develops the differences between the two velocity profiles increases and extends into
the wake.
Figure 5.13 shows the turbulence intensity through the centre of the domain with
ANSYS-CFX (red) and ANSYS-Fluent (blue) for the channel (solid line) and duct
(dashed line) simulations as well as the experimental study (o) [185]. It shows that
X
ANSYS-CFX predicts a higher turbulence intensity than ANSYS-Fluent for 2 < 2R <
8 and that ANSYS-CFX predicts a lower turbulence intensity for the duct simula-
tion. ANSYS-Fluent predicts near identical turbulence intensity for both channel and
duct simulations. Apart from these discrepancies ANSYS-CFX and ANSYS-Fluent
predicted similar turbulence intensity although both under-predict the intensity of the
X
experimental data [185] until well into the far wake at 2R < 17.
Figures 5.14 and 5.15 show the turbulence intensity behind the disc at various lo-
cations with ANSYS-CFX (red) and ANSYS-Fluent (blue) for the channel (solid line)
and duct (dashed line) simulations as well as the experimental study (o) [185]. Figure
x
Figure 5.11: Normalized velocity profile behind the disc at 2R = 4 with ANSYS-CFX
(red), ANSYS-Fluent (blue) and experimental data (o) [185] for the channel (solid
line) and duct (dashed line) simulations
5.14 shows the velocity profile four diameters behind the disc and shows that ANSYS-
Fluent predicts very similar turbulence intensity for both channel and duct simulations
within the roof boundary layer. ANSYS-CFX predicts similar turbulence intensity be-
x
low 2R = 1 although above this the two turbulence profiles differ with the channel
simulation having a higher intensity. Both ANSYS-Fluent and ANSYS-CFX predict
turbulent intensity peaks at the edge of the wake although neither predict the peak
at the centre of the wake which was observed in the experimental study [185]. Fig-
ure 5.15 shows the turbulence intensity at locations further behind the disc and shows
that further downstream the turbulence intensity peaks at the edge of the wake which
reduces to an almost uniform intensity. The variation between the channel and duct
simulations also diminish for the bottom half of the domain with a clear difference at
the roof of the domain. It shows that the turbulence intensity of the experimental data
x
[185] reduces quicker than the numerical data although at approximately 2R = 15 all
the data sets produce similar turbulence intensities.
x x
(a) 2R =7 (b) 2R = 11
x x
(c) 2R = 15 (d) 2R = 20
Figure 5.12: Normalized velocity profile behind the disc with ANSYS-CFX (red),
ANSYS-Fluent (blue) and experimental data (o) [185] for the channel (solid line) and
duct (dashed line) simulations
Figure 5.13: Turbulence intensity through the centre of the domain with ANSYS-CFX
(red), ANSYS-Fluent (blue) and experimental data (o) [185] for the channel (solid
line) and duct (dashed line) simulations
x
Figure 5.14: Turbulent intensity behind the disc at 2R = 4 with ANSYS-CFX (red),
ANSYS-Fluent (blue) and experimental data (o) [185] for the channel (solid line) and
duct (dashed line) simulations
x x
(a) 2R =7 (b) 2R = 11
x x
(c) 2R = 15 (d) 2R = 20
Figure 5.15: Turbulence intensity behind the disc with ANSYS-CFX (red), ANSYS-
Fluent (blue) and experimental data (o) [185] for the channel (solid line) and duct
(dashed line) simulations
5.3 Actuator disc study results
This study reproduced the EU FP5 project ’Model EXperiments In COntrolled COn-
ditions’ or ’MEXICO’ project as discussed in section 4.3 and detailed in the reference
[106]. In this section the actuator disc (AD) method, which was used to simulate the
turbine, is compared with the experimental MEXICO data [118][119] and a previous
numerical study [120].
Figure 5.16 shows the velocity profile through the centre of the domain at hub
height offset in the horizontal (z direction) by 1.84797m or z/R = 0.82132. It shows
the velocity for each wind speed 10m/s (blue), 15m/s (green) and 24m/s (red) for this
study (solid line), previous numerical study (o) [120] and the experimental data (.)
[106]. It shows that all three data sets show a similar velocity profile although the
numerical studies do not have as steep a drop as the experimental data [118][119]
just behind the turbine (discussed in §6.2.1). Both numerical studies under-predict the
velocity deficit within the wake for all but the slowest wind speed. The data produced
in this study predicts the highest velocity of all three data sets even before the turbine,
shown in Figure 5.17. The greatest variation in velocity between the data sets is within
the wake of the 15m/s wind speed.
Figures 5.18 and 5.19 show the normalized cross-sectional velocity profiles at var-
ious locations throughout the domain. They show the normalized data from the 10m/s
simulation although the data was normalized based on the inlet velocities used as de-
scribed in section 4.3. Figure 5.18 shows the velocity profile at x/R = 0.4 before the
turbine. It shows very good agreement between the two numerical studies with the
highest difference at the edge of the velocity stream.
Figure 5.19 shows the normalized velocity profiles at six locations downwind of
the turbine produced by this study (blue) and previous numerical study [120]. It shows
the characteristics of the velocity profiles of the two numerical studies are in good
Figure 5.16: Velocity profile through the centre of the domain at 1.84797m or z/R =
0.82132 for this study (solid line), experimental data (.) [106] and previous numerical
study (o) [120] for three inlet velocities
agreement for all locations. The velocity away from the wake and the wake expansion
is also in very good agreement for all locations. The discrepancy between the two
studies is at the velocity magnitude towards the edge and centre of the wake. This
study shows a higher velocity deficit at the edge of the wake then the previous study
[120] for all locations except the furthest downwind. At the centre of the wake this
study predicts little to no velocity deficit but the previous study [120] predicts their
highest velocity deficit here (discussed in §§6.2.2).
Figure 5.17: Velocity variation towards the inlet for this study (solid line), experimental
data (.) [106] and previous numerical study (o) [120]
Figure 5.18: Normalized velocity profile at x/2R = −0.4 with this study (blue) and
previous numerical study [120] (red o)
x
(a) x
=0 (b) 2R = 0.3
2R
x x
(c) 2R = 0.6 (d) 2R =1
x x
(e) 2R = 1.5 (f) 2R =2
Figure 5.19: Normalized velocity profiles of this study (blue) and previous numerical
study [120] (red o) at various locations
5.4 Two turbine study results
This study investigates the interactions between two wind turbines. It features two
10m NREL phase VI turbines [85] but with three blades, aligned with a spacing be-
tween them of three and six diameters. The domain, mesh and boundary conditions
used in this study were discussed in section 4.4. All the data presented in this section
were calculated using ANSYS-CFX with a mesh containing 22 million cells (5 million
nodes) and the simulations had converged to a root-mean-square residual of 1 × 10−6
unless otherwise stated. The data produced has been compared with two previous stud-
ies, an experimental [197] and a numerical [195]. The data used for comparison are
normalized based on the inlet velocity and turbine radius. Although four wind speeds
have been simulated only the 7m/s simulation was compared as it is closest to the other
studies [197][195] in terms of the tip speed ratio (5.5). The simulated data for the other
wind speed was however used in section 5.5.
Figure 5.20 shows the normalized velocity profiles downstream of two turbines
spaced 3 diameters apart where the solid line is the work described in this study, (o)
and (x) represent numerical data at 5% and 10% ambient turbulence intensity [195]
and (*) represents experimental data [197]. The velocity profiles are at 50m, 60m and
70m respectively downwind of the first turbine. There is good agreement between the
results of this study with previous studies at 50m and 60m. There is less agreement at
70m particularly at the centre of the wake although the edge of the wake remains in
good agreement.
a) 50m b) 60m c) 70m
Figure 5.20: Normalized velocity profiles at 50m, 60m and 70m downstream respec-
tively of the first turbine with this study (solid line), previous numerical data at 5% (o)
and 10% (x) ambient turbulence [195] and experimental data (*) [197]
Simulation Cells Nodes

number
1 4878311 22065456
2 4878455 22065941
3 4880487 22063055
4 4887175 22095389
Table 5.3: Mesh variation in the repeatability study
5.4.1 Repeatability study results
A repeatability study was carried out where the mesh generation and solution were
calculated multiple times using the same parameters. Table 5.3 shows the mesh density
of each simulation and that each mesh was similar with a variation of 8864 cells and
32334 nodes, which is a variation of less than 0.2%. Figure 5.21 shows the velocity
difference between simulation 1 and 4, which had the largest mesh variation, along the
yz plane through the centre of the domain and Figure 5.22 shows the region around the
disc showing the highest difference is at the disc location. It shows that the variation
in mesh density had little effect on the velocity field producing a maximum variation
of 1%.
Figure 5.21: Velocity difference along yz plane between two simulations used in the
repeatability study
Figure 5.22: Detailed view of the velocity difference
5.5 Rotor diameter study results
This study investigated the influence rotor diameter size has on wind turbine wakes and
the aerodynamic performance on a downwind turbine. The turbines in this study were
all based on the 10m NREL phase VI turbine [85] with three blades. Here two aligned
turbines with spacings of 30m and 60m are simulated. The first turbine features three
different diameters of 7.5m, 10m and 12.5m which were 30m away from the inlet while
the second turbine always featured a 10m diameter rotor. All the data presented in this
section was calculated using ANSYS-CFX with a mesh containing 22 million cells (5
million nodes) and the simulations had converged to a root-mean-square residual of
1 × 10−6 unless otherwise stated. Apart from the different rotor diameters all other
aspects are in line with the two turbine study as described in section 4.5.
The power of each turbine can be calculated using equation (5.4). In this study
an approximate power extracted by the turbine has been calculated using the average
velocity at the rotor and a power coefficient (C p ) of 0.4, which was used as the normal
range for a wind turbine is between 0.35-0.45 [205].
1
P = ρπr2 u3C p . (5.4)
2
Where ρ is the air density, r is the rotor radius, u is the wind speed which in this
situation is the average velocity at the rotor and C p is the power coefficient.
To investigate the effect rotor diameter has on turbine interactions a baseline needs
to be established. In this case it was the two 10m turbines simulated in section 5.4 and
unless otherwise stated, any comparison is made in regard to this configuration.
5.5.1 First turbine
Before the interactions between the turbines is analysed the performance of the first
turbine is considered. Figure 5.23 shows the velocity profiles at the first turbine for
each wind speed with the 7.5m turbine (blue), 10m turbine (green) and 12.5m turbine
(red). Table 5.4 shows the average velocity at the rotor of the first turbine. Figure 5.23
and Table 5.4 show that as the rotor diameter increases the average velocity measured
at the rotor decreased as expected as more momentum is extracted from the flow.
Table 5.5 shows the approximate power extracted by the first turbine, using equa-
tion (5.4), and shows that as the rotor diameter increases so does the power extracted,
as expected. Table 5.6 shows the difference in the amount of power extracted and per-
centage difference of the 7.5m and 12.5m turbines compared to the 10m turbine. Each
rotor has a rotor diameter of 7.5m, 10m and 12.5m and hence swept area correspond-
ing to 44.2m2 , 78.5m2 and 123m2 respectively. Tables 5.4 - 5.6 show that as the rotor
diameter increases the power extracted increases even though the average velocity at
the rotor decreases.
Figure 5.23: Velocity profiles at the first turbine for each wind speed of 7m/s, 10m/s,
15m/s and 20m/s respectively and the different colours represent the different turbine
configurations; 7.5m turbine (blue), 10m turbine (green) and 12.5m turbine (red)
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 6.30 5.91 5.63
10 9.42 9.17 8.98
15 14.2 14.2 14.1
20 18.9 18.9 18.9
Table 5.4: Velocity at the first turbine to 3 s.f.
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 2660 3890 5250
10 8870 14500 21200
15 30400 53500 82900
20 71900 127000 199000
Table 5.5: Approximate power extracted by the first turbine to 3 s.f.
5.5.2 60m spacing wake
This subsection details the results of the 12 simulations where the spacing between
the two turbines was 60m specifically looking at the wake and velocities. Figure 5.24
shows the velocity profiles at 70m downwind of the first turbine and/or 10m behind
the second turbine. Each turbine configuration is shown in a different colour where the
velocity behind the 7.5m turbine is blue, behind the 10m turbine is green and behind
the 12.5m turbine is red. It shows that the turbine behind the 7.5m turbine experiences
the smallest velocity deficit and the turbine behind the 12.5m turbine experiences the
highest velocity deficit. It shows that the wake behind the 12.5m turbine is widest and
the wake behind the 7.5m turbine is the narrowest. Figures 5.25 and 5.26 show the
velocity profiles through the centre of the domain offset by 4m and 2m respectively
showing the velocity behind the 7.5m (blue), 10m (green) and 12.5m (red) turbines.
They show the velocity profiles reducing due to the turbine and then slowly recovering
as expected and that the velocity recovers quicker after the second turbine compared
to the first with the velocity behind the 7.5m turbine constantly higher than the other
Rotor size [m]

Wind speed [m/s] 7.5 % deficit 12.5 Factor
7 1240 31.8 1360 1.35
10 5670 39.0 6700 1.46
15 23100 43.2 29400 1.55
20 55000 43.3 72200 1.57
Table 5.6: Power difference of the first turbines compared with the baseline to 3 s.f.
turbines.
5.5.3 60m spacing second turbine velocity
Figure 5.27 shows the normalized velocity profiles at the rotor of the second turbine
for each wind speed. It shows that the greatest velocity variation is along the blade
section of the rotor approximately 4m off centre and for the 10m/s wind speed. The
least variation is found around the centre of the rotor where the hub would be located
approximately less than 1m from the centre. Table 5.7 shows the average velocity at
the second rotor for each wind speed. Table 5.8 shows the difference in average veloc-
ity both in magnitude and percentage, of the 7.5m and 12.5m turbines compared to the
10m turbine. Tables 5.7 and 5.8 show the greatest velocity magnitude difference is at
10m/s and smallest is at 15m/s for all turbines. The greatest velocity percentage differ-
ence is at 7m/s and that increasing the wind speed decreases this percentage difference.
The greatest average velocity at the rotor for each wind speed was always experienced
by the turbine behind the 7.5m turbine by a minimum of 4% of the free stream wind
speed.
Figure 5.24: Velocity at 70m behind the first turbine for each wind speed of 7m/s
10m/s, 15m/s and 20m/s respectively and the different colours represent the different
turbine configurations; 7.5m turbine (blue), 10m turbine (green) and 12.5m turbine
(red)
Figure 5.25: Velocity profile through the centre of the domain offset by 4m with ve-
locity behind the 7.5m (blue), 10m (green) and 12.5m (red) turbines
Figure 5.27: Velocity profiles at the second turbine for each wind speed of 7m/s 10m/s,
configurations; 7.5m turbine (blue), 10m turbine (green) and 12.5m turbine (red)
5.5.4 60m spacing second turbine power
This subsection details the results of the 12 simulations where the spacing between
the two turbines was 60m, specifically the power. Figure 5.28 shows the approximate
power extracted by the second turbine, using equation (5.4), for each wind speed using
the average velocities from Table 5.7. It shows that the greatest power is extracted by
the turbine behind the 7.5m turbine and the greatest power deficit occurs at the turbine
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 4.65 4.00 3.54
10 8.07 7.35 6.88
15 12.9 12.4 12.1
20 17.2 16.6 16.2
Table 5.7: Average velocity at the second turbine to 3 s.f.
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 0.649 16.2 -0.460 -11.5
10 0.724 9.86 -0.463 -6.30
15 0.525 4.23 -0.277 -2.23
20 0.665 4.02 -0.322 -1.94
Table 5.8: Velocity difference at the second turbine compared to the baseline to 3 s.f.
behind the 12.5m turbine. Table 5.9 shows the difference in power in terms of magni-
tude and percentage of the turbines behind the 7.5m and 12.5m turbines compared to
the turbine behind the 10m turbine. The greatest power difference is at 20m/s where an
extra 10kW was extracted by the turbine behind the 7.5m turbine compared to the tur-
bine behind the 10m. It shows that for all wind speeds the greatest percentage power
variation is at 7m/s with the percentage decreasing as the wind speed increases. It
shows that turbine behind the 7.5m turbine extracts the most power, nearly 60% more
power than the turbine behind the 10m turbine and 2.28 times more power than the
turbine behind the 12.5m turbine at 7m/s. Even at 10m/s the turbine behind the 7.5m
turbine extracts 12.5% more power then the turbine behind the 10m turbine and 19%
more power than the turbine behind the 12.5m turbine.
However, as the configuration contains two turbines, to get a fuller picture of the
situation the turbines need to be considered together. Figure 5.29 shows the combined
power of the first and second turbines taken from Table 5.5 and Figure 5.28 for each
wind speed. Table 5.10 shows the difference in power in terms of magnitude and per-
centage of the turbines configurations with different diameters compared to the turbine
Figure 5.28: Power extracted by the second turbine (7.5m turbine (blue), 10m turbine
(green) and 12.5m turbine (red)) at different wind speeds
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 686 57.0 -370 -30.7
10 2430 32.6 -1330 -17.7
15 4770 13.2 -2360 -6.54
20 10700 12.5 -4890 -5.72
Table 5.9: Difference in power extracted by the second turbines compared to the base-
line to 3 s.f.
Figure 5.29: Combined power extracted by the first and second turbines at different
wind speeds
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 - 551 -10.8 990 19.4
10 - 3240 -14.7 5370 24.4
15 - 18400 -20.5 27000 30.2
20 - 44300 -20.8 67300 31.7
Table 5.10: Comparison in the combined power extracted by the 7.5m and 12.5m
configurations with the baseline configuration to 3 s.f.
layout containing just 10m turbines. Figure 5.29 and Table 5.10 show that the greatest
power variation both in terms of magnitude and percentage difference occurs at 20m/s
and the smallest variation occurs at 7m/s. It also shows that the 12.5m turbine config-
uration extracts the greatest amount of power and the 7.5m configuration extracts the
least.
Table 5.11 shows the power deficit between the first and second turbines. It shows
that all turbines behind the 10m and 12.5m turbines experienced a power deficit of
between 30-85% depending on the turbine size and wind speed. The turbines behind
the 12.5m turbine experienced the largest power deficit for each wind speed. The
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 28.8 69.1 84.1
10 -11.7 48.6 71.1
15 -34.3 32.6 59.3
20 -33.8 32.6 59.5
Table 5.11: Power deficit of the second turbine to 3 s.f.
turbine behind the 7.5m turbine has the least power deficit of 28% at 7m/s and for all
other wind speeds the second turbine extracts more power then the first turbine, up to
34%.
5.5.5 30m spacing wake
The subsection details the results of the 12 simulations where the spacing between the
two turbines was 30m, specifically looking at the wake and velocities. Figure 5.30
shows the velocity profiles 70m downwind of the first turbine and/or 40m behind the
second turbine. Each turbine configuration is shown in a different colour where the
velocity behind the 7.5m turbine is blue, behind the 10m turbine is green and behind
the 12.5m turbine is red. The majority of the velocity profiles have one deficit peak
representing the far wake region. Similar to the 60m situation the velocity behind
the 7.5m turbine experiences the smallest velocity deficit and the turbine behind the
12.5m turbine experiences the highest velocity deficit. It shows that the wake behind
the 12.5m turbine is widest and the wake behind the 7.5m turbine is the narrowest.
Figures 5.31 and 5.32 show the velocity profiles through the centre of the domain offset
by 4m and 2m respectively showing the velocity behind the 7.5m (blue), 10m (green)
and 12.5m (red) turbines. Here very similar velocity profiles are observed to those of
the 60m spacing situation, with the velocity recovering quicker after the second turbine
compared to the first and that the velocity behind the 7.5m turbine is constantly higher
than the other turbines.
Figure 5.30: Velocity profiles 70m behind the first turbine for each wind speed of
7m/s 10m/s, 15m/s and 20m/s respectively with the different colours representing the
different turbine configurations; 7.5m turbine (blue), 10m turbine (green) and 12.5m
turbine (red)
5.5.6 30m spacing second turbine velocity
Figure 5.33 shows the normalized velocity profiles at the rotor of the second turbine
for each wind speed. It shows that the greatest velocity variation is along the blade
section of the rotor approximately 4m off centre and for the 10m/s wind speed which
is similar to the 60m spacing situation. The smallest variation is found around the
centre of the rotor where the hub would be located approximately less than 1m from
centre. Table 5.12 shows the average velocity at the second rotor for each wind speed.
Table 5.13 shows the difference in average velocity both magnitude and percentage
of the 7.5m and 12.5m turbines compared to the 10m turbine. Tables 5.12 and 5.13
show the highest velocity magnitude difference is at 10m/s and smallest is at 15m/s
for all turbines. The greatest velocity percentage difference occurs when Vz = 7m/s.
The greatest average velocity at the turbine for each wind speed was experience by the
turbine behind the 7.5m turbine, which was between 3.5-18% higher than that behind
the 10m turbine and up to 35% higher than that behind the 12.5m turbine. Increasing
the wind speed decreases this percentage difference for all turbines.
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 4.54 3.84 3.38
10 7.98 7.29 6.91
15 12.9 12.4 12.1
20 17.1 16.5 16.2
Table 5.12: Average velocity at the second turbine with a spacing of 30m to 3 s.f.
5.5.7 30m spacing second turbine power
This subsection details the results of the 12 simulations where the spacing between the
two turbines was 30m, specifically at the power. Figure 5.34 shows the approximate
power extracted by the second turbine, using equation (5.4), for each wind speed using
the average velocities from Table 5.12. It shows that the greatest power is extracted
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 0.698 18.2 0.469 12.2
10 0.695 9.54 0.376 5.15
15 0.464 3.75 0.223 1.80
20 0.583 3.54 0.276 1.67
Table 5.13: Velocity difference at the second turbine compared to that behind the base-
line to 3 s.f.
Figure 5.33: Velocity profiles at the second turbine for each wind speed of 7m/s 10m/s,
behind the 7.5m turbine (blue), 10m turbine (green) and 12.5m turbine (red)
Figure 5.34: Power extracted by the second turbine (7.5m turbine (blue), 10m turbine
(green) and 12.5m turbine (red)) at different wind speeds
by the turbine behind the 7.5m turbine and the least power is extracted by the turbine
behind the 12.5m turbine. Table 5.14 compares the power data showing the difference
in power in terms of magnitude and percentage of the turbines behind the 7.5m and
12.5m turbines compared to the turbine behind the 10m turbine. It shows that for all
wind speeds the greatest percentage power variation is at wind speed of 7m/s. The
greatest power different was at 20m/s, where an extra 13kW was extracted by the
turbine behind the 7.5m turbine compared to the turbine behind the 12.5m. It shows
that the turbine behind the 7.5m turbine extracts nearly 65% more power then the
turbine behind the 10m turbine and 2.44 times more power than the turbine behind
the 12.5m turbine at 7m/s. Increasing the wind speed reduces the power difference
between the turbines to approximately 5%.
Similar to the 60m spacing situation to get the full picture both turbines need to
be considered together. Figure 5.35 shows the combined power of the first and second
turbines taken from Table 5.5 and Figure 5.34 for each wind speed. Table 5.15 shows
the difference in power in terms of magnitude and percentage between configurations
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 695 64.9 -346 -32.3
10 2300 31.4 -1070 -14.7
15 4170 11.7 -1900 -5.32
20 9290 11.0 -4170 -4.94
Table 5.14: Difference in power of the second turbines compared to 10m turbine to 3
s.f.
Figure 5.35: Combination of the power extracted by the first and second turbines at
different wind speeds
with different diameters compared to the baseline configuration containing just 10m
turbines. Figure 5.35 and Table 5.15 shows that the greatest power variation both
in terms of magnitude and percentage difference occurs at 20m/s and the smallest
variation occurs at 7m/s. The difference between each configuration is significant with
differences ranging from approximately 1.5 − 110kW between the 7.5m and 12.5m
turbine configurations depending on the wind speed. It shows that the 12.5m turbine
configuration extracts the greatest amount of power and the 7.5m configuration extracts
the least.
Table 5.16 shows the power deficit of the second turbines compared to the first. It
Rotor size [m]

Wind speed [m/s] 7.5 % 12.5 %
7 - 543 -10.9 1010 20.4
10 - 3380 -15.5 5620 25.8
15 - 19000 -21.2 27500 30.8
20 - 45700 -21.6 68000 32.2
configurations with the baseline configuration to 3 s.f.
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 33.5 72.5 86.2
10 -8.21 49.8 70.7
15 -31.1 33.3 59.3
20 -30.5 33.4 59.6
Table 5.16: Power deficit of second turbine compared to the first turbine to 3 s.f.
shows that all turbines behind the 10m and 12.5m turbines experienced a power deficit
of approximately 33 − 86% depending on the turbine size and wind speed. It shows
that turbines behind the 12.5m turbine experiences the largest power deficit whatever
the wind speed. The turbine behind the 7.5m turbine has the least power deficit of 33%
at 7m/s and for all other wind speeds the second turbine extracts more power then the
first turbine by up to 31%.
5.5.8 Power comparison between spacings
Table 5.17 shows the difference in power of the second turbine at difference spacings
and at various wind speeds. It shows that for all but two situations the 60m spacing
extracts more power than the 30m spacing. It shows that for at the lowest wind speed
of 7m/s there is a 6.5 − 13% difference in power of the second turbine. However, for
all other situations there is very little difference in power extraction with a maximum
difference of 3%.
Rotor size [m]

Wind speed [m/s] 7.5 % 10 % 12.5 %
7 124 6.55 133 11.0 109 13.1
10 311 3.14 172 2.30 - 82.1 1.34
15 990 2.42 386 1.07 - 76.7 0.227
20 2350 2.45 922 1.08 206 0.256
Table 5.17: Comparison of the power extracted by the turbines with a spacing of 60m
and 30m to 3 s.f.
Chapter 6
Discussion
This chapter discusses the results of the five studies presented in Chapter 5. It was
shown that the actuator disc performed well when compared with experimental data
in the far wake for each of the cases simulated with the near wake predictions being
poorer. Here quantitative and qualitative critical evaluation of the results of each study
will be preformed and common discrepancies between the models analysed using data
from previous chapters and independent sources evident in the literature. This penul-
timate chapter then closes with an examination of the rotor diameter study results,
including the possible affect these may have in regard to future developments.
6.1 Computational techniques
This section reviews the computational processes employed throughout the work de-
scribed in this thesis and critiques reasoning and associated consequences. All the
work described in this thesis was conducted using CFD methods using two primary
computational codes together with three different rotor models.
171
CHAPTER 6. DISCUSSION 172
6.1.1 RANS and LES
The RANS equations were used together with the k-ω SST turbulence model for each
simulation. The reasoning for solving the RANS equations as opposed to a LES
method was twofold, each relating to computational expense. LES is considered more
accurate as turbulence is resolved more precisely, particularly anisotropy. LES is also
transient requiring the solution to be resolved in both time and space. However, this
additional accuracy significantly increases the computational time. RANS is computa-
tionally less expensive as it allows steady state simulations to be conducted and is less
dependent on the particular discretization process. Considering the main rotor method
implemented in this work was the actuator disc, itself a time averaged representation, it
was considered that the RANS solution more appropriate. The justification is evident
in this thesis (§5.3 and §5.4) where no significant changes have been observed when
compared with modelling analogs (e.g. LES [195] and DES [120][121]).
6.1.2 Cell type
Subsections 4.2.5 and 5.2.3 presented and compared results of two different cell types
available. It was found that there was little benefit of using quadratic cells instead
of linear cells for the cases investigated here. The quadratic cells predicted a slightly
higher velocity deficit at the centre of the wake for all mesh densities, although no
quadratic mesh predicted a higher deficit then a finer linear mesh. The linear cells
showed a reduction in turbulence intensity after the peak (Figure 5.9) this is likely due
to the cells adding some numerical diffusion to certain regions of the model, as all
models used the same advection scheme (high resolution scheme), thus inaccurately
predicting the three-dimensional turbulent behaviour of the flow field. However, a
maximum variation of only 3% was observed, showing the linear cells were adequate
to predict the turbulence in the cases examined here. The major difference between the
cell types was in the computational time. The quadratic cells took at least 60% more
iterations to solve which translated into a minimum of 48% longer to solve (Table
5.2). Considering the results from these mesh densities, there was no clear additional
accuracy observed considering the additional computational expense. Whence there
was no reason to use quadratic cells instead of linear cells given the cases considered.
This is due to the relatively simple geometry and flow field which has no significant
recirculation or large scale turbulence structures.
Every CFD simulation requires boundary conditions which can have a significant im-
pact on the solution. The importance of the appropriate boundary conditions can not
be understated although their impact are not always known.
As part of the porous disc study (§§4.2.2 and §5.2) four boundary conditions were
investigated, two turbulence length scales and two different definitions at the domain’s
upper boundary (representing a duct and channel). It was observed that the duct pro-
duced a higher central velocity magnitude and marginally lower turbulence intensity
than the channel flow. This is to be expected due to the presence of the additional wall
at the top of the domain. This wall creates an additional boundary layer which restricts
and slows the flow near the wall. This deceleration along the wall forces the flow at
the centre to increase in velocity in order to maintain the same mass flow rate which
is visible in Figure 5.2. However, the influence of the approximately inviscid central
core (0.5 < y/2R < 2.5) is minimal, meaning that modelling an open channel flow as
a duct incurs minor errors, whilst reducing computational expense.
The resistance of the disc, in the porous disc study (§§4.2.2 and §5.2), was pre-
scribed using a predefined thrust coefficient (CT ). The differences between the cal-
culated CT values from the channel and duct flow can be attributed to the additional
boundary layer (i.e. a new wall boundary condition at the upper boundary) and subse-
quent small velocity increase, in the aforementioned inviscid core. The CT values were
calculated using a free stream velocity based on the inlet velocity which is reasonable
for the channel flow, as this velocity profile is fully developed. However, this is not
the case for the duct flow; taking this into account the CT values were recalculated, as
shown in Table 6.1, using a new free stream velocity of 0.34m/s obtained at the domain
origin in the absence of the disc. Considering these new CT values the duct situation
predicts the CT closest to the experimental data [185].
Study
Disc Experimental Numerical Channel Duct New
num- [185] [141] duct
ber
1 0.61 0.65 0.60 0.64 0.61
2 0.86 0.91 0.86 0.93 0.88
3 0.97 0.98 0.93 1.00 0.96
Table 6.1: New CT values
Similar results were found in the actuator disc study (§4.3 and §5.3) with the walls
defined as both free slip and nonslip. The different boundary conditions had little
influence on the results, as the boundary layer did not develop enough to influence the
velocity at the disc.
Two different turbulence definitions were investigated by defining the turbulence
length scale (§§4.2.2 and §5.2). Both definitions produced similar turbulence intensi-
ties at the inlet although there was a significant difference behind the disc. The second
definition referred to as inlet 2 (§4.2 and §5.2) predicted a better velocity and turbu-
lence intensity profile when compared to the experimental data [185]. This is due to a
reduction in the turbulent dissipation throughout the domain, prolonging the turbulence
generated at the inlet and disc which was overly dissipated using inlet 1 (§4.2).
6.1.4 ANSYS-CFX and ANSYS-Fluent
As part of this thesis two commercially available CFD codes were used, specifically
ANSYS-CFX and ANSYS-Fluent. The user customization of each code is slightly dif-
ferent, ANSYS-CFX uses CEL and ANSYS-Fluent uses UDFs, discussed in section
2.2. UDFs are considered more flexible although for many applications this is not a
consideration1 . UDFs need to be complied into ANSYS-Fluent and they require more
lines of code to be written so take longer to write and debug. In terms of implemen-
tation the two codes are very similar and although the author found CEL easier and
quicker, it is personal preference.
Subsection 5.2.4 presented results of both codes on analogous models of a porous
disc within a water channel. The two codes predicted very similar results, not enough
for either to be considered more numerically accurate. What differences did exist had
minimal affect on the overall flow field. ANSYS-Fluent predicted almost identical
velocity profiles for both channel and duct simulation away from the boundary but a
large difference at the boundary. This shows that the boundary layer did not develop
enough to change the velocity profile, leaving the inviscid core intact (Figure 5.11).
However, ANSYS-CFX predicted a slight difference between the duct and channel
velocity profiles and hence a larger boundary layer, which had an affect on the velocity
profile although this was minimal at the centre of the domain. ANSYS-CFX predicted
a slightly higher peak in turbulence intensity compared to ANSYS-Fluent, this may be
due to boundary layer prediction as the peak is due to the merger of the shear layer
between the wake and free stream velocity.
The main difference between the two codes was in the time it took for the simu-
lations to converge. ANSYS-CFX took a few hundred iterations and approximately
one hour to converge whereas ANSYS-Fluent took several thousand iterations and a
1 The author has not found a situation which is beyond CEL
few hours to converge2 . It appears that ANSYS-Fluent processed each iteration much
quicker than ANSYS-CFX however, it took ANSYS-Fluent considerably more iter-
ations to converge. The reason for this difference is that ANSYS-CFX has a cou-
pled solver which solves for all variables at once by creating a global solution matrix.
However, the ANSYS-Fluent simulations were conducted using a segregated solver
which solves for each variable sequentially [54]. Coupled solvers converge signif-
icantly quicker although they require more memory (1.5 − 2 times) as they need to
store data on all variable not just one [54].
Considering that there was no numerical difference between the two codes and the
additional computational expense ANSYS-Fluent required, ANSYS-CFX is the more
efficient.
6.2 Rotor models
Within this thesis three different rotor models have been described and implemented,
which have been compared with other numerical models and experimental data. This
section will attempt to critically evaluate the differences between the models and any
discrepancies found between these and analogous experimental data evident in the
literature.
6.2.1 Near wake discrepancy
Throughout this thesis, discrepancies have been indicated when comparing numeri-
cally simulated results with experiments especially in the near wake. Although good
agreement has been demonstrated with comparisons made with other numerical stud-
ies. Some of the salient reasons for such apparent inconsistencies (or otherwise) will
be discussed in a little more detail throughout this section.
2 Simulations were conducted using a dual core 2.27Ghz Intel Xeon desk top computer.
Porous disc
The porous disc study predicted a greater velocity deficit and a turbulence intensity
peak at far lower magnitude and further downstream of the disc than observed experi-
mentally [185]. This is due to the definition of the discs within the model as opposed to
the physical discs. The discs within the experimental study [185] extracted momentum
from the flow by converting the velocity into small scale turbulence and thus creating
a high level of turbulence behind the disc and produced a variable three-dimensional
momentum loss. However, the porous discs used in the simulations were defined with
an isotropic 1D momentum loss which extracted momentum from the flow explicitly,
reducing the velocity with no added turbulence. This explains the relatively high levels
of turbulence behind the rotor for the experimental data [185] and the lack of this peak
in the modelled results. The turbulence intensity, of the simulations peaked, further
downstream than observed experimentally (Figure 5.4), this is due to the merger of
the shear layers at the edge of the wake created by the velocity deficit; this delayed
turbulence also explains the velocity deficit. The higher turbulence levels in the exper-
imental study increased the velocity recovery, resulting in a greater velocity deficit in
the simulations.
Application of the one-dimensional momentum source as a predefined constant
unidirectional loss introduces an inherent flaw into the porous disc modelling proce-
dure. Allowing for this, these models performed well and predicted some characteris-
tics of the velocity profile and turbulence levels remarkably well. The model generally
predicted both the velocity and turbulence intensity magnitudes lower than the exper-
imental data [185] in the near wake region, with the discrepancy reducing as the flow
moves downstream.
Actuator disc
In the actuator disc study, the results obtained were compared with another actuator
disc model and experimental data. All three data sets showed similar velocity profiles,
although both actuator disc models under-predicted the magnitude of the velocity drop
just behind the turbine. This inconsistency is not just predicted by the actuator disc but
has been observed in other numerical studies which have tried to replicate the experi-
ment [118][119]. Scheper et al. [119] noted that no study was able to predict both the
loads and velocities. Other researchers such as Réthoré et al. [120] noted that either
the same force distribution or velocity profile may be predicted but not both, to date
there is little physical interpretation for this phenomenon. However, a possible reason
maybe that the force exerted by the turbine is greater than the measured values along
the blade suggests. This being evidenced by predicting a closer wake comparison by
increasing force along the blade. The data produced in this study predicts the highest
velocity before the turbine which is due to the inlet geometry and defined inlet velocity,
which had already been reduced compared to the experimental data.
This discrepancy is due to the near wake being almost completely characterized by
the turbine’s geometry which results in a transient flow region, not fully predicted by
a time averaged model like the actuator disc. This is particularly well demonstrated in
previous studies of the MEXICO data [118][119] where the transition between aero-
foil sections created unexpected velocity fluctuations which were not predicted by any
model but affected the near wake.
6.2.2 Actuator disc
Even though the actuator disc has been the main model used in this thesis there are
various types within the literature (§2.2.2).
Normalized mean velocity

x/R This study Réthoré et al. [120] Difference (%) Difference (m/s)
-0.4 0.919 0.955 - 3.78 - 0.351
0.0 0.873 0.841 3.80 0.353
0.3 0.818 0.710 15.2 1.74
0.6 0.818 0.764 7.02 0.652
1.0 0.804 0.739 8.81 0.819
1.5 0.809 0.735 10.0 0.930
2.0 0.826 0.688 20.2 1.88
Table 6.2: Average velocity difference between this study and Réthoré et al. [120]
Within the actuator disc study two different discs were compared in terms of veloc-
ity profiles. Very good agreement was found along the centre line and wake expansion.
The greater discrepancy was observed at the cross-sections behind the turbine at the
wake centre. The actuator disc in this study predicted little velocity deficit at the cen-
tre of the wake, whereas Réthoré et al. [120] predicted their highest velocity deficit
at this position. This discrepancy is due to the description of the model momentum
sink, specifically the hub and nacelle which was excluded in this study but was pre-
sumably included by Réthoré et al. [120]. This explains the difference between the
velocity profiles and why Réthoré et al. [120] predicted almost zero velocity at the
centre just behind the turbine. Further evidence of this is shown in Table 6.2, demon-
strating that Réthoré et al. [120] predicted a lower velocity at the rotor (approximately
4%) and hence greater momentum loss, through the inclusion of the hub. Despite the
higher velocity deficit just behind the rotor the average velocity was less than 10% for
0.6 < x/R < 1. The discrepancy for x/R ≥ 1 increases as the mesh used in this study
was only refined up to x/R = 1, so this was expected.
6.2.3 Actuator line
The actuator line method (§3.5.2) is a transient model which represents the turbine as
lines instead of a disc. As part of this thesis the AD and AL methods were implemented
(a) (b)
Figure 6.1: Force distribution of (a) actuator disc and (b) actuator line methods
within ANSYS-CFX. Figure 6.1 shows the force distribution of the two methods. The
ALM is able to predict the tip and root vortices shown in Figure 6.2 which the ADM
cannot. This additional accuracy is at the expense of computational time with the ADM
taking a matter of hours but the ALM taking days.
Comparison between the ADM and ALM was conducted as part of the two turbine
study (§4.4). In this instance the ALM was conducted by another author [120]. The
results (§5.4) showed good agreement between the models in terms of the velocity
profile behind the second turbine, showing the additional computational expense is
unnecessary for turbine interactions. The ALM shows little advantage over the ADM
in the far wake as the modelled vortices have broken up as part of the transition to the
far wake.
Considering turbine interactions, specifically two turbines with 30m spacing, the
velocity profiles at 50m and 70m behind the first turbine were compared (§4.4); at 50m
the results were in very good agreement. There is less agreement at 70m particularly
(a) (b)
Figure 6.2: Isosurface of the vorticity (a) actuator disc and (b) actuator line methods
at the centre of the wake although the edge of the wake remains in good agreement.
It appears that at this distance behind the turbine, the wake has recovered more in this
study, than in the numerical study detailed the literature [195]. There are two reasons
for this discrepancy, being far behind the wake. Firstly, the turbulence which can not
be compared may be different. Secondly, it could be due to the differences between the
models, specifically the lack of tip vortices within this study, which may have caused
the wake to destabilize and recover; this has been shown before to influence the wake
[198].
6.2.4 Ground affect and shear
The inlets used in the work described in this thesis are nearly all defined with an uni-
form velocity profile. For the turbines interaction studies ground effects were not con-
sidered. This was justified (§5.4) as good agreement was found with both experimental
and numerical studies which had considered these effects, and showed they had no sig-
nificant influence on the velocity field. Zahle & Sørensen [156] showed that including
the ground effects without a shear inflow increased the power extracted by 0.4% and
adding a shear inflow increased this further to 5.5%.
6.3 Scaling
Scaling is a significant part of fluid mechanics as there is no point in building a full

scale plane just to find out if it can fly. Instead smaller, less expensive, scale models
or prototypes are created which should obey what are known as scaling laws, which if
satisfied the condition of similarity exists between models.
Scaling is accomplished through a branch of fluid mechanics known as dimen-
sional analysis. A key concept in dimensional analysis involves reducing the number
of flow variables, commonly to non-dimensional variables such as the Reynolds num-
ber (§1.5.2) [18]. There are several methods of reducing the dimensional variables
of the flow, the Buckingham pi-theorem [206] being the most commonly used. Mod-
els are considered completely similar if all relevant dimensionless parameters have
the same values. This is not always possible, instead particular similarity is normally
considered such as:
• geometric similarity where all body dimensions have the same length-scale ratio;
• kinematic similarity where models have the same length-scale and time-scale
ratios;
• dynamic similarity where models have the same length-scale, time-scale and
force-scale ratios.
For incompressible flow with no free surface two models have dynamic similarity if
the Reynolds numbers are equal [18].
6.3.1 Reynolds number
To maintain the Reynolds number when scaling, the velocity is increased to compen-
sate for the reduction in length. In practice it is not always possible to maintain the
Reynolds number when scaling a model. Increasing the velocity can create additional
effects such as cavitation or compressibility. One possibility would be to change the
viscosity as well as the velocity but this is not normally practical.
Scaling problems have been a factor for two studies in this thesis, specifically the
porous disc and interaction studies. In the porous disc study geometrical similarity was
not possible as it would have resulted in unrealistic swirl and hence cavitation in the
flow due to the high RPM required to maintain the required tip speed ratio [185].
There were similar scaling issues in the two turbine study. All three studies were of
course geometrically similar featuring a scaled turbine although the Reynolds number
were different. The three studies had Reynolds numbers of 660 thousand, 4.5 million
and 40 million for the experimental study [197], for this study and numerical study
[195] respectively. Similar Reynolds numbers would have caused compressible effects
in the experimental study (inlet velocity ≈ 710m/s and Mach number ≈ 2 ). If the
same Reynolds number was maintained in this work it would have caused unrealistic
operating conditions with an inlet velocity of 64m/s which is 2.5 times the cut off
point. This increase in velocity requires the rotor velocity to increase to maintain the
tip speed ratio. The rotor RPM needs to increase to 343RPM which pushes the relative
velocity at the blade tip to cause compressible effects (Ma ≈ 0.6).
However, there is a logarithmic relationship between the rotor diameter and the
Reynolds number of the three studies, shown in Figure 6.3 which explains the good
agreement shown in Figure 5.20.
Figure 6.3: Log relationship between the Reynolds number and rotor diameter
6.4 Rotor diameter
This section discusses the results detailed in section 5.5 which investigated the influ-
ence rotor diameter size has on wind turbine wakes and the aerodynamic performance
on a downwind turbine.
To begin with the first turbine was analysed in isolation, here it was found that as
the diameter increased the average velocity at the rotor decreased. This decrease did
not significantly influence the power which increases as the rotor diameter increases
due to increase in swept area being far greater than the slight drop in velocity. The
three turbines have a diameter of 7.5m, 10m and 12.5m which corresponds to a swept
area of 44.2m2 , 78.5m2 and 123m2 respectively which is an increase of 43.7% and
36.2% compared to a decrease in velocity of less than 7%. To investigate the affect the
rotor diameter has on turbine interactions a baseline needed to be established. In this
study the two 10m turbines (§4.4 and §5.4) were used for such purposes.
6.4.1 Wake
A downwind turbine experiences a velocity deficit and hence a power deficit compared
to the first. It was shown that a turbine behind the 7.5m turbine experiences the smallest
velocity deficit and narrowest wake. It had a minimum of 4% increase in velocity
compared to the other turbines if spaced 60m apart and 3.5% when spaced 30m apart.
Whereas, the turbine behind the 12.5m turbine experiences the highest velocity deficit
and widest wake. This was expected as the 12.5m turbine extracts the most power and
hence momentum reducing the velocity the most and as it is the largest it produces the
widest wake. Figure 6.4 visualizes this showing the velocity contours on the yz−plane.
It shows the interactions between the wakes and the difference in wake width between
the three configurations. Figure 6.5 shows the streamlines through the domain and
that there is more swirl in the flow after the second turbine which improves velocity
recovery.
6.4.2 Power
It was shown (§§5.5.4 and §§5.5.7) that this velocity increase does translate into a
power increase at the second turbine of a minimum of 12.5% if spaced 60m apart
and 11% when spaced 30m apart. This is expected as the only variation in the power
equation is the velocity component which is cubed.
Combining the power of both turbines shows that despite the increase in power of
the second turbine, the 7.5m turbine configuration extracts the lowest power overall
and the 12.5m turbine configuration extracts the most power. This is because although
the second turbine in the 7.5m configuration generates at least 11% more then the
other turbines the first turbine extracts at least 30% less power. Considering the three
configurations in this study it can be seen that rotor diameter can have a significant
affect on the power extracted by as much as 40%.
(a)
(b)
(c)
Figure 6.4: Velocity contours on the YZ plane for the (a) 7.5m configuration, (b) 10m
configuration and the (c) 12.5m configuration
(a)
(b)
Figure 6.5: Streamlines viewed along (a) the yz−plane and (b) the xy−plane
Comparing the results of the different spacing situations showed that the 60m spac-
ing extracted more power than 30m spacing. This is to be expected as the 60m spacing
allows the velocity deficit to recover for longer. However, the additional distance only
produced a maximum power difference of 3% for all wind speeds except for 7m/s
which was 13%; this slight increase in power does not justify the additional spacing in
terms of the energy density.
As previously discussed the second turbine experiences the largest velocity deficit
and hence power deficit of turbines within a wind farm [13][187][188]. This work
shows that the second turbine of the baseline configuration had a power loss of 33-
70% depending on wind speed and spacing. Increasing the first turbine size increases
the power loss of the second turbine to 60-84% depending on wind speed and spacing.
However, the 7.5m turbine configuration not only reduced the power loss, the second
turbine actually extracted up to 30% more power. This shows that decreasing the rotor
diameter can improve the power extracted by the second turbine and a 25% decrease
in rotor size results in the second turbine producing more power than the first.
Considering this a configuration with a 9m diameter rotor 60m in front of a 10m
turbine was simulated. Table 6.3 shows the power extracted by a 9m turbine con-
figuration. It shows that for this configuration the first turbine produces more power
than the second as normal. This suggests that a first turbine with a diameter between
7.5−9m may produce a configuration where both turbines can extract the same amount
of power. Comparing the 9m configuration to other configurations finds that the first
turbine does extract less power and the second turbines extracts more power than the
baseline. The 9m turbine has a swept area which is 19% smaller than a 10m turbine
which is reflected in the power extracted in the first turbine which is approximately
19% less than the 10m turbine. The second turbine extracts approximately 4-19%
more power than the 10m turbine. This results in the power extracted by the 9m con-
figuration to be 5-10% less than the 10m turbine configuration.
9m diameter
Wind speed [m/s] 1st 2nd Total
7 3400 1430 4830
10 12200 8350 20600
15 43500 37800 81300
20 103000 89200 192000
Table 6.3: Power of 9m configuration
Figure 6.6 shows the power extracted by each turbine at different wind speeds. It
visually shows the large variation in the power extracted by the first turbine due to the
rotor diameter and small variation in the second turbine due to the velocity variation.
It shows visually the power loss of the second turbine particularly behind the 12.5m
turbine which extracts the most overall power.
6.4.3 Power Coefficient
The results presented in sections 5.5.4 and 5.5.7 featured a fixed power coefficient (CP ).
CP was fixed to show the direct effect of velocity on the power however, in reality the
CP is not fixed but varies with wind speed and tip speed ratio. This section presents
modified results from section 5.5.4 featuring a variable CP . The CP was calculated
using equation (6.1). It varies based on the axial interference factor (a) which is based
on the ratio between the velocity at the rotor and the free stream velocity.
CP = 4a(1 − a)2 (6.1)

Ur
a = 1− (6.2)
U0
Table 6.4 shows the new CP values, for all except two situations the CP value is
lower than the 0.4 as employed in section 5.5. The 7.5m configuration experiences the
highest reduction in CP due to the increase in the velocity at the turbine. This reduction
Figure 6.6: Power of each turbine at different wind speeds. Blue is the 7.5m turbine,
green is the 10m turbine, red is the 12.5m turbine and light blue is the 9m turbine.
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 0.323 0.444 0.507
10 0.205 0.2795 0.332
15 0.189 0.199 0.208
20 0.192 0.199 0.195
Table 6.4: New power coefficient values to 3 s.f.
in CP relates to a reduction in extracted power, shown in table 6.5, The exceptions

being two cases where a higher CP was calculated. Table 6.6 shows that a variable
CP exaggerates the differences between the turbines, with the 7.5m turbine extracting
60% less than the baseline (2 − 18% more than with a fixed CP ) and 12.5m turbine
extracting 1.53 − 1.74 times the baseline (3-20% higher than with a fixed CP ).
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 2140 4320 6650
10 4550 10100 17600
15 14400 26700 43100
20 34500 63300 97200
Table 6.5: Approximate power extract by the first turbine with a variable CP to 3 s.f.
Variable 7.5 m turbine 12.5 m turbine

Wind speed [m/s] Magnitude % deficit Magnitude factor
7 - 2180 50.4 2330 1.54
10 - 5580 55.1 7500 1.74
15 - 12300 46.0 16400 1.61
20 - 28700 45.4 34000 1.54
Table 6.6: Power difference of the first turbines compared with the baseline using a
variable CP to 3 s.f.
When considering the second turbine there is an issue with the definition of the CP
value, as it involves the power in the wind. As the first turbine has extracted power from
the wind the free stream velocity is not applicable for the second turbine. Considering
this a new CP value which is referred to CP2 was, calculated using a velocity one
diameter in front of the turbine (referred to as the local free stream velocity), before
the velocity is influenced by the turbine as shown in Figures 5.25 and 5.26. Table 6.7
shows the CP2 values for each turbine and that for each situation (except the 12.5m
configuration at 7m/s) the CP2 < CP . This translates into a power deficit for each
situation, save the aforementioned exception (Table 6.8).
Rotor size [m]

Wind speed [m/s] 7.5 10 12.5
7 0.2530991 0.4074242 0.5405362
10 0.1843507 0.2514721 0.3081106
15 0.1879297 0.1950545 0.2045937
20 0.1907710 0.195539 0.1906686
Table 6.7: CP values calculated for the second turbine referred to as CP2
Rotor size [m]

7.5 10 12.5
Wind speed [m/s]
CP CP2 CP CP2 CP CP2
7 1526 1196 1335 1226 1055 1126
10 5084 4566 5207 4698 5103 4734
15 19343 19192 17972 17593 17502 17246
20 46212 45881 42617 41786 39376 38414
Table 6.8: Approximate power extract by the second turbine for different CP values
Table 6.9 shows a comparison between the power extracted by the second turbines
using the variable CP values. It shows similar characteristics to the fixed CP results.
Using a variable CP the results are much closer with a maximum variation of 10%
from the baseline configuration (except 12.5m configuration at 7m/s). These results
are much closer now due to the variation in the power coefficient of the device.
Table 6.10 shows the power deficit of the second turbine compared to the first.
Here, very similar predictions to the fixed CP value data (Table 5.11 and 5.16) are
observed, with the 12.5m configuration having the largest power deficit of 60-83% and
the 7.5m configuration extracting more power from the second turbine.
Rotor size [m]

7.5 12.5
Wind speed
CP CP2 CP CP2
[m/s]
Total % Total % Total % Total %
7 190 14.3 -29 -2.45 -279 20.9 -99 -8.10
10 -122 -2.36 -131 -2.81 -103 1.98 36 0.78
15 1370 7.63 1598 9.09 -470 2.62 -346 -1.97
20 3595 8.44 4095 9.80 -3240 7.60 -3371 -8.07
Table 6.9: Difference in power of the second turbines compared to 10m turbine
Rotor size [m]

Wind speed 7.5 10 12.5
[m/s] CP CP2 CP CP2 CP CP2
7 28.8 44.2 69.1 71.6 84.1 83.1
10 -11.7 -0.334 48.6 53.6 71.1 73.1
15 -34.3 -33.3 32.6 34.0 59.3 59.9
20 -33.8 -32.9 32.6 33.9 59.5 60.5
Table 6.10: Power deficit of the second turbine
Table 6.11 combines the power of the two turbines to give their overall power
extracted. It shows little difference between the power extracted given the two different
power coefficients. The maximum variation is 10% for the 7.5m configuration at 7m/s,
all other power values have a variation of below 5%. This shows that using the free
stream velocity to define the CP is reasonable for both turbines even though there is
less power for the second to extract. Table 6.12 shows the comparison between the
different configurations. Furthermore, similar results to the fixed CP results with the
12.5m configuration extracting the most power and the 7.5m configuration extracting
the least are also shown. Comparing the results for different power coefficient values
shows that using a variable coefficient of performance predicts higher power extraction
at low values for the 10m and 12.5m configurations but lower power for all other
situations. A variable CP also predicts a greater variation between the configurations
of up to 50%, due to a larger variation predicted for the first turbine.
Rotor size [m]

7.5 10 12.5
Wind speed [m/s]
CP CP2 CP CP2 CP CP2
7 3669 3339 5653 5544 7709 7779
10 9635 9117 15340 14831 22738 22369
15 33744 33593 44642 44262 60553 60297
20 80748 80416 105875 105044 136630 135669
Table 6.11: Combined power of the first and second turbines for two variable CP values
Rotor size [m]

7.5 12.5
Wind speed
CP CP2 CP CP2
[m/s]
Total % Total % Total % Total %
7 -1984 -35.1 -2205 -39.8 2055 36.4 2235 40.3
10 -5705 -37.2 -5714 -38.5 7397 48.2 7538 50.8
15 -10897 -24.4 -10669 -24.1 15911 35.6 16035 36.2
20 -25127 -23.7 -24627 -23.4 30755 29.0 30625 29.2
configurations with the power extracted by the 10m configuration for two variable CP
values
6.4.4 Overall performance
There are different ways to measure the performance of a wind-farm. The power pro-
duction is normally the most important and in this case the 12.5m configuration clearly
extracts the most. A wind-farm’s performance is also measured in terms of power ex-
tracted related to the installed capacity. Considering these cases in that regard and
considering the first turbine to be operating at its installed capacity, the 7.5m config-
uration performs best. It performs at least 8% closer to its installed capacity then the
other configurations over all wind speeds using a fixed CP and 3% closer using a vari-
able CP . This is because it features the least amount of wake interactions and each
turbines performs closest to its optimum. Considering these two measures the 12.5m
configuration extracts the most power but the 7.5m configuration performs closest to
installation.
6.4.5 Wind-farm
It has been shown that reducing the diameter of the first turbine can improve the second
turbines performance by as much as 65% at low wind speeds and 10% at high speeds
(§5.5.4). However, decreasing the diameter of the first turbine can also significantly
decreased the pairs overall power. However, the velocity behind the second turbine is
significantly higher (up to 10%) which would improve the performance of additional
turbines, so would improve the performance of more than one turbine. This could
allow the smaller 7.5m configuration to improve the performance of more turbines and
make up for its own loss in power. Unfortunately, only two turbines were considered
in this study but considering a full wind-farm or a row could have benefits not found
for two turbines. It was shown that using larger turbines only effected the final velocity
behind the second turbine by 4% compared to the baseline and it was shown that the
36% increase in swept area of the larger turbine more than makes up for the reduced
velocity. It was also shown that the presence of the second turbine increased the wake
recovery. To this end there may be scope of using a range of turbines to extract more
power, with some large turbines to extract the energy and smaller ones which promote
wake recovery while still extracting some power.
Chapter 7
Conclusions & Recommendations
This final chapter describes the major conclusions of the five studies conducted as part
of this thesis, highlighting the novel findings. The chapter concludes with recommen-
dations for further work which could be conducted to extend the work presented in this
thesis.
The aims and objectives of this thesis as described in Section 1.2 have been satis-
fied. Objective 1 (§1.2) was satisfied in Section 2.2 where three rotor models were de-
scribed and implemented within ANSYS-CFX and ANSYS-Fluent. Objective 2 (§1.2)
was satisfied in sections 5.1, 5.2 and 5.3 where the methods employed in this thesis
were benchmarked with previous theoretical, experimental and numerical work. Ob-
jective 3 (§1.2) was satisfied in section 5.4 where the aerodynamic interactions of two
similar axially aligned turbines was predicted. Objective 4 and 5 (§1.2) were satisfied
in subsections 5.5.4, 5.5.7 and 5.5.8 where the performance and power extraction of
two aligned turbines with different rotor diameters was presented. Objective 6 (§1.2)
was satisfied in Section 6.4 where possible benefits of using different sized turbines
was discussed. Objective 7 (§1.2) was satisfied through the work presented at two con-
ferences [1][2] and published within peer-reviewed journals [3]. Through these objec-
tives the aims of the thesis were satisfied and a number of conclusions were reached
196
CHAPTER 7. CONCLUSIONS & RECOMMENDATIONS 197
which are detailed in this chapter.
Novel contributions to knowledge

The novel contributions to knowledge provided by this thesis are concerned with
the aerodynamic interaction between two turbines with different rotor diameters and
how this affects performance. To the authors knowledge, this has not been investigated
before either numerically or experimentally in any detail before. The work described
throughout this thesis found that changing the rotor diameter of the first turbine1 has a
significant affect on both the downstream turbine and the performance of the pair as a
whole; resulted in:
• a velocity increase at the second turbine by as much as 18% thereby resulting in

a 65% increase in power extraction;
• conversely, the velocity can decrease as much as 12% thereby resulting in a 32%
decrease in power extraction.
However, it was found that the performance of the second turbine was not as significant
as the first turbine in terms of overall power extraction.
• A reduction in the first turbine’s rotor diameter decreases the overall power ex-
tracted by the pair by between 10-20%.
• An increase in the first turbine’s rotor diameter increased the overall power ex-
tracted by the pair by 20-32%.
Other possible novelty within this thesis may be to attributed implementation of the
actuator disc and actuator line methods within commercially available CFD software,
specifically ANSYS-CFX and ANSYS-Fluent. The work described through out this
thesis, is certainly one of the first independent attempts to conduct simulations of tur-
bine interactions using ANSYS-CFX and most probably the first to use actuator meth-
ods which are entirely calculated by this code2 . Hence, the thesis contains the most
detailed and comprehensive description of the implementation process of the actuator
methods within ANSYS-CFX and ANSYS-Fluent.
1 (by ±25% in this case)
2 Before starting this project the author could not find any literature in regard to the actuator methods
within commercially available CFD codes. However, recently some preliminary work in this area has
come to light [87]. Here, the actuator disc and line methods were implemented in ANSYS-CFX; a single
turbine using prescribed forces from an external program was simulated.
7.1 Conclusions
A part from the key contributions to knowledge there are a number of other conclu-
sions. The first part of this thesis involved implementation and benchmarking the rotor
methods detailed in section 2.2. The benchmarking was conducted using previous the-
oretical, numerical and experimental data for turbines in water and air. The second part
involved investigating the aerodynamic interactions of two turbines and how the rotor
diameter of the first turbine affects these interactions.
7.1.1 Porous disk in water
Porous discs have been used both experimentally and numerically to represent turbines.
This study compared CFD results of a porous disc in water with that of experimental
[185] and numerical [141] studies and found:
• the modelling method was sufficient to predict the far field velocity characteris-
tics of a porous disc;
• the CFD models was unable to match the velocity deficit and turbulence magni-
tude just behind the disc of the experimental study [185], although the agreement
improved further downstream. This discrepancy was attributed to small scale
turbulence present in the experiments and the momentum extraction method em-
ployed by the models.
The main discrepancies observed were attributed to the difference between the CFD
and experimental setup in terms of the momentum source. In the CFD simulations
momentum was extracted explicitly rather than converting it into small scale turbulence
as in the experimental case [185]. Moreover, only an unidirectional momentum source
was used, which did not account for the three dimensional effects of the real discs used
in the experimental study [185].
7.1.2 Influence of the boundary conditions
The influence of the boundary conditions were investigated using two different turbu-
lence definitions as well as two different definitions of the top of the domain represent-
ing a duct and channel.
• The two turbulence boundary conditions produced very similar inlet profiles but
different profiles downstream of the disc, showing the importance of the defini-
tion of turbulence throughout the domain.
• It was found that the channel and duct simulations predicted very similar results
with the duct predicting a slightly higher velocity magnitude for the majority of
the domain which was discussed in subsection 6.2.1.
It demonstrated the usefulness of the duct boundary conditions (for computational

ease) for representing a porous disc in open channel flow when simulating far field
effects; given the particular velocity profile applied at the inlet.
7.1.3 Actuator disc in air
The actuator disc method is a very popular method for modelling turbines, however the
majority of CFD studies use in-house software. Previous experimental and numerical
studies were compared with CFD results conducted as part of this thesis. It featured an
actuator disc in air implemented within the commercially available CFD code ANSYS-
CFX and found:
• good agreement in terms of velocity profile compared with previous experimen-

tal and numerical studies. The main discrepancy was just behind the rotor and
due to the force defined by the CFD model being too small to accurately predict
the velocity deficit measured in the experimental data. This is likely the cause
as it was found by other studies and my inlet velocity was shown to be too high
compared to other data sets.
• For this particular experiment it is possible to produce the correct load profile
along the blade or the correct velocity deficit but not both, which was also found
by others.
• Good agreement with previous numerical results in terms of average velocity

profile behind the turbine and hence momentum loss. Difference in velocity
profile is due to the distribution of the force along the turbine, specifically at the
hub as the nacelle was not considered in the work in this thesis.
This study showed it is possible to implement the actuator disc methods within com-
mercially available codes such as ANSYS-CFX. Therefore negating the need to use
specially developed codes to implement the actuator methods.
7.1.4 Two turbine interactions
Within the literature there is no definitive two turbine study or benchmark, as such
CFD results of two actuator discs were compared with normalized data from a previous
experimental study of two model turbines and a previous numerical study of two full
size turbines using the actuator line. Although each study featured different turbines
both in terms of size and blade geometry and at different wind speeds it was found
that:
• there was a good agreement between all data sets in terms of normalized velocity
profiles. This shows it is possible to compare normalized data of different turbine
studies despite the different rotor diameters and blade geometry.
• Equally, the results of this study were predicted using an actuator disc and RANS
to get a time averaged prediction of the turbines. In comparison the experimental
study used model turbines and the numerical study used the actuator line method,
both of which would have had tip and root vortices which the actuator disc can
not predict. This shows that the influence of these vortices is not significant for
the interactions between two turbines. This is compelling as the actuator disc
method requires significantly less computational resources than the actuator line
method.
• The second turbine experiences a power deficit of approximately 32-72% de-

pending on the wind speed.
7.1.5 Influence of rotor diameter
It is known that the second turbine experiences the largest velocity deficit and hence
power deficient within a wind-farm. Therefore improving its performance should im-
prove the overall performance of the wind-farm. There are a number of approaches
investigated within the literature such as yawing the first turbine or reducing its RPM.
This study is the first to investigate the affect of the first turbine’s rotor diameter has
on performance. The key conclusions have already been detailed in the novel contri-
butions to knowledge section at the start of the chapter, the other conclusions found
are:
• that increasing the rotor diameter decreases the average velocity at the rotor how-
ever, the power increases as the increase in swept area is greater than the decrease
in velocity;
• the second turbine improves the velocity recovery;
• reducing the first turbine diameter increases the power extracted by the second
turbine to the point that it extracts up to 34% more power;
• a variable CP exaggerates the differences between the first turbines but predicts
similar power extraction for the second turbines;
• two variable CP values predict similar power showing that using the free stream
velocity to define the CP is reasonable for both turbines even though there is less
power for the second to extract;
• the spacing between the turbines has little affect on the overall power extracted,
with the maximum difference of 3% found for all but the lowest wind speeds.
7.2 Recommended work
This thesis has shown that the rotor diameter of the first turbine has a significant affect
on the performance of a downstream turbine. This is the first study into the influ-
ence of rotor diameter size and was limited to two axially aligned turbines, which was
investigated as it represents the worst case scenario.
There are multiple avenues this research could take, four possible investigations
are detailed below.
• The influence rotor diameter has on multiple turbines. It was found that reducing
the rotor diameter increased the velocity throughout the domain. This suggests
that more than one turbine would be effected by varying the rotor diameter.
• The influence of the rotor diameter within a turbine array. It was found that
the second turbine increased wake recovery. It may be possible to use different
sized turbines to improve power extraction and wake recovery by using smaller,
cheaper turbine to improve the velocity for larger turbines within an array.
• The interactions of offset turbines could be investigated. The different rotor di-
ameters produced different sized wakes, which would influence the interactions
of offset or misaligned turbines. Using different sized turbines allows the tur-
bines to be offset both vertically and horizontally, which has not been conducted
before and would likely have an affect on the turbine interactions.
• The work presented conducted CFD simulation which solved the RANS equa-
tions using uniform inlets with no ground or shear affects. Investigations could
be conducted to determine these influences on the interaction, specifically con-
sidering different hub heights.
References
[1] Johnson BMC. Analytical and numerical study of a constantly loaded actuator
disk using CFX. In: 8th PhD Seminar on Wind Energy in Europe. Zurich,
Switzerland; 2012. p. 5.
[2] Johnson BMC. Implementation of Actuator Methods in ANSYS-CFX. In: 9th

PhD Seminar on Wind Energy in Europe. Gotland, Sweden; 2013. p. 5.
[3] Johnson B, Whitty J, Howe J, Francis J. Computational Actuator disk models

for wind and tidal applications. Journal of Renewable Energy. 2014;p. In press.
[4] Whitty J, Haydock T, Johnson B, Howe J. On the Deflexion of Anisotropic

Structural Composite Aerodynamic Components. Journal of Wind Energy.
2014;(Article ID 987414):13.
[5] DIRECTIVE 2009/28/EC OF THE EUROPEAN PARLIAMENT AND OF

THE COUNCIL of 23 April 2009, Official Journal of the European Union;
2009.
[6] Sawyer S, Rave K. Global wind report: Annual market update 2012. Fried L,
Sawyer S, Shukla S, Qiao L, editors. GWEC; 2012.
[7] Johnston A, editor. Tracking Clean Energy Progress 2013: IEA input to the
Clean Energy Ministerial. IEA publications; 2013.
204
REFERENCES 205
[8] Chowdhury S, Zhang J, Messac A, Castillo L. Unrestricted wind farm lay-

out optimization (UWFLO): Investigating key factors influencing the maximum
power generation. Renewable Energy. 2012;38(1):16–30.
[9] Katic I, Højstrup J, Jensen NO. A simple model for cluster efficiency. In:
EWEC. Italy; 1986. p. 407–410.
[10] Burton T, Sharpe D, Jenkins N, Bossanyi E. Wind Energy HandBook. John

Wiley and Sons; 2001.
[11] Manwell JF, McGowan JG, Rogers AL. Wind Energy Explained: Theory, De-
sign and Application. John Wiley and Sons; 2009.
[12] Sørensen T, Nielsen P, Thøgersen ML. Recalibrating wind turbine wake model
parameters validating the wake model performance for large offshore wind
farms. In: EWEC. Athens, Greece; 2006. p. 6.
[13] Barthelmie RJ, Frandsen ST, Hansen K, Schepers JG, Rados K, Schlez W, et al.
Modelling the impact of wakes on power output at Nysted and Horns Rev. In:
EWEC. Marseille, France; 2009. p. 10.
[14] Barthelmie RJ, Frandsen ST, Rathmann O, Hansen K, Politis ES, Prospathopou-
los J, et al. Flow and wakes in large wind farms in complex terrain and offshore.
In: EWEC. Brussels, Belgium; 2008. p. 10.
[15] Barthelmie RJ, Frandsen ST, Nielsen NM, Pryor SC, Réthoré PE, Jørgensen
HE. Modelling and measurements of power losses and turbulence intensity
in wind turbine wakes at middelgrunden offshore wind farm. Wind Energy.
2007;10:217–228.
[16] Réthoré PE. State of the Art in Wind Farm Layout Optimization; 2010.
REFERENCES 206
[17] Vermeer LJ, Sørensen JN, Crespo A. Wind turbine wake aerodynamics.
Progress in Aerospace Sciences. 2003;39:467–510.
[18] White FM. Fluid Mechanics (Sixth Edition). McGraw-Hill; 2009.
[19] Newton I. The Principia: Mathematical Principles of Natural Philosophy. Snow-

ballpublishing; 2011.
[20] Clausius R. Über die bewegende Kraft der Wärme. Poggendorff’s Annalen der
Physick. 1850;79:368–397,500–524.
[21] Clausius R. On the Moving Force of Heat, and the Laws regarding the Nature of
Heat itself which are deducible therefrom. Philosophical Magazine and Journal
of Science. 1851;2:1–21, 102–119.
[22] Snel H. Review of Aerodynamics for Wind Turbines. Wind Energy.

2003;6:203–211.
[23] Bechmann A. Large-Eddy Simulation of Atmospheric Flow over Complex Ter-

rain. Technical University of Denmark; 2006.
[24] Wilcox DC. Turbulence Modeling in CFD. DCW; 2006.
[25] Stull RB. An Introduction to Boundary Layer Meteorology. Kluwer Academic

Publishers; 1988.
[26] Reynolds O. An Experimental Investigation of the Circumstances Which Deter-

mine Whether the Motion of Water Shall Be Direct or Sinuous, and of the Law
of Resistance in Parallel Channels. Philosophical Transactions of the Royal
Society of London. 1883;174:935–982.
REFERENCES 207
[27] Reynolds O. On the Dynamical Theory of Incompressible Viscous Fluids and

the Determination of the Criterion. Philosophical Transactions of the Royal
Society of London (A). 1895;186:123–164.
[28] Jones WP, Launder BE. The prediction of laminarization with a two-equation
model of turbulence. International Journal of Heat and Mass Transfer.
1972;15(2):301–314.
[29] Menter FR. Zonal two-equation k-ω models for aerodynamic flows. AIAA.
1993;(93-2906).
[30] Boussinesq MJ. Théorie de l’écoulement tourbillonnant et tumultueux des liq-

uides. Gauthier-Villars et fils; 1897.
[31] Réthoré PE. Wind turbine wake in atmospheric turbulence. Aalborg University;
2009.
[32] Smagorinsky J. General Circulation Experiments with the Primitive Equations.

Monthly Weather Review. 1963;91:99–164.
[33] Lilly DK. On the application of the eddy viscosity concept in the Inertial sub-
range of turbulence. NCAR; 1966. Manuscript 123.
[34] Lilly DK. The representation of small-scale turbulence in numerical simulation

experiments. NCAR; 1966. Manuscript 281.
[35] Blazek J. Computational Fluid Dynamics: Principles and Applications. Else-

vier; 2001.
[36] Piomelli U. Large-Eddy and Direct Simulation of Turbulent Flows; 2014.

REFERENCES 208
[37] Schmitt FG. About Boussinesq’s turbulent viscosity hypothesis: historical re-
marks and a direct evaluation of its validity. Comptes Rendus Mécanique.
2007;335:617–627.
[38] Prandtl L. Über die ausgebildete Turbulenz. ZAMM. 1925;5:136–139.
[39] Prandtl L. Über ein neues Formelsystem für die ausgebildete Turbulenz. Nachr
Akad Wiss Gttingen, Math Phys. 1945;p. 6–19.
[40] Baldwin BS, Barth TJ. A One-Equation Turbulence Transport Model for High
Reynolds Number Wall-Bounded Flows. NASA; 1990. TM 102847.
[41] Spalart PR, Allmaras SR. A One-Equation Turbulence Model for Aerodynamic
Flows. AIAA. 1992;(92-0439).
[42] Wilcox DC. Comparison of two-equation turbulence models for boundary layers
with pressure gradient. AIAA Journal. 1993;31:1414–1421.
[43] Bredberg J. On Two-equation Eddy-Viscosity Models. Chalmers University of

Technology; 2001.
[44] Kolmogorov AN. Equations of Turbulent Motion of an Incompressible Fluid.

Izv Akad Nauk SSSR. 1942;6:56–58.
[45] Launder BE, Sharma BI. Application of the Energy Dissipation Model of Tur-
bulence to the Calculation of Flow Near a Spinning Disc. Letters in Heat and
Mass Transfer. 1974;1(2):131–138.
[46] Reynolds WC. Computation of tubulent flows. Annual Review Fluid Mechan-
ics. 1976;8:183–208.
[47] Wilcox DC. Reassessment of the Scale Determining Equations for Advanced
Turbulence Models. AIAA Journal. 1988;26(11):1311–1320.
REFERENCES 209
[48] Menter FR. Turbulence Modeling for Engineering Flows. ANSYS inc.; 2011.
[49] Menter FR. Two-equation eddy-viscosity turbulence models for engineering

applications. AIAA Journal. 1994;32(8):1598–1605.
[50] Menter FR, Kuntz M, Langtry R. Ten Years of Industrial Experience with the
SST Turbulence Model. In: Turbulence, Heat and Mass Transfer. vol. 4; 2003.
p. 625–632.
[51] Johnson DA, King LS. A new turbulence closure model for boundary layer
flows with strong adverse pressure gradients and seperation. AIAA. 1984;(1984-
0175).
[52] Kalitzin G, Medic G, Iaccarino G, Durbin P. Near-wall behavior of {RANS} tur-

bulence models and implications for wall functions. Journal of Computational
Physics. 2005;204(1):265 – 291.
[53] Rotta JC. Statistische Theorie nichthomogener Turbulenz. Zeitschrift fur

Physik. 1951;129:547–572.
[54] ANSYS FLUENT Solver User Guide. ANSYS Inc.; 2012.
[55] Courant R, Friedrichs KO, Lewy H. Uber die partiellen Diflerenzengleichun-

gen. der mathematischen Physik Math. 1928;100:32–74. Transl.: On the Partial
Difference Equations of Mathematical Physics. IBM Journal, 11 (1967), pp.
215-234.
[56] Zou GW, Liu SL, Chow WK, Gao Y. Large eddy simulation ofTurbulent flows.
International Journal of Architectural Science. 2006;7(1):26–34.
[57] Ferziger JH. Direct and Large Eddy Simulation of Turbulence. In: Vincent A,
editor. Numerical Methods in Fluid Mechanics. vol. 16. CRM; 1998. p. 53–73.
REFERENCES 210
[58] Iliescu T. Large Eddy Simulation for Turbulent Flows. University of Pittsburgh;
2000.
[59] Deck S. Zonal detached-eddy simulation of the flow around a high-lift configu-
ration with deployed slat and flap. AIAA. 2005;43:2372–2384.
[60] Spalart PR, Jou WH, Strelets M, Allmaras SR. Comments on the feasibility of
LES for wing and on a hybrid RANS/LES approach. In: 1st ASOSR CONFER-
ENCE on DNS/LES; 1997. p. 137–147.
[61] Spalart PR. Detached-Eddy Simulation. Annual Review of Fluid Mechanics.

2009;41(1):181–202.
[62] Moin P, Mahesh K. DIRECT NUMERICAL SIMULATION: A Tool in Turbu-

lence Research. Annual Review of Fluid Mechanics. 1998;30(1):539–578.
[63] Turner MJ, Clough RWaHC, J TL. Stiffness and Deflection Analysis of Com-
plex Structures. J of Aero Sci. 1956;23(9):805–823.
[64] McDonald PW. The Cornputatdon of Transonic Flow through Two-Dimensional

Gas Turbine Cascades. ASME; 1971. 71-GT-89.
[65] ANSYS CFX Solver Theory Guide. ANSYS Inc.; 2010.
[66] ANSYS FLUENT Solver Theory Guide. ANSYS Inc.; 2012.
[67] Sanderse B, van der Pijl SP, Koren B. Review of CFD for wind-turbine wake
aerodynamics. Wind Energy. 2011;14:799–819.
[68] Wu Y, Porté-Agel F. Large-eddy simulation of wind-turbine wakes: Evaluation

of turbine parameterizations. Boundary Layer Meteorology. 2011;138.
[69] Aubrun S, Loyer S, Hancock PE, Hayden P. Wind turbine wake proper-
ties: Comparison between a non-rotating simplified wind turbine model and
REFERENCES 211
a rotating model. Journal of Wind Engineering and Industrial Aerodynamics.

2013;120:1–8.
[70] Calaf M, Meneveau C, Meyers J. Large eddy simulation study of fully devel-
oped wind-turbine array boundary layers. Physics of Fluids. 2010;22:22–38.
[71] Port-Agel F, Wu YT, Chen CH. A Numerical Study of the Effects of Wind
Direction on Turbine Wakes and Power Losses in a Large Wind Farm. Energies.
2013;6(10):5297–5313.
[72] Castellani F, Vignaroli A. An application of the actuator disc model for wind
turbine wakes calculations. Applied Energy. 2013;101:432–440.
[73] Sørensen JN, Kock CW. A model for unsteady rotor aerodynamics. Journal of
Wind Engineering and Industrial Aerodynamics. 1995;58(3):259–275.
[74] Troldborg N. Actuator line modeling of wind turbine wakes. Technical Univer-
sity of Denmark; 2008.
[75] Mikkelsen R. Actuator disc methods applied to wind turbines. Technical Uni-
versity of Denmark; 2002.
[76] Porté-Agel F, Lu H, Wu Y. A large-eddy simulation framework for wind energy

applications. In: The Fifth International Symposium on Computational Wind
Engineering; 2010. p. 21.
[77] Porté-Agel F, Wu Y, Lu H, Conzemius RJ. Large-eddy simulation of atmo-

spheric boundary layer flow through wind turbines and wind farms. Journal of
Wind Engineering and Industrial Aerodynamics. 2011;99(4):154–168.
[78] Shen WZ, Mikkelsen R, Sørensen J, Bak C. Tip loss corrections for wind turbine
computations. Wind Energy. 2005;8(4):457–475.
REFERENCES 212
[79] Branlard E. Wind turbine tip-loss corrections. Technical University of Denmark;

2011.
[80] Chicacausa AM. Development of a model for computing tip loss corrections in
wind turbines. University of Zaragoza; 2013.
[81] Prandtl L. Applications of modern hydrodynamics to aeronautics. NACA; 1921.

report No. 116.
[82] Xu G, Sankar L. Computational studies of horizontal axis wind turbines. In:

ASME Wind Energy Symposium; 1999. p. 6.
[83] Lindenburg C. Investigation into rotor blade aerodynamics. ECN; 2003. ECN-
C-03-025.
[84] Sant T. mproving BEM-based aerodynamics models in Wind turbine design

codes. Delft University Wind Energy Research Institute; 2007.
[85] Hand MM, Simms DA, Fingersh LJ, Jager DW, Cotrell JR, Schreck S, et al. Un-
steady Aerodynamics Experiment Phase VI: Wind Tunnel Test Configurations
and Available Data Campaigns. NREL/TP-500-29955; 2001.
[86] Bracewell R. The Fourier Transform and Its Applications. Third edition ed.
New York: McGraw-Hill; 2000.
[87] Keck RE. A numerical investigation of nacelle anemometry for a HAWT using
actuator disc and line models in CFX. Renewable Energy. 2012;48(0):72–84.
[88] ANSYS FLUENT 14.5 UDF Manual. ANSYS Inc.; 2012.
[89] Sørensen JN, Shen WZ. Numerical modeling of wind turbine wakes. Journal of
Fluids Engineering. 2002;124:393–399.
REFERENCES 213
[90] Snel H. Review of the present status of rotor aerodynamics. Wind Energy.
1998;1:46–69.
[91] Sanderse B. Aerodynamics of wind turbine wakes Literature review. ECN Wind
Energy; 2009. ECN-E–09-016.
[92] Sørensen JN. 2.08 - Aerodynamic Analysis of Wind Turbines. In: Sayigh A,
editor. Comprehensive Renewable Energy. Oxford: Elsevier; 2012. p. 225–241.
[93] Rankine WJM. On the Mechanical Principles of the Action of Propellers. Trans-
actions of the Institution of Naval Architects. 1865;6:13–30.
[94] Froude RE. On the part played in propulsion by differences of fluid pressure.
Transactions of the Institute of Naval Architects. 1889;30:390–405.
[95] Hansen MOL. Aerodynamics of Wind Turbines (Second Edition). Earthscan;

2008.
[96] Betz A. Das Maximum der theoretisch möglichen Ausnützung des Windes
durch Windmotoren. Zeitschrift für das gesamte Turbinenwesen. 1920;26:307–
309.
[97] Bergey KH. The Lanchester-Betz limit. Journal of Energy. 1979;3:382–384.
[98] van Kuik GAM. The Lanchester-Betz-Joukowsky limit. Wind Energ.

2007;10:289–291.
[99] Okulov VL, van Kuik GAM. The Betz-Joukowsky limit: on the contribution to
rotor aerodynamics by the British, German and Russian scientific schools. Wind
Energy. 2012;15(2):335–344.
[100] Lanchester FW. A contribution to the theory of propulsion and the screw pro-
peller. Journal of the American Society for Naval Engineers. 1915;27:509–510.
REFERENCES 214
[101] Joukowsky NE. Windmill of the NEJ type. In: Transactions of the Central In-
stitute for Aero-Hydrodynamics of Moscow; 1920. And in: Joukowsky NE.
Collected Papers Vol VI . The Joukowsky Institute for AeroHydrodynamics:
Moscow, Russia, 1937; VI : 405-409 (in Russian).
[102] Crespo A, Hernández J, Frandsen S. Survey of Modeling Methods for Wind

Turbine Wakes and Wind Farms. Wind Energy. 1999;2:1–24.
[103] Butterfield CP, Musial WP, Simms DA. Combined Experiment Phase I Final
Report. NREL/TP-257-4655; 1992.
[104] Simms DA, Hand MM, Fingersh LJ, Jager DW. Unsteady Aerodynamics Exper-
iment Phases II-IV Test Configurations and Available Data Campaigns. NREL;
1999. TP-500-25950.
[105] Simms DA, Schreck S, Hand MM, Fingersh LJ. NREL Unsteady Aerodynamics
Experiment in the NASA-Ames Wind Tunnel: A Comparion of Predictions to
Measurements. NREL; 2001. TP-500-29494.
[106] Schepers JG, Snel H. Model experiments in controlled conditions, final report.
ECN; 2007. ECN-E-07-042.
[107] Stoevesandt B, Herraez Hernandez I, Plischka H, Peinke J. 3D - Rotational

Effects - A Lesson to Learn from the MEXCIO Experiment. In: EWEA; 2012.
p. 6.
[108] Bechmann A, Sørensen NN, Zahle F. CFD simulations of the MEXICO rotor.
Wind Energy. 2011;14(5):677–689.
[109] Sørensen NN, Bechmann A, Réthoré PE, Zahle F. Near wake Reynolds-
averaged Navier-Stokes predictions of the wake behind the MEXICO rotor in
axial and yawed flow conditions. Wind Energy. 2012;17:75–86.
REFERENCES 215
[110] Lutz T, Meister K, Krämer E. Near wake studies of the mexico rotor. In: Euro-
pean Wind Energy Conference (EWEC); 2011. p. 161–165.
[111] Carrión M, Steijl R, Woodgate M, Barakos G, Munduate X, Gomez-Iradi S.

Computational fluid dynamics analysis of the wake behind the MEXICO rotor
in axial flow conditions. Wind Energy. 2014;p. 1–23.
[112] Grasso F, van Garrel A. Near Wake Simulation of Mexico rotor in Axial and
Yawed Flow Conditions with Lifting Line Free Wake Code. In: Wake Confer-
ence; 2011. p. 9.
[113] Micallef D, Kloosterman M, Ferreira CS, Sant T, van Bussel G. Validating

BEM, Direct and Inverse Free Wake Models with the MEXICO Experiment.
48th AIAA Aerospace Sciences Meeting; 2010. AIAA 2010-462.
[114] Breton S, Sibuet C, Masson C. Using the Actuator Surface Method to Model
the Three-Bladed MEXICO Wind Turbine. In: 48th AIAA Aerospace Sciences
meeting; 2010. p. 11.
[115] Shen W, Sørensen J, Yang H. Actuator Line/Navier-Stokes Computations for

Flows past the Yawed MEXICO Rotor. In: Wake Conference; 2011. p. 6.
[116] Shen WZ, Zhu WJ, Sørensen J. Actuator line/Navier-Stokes computations for
the MEXICO rotor: comparison with detailed measurements. Wind Energy.
2012;15(5):811–825.
[117] Schepers JG, Pascal L, Snel H. First results from Mexnext: Analysis of detailed
aerodynamic measurements on a 4.5 m diameter rotor placed in the large Ger-
man Dutch Wind Tunnel DNW. In: European Wind Energy Conference and
Exhibition (EWEC); 2010. p. 9.
REFERENCES 216
[118] Schepers JG, Boorsma K, Cho T, Schaffarczy SGIG, Jeromin A, Shen WZ,
et al. Final report of IEA Task 29, Mexnext (Phase 1), Analysis of Mexico wind
tunnel measurements. Energy Research Center of the Netherlands, ECN; 2011.
ECN-E-12-004.
[119] Schepers JG, Boorsma K, Munduate X. Final Results from Mexnext-I: Analysis
of detailed aerodynamic measurements on a 4.5 m diameter rotor placed in the
large German Dutch Wind Tunnel DNW. In: The Science of Making Torque
From Wind; 2012. p. 24.
[120] Réthoré PE, Zahle F, Sørensen N, Bechmann A. CFD Simulations of the MEX-
ICO Wind Tunnel and Wind Turbine. In: EWEA Annual Event; 2011. p. 8.
[121] Réthoré PE, Zahle F, Sørensen N, Bechmann A. MEXICO wind tunnel and
wind turbine modelled in CFD. In: 41st AIAA Fluid Dynamics Conference and
Exhibit; 2011. p. 10.
[122] Glauert H. Aerodynamic Theory. vol. 4. Durand WF, editor; 1935.
[123] Lissaman PBS. Energy effectiveness of arbitrary arrays of wind turbines. In:
AIAA, Aerospace Sciences Meeting; 1979. p. 7.
[124] Ainslie JF. Development of an eddy viscosity model for wind turbine wakes.
In: 7th EWEA Wind Energy Conference; 1985. p. 61–66.
[125] Ainslie JF. Calculating the flowfield in the wake of wind turbines. Journal of
Wind Engineering and Industrial Aerodynamics. 1988;27:213–224.
[126] Adams BM, Quarton DC. Dynamic loads in wind farms. Hassan G, editor. Joule
Project; 1996.
REFERENCES 217
[127] Crespo A, Manuel F, Moreno D, Fraga E, Hernández J. Numerical analysis of

wind turbine wakes. Proceedings Delphi Workshop On Wind Energy Applica-
tions. 1985;p. 15–25.
[128] Crespo A, Hernández J. Numerical modelling of the flow field in a wind tur-
bine wake. 3rd Joint ASCE/ASME Mechanics Conference Forum on Turbulent
Flows. 1989;p. 121–127.
[129] Crespo A, Manuel F, Hernández J. Numerical modelling of wind turbine wakes.

In: EWEC. Madrid; 1990. p. 166–170.
[130] Crespo A, Chacón L, Hernández J, Manuel F, Grau JC. UPMPARK: a parabolic

3D code to model wind farms. EWEC. 1994;p. 454–459.
[131] Gómez-Elvira R, Crespo A, Migoya E, Manuel F, Hernández J. Anisotropy of

turbulence in wind turbine wakes. Journal of Wind Engineering and Industrial
Aerodynamics. 2005;93:797–814.
[132] Zahle F. Wind Turbine Aerodynamics Using an Incompressible Overset Grid

Method. Imperial College London; 2006.
[133] Michelsen JA. Basis3D - a platform for development of multiblock PDE solvers.
Technical University of Denmark; 1992. AFM 92-05.
[134] Michelsen JA. Block structured multigrid solution of 2D and 3D elliptic PDEs.
Technical University of Denmark; 1994.
[135] Sørensen NN. General purpose flow solver applied to flow over hills. Risoe
National Laboratory; 1995.
[136] Ivanell S. Numerical Computations of Wind Turbine Wakes. KTH Royal Insti-
tute of Technology; 2009.
REFERENCES 218
[137] Laursen J, Enevoldsen P, Hjort S. 3D CFD Quantification of the Performance

of a Multi-Megawatt Wind Turbine. Journal of Physics: Conference Series.
2007;75(1):1–13.
[138] Yelmule MM, EswaraRao Anjuri VSJ. CFD predictions of NREL phase VI rotor
experiments in NASA/AMES wind tunnel. International Journal of Renewable
Energy Research. 2013;3(2):261–269.
[139] Choi NJ, Nam SH, Jeong JH, Kim KC. Numerical study on the horizontal
axis turbines arrangement in a wind farm: Effect of separation distance on the
turbine aerodynamic power output. Journal of Wind Engineering and Industrial
Aerodynamics. 2013;117(0):11–17.
[140] Harrison ME, Batten WMJ, Myers LE, Bahaj AS. A comparison between CFD
simulations and experiments for predicting the far wake of horizontal axis tidal
turbines. In: 8th European Wave and Tidal Energy Conferences. Uppsala, Swe-
den; 2009. p. 10.
[141] Harrison ME, Batten WMJ, Myers LE, Bahaj AS. Comparison between CFD
simulations and experiments for predicting the far wake of horizontal axis tidal
turbines. IET Renewable Power Generation. 2010;4:613–627.
[142] Mahu R, Popescu F. NREL Phase VI rotor modeling and simulation using
ANSYS FLUENT 12.1. Mathematical Modeling In Civil Engineering. 2011;p.
185–194.
[143] Mahu R, Popescu F, Frunzulica F, Dumitrache A. 3D CFD Modeling

and Simulation of NREL Phase VI Rotor. AIP Conference Proceedings.
2011;1389(1):1520–1523.
REFERENCES 219
[144] Qunfeng L, Jin C, Jiangtao C, Ning Q, Danao LAM. Study of CFD simula-
tion of a 3-D wind turbine. In: Materials for Renewable Energy Environment
(ICMREE). vol. 1; 2011. p. 596–600.
[145] WuBow S, Sitzki L, Hahm T. 3D simulation of the turbulent wake behind a

wind turbine. Journal of Physics Conference Series 75; 2007. 012033.
[146] Makridis A, Chick J. CFD Modeling of the wake interactions of two wind
turbines on a Gaussian hill. In: EACWE 5; 2009. p. 9.
[147] Makridis A, Chick J. Validation of a {CFD} model of wind turbine wakes

with terrain effects. Journal of Wind Engineering and Industrial Aerodynamics.
2013;123, Part A(0):12–29.
[148] Stovall T, Pawlas G, Moriarty P. Wind farm wake simulations in OpenFOAM.

48th AIAA Aerospace Sciences Meeting; 2010. AIAA 2010-825.
[149] Réthoré PE, Sørensen NN, Bechmann A, Zahle F. Study of the atmospheric
wake turbulence of a CFD actuator disc model. In: EWEC. Marseille, France;
2009. p. 9.
[150] El Kasmi A, Masson C. An extended k-ε model for turbulent flow through
horizontal-axis wind turbines. Journal of Wind Engineering and Industrial Aero-
dynamics. 2008;96:103–122.
[151] Cabezón D, Migoya E, Crespo A. Comparison of turbulence models for the

computational fluid dynamics simulation of wind turbine wakes in the atmo-
spheric boundary layer. Wind Energy. 2011;14(7):909–921.
[152] Rados KG, Prospathopoulos JM, Stefanatos NC, Politis ES, Chaviaropoulos PK,
Zervos A. CFD modeling issues of wind turbine wakes under stable atmospheric
conditions. In: EWEC. Marseille, France; 2009. p. 8.
REFERENCES 220
[153] Johansen J, Sørensen NN, Michelsen JA, Schreck S. Detached-eddy simulation

of flow around the NREL phase VI blade. Wind Energy. 2002;5:185–197.
[154] Troldborg N, Sørensen JN, Mikkelsen R. Actuator Line Simulation of Wake of

Wind Turbine Operating in Turbulent Inflow. The Science of Making Torque
from Wind Journal of Physics: Conference Series. 2007;75.
[155] Zahle F, Sørensen NN. Overset grid flow simulation on a modern wind turbine.
AIAA; 2008. Paper 20086727.
[156] Zahle F, Sørensen NN. In: Upwind Navier-Stokes Rotor Flow Simulations
with ground proximity and shear. European Wind Energy Association (EWEA);
2010. p. 255–263.
[157] Mann J. The spatial structure of neutral atmospheric surface-layer turbulence.

Journal of Fluid Mechanics. 1994;272:141–168.
[158] Mann J. Wind field simulation. Probabilistic Engineering Mechanics 13.

1998;13(4):269–282.
[159] Mikkelsen R, Sørensen JN, Troldborg N. Prescribed wind shear modelling with
the actuator line technique. In: EWEC. Milan, Italy; 2007. p. 5.
[160] Meyers J, Meneveau C. Large eddy simulations of large wind-turbine arrays

in the atmospheric boundary layer. 48th AIAA Aerospace Sciences Meeting;
2010. AIAA 2010-827.
[161] Sørensen JN, Shen WZ, Munduate X. Analysis of wake states by a full-field
actuator disc model. Wind Energy. 1998;1:73–88.
[162] Ivanell S, Sørensen JN, Mikkelsen R, Henningson D. Analysis of numerically

generated wake structures. Wind Energy. 2009;12:63–80.
REFERENCES 221
[163] Whale J, Anderson CG, Bareiss R, Wagner S. An experimental and numerical

study of the vortex structure in the wake of a wind turbine. Journal of Wind
Engineering and Industrial Aerodynamics. 2000;84:1–21.
[164] Prospathopoulos JM, Politis ES, Chaviaropoulos PK. Modelling wind turbine
wakes in complex terrain. In: EWEC. Brussels, Belgium; 2008. .
[165] Antoniadis AF, Drikakis D, Zhong B, Barakos G, Steijl R, Biava M, et al. As-
sessment of CFD methods against experimental flow measurements for heli-
copter flows. Aerospace Science and Technology. 2012;19(1):86–100.
[166] Wald QR. The aerodynamics of propellers. Progress in Aerospace Sciences.

2006;42(2):85–128.
[167] Réthoré PE, Zahle F, Sørensen NN. Comparison of an actuator disc model with
a full Rotor CFD model under uniform and shear inflow condition. In: 4th PhD
seminar on wind energy in Europe; 2008. p. 123–126.
[168] Réthoré PE, Zahle F, Sørensen NN. Validation of an actuator disc model. In:
EWEC. Warsaw, Poland; 2010. p. 10.
[169] Conway JT. Analytical solution for the actuator disk with variable radial distri-
bution of load. Journal of Fluid Mechanics. 1995;297:327–355.
[170] Réthoré PE, Sørensen NN. Actuator disc model using a modified Rhie-
Chow/SIMPLE pressure correction algorithm. Comparison with analytical so-
lutions. In: EWEC. Brussels, Belgium; 2008. p. 11.
[171] Sørensen JN, Mikkelsen R, Troldborg N. Simulation and Modelling of Turbu-

lence in WInd Farms. In: EWEC. Milan, Italy; 2007. p. 10.
REFERENCES 222
[172] Ivanell S, Mikkelsen R, Sørensen JN, Henningson D. Three dimensional ac-

tuator disc modelling of wind farm wake interaction. In: EWEC. Brussels,
Belgium; 2008. p. 20.
[173] Alinot C, Masson C. Aerodynamic simulations of wind turbines operating in

atmospheric boundary layer with various thermal stratifications. ASME 2002
Wind Energy Symposium. 2002;p. 206–215.
[174] Mikkelsen R, Sørensen JN, Øye S, Troldborg N. Analysis of Power Enhance-

ment for a Row of Wind Turbines Using the Actuator Line Technique. Journal
of Physics: Conference Series 75. 2007;(012044).
[175] Réthoré PE, Troldborg N, Zahle F, Sørensen NN. Comparison of the Near Wake
of Different Kinds of Wind Turbine CFD Models. In: Wake Conference. Got-
land; 2011. p. 7.
[176] Shen WZ, Sørensen JN, Zhang JH. Actuator surface model for wind turbine
flow computations. In: EWEC. Milan, Italy; 2007. p. 8.
[177] Shen WZ, Zhang JH, Sørensen JN. The Actuator Surface Model: A New
Navier-Stokes Based Model for Rotor Computations. Journal of Solar Energy
Engineering. 2009;131(1):9.
[178] Dobrev I, Massouh F, Maalouf B. Lifting surface method for prediction of rotor
vortex wake. 19ème Congrès Français de Mécanique, Marseille. 2009;p. 24–28.
[179] Sibuet Watters C, Breton SP, Masson C. Application of the actuator surface
concept to wind turbine rotor aerodynamics. Wind Energy. 2010;13:433–447.
[180] Dobrev I, Massouh F, Rapin M. Actuator surface hybrid model. Journal of

Physics: Conference Series. 2007;75(1):1–7.
REFERENCES 223
[181] Sørensen NN, Hansen MOL. Rotor performance predictions using a Navier-
Stokes method. AIAA; 1998. Paper 98-0025.
[182] Zahle F, Sørensen NN, Johansen J. Wind turbine rotor-tower interaction using
an incompressible overset grid method. Wind Energy. 2009;12(6):594–619.
[183] Miller W, Sitaraman J, Mavriplis D. Large Scale Simulations of Wind Farm

Aerodynamics in Turbulent Atmospheric Inflow Conditions; 2013.
[184] Bahaj AS, Myers LE, Thompson G. Characterising the wake of horizontal axis
marine current turbines. In: 7th European Wave and Tidal Energy Conference.
Porto, Portugal; 2007. p. 9.
[185] Myers LE, Bahaj AS. Experimental analysis of the flow field around horizonta
laxis tidal turbines by use of scale mesh disk rotor simulators. Ocean Engineer-
ing. 2010;37:218–227.
[186] ANSYS Modeling and Meshing Guide. ANSYS Inc.; 2005.
[187] Beyer HG, Pahlke T, Schmidt W, Waldl HP, de Witt U. Wake effects in a
linear wind farm. Journal of Wind Engineering and Industrial Aerodynamics.
1994;51(3):303–318.
[188] Dahlberg J, Thor S. Power performance and wake effects in the closely spaced
Lillgrund offshore wind farm. In: European Offshore Wind Conference; 2009.
p. 10.
[189] Wu YT, Porté-Agel F. Simulation of Turbulent Flow Inside and Above Wind
Farms: Model Validation and Layout Effects. Boundary-Layer Meteorology.
2013;146:181–205.
[190] Chamorro L, Arndt R, Sotiropoulos F. Turbulent Flow Properties Around a

Staggered Wind Farm. Boundary-Layer Meteorology. 2011;141(3):349–367.
REFERENCES 224
[191] Chamorro LP, Porté-Agel F. Turbulent Flow Inside and Above a Wind Farm: A
Wind-Tunnel Study. Energies. 2011;4(11):1916–1936.
[192] Nilsson K. Numerical computations of wind turbine wakes and wake interaction
: Optimization and control; 2012. QC 20130111.
[193] Barthelmie RJ, Hansen K, Frandsen ST, Rathmann O, Schepers JG, Schlez W,
et al. Modelling and measuring flow and wind turbine wakes in large wind farms
offshore. Wind Energy. 2009;12(5):431–444.
[194] Yang X, Kang S, Sotiropoulos F. Computational study and modeling of

turbine spacing effects in infinite aligned wind farms. Physics of Fluids.
2012;24(11):115107 – 115107–28.
[195] Troldborg N, Larsen GC, Madsen HA, Hansen KS, Sørensen JN, Mikkelsen R.
Numerical simulations of wake interaction between two wind turbines at various
inflow conditions. Wind Energy. 2011;14(7):859–876.
[196] Adaramola M, Krogstad PA. Wind tunnel simulation of wake effects on wind
turbine performance. In: EWEC; 2010. p. 64–67.
[197] Bartl J, Pierella F, Sætrana L. Wake Measurements Behind an Array of Two

Model Wind Turbines. Energy Procedia. 2012;24:305–312.
[198] Fletcher TM, Brown RE. Simulation of wind turbine wake interaction using the
vorticity transport model. Wind Energy. 2010;13(7):587–602.
[199] Ozbay A. Experimental investigations on the wake interferences of multiple

wind turbines. Iowa State University; 2012.
[200] Turnock SR, Phillips AB, Banks J, Nicholls-Lee R. Modelling tidal current
turbine wakes using a coupled RANS-BEMT approach as a tool for analysing
REFERENCES 225
power capture of arrays of turbines. Ocean Engineering. 2011;38(11):1300–

1307.
[201] O’Doherty DM, Mason-Jones A, O’Doherty T, Byrne CBB. Considerations of

improved tidal stream turbine performance using double rows of contra-rotating
blades. In: 8th European Wave and Tidal Energy (EWTEC); 2009. p. 9.
[202] Steinbuch M, de Boer WW, Bosgra OH, Peters SAWM, Ploeg J. Optimal control
of wind power plants. Journal of Wind Engineering and Industrial Aerodynam-
ics. 1988;27(1-3):237–246.
[203] Johnson KE, Thomas N. Wind farm control: Addressing the aerodynamic in-
teraction among wind turbines. In: American Control Conference; 2009. p.
2104–2109.
[204] Frandsen S, Barthelmie R, Pryor S, Rathmann O, Larsen S, Højstrup J, et al.

Analytical modelling of wind speed deficit in large offshore wind farms. Wind
Energy. 2006;9(1-2):39–53.
[205] Wind Turbine Power Calculations. The Royal Academy of Engineering;.
[206] Buckingham E. On Physically Similar Systems: Illustrations of the Use of

Dimensional Equations. Physics Review. 1914;4(4):345–376.
Appendix A
Reynolds number
This appendix details two methods which can be used to proof that the Reynolds num-
ber, equation (A.2), is a ratio of the inertia force to the viscous force.
Inertia f orce
Re = (A.1)
Viscous f orce
ρul
= (A.2)
µ
Equations (A.3) and (A.4) define the inertia and viscous force respectively.
Inertia f orce = ma = ρAla (A.3)
Viscous f orce = µAγ̇ (A.4)
Where m is the mass, a is the acceleration, A is the area, l is a length, ρ is the

density, u is the velocity, µ is the dynamic viscosity and γ̇ is the shear rate.
Combining these definitions produces equation (A.5).
226
APPENDIX A. REYNOLDS NUMBER 227
ρAl a
= (A.5)
µA γ̇
ρl a
= (A.6)
µ γ̇
Comparing equations (A.5) and (A.2) shows that if the Reynolds number is a ratio
a
of the inertia force to the viscous force then γ̇ = u. Equation (A.7) shows that this is
consistent using dimensional analysis.
a [ms−2 ]
= −1 = [ms−1 ] = u (A.7)
γ̇ [s ]
Proof 1:
Considering strain
Dε
a = ×u (A.8)
Dt
Dε
γ̇ = (A.9)
Dt
Combining theses variables as shown in equation (A.10), then expanding and simpli-
fying shows that:
a
= aγ̇−1 (A.10)
γ̇
−1
Dε Dε
= ×u× (A.11)
Dt Dt
= u (A.12)
D
where Dt is the substantial derivative, a is the acceleration, γ̇ is the shear rate, ε is
APPENDIX A. REYNOLDS NUMBER 228
the strain, u is the velocity and t is the time.
Proof 2:
This proof is provided online on the NASA website1 .
du
Inertia f orce = ρu (A.13)
dx
d2u
Viscous f orce = µ 2 (A.14)
dx
Inertia f orce
Re = (A.15)
Viscous f orce
ρu du
dx
= 2 (A.16)
µ ddxu2
The velocity gradient is proportional to the velocity divided by a length scale l. Sim-
ilarly, the second derivative is proportional to the velocity divided by the length scale
squared.
ρu du
dx ρu ul ρuul 2
= = (A.17)
2
µ ddxu2 µ lu2 µul
ρul
= (A.18)
µ
Where ρ is the density, u is the speed, µ is the dynamic viscosity, γ is the shear
µ
stress and ν = ρ is the kinematic viscosity.
1 https://www.grc.nasa.gov/www/k-12/airplane/reynolds.html

Johnson Ben Final E-Thesis (Master Copy)

Uploaded by

Copyright:

Available Formats

Johnson Ben Final E-Thesis (Master Copy)

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Johnson Ben Final E-Thesis (Master Copy)

Uploaded by

Copyright:

Available Formats

Computational Fluid Dynamics (CFD) Modelling

of Renewable Energy Turbine Wake Interactions

Benjamin Michael Carver Johnson

Concurrent registration for two or more academic awards

Material submitted for another award

Signature of Candidate _______________________________________________________

Type of Award PhD

School School of CEPS

7 Conclusions & Recommendations 196

A Reynolds number 226

4.1 Mesh densities using linear and quadratic cells . . . . . . . . . . . . 111

5.1 Thrust coefficient results to 3 s.f. . . . . . . . . . . . . . . . . . . . 130

6.1 New CT values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174

1.1 Types of wind turbine . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.1 Difference in grid resolution using a wall function . . . . . . . . . . 57

3.1 1D momentum theory . . . . . . . . . . . . . . . . . . . . . . . . . . 82

4.1 1D momentum theory . . . . . . . . . . . . . . . . . . . . . . . . . 104

This thesis presents Computational Fluid Dynamics (CFD) simulations of renewable

a Axial induction factor [−]

a0 Angular or tangential induction factor [−]

a1 Constant (fraction of the rate of turbulence production) [−]

A Boundary layer model constant [−]

Ar Rotor swept area [m2 ]

b Distance between a grid point and source point [m]

B Number of blades [−]

c Chord length [m]

cp Specific heat capacity [Nm/K]

C Closure coefficient [−]

C1,2 Prescribed matrix [−]

CD Drag coefficient [−]

Ci jk Turbulent transport tensor [−]

CL Lift coefficient [−]

Cp Power coefficient [−]

CS Smagorinsky constant [−]

CT Thrust coefficient [−]

CDkw Cross diffusion term [−]

D Drag force [N]

DAD Actuator Disc [N/m3 ]

DADUT Actuator Disc with Uniform Thrust definition term [N/m2 ]

DADV T Actuator Disc with Variable Thrust definition term [N/m]

DADL Actuator Disc with Lift definition term [N/m]

DAL Actuator Line definition term [N/m]

Dn Disc definition term [N/m]

d Rotor thickness/depth [m]

eL Lift unit vector [−]

eD Drag unit vector [−]

f Base term used in the definition of the actuator methods [N/m2 ]

F1 ,F2 Blending function [−]

FP Prandtl’s correction factor [−]

FG Glauert’s correction factor [−]

FS Shen et al. correction factor [−]

G Filter function [−]

H(x) Step function [−]

i, j Represents components in the i, jth direction [−]

k Specific turbulence kinetic energy [m2 /s2 ]

k Thermal diffusivity [m2 /s]

K Resistance coefficient [−]

l Characteristic length [m]

L Lift force [N]

Ln Actuator line area [m2 ]