add PowerCtGeneric with power and ct based on generic analytical/empirical formulation

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]

py_wake/utils/generic_power_ct_curves.py

+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()

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)