Skip to content
Snippets Groups Projects
Commit a5882550 authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

Resolve "Include height in fuga"

parent a619ef52
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
py_wake/aep/_aep.py
+ 5
− 4
View file @ a5882550
......@@ -32,10 +32,10 @@ class AEP():
self.WD_ilk, self.WS_ilk, self.TI_ilk, self.P_lk = self.site.local_wind(x_i, y_i, wd, ws)
# Calculate down-wind and cross-wind distances
dw_iil, cw_iil, dw_order_l = self.site.wt2wt_distances(x_i, y_i, h_i, self.WD_ilk.mean(2))
dw_iil, cw_iil, dh_iil, dw_order_l = self.site.wt2wt_distances(x_i, y_i, h_i, self.WD_ilk.mean(2))
self.WS_eff_ilk, self.TI_eff_ilk, self.power_ilk, self.ct_ilk =\
self.wake_model.calc_wake(self.WS_ilk, self.TI_ilk, dw_iil, cw_iil, dw_order_l, type_i)
self.wake_model.calc_wake(self.WS_ilk, self.TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_l, type_i)
def calculate_AEP(self, x_i, y_i, h_i=None, type_i=None):
self._run_wake_model(x_i, y_i, h_i, type_i)
......@@ -65,8 +65,9 @@ class AEP():
h_j = np.zeros_like(x_j) + h
_, WS_jlk, _, P_lk = self.site.local_wind(x_j, y_j, wd, ws)
dw_ijl, cw_ijl, _ = self.site.distances(x_i, y_i, h_i, x_j, y_j, h_j, self.WD_ilk.mean(2))
WS_eff_jlk = self.wake_model.wake_map(self.WS_ilk, self.WS_eff_ilk, dw_ijl, cw_ijl, self.ct_ilk, type_i, WS_jlk)
dw_ijl, cw_ijl, dh_ijl, _ = self.site.distances(x_i, y_i, h_i, x_j, y_j, h_j, self.WD_ilk.mean(2))
WS_eff_jlk = self.wake_model.wake_map(self.WS_ilk, self.WS_eff_ilk, dw_ijl,
cw_ijl, dh_ijl, self.ct_ilk, type_i, WS_jlk)
return X_j, Y_j, WS_eff_jlk, P_lk
......
py_wake/examples/data/iea37/_iea37.py
+ 1
− 1
View file @ a5882550
......@@ -32,7 +32,7 @@ def main():
site = IEA37_Site(16)
x, y = site.initial_position.T
dw, cw, dw_order = site.wt2wt_distances(x, y, 70, np.array([[0]]))
dw, cw, dh, dw_order = site.wt2wt_distances(x, y, 70, np.array([[0]]))
print(dw.shape)
WD_ilk, WS_ilk, TI_ilk, P_lk = site.local_wind(x, y)
print(WS_ilk.shape)
......
py_wake/site/_site.py
+ 9
− 6
View file @ a5882550
......@@ -131,7 +131,7 @@ class UniformSite(Site):
def distances(self, src_x_i, src_y_i, src_h_i, dst_x_j, dst_y_j, dst_h_j, wd_il):
wd_l = np.mean(wd_il, 0)
dx_ij, dy_ij, dh_ij = [np.subtract(*np.meshgrid(src_i, dst_j, indexing='ij'))
dx_ij, dy_ij, dh_ij = [np.subtract(*np.meshgrid(dst_j, src_i, indexing='ij')).T
for src_i, dst_j in [(src_x_i, dst_x_j),
(src_y_i, dst_y_j),
(src_h_i, dst_h_j)]]
......@@ -141,12 +141,15 @@ class UniformSite(Site):
cos_l = np.cos(theta_l)
sin_l = np.sin(theta_l)
dw_il = (cos_l[na, :] * src_x_i[:, na] + sin_l[na] * src_y_i[:, na])
dw_ijl = (cos_l[na, na, :] * dx_ij[:, :, na] + sin_l[na, na, :] * dy_ij[:, :, na])
cw_ijl = np.sqrt((-sin_l[na, na, :] * dx_ij[:, :, na] +
cos_l[na, na, :] * dy_ij[:, :, na])**2 +
dh_ij[:, :, na]**2)
dw_ijl = (-cos_l[na, na, :] * dx_ij[:, :, na] - sin_l[na, na, :] * dy_ij[:, :, na])
cw_ijl = np.abs(sin_l[na, na, :] * dx_ij[:, :, na] +
-cos_l[na, na, :] * dy_ij[:, :, na])
dh_ijl = np.zeros_like(dw_ijl)
dh_ijl[:, :, :] = dh_ij[:, :, na]
dw_order_indices_l = np.argsort(-dw_il, 0).astype(np.int).T
return dw_ijl, cw_ijl, dw_order_indices_l
return dw_ijl, cw_ijl, dh_ijl, dw_order_indices_l
def elevation(self, x_i, y_i):
return np.zeros_like(x_i)
......
py_wake/tests/test_wake_models/test_fuga.py
+ 21
− 11
View file @ a5882550
......@@ -22,27 +22,37 @@ def test_fuga():
x_j = np.linspace(-1500, 1500, 500)
y_j = np.linspace(-1500, 1500, 300)
X, Y, Z2 = aep.wake_map(x_j, y_j, 70, wt_x, wt_y, h_i=70, wd=[30], ws=[10])
npt.assert_array_almost_equal(
Z2[140, 100:400:10],
[10.0547, 10.0519, 10.0718, 10.0093, 9.6786, 7.8589, 6.8539, 9.2199,
9.9837, 10.036, 10.0796, 10.0469, 10.0439, 9.1866, 7.2552, 9.1518,
10.0449, 10.0261, 10.0353, 9.9256, 9.319, 8.0062, 6.789, 8.3578,
9.9393, 10.0332, 10.0191, 10.0186, 10.0191, 10.0139], 4)
_, _, Z70 = aep.wake_map(x_j, y_j, 70, wt_x, wt_y, h_i=70, wd=[30], ws=[10])
X, Y, Z73 = aep.wake_map(x_j, y_j, 73, wt_x, wt_y, h_i=70, wd=[30], ws=[10])
if 0:
import matplotlib.pyplot as plt
c = plt.contourf(X, Y, Z2, np.arange(6, 10.5, .1))
c = plt.contourf(X, Y, Z70, np.arange(6, 10.5, .1))
plt.colorbar(c)
plt.plot(X[0], Y[140])
wts.plot(wt_x, wt_y)
plt.figure()
plt.plot(X[0], Z2[140, :])
plt.plot(X[0, 100:400:10], Z2[140, 100:400:10], '.')
plt.plot(X[0], Z70[140, :], label="Z=70m")
plt.plot(X[0], Z73[140, :], label="Z=73m")
plt.plot(X[0, 100:400:10], Z70[140, 100:400:10], '.')
plt.legend()
plt.show()
npt.assert_array_almost_equal(
Z70[140, 100:400:10],
[10.0547, 10.0519, 10.0718, 10.0093, 9.6786, 7.8589, 6.8539, 9.2199,
9.9837, 10.036, 10.0796, 10.0469, 10.0439, 9.1866, 7.2552, 9.1518,
10.0449, 10.0261, 10.0353, 9.9256, 9.319, 8.0062, 6.789, 8.3578,
9.9393, 10.0332, 10.0191, 10.0186, 10.0191, 10.0139], 4)
npt.assert_array_almost_equal(
Z73[140, 100:400:10],
[10.0542, 10.0514, 10.0706, 10.0075, 9.6778, 7.9006, 6.9218, 9.228,
9.9808, 10.0354, 10.0786, 10.0464, 10.0414, 9.1973, 7.3099, 9.1629,
10.0432, 10.0257, 10.0344, 9.9236, 9.3274, 8.0502, 6.8512, 8.3813,
9.9379, 10.0325, 10.0188, 10.0183, 10.019, 10.0138], 4)
def cmp_fuga_with_colonel():
# move turbine 1 600 300
......
py_wake/tests/test_wake_models/test_noj.py
+ 7
− 6
View file @ a5882550
......@@ -24,9 +24,9 @@ def test_NOJ_Nibe_result():
WS_ilk = np.array([[[8.1]], [[8.1]], [[8.1]]])
TI_ilk = np.zeros_like(WS_ilk)
site = UniformSite([1], 0.1)
dw_iil, cw_iil, dw_order_indices_l = site.wt2wt_distances(
dw_iil, cw_iil, dh_iil, dw_order_indices_l = site.wt2wt_distances(
x_i=[0, 0, 0], y_i=[0, -40, -100], h_i=[50, 50, 50], wd_il=[[0]])
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dw_order_indices_l, [0, 1, 1])[0]
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_indices_l, [0, 1, 1])[0]
npt.assert_array_almost_equal(WS_eff_ilk[:, 0, 0], [8.1, 4.35, 5.7])
......@@ -58,8 +58,8 @@ def test_NOJ_two_turbines_in_row(wdir, x, y):
WS_ilk = np.array([[[8.1]], [[8.1]], [[8.1]]])
TI_ilk = np.zeros_like(WS_ilk)
site = UniformSite([1], 0.1)
dw_iil, cw_iil, dw_order_indices_l = site.wt2wt_distances(x_i=x, y_i=y, h_i=[50, 50, 50], wd_il=[[wdir]])
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dw_order_indices_l, [0, 0, 0])[0]
dw_iil, cw_iil, dh_iil, dw_order_indices_l = site.wt2wt_distances(x_i=x, y_i=y, h_i=[50, 50, 50], wd_il=[[wdir]])
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_indices_l, [0, 0, 0])[0]
ws_wt3 = 8.1 - np.sqrt((8.1 * 2 / 3 * (20 / 26)**2)**2 + (8.1 * 2 / 3 * (20 / 30)**2)**2)
npt.assert_array_almost_equal(WS_eff_ilk[:, 0, 0], [8.1, 4.35, ws_wt3])
......@@ -74,8 +74,9 @@ def test_NOJ_6_turbines_in_row():
WD_ilk, WS_ilk, _, _ = site.local_wind(x, y, [0], [11])
TI_ilk = np.zeros_like(WS_ilk)
site = UniformSite([1], 0.1)
dw_iil, cw_iil, dw_order_indices_l = site.wt2wt_distances(x_i=x, y_i=y, h_i=[50] * n_wt, wd_il=WD_ilk.mean(2))
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dw_order_indices_l, [0] * n_wt)[0]
dw_iil, cw_iil, dh_iil, dw_order_indices_l = site.wt2wt_distances(
x_i=x, y_i=y, h_i=[50] * n_wt, wd_il=WD_ilk.mean(2))
WS_eff_ilk = wake_model.calc_wake(WS_ilk, TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_indices_l, [0] * n_wt)[0]
np.testing.assert_array_almost_equal(
WS_eff_ilk[1:, 0, 0], 11 - np.sqrt(np.cumsum(((11 * 2 / 3 * 20**2)**2) / (20 + 8 * np.arange(1, 6))**4)))
......
py_wake/wake_model.py
+ 36
− 5
View file @ a5882550
......@@ -6,12 +6,34 @@ import numpy as np
class WakeModel(ABC):
# args4deficit = ['WS_lk', 'WS_eff_lk', 'D_src_l', 'D_dst_jl', 'dw_jl', 'cw_jl', 'ct_lk']
"""
Base class for wake models
Make a subclass and implement calc_deficit and calc_effective_WS
Implementations of linear and squared sum available through inherritance
Prefixs:
i: turbine
j: downstream points/turbines
k: wind speed
l: wind direction
m: turbine and wind direction (il.flatten())
n: from_turbine, to_turbine and wind direction (iil.flatten())
Arguments available for calc_deficit (specifiy in args4deficit):
- WS_lk: Local wind speed without wake effects
- WS_eff_lk: Local wind speed with wake effects
- D_src_l: Diameter of source turbine
- D_dst_jl: Diameter of destination turbine
- dw_jl: Downwind distance from turbine i to point/turbine j
- hcw_jl: Horizontal cross wind distance from turbine i to point/turbine j
- cw_jl: Cross wind(horizontal and vertical) distance from turbine i to point/turbine j
- ct_lk: Thrust coefficient
"""
def __init__(self, windTurbines):
self.windTurbines = windTurbines
def calc_wake(self, WS_ilk, TI_ilk, dw_iil, cw_iil, dw_order_indices_l, types_i):
def calc_wake(self, WS_ilk, TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_indices_l, types_i):
I, L = dw_iil.shape[1:]
i1, i2, _ = np.where((np.abs(dw_iil) + np.abs(cw_iil) + np.eye(I)[:, :, na]) == 0)
if len(i1):
......@@ -28,10 +50,12 @@ class WakeModel(ABC):
WS_eff_mk = WS_mk.copy()
dw_n = dw_iil.flatten()
cw_n = cw_iil.flatten()
dh_n = dh_iil.flatten()
power_ilk = np.zeros((I, L, K))
ct_ilk = np.zeros((I, L, K))
types_i = np.asarray(types_i)
D_i = self.windTurbines.diameter(types_i)
H_i = self.windTurbines.hub_height(types_i)
i_wd_l = np.arange(L)
for j in range(I):
......@@ -56,7 +80,10 @@ class WakeModel(ABC):
'D_src_l': lambda: D_i[i_wt_l],
'D_dst_jl': lambda: D_i[dw_order_indices_l[:, j + 1:]].T,
'dw_jl': lambda: dw_n[n_dw],
'cw_jl': lambda: cw_n[n_dw],
'cw_jl': lambda: np.sqrt(cw_n[n_dw]**2 + dh_n[n_dw]**2),
'hcw_jl': lambda: cw_n[n_dw],
'dh_jl': lambda: dh_n[n_dw],
'h_l': lambda: H_i[i_wt_l],
'ct_lk': lambda: ct_lk}
args = {k: arg_funcs[k]() for k in self.args4deficit}
deficit_nk[n_dw] = self.calc_deficit(**args)
......@@ -68,8 +95,9 @@ class WakeModel(ABC):
WS_eff_ilk = WS_eff_mk.reshape((I, L, K))
return WS_eff_ilk, TI_ilk, power_ilk, ct_ilk
def wake_map(self, WS_ilk, WS_eff_ilk, dw_ijl, cw_ijl, ct_ilk, types_i, WS_jlk):
def wake_map(self, WS_ilk, WS_eff_ilk, dw_ijl, cw_ijl, dh_ijl, ct_ilk, types_i, WS_jlk):
D_i = self.windTurbines.diameter(types_i)
H_i = self.windTurbines.hub_height(types_i)
I, J, L = dw_ijl.shape
K = WS_ilk.shape[2]
......@@ -85,7 +113,10 @@ class WakeModel(ABC):
'D_src_l': lambda: D_i[i][na],
'D_dst_jl': lambda: None,
'dw_jl': lambda: dw_ijl[i, :, l][m][:, na],
'cw_jl': lambda: cw_ijl[i, :, l][m][:, na],
'cw_jl': lambda: np.sqrt(cw_ijl[i, :, l][m]**2 + dh_ijl[i, :, l][m]**2)[:, na],
'hcw_jl': lambda: cw_ijl[i, :, l][m][:, na],
'dh_jl': lambda: dh_ijl[i, :, l][m][:, na],
'h_l': lambda: H_i[i][na],
'ct_lk': lambda: ct_ilk[i, l][na]}
args = {k: arg_funcs[k]() for k in self.args4deficit}
......
py_wake/wake_models/fuga.py
+ 3
− 3
View file @ a5882550
......@@ -8,7 +8,7 @@ from numpy import newaxis as na
class FugaWakeModel(WakeModel, LinearSum):
ams = 5
invL = 0
args4deficit = ['WS_lk', 'WS_eff_lk', 'dw_jl', 'cw_jl', 'ct_lk']
args4deficit = ['WS_lk', 'WS_eff_lk', 'dw_jl', 'hcw_jl', 'dh_jl', 'h_l', 'ct_lk']
def __init__(self, LUT_path, windTurbines):
WakeModel.__init__(self, windTurbines)
......@@ -68,8 +68,8 @@ class FugaWakeModel(WakeModel, LinearSum):
z = np.maximum(np.minimum(z, self.z[-1]), self.z[0])
return self.lut_interpolator((x, y, z))
def calc_deficit(self, WS_lk, WS_eff_lk, dw_jl, cw_jl, ct_lk):
mdu_jl = self.interpolate(dw_jl, cw_jl, 70)
def calc_deficit(self, WS_lk, WS_eff_lk, dw_jl, hcw_jl, dh_jl, h_l, ct_lk):
mdu_jl = self.interpolate(dw_jl, hcw_jl, h_l + dh_jl)
deficit_jlk = mdu_jl[:, :, na] * (ct_lk * WS_eff_lk**2 / WS_lk)
return deficit_jlk
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment