From fd9e885b2595b7434c12d054355b4c278ef9108c Mon Sep 17 00:00:00 2001
From: mmpe <mmpe@dtu.dk>
Date: Wed, 18 Aug 2021 09:09:46 +0200
Subject: [PATCH] replace nummerical ct calculation with cubic solution from
plaa
---
py_wake/examples/data/dtu10mw/_dtu10mw.py | 3 ++-
.../test_generic_wind_turbines.py | 14 +++++++++-----
py_wake/utils/generic_power_ct_curves.py | 11 +++++++++--
py_wake/wind_turbines/generic_wind_turbines.py | 4 +++-
4 files changed, 23 insertions(+), 9 deletions(-)
diff --git a/py_wake/examples/data/dtu10mw/_dtu10mw.py b/py_wake/examples/data/dtu10mw/_dtu10mw.py
index 9e5d4e804..1888e62cc 100644
--- a/py_wake/examples/data/dtu10mw/_dtu10mw.py
+++ b/py_wake/examples/data/dtu10mw/_dtu10mw.py
@@ -64,7 +64,8 @@ class DTU10MW(WindTurbine):
'DTU10MW',
diameter=178.3,
hub_height=119,
- powerCtFunction=PowerCtTabular(u, p * 1000, 'w', ct_curve[:, 1], method=method))
+ powerCtFunction=PowerCtTabular(u, p * 1000, 'w', ct_curve[:, 1], ws_cutin=4, ws_cutout=25,
+ ct_idle=0.059, method=method))
DTU10WM_RWT = DTU10MW
diff --git a/py_wake/tests/test_windturbines/test_generic_wind_turbines.py b/py_wake/tests/test_windturbines/test_generic_wind_turbines.py
index 0dec1eb35..7c08f0ffd 100644
--- a/py_wake/tests/test_windturbines/test_generic_wind_turbines.py
+++ b/py_wake/tests/test_windturbines/test_generic_wind_turbines.py
@@ -12,12 +12,14 @@ from py_wake.site.xrsite import XRSite
def test_GenericWindTurbine():
- for ref, ti, p_tol, ct_tol in [(V80(), .1, 0.03, .16),
- (WindTurbine.from_WAsP_wtg(wtg_path + "Vestas V112-3.0 MW.wtg"), .05, 0.035, .07),
- (DTU10MW(), .05, 0.06, .13)]:
+ for ref, ti, cut_in, cut_out, p_tol, ct_tol in [(V80(), .1, 4, None, 0.03, .14),
+ (WindTurbine.from_WAsP_wtg(wtg_path +
+ "Vestas V112-3.0 MW.wtg"), .05, 3, 25, 0.035, .07),
+ (DTU10MW(), .05, 4, 25, 0.06, .12)]:
+
power_norm = ref.power(np.arange(10, 20)).max()
wt = GenericWindTurbine('Generic', ref.diameter(), ref.hub_height(), power_norm / 1e3,
- turbulence_intensity=ti, ws_cutin=None)
+ turbulence_intensity=ti, ws_cutin=cut_in, ws_cutout=cut_out)
if 0:
u = np.arange(0, 30, .1)
@@ -27,10 +29,12 @@ def test_GenericWindTurbine():
plt.plot(u, ref.power(u) / 1e6, label=ref.name())
plt.ylabel('Power [MW]')
- plt.legend()
+ plt.legend(loc='center left')
ax = plt.twinx()
ax.plot(u, ct, '--')
ax.plot(u, ref.ct(u), '--')
+
+ ax.set_ylim([0, 1])
plt.ylabel('Ct')
plt.show()
diff --git a/py_wake/utils/generic_power_ct_curves.py b/py_wake/utils/generic_power_ct_curves.py
index 92fd866ee..d78bcb57d 100644
--- a/py_wake/utils/generic_power_ct_curves.py
+++ b/py_wake/utils/generic_power_ct_curves.py
@@ -63,17 +63,24 @@ def standard_power_ct_curve(power_norm, diameter, turbulence_intensity=.1,
p_gear = p_gen / (1 - generator_loss)
p_aero = (p_gear + gear_loss_const * power_norm) / (1 - gear_loss_var)
cp = p_aero * 1000 / (.5 * rho * area * wsp_lst ** 3)
+ # cp is too high at low ws due to constant loss, so limit to 16/27
cp = np.minimum(cp, 16 / 27)
def cp2ct(cp):
- a = np.array([newton(lambda a, cp=cp:4 * a * (1 - a)**2 - cp, .1) for cp in cp])
+ # solve cp = 4 * a * (1 - a)**2 for a
+ y = 27.0 / 16.0 * cp
+ a = 2.0 / 3.0 * (1 - np.cos(np.arctan2(2 * np.sqrt(y * (1.0 - y)), 1 - 2 * y) / 3.0))
return 4 * a * (1 - a)
ct_lst = cp2ct(cp)
+ # scale ct, such that the constant region (~cut-in to rated) equals <constant_ct>
+ # First part (~0-2m/s) is constant at 8/9 due to cp limit and must be disregarded
ct_below_lim = np.where(ct_lst < 8 / 9 - 1e-6)[0][0]
+ # find index of most constant ct after the disregarded 8/9 region
constant_ct_idx = ct_below_lim + np.argmin(np.abs(np.diff(ct_lst[ct_below_lim:])))
- if ct_lst[constant_ct_idx] < cp2ct([max_cp]):
+ if ct_lst[constant_ct_idx] < cp2ct(max_cp):
+ # if TI is high, then there is no constant region
constant_ct_idx = 0
f = constant_ct / ct_lst[constant_ct_idx]
ct_lst = np.minimum(8 / 9, ct_lst * f)
diff --git a/py_wake/wind_turbines/generic_wind_turbines.py b/py_wake/wind_turbines/generic_wind_turbines.py
index bc310fa27..c85c80883 100644
--- a/py_wake/wind_turbines/generic_wind_turbines.py
+++ b/py_wake/wind_turbines/generic_wind_turbines.py
@@ -9,7 +9,7 @@ class GenericWindTurbine(WindTurbine):
air_density=1.225, max_cp=.49, constant_ct=.8,
gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03,
ws_lst=np.arange(.1, 30, .1), ws_cutin=None,
- ws_cutout=None, power_idle=0, ct_idle=0, method='linear',
+ ws_cutout=None, power_idle=0, ct_idle=None, method='linear',
additional_models=[SimpleYawModel()]):
"""Wind turbine with generic standard power curve based on max_cp, rated power and losses.
Ct is computed from the basic 1d momentum theory
@@ -64,6 +64,8 @@ class GenericWindTurbine(WindTurbine):
u, p, ct = [v[u >= ws_cutin] for v in [u, p, ct]]
if ws_cutout is not None:
u, p, ct = [v[u <= ws_cutout] for v in [u, p, ct]]
+ if ct_idle is None:
+ ct_idle = ct[-1]
powerCtFunction = PowerCtTabular(u, p * 1000, 'w', ct, ws_cutin=ws_cutin, ws_cutout=ws_cutout,
power_idle=power_idle, ct_idle=ct_idle, method=method,
additional_models=additional_models)
--
GitLab