diff --git a/py_wake/deflection_models/jimenez.py b/py_wake/deflection_models/jimenez.py
index bd7fe01f4b6af47460621faf8df94a03af36849d..a8f4b525db43e1082c5fabf364c13f68e8206255 100644
--- a/py_wake/deflection_models/jimenez.py
+++ b/py_wake/deflection_models/jimenez.py
@@ -9,22 +9,25 @@ class JimenezWakeDeflection(DeflectionModel):
the wake deflection of a wind turbine in yaw. Wind Energ., 13: 559-572. doi:10.1002/we.380
"""
- args4deflection = ['D_src_il', 'yaw_ilk', 'ct_ilk']
+ args4deflection = ['D_src_il', 'yaw_ilk', 'ct_ilk', 'tilt_ilk']
def __init__(self, N=20, beta=.1):
self.beta = beta
self.N = N
- def calc_deflection(self, dw_ijl, hcw_ijl, dh_ijl, D_src_il, yaw_ilk, ct_ilk, **kwargs):
+ def calc_deflection(self, dw_ijl, hcw_ijl, dh_ijl, D_src_il, yaw_ilk, tilt_ilk, ct_ilk, **kwargs):
dw_lst = (np.logspace(0, 1.1, self.N) - 1) / (10**1.1 - 1)
dw_ijxl = dw_ijl[:, :, na] * dw_lst[na, na, :, na]
- theta_ilk = np.deg2rad(yaw_ilk)
+ theta_yaw_ilk, theta_tilt_ilk = np.deg2rad(yaw_ilk), np.deg2rad(-tilt_ilk)
+ theta_ilk = np.sqrt(theta_yaw_ilk**2 + theta_tilt_ilk**2)
+ theta_deflection_ilk = np.arctan2(theta_tilt_ilk, theta_yaw_ilk)
denominator_ilk = np.cos(theta_ilk)**2 * np.sin(theta_ilk) * (ct_ilk / 2)
nominator_ijxl = (1 + (self.beta / D_src_il)[:, na, na, :] * dw_ijxl)**2
alpha = denominator_ilk[:, na, na] / nominator_ijxl[..., na]
deflection_ijlk = np.trapz(np.sin(alpha), dw_ijxl[..., na], axis=2)
- self.hcw_ijlk = hcw_ijl[..., na] + deflection_ijlk
- return dw_ijl[..., na], self.hcw_ijlk, dh_ijl[..., na]
+ self.hcw_ijlk = hcw_ijl[..., na] + deflection_ijlk * np.cos(theta_deflection_ilk[..., na])
+ self.dh_ijlk = dh_ijl[..., na] + deflection_ijlk * np.sin(theta_deflection_ilk[..., na])
+ return dw_ijl[..., na], self.hcw_ijlk, self.dh_ijlk
def main():
diff --git a/py_wake/flow_map.py b/py_wake/flow_map.py
index 4d2915e0072cc7dbce25120374a4e4706f9acd37..d038bba728bc9be60b2ffc503615170f3518f964 100644
--- a/py_wake/flow_map.py
+++ b/py_wake/flow_map.py
@@ -174,14 +174,15 @@ class FlowMap(FlowBox):
def plot_windturbines(self, normalize_with=1, ax=None):
fm = self.windFarmModel
yaw = self.simulationResult.yaw.sel(wd=self.wd[0]).mean(['ws']).data
+ tilt = self.simulationResult.tilt.sel(wd=self.wd[0]).mean(['ws']).data
if self.plane[0] == "YZ":
x_i, y_i = self.simulationResult.x.values, self.simulationResult.y.values
h_i = self.simulationResult.h.values
z_i = self.simulationResult.windFarmModel.site.elevation(x_i, y_i)
- fm.windTurbines.plot_yz(y_i, z_i, h_i, wd=self.wd, yaw=yaw, normalize_with=normalize_with, ax=ax)
+ fm.windTurbines.plot_yz(y_i, z_i, h_i, wd=self.wd, yaw=yaw, tilt=tilt, normalize_with=normalize_with, ax=ax)
else: # self.plane[0] == "XY":
fm.windTurbines.plot_xy(self.simulationResult.x, self.simulationResult.y, self.simulationResult.type.data,
- wd=self.wd, yaw=yaw, normalize_with=normalize_with, ax=ax)
+ wd=self.wd, yaw=yaw, tilt=tilt, normalize_with=normalize_with, ax=ax)
def plot_wake_map(self, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
normalize_with=1, ax=None):
diff --git a/py_wake/tests/test_deflection_models/test_jimenez.py b/py_wake/tests/test_deflection_models/test_jimenez.py
new file mode 100644
index 0000000000000000000000000000000000000000..62ed1913070ed076769f47b9e345fbb2aa7d4bab
--- /dev/null
+++ b/py_wake/tests/test_deflection_models/test_jimenez.py
@@ -0,0 +1,41 @@
+from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussian
+from py_wake.examples.data.hornsrev1 import V80
+from py_wake.examples.data.iea37._iea37 import IEA37Site
+from py_wake.deflection_models.jimenez import JimenezWakeDeflection
+from py_wake.flow_map import YZGrid, XYGrid, Points
+import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+from py_wake.tests import npt
+
+
+@pytest.mark.parametrize('yaw,tilt,cx,cz', [
+ (30, 0, -30, 70),
+ (0, 30, 0, 100),
+ (29.5805, 5, -30, 75)])
+def test_yaw_tilt(yaw, tilt, cx, cz):
+ site = IEA37Site(16)
+ x, y = [0], [0]
+ windTurbines = V80()
+
+ D = windTurbines.diameter()
+ wfm = IEA37SimpleBastankhahGaussian(site, windTurbines, deflectionModel=JimenezWakeDeflection())
+ x_lst = np.linspace(-200, 1000, 100)
+ y_lst = np.linspace(-100, 100, 201)
+ z_lst = np.arange(10, 200, 1)
+
+ fm_yz = wfm(x, y, wd=270, ws=10, yaw=yaw, tilt=tilt).flow_map(YZGrid(x=5 * D, y=y_lst, z=z_lst))
+ fm_xz = wfm(x, y, wd=180, ws=10, yaw=yaw, tilt=tilt).flow_map(YZGrid(x=0, y=x_lst, z=z_lst))
+ fm_xy = wfm(x, y, wd=270, ws=10, yaw=yaw, tilt=tilt).flow_map(XYGrid(x=x_lst))
+ if 0:
+ axes = plt.subplots(3, 1)[1].flatten()
+ fm_xy.plot_wake_map(ax=axes[0])
+ axes[0].axvline(5 * D)
+ fm_xz.plot_wake_map(ax=axes[1])
+ axes[1].axvline(5 * D)
+ fm_yz.plot_wake_map(ax=axes[2])
+ plt.show()
+
+ wake_center = fm_yz.WS_eff.where(fm_yz.WS_eff == fm_yz.WS_eff.min(), drop=True).squeeze()
+ npt.assert_allclose(wake_center.y, cx, atol=.1)
+ npt.assert_allclose(wake_center.h, cz, atol=.24)
diff --git a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
index 8c6ac2c4b7f2a6e08015bc13544747fbecbe8cd3..ee131fbea11837ab56b6731c5426d9713f2de247 100644
--- a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
+++ b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
@@ -7,7 +7,7 @@ from py_wake.site._site import UniformSite
from py_wake.tests import npt
from py_wake.examples.data.hornsrev1 import HornsrevV80, Hornsrev1Site, wt_x, wt_y, V80
from py_wake.tests.test_files.fuga import LUT_path_2MW_z0_0_03
-from py_wake.flow_map import HorizontalGrid
+from py_wake.flow_map import HorizontalGrid, XYGrid
from py_wake.wind_farm_models.engineering_models import All2AllIterative, PropagateDownwind
from py_wake.deficit_models.noj import NOJDeficit
from py_wake.superposition_models import SquaredSum, WeightedSum
@@ -173,7 +173,10 @@ def test_All2AllIterativeDeflection():
wake_deficitModel=NOJDeficit(),
superpositionModel=SquaredSum(),
deflectionModel=JimenezWakeDeflection())
- wf_model([0], [0])
+ sim_res = wf_model([0, 500], [0, 0], wd=270, ws=10, yaw=30)
+ if 0:
+ sim_res.flow_map(XYGrid(x=np.linspace(-200, 1000, 100))).plot_wake_map()
+ plt.show()
def test_dAEP_2wt():
diff --git a/py_wake/tests/test_windturbines/test_power_ct_curves.py b/py_wake/tests/test_windturbines/test_power_ct_curves.py
index 9b2f7f9e0061f10e9e86a8dbdf95fa11cc1255a0..b4bf0f388ca37a667931109514479860cb8e2d3b 100644
--- a/py_wake/tests/test_windturbines/test_power_ct_curves.py
+++ b/py_wake/tests/test_windturbines/test_power_ct_curves.py
@@ -30,7 +30,7 @@ def test_TabularPowerCtCurve(method, unit, p_scale, p_ref, ct_ref):
u_ct, ct = hornsrev1.ct_curve.T
npt.assert_array_equal(u_p, u_ct)
curve = PowerCtTabular(ws=u_p, power=p, power_unit=unit, ct=ct, ws_cutin=4, ws_cutout=25, method=method)
- npt.assert_array_equal(curve.optional_inputs, ['Air_density', 'yaw'])
+ npt.assert_array_equal(curve.optional_inputs, ['Air_density', 'tilt', 'yaw'])
npt.assert_array_equal(curve.required_inputs, [])
u = np.arange(0, 30, .1)
@@ -61,7 +61,7 @@ def test_MultiPowerCtCurve():
curve = PowerCtFunctionList('mode', [PowerCtTabular(ws=u_p, power=p, power_unit='w', ct=ct),
PowerCtTabular(ws=u_p, power=p * 1.1, power_unit='w', ct=ct + .1)])
- npt.assert_array_equal(sorted(curve.optional_inputs), ['Air_density', 'yaw'])
+ npt.assert_array_equal(sorted(curve.optional_inputs), ['Air_density', 'tilt', 'yaw'])
npt.assert_array_equal(list(curve.required_inputs), ['mode'])
u = np.arange(0, 30, .1)
@@ -98,7 +98,7 @@ def test_MultiMultiPowerCtCurve_subset():
PowerCtTabular(ws=u_p, power=p + 6, power_unit='w', ct=ct)]),
])
- npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'Air_density'])
+ npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'tilt', 'Air_density'])
npt.assert_array_equal(sorted(curve.required_inputs)[::-1], ['type', 'mode'])
u = np.zeros((2, 3, 4)) + np.arange(3, 7)[na, na, :]
@@ -125,7 +125,7 @@ def test_2d_tabular():
curve = PowerCtNDTabular(['ws', 'boost'], [u_p, [0, 1]],
np.array([p_c, 2 * p_c]).T, 'w',
np.array([ct_c, ct_c]).T)
- npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'Air_density'])
+ npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'tilt', 'Air_density'])
npt.assert_array_equal(sorted(curve.required_inputs)[::-1], ['boost'])
u = np.linspace(3, 25, 10)
@@ -175,7 +175,7 @@ def test_2d_tabular_default_value():
def test_FunctionalPowerCtCurve():
curve = CubePowerSimpleCt()
- npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'Air_density'])
+ npt.assert_array_equal(sorted(curve.optional_inputs)[::-1], ['yaw', 'tilt', 'Air_density'])
npt.assert_array_equal(curve.required_inputs, [])
u = np.arange(0, 30, .1)
p, ct = curve(u, Air_density=1.3)
diff --git a/py_wake/tests/test_windturbines/test_power_ct_wind_turbines.py b/py_wake/tests/test_windturbines/test_power_ct_wind_turbines.py
index 1f5373fafffa903929be52c78e45aad69621119a..a7d03f79f8572fe911b3fe6c7acd51516ccebdce 100644
--- a/py_wake/tests/test_windturbines/test_power_ct_wind_turbines.py
+++ b/py_wake/tests/test_windturbines/test_power_ct_wind_turbines.py
@@ -58,7 +58,7 @@ def test_PowerCtTabular(method, unit, p_scale, p_ref, ct_ref):
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, [])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
sim_res = wfm([0], [0], ws=u, wd=[0]).squeeze()
p = sim_res.Power.values
ct = sim_res.CT.values
@@ -89,7 +89,7 @@ def test_MultiPowerCtCurve():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['mode'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
sim_res = wfm([0], [0], ws=u, wd=0, mode=0)
p0, ct0 = sim_res.Power.squeeze().values, sim_res.CT.squeeze().values,
@@ -126,7 +126,7 @@ def test_MultiMultiPowerCtCurve_subset():
wfm = get_wfm(curves)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['mode', 'mytype'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
u = np.zeros((2, 3, 4)) + np.arange(3, 7)[na, na, :]
type_2 = np.array([0, 1])
@@ -157,7 +157,7 @@ def test_2d_tabular():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['boost'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
sim_res = wfm([0], [0], wd=0, ws=u, boost=0)
p = sim_res.Power.values
@@ -198,7 +198,7 @@ def test_PowerCtXr():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['boost'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
sim_res = wfm([0], [0], wd=0, ws=u, boost=0)
p = sim_res.Power.values
@@ -232,7 +232,7 @@ def test_FunctionalPowerCtCurve():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, [])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
u = np.arange(0, 30, .1)
sim_res = wfm([0], [0], wd=0, ws=u)
@@ -261,7 +261,7 @@ def test_continuous():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, [])
- npt.assert_array_equal(oi, ['Air_density', 'boost', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'boost', 'tilt', 'yaw'])
u = np.arange(0, 30, .1)
sim_res = wfm([0], [0], wd=0, ws=u)
@@ -294,7 +294,7 @@ def test_unused_input():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, [])
- npt.assert_array_equal(oi, ['Air_density', 'boost', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'boost', 'tilt', 'yaw'])
with pytest.raises(TypeError, match=r"got unexpected keyword argument\(s\): 'mode'"):
wfm([0], [0], boost=1, mode=1)
@@ -304,7 +304,7 @@ def test_missing_input():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['boost'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
with pytest.raises(KeyError, match="Argument, boost, required to calculate power and ct not found"):
wfm([0], [0])
@@ -319,7 +319,7 @@ def test_missing_input_PowerCtFunctionList():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['mode'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
with pytest.raises(KeyError, match="Argument, mode, required to calculate power and ct not found"):
wfm([0], [0])
@@ -332,7 +332,7 @@ def test_DensityScaleAndSimpleYawModel():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, [])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
u = np.arange(4, 25, 1.1)
yaw = 30
@@ -378,7 +378,7 @@ def test_WSFromLocalWind():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['WD'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
u = np.arange(4, 25, .1)
sim_res = wfm([0], [0], wd=[2], ws=u)
p = sim_res.Power.values.squeeze()
@@ -398,7 +398,7 @@ def test_TIFromWFM():
wfm = get_wfm(curve)
ri, oi = wfm.windTurbines.function_inputs
npt.assert_array_equal(ri, ['TI_eff'])
- npt.assert_array_equal(oi, ['Air_density', 'yaw'])
+ npt.assert_array_equal(oi, ['Air_density', 'tilt', 'yaw'])
u = np.arange(4, 25, .1)
with pytest.raises(KeyError, match='Argument, TI_eff, needed to calculate power and ct requires a TurbulenceModel'):
wfm([0], [0], wd=[2], ws=u)
diff --git a/py_wake/utils/model_utils.py b/py_wake/utils/model_utils.py
index a2f593c9391928f665f6d25567dbdb038028759a..702767f8baf89ca76d2562d135dbe6d906dc2ea0 100644
--- a/py_wake/utils/model_utils.py
+++ b/py_wake/utils/model_utils.py
@@ -113,7 +113,7 @@ def get_signature(cls, kwargs={}, indent_level=0):
return "%s(%s)" % (cls.__name__, arg_str)
-def get_model_input(wfm, x, y, ws=10, wd=270, yaw=[[[0]]]):
+def get_model_input(wfm, x, y, ws=10, wd=270, yaw=[[[0]]], tilt=[[[0]]]):
ws, wd = [np.atleast_1d(v) for v in [ws, wd]]
x, y = map(np.asarray, [x, y])
wfm.site.distance.setup(src_x_i=[0], src_y_i=[0], src_h_i=[0],
@@ -124,6 +124,7 @@ def get_model_input(wfm, x, y, ws=10, wd=270, yaw=[[[0]]]):
args = {'dw_ijl': dw_ijl, 'hcw_ijl': hcw_ijl, 'dh_ijl': dh_ijl,
'D_src_il': np.atleast_1d(wfm.windTurbines.diameter())[na]}
args.update({k: sim_res[n].ilk() for k, n in [('yaw_ilk', 'yaw'),
+ ('tilt_ilk', 'tilt'),
('WS_ilk', 'WS'),
('WS_eff_ilk', 'WS_eff'),
('ct_ilk', 'CT')]})
diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index 98e2dd75199cd57ba9be4fe8d15d4afb8d40d054..7e975e1f3964226fd721ad39b1ca37e952690b39 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -150,7 +150,8 @@ class EngineeringWindFarmModel(WindFarmModel):
deficit, blockage = self._add_blockage(deficit, dw_ijlk, **kwargs)
return deficit, uc, sigma_sqr, blockage
- def calc_wt_interaction(self, x_i, y_i, h_i=None, type_i=0, wd=None, ws=None, time=False, yaw_ilk=None, **kwargs):
+ def calc_wt_interaction(self, x_i, y_i, h_i=None, type_i=0, wd=None, ws=None, time=False,
+ yaw_ilk=None, tilt_ilk=None, **kwargs):
"""See WindFarmModel.calc_wt_interaction"""
h_i, D_i = self.windTurbines.get_defaults(len(x_i), type_i, h_i)
x_i, y_i, type_i = [np.asarray(v) for v in [x_i, y_i, type_i]]
@@ -200,8 +201,8 @@ class EngineeringWindFarmModel(WindFarmModel):
def add_arg(name, optional):
if name in wt_kwargs: # custom WindFarmModel.__call__ arguments
return
- elif name in {'yaw', 'type'}: # fixed WindFarmModel.__call__ arguments
- wt_kwargs[name] = {'yaw': yaw_ilk, 'type': type_i}[name]
+ elif name in {'yaw', 'tilt', 'type'}: # fixed WindFarmModel.__call__ arguments
+ wt_kwargs[name] = {'yaw': yaw_ilk, 'tilt': tilt_ilk, 'type': type_i}[name]
elif name in lw:
wt_kwargs[name] = lw[name].values
elif name in self.site.ds:
@@ -223,10 +224,13 @@ class EngineeringWindFarmModel(WindFarmModel):
if yaw_ilk is None:
yaw_ilk = np.zeros((I, L, K))
+ if tilt_ilk is None:
+ tilt_ilk = np.zeros((I, L, K))
kwargs = {'localWind': lw,
'WS_eff_ilk': WS_eff_ilk, 'TI_eff_ilk': TI_eff_ilk,
- 'x_i': x_i, 'y_i': y_i, 'h_i': h_i, 'D_i': D_i, 'yaw_ilk': yaw_ilk,
+ 'x_i': x_i, 'y_i': y_i, 'h_i': h_i, 'D_i': D_i,
+ 'yaw_ilk': yaw_ilk, 'tilt_ilk': tilt_ilk,
'I': I, 'L': L, 'K': K, **wt_kwargs}
WS_eff_ilk, TI_eff_ilk, ct_ilk = self._calc_wt_interaction(**kwargs)
if 'TI_eff' in wt_kwargs:
@@ -281,6 +285,7 @@ class EngineeringWindFarmModel(WindFarmModel):
'TI_ilk': get_ilk('TI'),
'TI_eff_ilk': get_ilk('TI_eff'),
'yaw_ilk': get_ilk('yaw'),
+ 'tilt_ilk': get_ilk('tilt'),
'D_src_il': lambda: wt_d_i[:, na],
'D_dst_ijl': lambda: np.zeros_like(dh_ijl),
'h_il': lambda: wt_h_i.data[:, na],
@@ -291,7 +296,7 @@ class EngineeringWindFarmModel(WindFarmModel):
** {k: arg_funcs[k]() for k in self.deflectionModel.args4deflection})
else:
dw_ijlk, hcw_ijlk, dh_ijlk = dw_ijl[..., na], hcw_ijl[..., na], dh_ijl[..., na]
- arg_funcs.update({'cw_ijlk': lambda: np.hypot(dh_ijl[..., na], hcw_ijlk),
+ arg_funcs.update({'cw_ijlk': lambda: np.hypot(dh_ijlk, hcw_ijlk),
'dw_ijlk': lambda: dw_ijlk, 'hcw_ijlk': lambda: hcw_ijlk, 'dh_ijlk': lambda: dh_ijlk})
args = {k: arg_funcs[k]() for k in self.args4deficit if k != 'dw_ijlk'}
@@ -414,7 +419,7 @@ class PropagateDownwind(EngineeringWindFarmModel):
def _calc_wt_interaction(self, localWind,
WS_eff_ilk, TI_eff_ilk,
- x_i, y_i, h_i, D_i, yaw_ilk,
+ x_i, y_i, h_i, D_i, yaw_ilk, tilt_ilk,
I, L, K, **kwargs):
"""
Additional suffixes:
@@ -441,6 +446,7 @@ class PropagateDownwind(EngineeringWindFarmModel):
WS_eff_mk = []
TI_eff_mk = []
yaw_mk = ilk2mk(yaw_ilk)
+ tilt_mk = ilk2mk(tilt_ilk)
ct_jlk = []
if self.turbulenceModel:
@@ -517,6 +523,7 @@ class PropagateDownwind(EngineeringWindFarmModel):
'TI_eff_ilk': lambda: TI_eff_mk[-1][na],
'D_src_il': lambda: D_i[i_wt_l][na],
'yaw_ilk': lambda: yaw_mk[m][na],
+ 'tilt_ilk': lambda: tilt_mk[m][na],
'D_dst_ijl': lambda: D_i[dw_order_indices_dl[:, j + 1:]].T[na],
'h_il': lambda: h_i[i_wt_l][na],
'ct_ilk': lambda: ct_lk[na],
@@ -621,8 +628,7 @@ class All2AllIterative(EngineeringWindFarmModel):
def _calc_wt_interaction(self, localWind,
WS_eff_ilk, TI_eff_ilk,
- x_i, y_i, h_i, D_i, yaw_ilk,
- # dw_iil, hcw_iil, cw_iil, dh_iil, dw_order_indices_dl,
+ x_i, y_i, h_i, D_i, yaw_ilk, tilt_ilk,
I, L, K, **kwargs):
lw = localWind
WS_eff_ilk_last = WS_eff_ilk.copy()
@@ -634,6 +640,7 @@ class All2AllIterative(EngineeringWindFarmModel):
'TI_ilk': lw.TI.ilk((I, L, K)),
'TI_eff_ilk': lw.TI.ilk((I, L, K)),
'yaw_ilk': yaw_ilk,
+ 'tilt_ilk': tilt_ilk,
'D_src_il': D_src_il,
'D_dst_ijl': D_src_il[na],
'cw_ijlk': np.sqrt(hcw_iil**2 + dh_iil**2)[..., na],
@@ -649,7 +656,7 @@ class All2AllIterative(EngineeringWindFarmModel):
args['WS_eff_ilk'] = WS_eff_ilk
if self.deflectionModel:
dw_ijlk, hcw_ijlk, dh_ijlk = self.deflectionModel.calc_deflection(
- dw_ijl=dw_iil, hcw_ijl=dw_iil, dh_ijl=dh_iil, **args)
+ dw_ijl=dw_iil, hcw_ijl=hcw_iil, dh_ijl=dh_iil, **args)
args.update({'dw_ijlk': dw_ijlk, 'hcw_ijlk': hcw_ijlk, 'dh_ijlk': dh_ijlk,
'cw_ijlk': np.hypot(dh_iil[..., na], hcw_ijlk)})
else:
diff --git a/py_wake/wind_farm_models/wind_farm_model.py b/py_wake/wind_farm_models/wind_farm_model.py
index 41a9ada8eefbd51e39f43941b1b92ccca3459dbd..26963921b2436437fc3fc4d4a3ac1b06d2d2c2f9 100644
--- a/py_wake/wind_farm_models/wind_farm_model.py
+++ b/py_wake/wind_farm_models/wind_farm_model.py
@@ -19,7 +19,8 @@ class WindFarmModel(ABC):
self.site = site
self.windTurbines = windTurbines
- def __call__(self, x, y, h=None, type=0, wd=None, ws=None, yaw=None, time=False, verbose=False, **kwargs):
+ def __call__(self, x, y, h=None, type=0, wd=None, ws=None, yaw=None,
+ tilt=None, time=False, verbose=False, **kwargs):
"""Run the wind farm simulation
Parameters
@@ -54,18 +55,20 @@ class WindFarmModel(ABC):
raise ValueError(
'Custom *yaw*-keyword arguments not allowed to avoid confusion with the default "yaw" keyword')
yaw_ilk = fix_shape(yaw, (I, L, K), allow_None=True)
+ tilt_ilk = fix_shape(tilt, (I, L, K), allow_None=True)
if len(x) == 0:
lw = UniformSite([1], 0.1).local_wind(x_i=[], y_i=[], h_i=[], wd=wd, ws=ws)
z = xr.DataArray(np.zeros((0, len(lw.wd), len(lw.ws))), coords=[('wt', []), ('wd', lw.wd), ('ws', lw.ws)])
- return SimulationResult(self, lw, [], yaw, z, z, z, z, kwargs)
- res = self.calc_wt_interaction(x_i=np.asarray(x), y_i=np.asarray(y), h_i=h, type_i=type, yaw_ilk=yaw_ilk,
+ return SimulationResult(self, lw, [], yaw, tilt, z, z, z, z, kwargs)
+ res = self.calc_wt_interaction(x_i=np.asarray(x), y_i=np.asarray(y), h_i=h, type_i=type,
+ yaw_ilk=yaw_ilk, tilt_ilk=tilt_ilk,
wd=wd, ws=ws, time=time, **kwargs)
WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk, localWind, wt_inputs = res
return SimulationResult(self, localWind=localWind,
- type_i=np.zeros(len(x), dtype=int) + type, yaw_ilk=yaw_ilk,
-
+ type_i=np.zeros(len(x), dtype=int) + type,
+ yaw_ilk=yaw_ilk, tilt_ilk=tilt_ilk,
WS_eff_ilk=WS_eff_ilk, TI_eff_ilk=TI_eff_ilk,
power_ilk=power_ilk, ct_ilk=ct_ilk, wt_inputs=wt_inputs)
@@ -146,7 +149,7 @@ class SimulationResult(xr.Dataset):
"""Simulation result returned when calling a WindFarmModel object"""
__slots__ = ('windFarmModel', 'localWind', 'wt_inputs')
- def __init__(self, windFarmModel, localWind, type_i, yaw_ilk,
+ def __init__(self, windFarmModel, localWind, type_i, yaw_ilk, tilt_ilk,
WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk, wt_inputs):
self.windFarmModel = windFarmModel
lw = localWind
@@ -187,6 +190,11 @@ class SimulationResult(xr.Dataset):
else:
self.add_ilk('yaw', yaw_ilk)
self['yaw'].attrs['Description'] = 'Yaw misalignment [deg]'
+ if tilt_ilk is None:
+ self['tilt'] = self.Power * 0
+ else:
+ self.add_ilk('tilt', tilt_ilk)
+ self['tilt'].attrs['Description'] = 'Rotor tilt [deg]'
# for backward compatibility
for k in ['WD', 'WS', 'TI', 'P', 'WS_eff', 'TI_eff']:
@@ -245,7 +253,7 @@ class SimulationResult(xr.Dataset):
self.Weibull_k.ilk())
aep = weighted_power * self.Sector_frequency.ilk() * hours_pr_year * 1e-9
ws = (self.ws.values[1:] + self.ws.values[:-1]) / 2
- return xr.DataArray(aep, [('wt', self.wt), ('wd', self.wd), ('ws', ws)])
+ return xr.DataArray(aep, [('wt', self.wt.values), ('wd', self.wd.values), ('ws', ws)])
else:
weighted_power = power_ilk * self.P.ilk() / norm
if 'time' in self and weighted_power.shape[2] == 1:
@@ -423,7 +431,7 @@ class SimulationResult(xr.Dataset):
ds = xr.load_dataset(filename)
lw = LocalWind(ds.x, ds.y, ds.h, ds.wd, ds.ws, time=False, wd_bin_size=ds.attrs['wd_bin_size'],
WD=ds.WD, WS=ds.WS, TI=ds.TI, P=ds.P)
- sim_res = SimulationResult(wfm, lw, type_i=ds.type.values, yaw_ilk=ds.yaw,
+ sim_res = SimulationResult(wfm, lw, type_i=ds.type.values, yaw_ilk=ds.yaw, tilt_ilk=ds.tilt,
WS_eff_ilk=ds.WS_eff.ilk(), TI_eff_ilk=ds.TI_eff.ilk(), power_ilk=ds.Power, ct_ilk=ds.CT,
wt_inputs={})
diff --git a/py_wake/wind_turbines/_wind_turbines.py b/py_wake/wind_turbines/_wind_turbines.py
index 35fdb9d70bbb00c7d5235043c4566d9aa0130685..09d3f8220db690ea41b9155559aa32009a7eae51 100644
--- a/py_wake/wind_turbines/_wind_turbines.py
+++ b/py_wake/wind_turbines/_wind_turbines.py
@@ -128,7 +128,7 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
def enable_autograd(self):
self.powerCtFunction.enable_autograd()
- def plot_xy(self, x, y, types=None, wd=None, yaw=0, normalize_with=1, ax=None):
+ def plot_xy(self, x, y, types=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None):
"""Plot wind farm layout including type name and diameter
Parameters
@@ -158,6 +158,7 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
assert len(x) == len(y)
types = (np.zeros_like(x) + types).astype(int) # ensure same length as x
yaw = np.zeros_like(x) + yaw
+ tilt = np.zeros_like(x) + tilt
x, y, D = [np.asarray(v) / normalize_with for v in [x, y, self.diameter(types)]]
R = D / 2
@@ -168,8 +169,9 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
ax.plot(x_, y_, 'None', )
else:
for wd_ in np.atleast_1d(wd):
- c, s = np.cos(np.deg2rad(90 + wd_ - yaw_)), np.sin(np.deg2rad(90 + wd_ - yaw_))
- ax.plot([x_ - s * r, x_ + s * r], [y_ - c * r, y_ + c * r], lw=1, color=colors[t])
+ circle = Ellipse((x_, y), 2 * r * np.sin(np.deg2rad(tilt)), 2 * r,
+ angle=90 - wd_ + yaw_, ec=colors[t], fc="None")
+ ax.add_artist(circle)
for t, m, c in zip(np.unique(types), markers, colors):
# ax.plot(np.asarray(x)[types == t], np.asarray(y)[types == t], '%sk' % m, label=self._names[int(t)])
@@ -180,7 +182,7 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
ax.legend(loc=1)
ax.axis('equal')
- def plot_yz(self, y, z=None, h=None, types=None, wd=270, yaw=0, normalize_with=1, ax=None):
+ def plot_yz(self, y, z=None, h=None, types=None, wd=270, yaw=0, tilt=0, normalize_with=1, ax=None):
"""Plot wind farm layout in yz-plane including type name and diameter
Parameters
@@ -216,10 +218,12 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
from matplotlib.patches import Circle
yaw = np.zeros_like(y) + yaw
+ tilt = np.zeros_like(y) + tilt
y, z, h, D = [v / normalize_with for v in [y, z, h, self.diameter(types)]]
for i, (y_, z_, h_, d, t, yaw_) in enumerate(
zip(y, z, h, D, types, yaw)):
- circle = Ellipse((y_, h_ + z_), d * np.sin(np.deg2rad(wd - yaw_)), d, ec=colors[t], fc="None")
+ circle = Ellipse((y_, h_ + z_), d * np.sin(np.deg2rad(wd - yaw_)),
+ d, angle=-tilt, ec=colors[t], fc="None")
ax.add_artist(circle)
ax.plot([y_, y_], [z_, z_ + h_], 'k')
ax.plot(y_, h_, 'None')
@@ -232,8 +236,8 @@ Use WindTurbines(names, diameters, hub_heights, power_ct_funcs) instead""", Depr
ax.legend(loc=1)
ax.axis('equal')
- def plot(self, x, y, type=None, wd=None, yaw=0, normalize_with=1, ax=None):
- return self.plot_xy(x, y, type, wd, yaw, normalize_with, ax)
+ def plot(self, x, y, type=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None):
+ return self.plot_xy(x, y, type, wd, yaw, tilt, normalize_with, ax)
@staticmethod
def from_WindTurbine_lst(wt_lst):
diff --git a/py_wake/wind_turbines/power_ct_functions.py b/py_wake/wind_turbines/power_ct_functions.py
index fcca9caa6c4276df914083caf59a2699cdfe72ce..c6cd9aee40194593112db2460a381bba7357266f 100644
--- a/py_wake/wind_turbines/power_ct_functions.py
+++ b/py_wake/wind_turbines/power_ct_functions.py
@@ -98,11 +98,18 @@ class SimpleYawModel(AdditionalModel):
"""Simple model that replace ws with cos(yaw)*ws and scales the CT output with cos(yaw)**2"""
def __init__(self):
- AdditionalModel.__init__(self, input_keys=['ws', 'yaw'], optional_inputs=['yaw'], output_keys=['power', 'ct'])
-
- def __call__(self, f, ws, yaw=None, **kwargs):
- if yaw is not None:
- co = np.cos(np.deg2rad(fix_shape(yaw, ws, True)))
+ AdditionalModel.__init__(
+ self, input_keys=[
+ 'ws', 'yaw', 'tilt'], optional_inputs=['yaw', 'tilt'],
+ output_keys=['power', 'ct'])
+
+ def __call__(self, f, ws, yaw=None, tilt=None, **kwargs):
+ if yaw is not None or tilt is not None:
+ co = 1
+ if yaw is not None:
+ co *= np.cos(np.deg2rad(fix_shape(yaw, ws, True)))
+ if tilt is not None:
+ co *= np.cos(np.deg2rad(fix_shape(tilt, ws, True)))
power_ct_arr = f(ws * co, **kwargs) # calculate for reduced ws (ws projection on rotor)
if kwargs['run_only'] == 1: # ct
# multiply ct by cos(yaw)**2 to compensate for reduced thrust