replace nummerical ct calculation with cubic solution from plaa

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

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

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)

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)