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 5cf38b2163c28a0c5c25589ce3b5bb30ee97abbb..8a162c3671098baa75d222089136d8bcb529d0b2 100644
--- a/py_wake/tests/test_windturbines/test_power_ct_curves.py
+++ b/py_wake/tests/test_windturbines/test_power_ct_curves.py
@@ -3,10 +3,13 @@ import pytest
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
-from py_wake.examples.data import hornsrev1
+from py_wake.examples.data import hornsrev1, wtg_path
from py_wake.tests import npt
from py_wake.wind_turbines.power_ct_functions import CubePowerSimpleCt, PowerCtNDTabular, DensityScale, \
- PowerCtTabular, PowerCtFunction, PowerCtFunctionList, PowerCtXr
+ PowerCtTabular, PowerCtFunction, PowerCtFunctionList, PowerCtXr, PowerCtGeneric
+from py_wake.examples.data.hornsrev1 import V80
+from py_wake.wind_turbines._wind_turbines import WindTurbine
+from py_wake.examples.data.dtu10mw import DTU10MW
@pytest.mark.parametrize('method,unit,p_scale,p_ref,ct_ref', [
@@ -168,6 +171,70 @@ def test_FunctionalPowerCtCurve():
0.245, 0.092, 0.031, 0.03]) * 1.3 / 1.225, 3)
+def test_PowerCtGeneric():
+ 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)]:
+ power_norm = ref.power(15)
+ curve = PowerCtGeneric(
+ power_norm=power_norm / 1000,
+ diameter=ref.diameter(),
+ turbulence_intensity=ti,
+ ws_cutin=None,
+ max_cp=.49,
+ constant_ct=.8)
+
+ if 0:
+ u = np.arange(0, 30, .1)
+ p, ct = curve(u)
+ plt.plot(u, p / 1e6, label='Generic')
+
+ plt.plot(u, ref.power(u) / 1e6, label=ref.name())
+
+ plt.ylabel('Power [MW]')
+ plt.legend()
+ ax = plt.twinx()
+ ax.plot(u, ct, '--')
+ ax.plot(u, ref.ct(u), '--')
+ plt.ylabel('Ct')
+ plt.show()
+
+ u = np.arange(5, 25)
+ p, ct = curve(u)
+ p_ref, ct_ref = ref.power_ct(u)
+ # print(np.abs(p_ref - p).max() / power_norm)
+ npt.assert_allclose(p, p_ref, atol=power_norm * p_tol)
+ # print(np.abs(ct_ref - ct).max())
+ npt.assert_allclose(ct, ct_ref, atol=ct_tol)
+
+
+def test_PowerCtGeneric2():
+ ref = V80()
+ power_norm = ref.power(15)
+ curve = PowerCtGeneric(
+ power_norm=power_norm / 1000,
+ diameter=ref.diameter(),
+ turbulence_intensity=0,
+ ws_cutin=3,
+ max_cp=.49,
+ constant_ct=.8)
+
+ if 0:
+ u = np.arange(0, 30, .1)
+ p, ct = curve(u)
+ plt.plot(u, p / 1e6, label='Generic')
+
+ plt.plot(u, ref.power(u) / 1e6, label=ref.name())
+
+ plt.ylabel('Power [MW]')
+ plt.legend()
+ ax = plt.twinx()
+ ax.plot(u, ct, '--')
+ ax.plot(u, ref.ct(u), '--')
+ plt.ylabel('Ct')
+ plt.show()
+
+
def get_continuous_curve(key, optional):
u_p, p = np.asarray(hornsrev1.power_curve).T.copy()
ct = hornsrev1.ct_curve[:, 1]
diff --git a/py_wake/utils/generic_power_ct_curves.py b/py_wake/utils/generic_power_ct_curves.py
new file mode 100644
index 0000000000000000000000000000000000000000..f25f7f738ddba6ebf70879b76e5207b66d29fbc6
--- /dev/null
+++ b/py_wake/utils/generic_power_ct_curves.py
@@ -0,0 +1,90 @@
+import numpy as np
+from numpy import newaxis as na
+from scipy.optimize import optimize
+from scipy.optimize.zeros import newton
+
+
+def standard_power_ct_curve(power_norm, diameter, turbulence_intensity=.1,
+
+ rho=1.225, max_cp=.49, constant_ct=.8,
+ gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03,
+ wsp_lst=np.arange(3, 25, .1)):
+ """Generate standard power curve, extracted from WETB (original extracted from excel sheet made by Kenneth Thomsen)
+
+ Parameters
+ ----------
+ power_norm : int or float
+ Nominal power [kW]
+ diameter : int or float
+ Rotor diameter [m]
+ turbulence_intensity : float
+ Turbulence intensity [%]
+ rho : float, optional
+ Density of air [kg/m^3], default is 1.225
+ max_cp : float, optional
+ Maximum power coefficient
+ constant_ct : float, optional
+ Ct value in constant-ct region
+ gear_loss_const : float, optional
+ Constant gear loss [% of power_norm in kW]
+ gear_loss_var : float, optional
+ Variable gear loss [%]
+ generator_loss : float, optional
+ Generator loss [%]
+ converter_loss : float, optional
+ converter loss [%]
+ """
+
+ area = (diameter / 2) ** 2 * np.pi
+
+ sigma_lst = wsp_lst * turbulence_intensity
+
+ p_aero = .5 * rho * area * wsp_lst ** 3 * max_cp / 1000
+
+ # calc power - gear, generator and conv loss
+ gear_loss = gear_loss_const * power_norm + gear_loss_var * p_aero
+ p_gear = p_aero - gear_loss
+ p_gear[p_gear < 0] = 0
+ p_generator_loss = generator_loss * p_gear
+ p_gen = p_gear - p_generator_loss
+ p_converter_loss = converter_loss * p_gen
+ p_raw = p_gen - p_converter_loss
+ p_raw[p_raw > power_norm] = power_norm
+
+ if turbulence_intensity > 0:
+ sigma = sigma_lst[:, na]
+ ndist = 1 / np.sqrt(2 * np.pi * sigma ** 2) * np.exp(-(wsp_lst - wsp_lst[:, na]) ** 2 / (2 * sigma ** 2))
+ power_lst = (ndist * p_raw).sum(1) / ndist.sum(1)
+ else:
+ power_lst = p_raw
+
+ p_gen = power_lst / (1 - converter_loss)
+ 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)
+ constant_ct_idx = np.where(np.abs(np.diff(cp)) < 1e-4)[0][0]
+ 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])
+ return 4 * a * (1 - a)
+
+ ct_lst = cp2ct(cp)
+ f = constant_ct / ct_lst[constant_ct_idx]
+ ct_lst = np.minimum(8 / 9, ct_lst * f)
+ return wsp_lst, power_lst, ct_lst
+
+
+def main():
+ if __name__ == '__main__':
+
+ import matplotlib.pyplot as plt
+ u = np.linspace(0., 30, 500)
+ u, p, ct = standard_power_ct_curve(10000, 178.3, 0.1)
+ plt.plot(u, p)
+ ax = plt.twinx()
+ ax.plot(u, ct)
+ plt.show()
+
+
+main()
diff --git a/py_wake/wind_turbines/power_ct_functions.py b/py_wake/wind_turbines/power_ct_functions.py
index 6db35ee3b80b2da4672ce765635473012569c9ec..b83a83b6d0b8c9c4fa7f985ae8e38396f1a74243 100644
--- a/py_wake/wind_turbines/power_ct_functions.py
+++ b/py_wake/wind_turbines/power_ct_functions.py
@@ -8,6 +8,7 @@ from py_wake.utils.gradients import fd
from scipy.interpolate.fitpack2 import UnivariateSpline
from autograd.numpy.numpy_boxes import ArrayBox
from autograd.builtins import SequenceBox
+from py_wake.utils.generic_power_ct_curves import standard_power_ct_curve
class PowerCtModel():
@@ -502,3 +503,61 @@ class CubePowerSimpleCt(PowerCtFunction):
self.dct_rated2cutout(ws),
0) # constant ct
return np.asarray([dp, dct])
+
+
+class PowerCtGeneric(PowerCtTabular):
+ def __init__(self, power_norm, diameter, turbulence_intensity=.1,
+ rho=1.225, max_cp=.49, constant_ct=.8,
+ gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03,
+ wsp_lst=np.arange(.1, 25, .1), ws_cutin=None,
+ ws_cutout=None, power_idle=0, ct_idle=0, method='linear', additional_models=default_additional_models):
+ """Generate generic standard power curve based on max_cp, rated power and losses.
+ Ct is computed from the basic 1d momentum theory
+
+ Parameters
+ ----------
+ power_norm : int or float
+ Nominal power [kW]
+ diameter : int or float
+ Rotor diameter [m]
+ turbulence_intensity : float
+ Turbulence intensity [%]
+ rho : float optional
+ Density of air [kg/m^3], defualt is 1.225
+ max_cp : float
+ Maximum power coefficient
+ constant_ct : float, optional
+ Ct value in constant-ct region
+ gear_loss_const : float
+ Constant gear loss [%]
+ gear_loss_var : float
+ Variable gear loss [%]
+ generator_loss : float
+ Generator loss [%]
+ converter_loss : float
+ converter loss [%]
+ ws_cutin : number or None, optional
+ if number, then the range [0,ws_cutin[ will be set to power_idle and ct_idle
+ ws_cutout : number or None, optional
+ if number, then the range ]ws_cutout,100] will be set to power_idle and ct_idle
+ power_idle : number, optional
+ see ws_cutin and ws_cutout
+ ct_idle : number, optional
+ see ws_cutin and ws_cutout
+ method : {'linear', 'phip','spline}
+ Interpolation method:\n
+ - linear: fast, discontinous gradients\n
+ - pchip: smooth\n
+ - spline: smooth, closer to linear, small overshoots in transition to/from constant plateaus)
+ additional_models : list, optional
+ list of additional models.
+ """
+
+ u, p, ct = standard_power_ct_curve(power_norm, diameter, turbulence_intensity, rho, max_cp, constant_ct,
+ gear_loss_const, gear_loss_var, generator_loss, converter_loss, wsp_lst)
+
+ if ws_cutin is not None:
+ ct[u < ws_cutin] = ct_idle
+ PowerCtTabular.__init__(self, 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)