Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
import numpy as np
from py_wake.deficit_models.deficit_model import DeficitModel
from py_wake.wind_farm_models.engineering_models import PropagateDownwind
from py_wake.superposition_models import LinearSum
from py_wake.rotor_avg_models.rotor_avg_model import RotorCenter
na = np.newaxis
def my_power(term, factor):
with np.warnings.catch_warnings():
# if term is 0, exp(log(0))=0 as expected for a positive factor
np.warnings.filterwarnings('ignore', r'divide by zero encountered in log')
return np.exp(factor * np.log(term))
def get_r96(D, CT, TI):
"""Computes the wake radius at 9.6D downstream location of a turbine from empirical relation
.. math::
R_{9.6D} = a_1 \\exp (a_2 C_T^2 + a_3 C_T + a_4) (b_1 TI + b_2) D
Inputs
----------
D: float
Wind turbine diameter
CT: float
Outputs WindTurbine object's thrust coefficient
TI: float
Ambient turbulence intensity
pars: list
GCL Model parameters [a1, a2, a3, a4, b1, b2]
Returns
-------
R96: float
Wake radius at 9.6D downstream location
"""
a1, a2, a3, a4, b1, b2 = [0.435449861, 0.797853685, -0.124807893, 0.136821858, 15.6298, 1.0]
R96 = a1 * (np.exp(a2 * CT * CT + a3 * CT + a4)) * (b1 * TI + b2) * D
return R96
def get_Rw(x, R, TI, CT):
"""Computes the wake radius at a location.
[1]-eq.3
.. math::
R_w = \\left(\\frac{105 c_1^2 }{2 \\pi}\\right)^{0.2} (C_T A (x + x_0))^{1/3}
with A, the area, and x_0 and c_1 defined as
.. math::
x_0 = \\frac{9.6 D}{\\left(\\frac{2 R_96}{k D} \\right)^3 - 1}
c_1 = \\left(\\frac{k D}{2}\\right)^{5/2}
\\left(\\frac{105}{2 \\pi} \\right)^{-1/2}
(C_T A x_0)^{-5/6}
with k and m defined as
.. math::
k = \\sqrt{\\frac{m + 1}{2}}
m = \\frac{1}{\\sqrt{1 - C_T}}
Inputs
----------
x: float or ndarray
Distance between turbines and wake location in the wind direction
R: float
Wind turbine radius
TI: float
Ambient turbulence intensity
CT: float
Outputs WindTurbine object's thrust coefficient
Returns
-------
Rw: float or ndarray
Wake radius at a location
"""
D = 2.0 * R
Area = np.pi * R * R
m = 1.0 / (np.sqrt(1.0 - CT))
k = np.sqrt((m + 1.0) / 2.0)
R96 = get_r96(D, CT, TI)
x0 = (9.6 * D) / (my_power(2.0 * R96 / (k * D), 3.0) - 1.0)
xx0 = x + x0
term1 = my_power(k * D / 2.0, 2.5)
term2 = my_power(105.0 / (2.0 * np.pi), -0.5)
term3 = my_power(CT * Area * x0, -5.0 / 6.0)
c1 = term1 * term2 * term3
Rw = my_power(105.0 * c1 * c1 / (2.0 * np.pi), 0.2) * my_power(CT * Area * xx0, 1.0 / 3.0)
Rw = np.where(x + x0 <= 0., 0., Rw)
return Rw, xx0, c1
def get_dU(x, r, R, CT, TI):
"""Computes the wake velocity deficit at a location
Inputs
----------
x: float
Distance between turbines and wake location in the wind direction
r: float
Radial distance between the turbine and the location
R: float
Wake producing turbine's radius [m]
CT: float
Outputs WindTurbine object's thrust coefficient
TI: float
Ambient turbulence intensity [-]
order: int, optional
Returns
-------
dU: float
Wake velocity deficit at a location
"""
CT = np.maximum(CT, np.finfo(np.float).eps)
Area = np.pi * R * R
Rw, xx0, c1 = get_Rw(x, R, TI, CT)
c1s = c1 * c1
term10 = (1 / 9)
term20 = my_power(CT * Area / (xx0 * xx0), 1. / 3.)
term310 = my_power(r, 1.5)
term320 = 1.0 / np.sqrt(3. * c1s * CT * Area * xx0)
term30 = term310 * term320
term41 = my_power(35. / (2. * np.pi), .3)
term42 = my_power(3. * c1s, -0.2)
term40 = term41 * term42
t4 = term30 - term40
dU1 = -term10 * term20 * t4 * t4
dU = dU1
dU = np.where((Rw < r) | (x <= 0), 0, dU)
return dU
class GCLDeficitModel(DeficitModel):
args4deficit = ['WS_ilk', 'D_src_il', 'dw_ijlk', 'cw_ijlk', 'ct_ilk', 'TI_ilk']
def wake_radius(self, dw_ijlk, D_src_il, TI_ilk, ct_ilk, **_):
return get_Rw(x=dw_ijlk, R=(D_src_il / 2)[:, na, :, na], TI=TI_ilk[:, na], CT=ct_ilk[:, na])[0]
def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, TI_ilk, **_):
eps = 1e-10
dw_ijlk_gt0 = np.maximum(dw_ijlk, eps)
R_src_il = D_src_il / 2.
dU = -get_dU(x=dw_ijlk_gt0, r=cw_ijlk, R=R_src_il[:, na, :, na],
CT=ct_ilk[:, na], TI=TI_ilk[:, na])
return WS_ilk[:, na] * dU * (dw_ijlk > eps)
class GCL(PropagateDownwind):
def __init__(self, site, windTurbines, rotorAvgModel=RotorCenter(), superpositionModel=LinearSum(),
deflectionModel=None, turbulenceModel=None):
PropagateDownwind.__init__(self, site, windTurbines, wake_deficitModel=GCLDeficitModel(),
rotorAvgModel=rotorAvgModel, superpositionModel=superpositionModel,
deflectionModel=deflectionModel, turbulenceModel=turbulenceModel)
def main():
if __name__ == '__main__':
from py_wake.examples.data.iea37._iea37 import IEA37Site
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
import matplotlib.pyplot as plt
# setup site, turbines and wind farm model
site = IEA37Site(16)
x, y = site.initial_position.T
windTurbines = IEA37_WindTurbines()
wf_model = GCL(site, windTurbines)
plt.figure()
print(wf_model)
# run wind farm simulation
sim_res = wf_model(x, y)
# calculate AEP
aep = sim_res.aep()
# plot wake map
flow_map = sim_res.flow_map(wd=30, ws=9.8)
flow_map.plot_wake_map()
flow_map.plot_windturbines()
plt.title('AEP: %.2f GWh' % aep)
plt.show()
main()