add BastankhahGaussian

py_wake/aep/_aep.py

@@ -37,8 +37,8 @@ class AEP():
         self.WS_eff_ilk, self.TI_eff_ilk, self.power_ilk, self.ct_ilk =\
             self.wake_model.calc_wake(self.WS_ilk, self.TI_ilk, dw_iil, cw_iil, dh_iil, dw_order_l, type_i)
 
-    def calculate_AEP(self, x_i, y_i, h_i=None, type_i=None):
-        self._run_wake_model(x_i, y_i, h_i, type_i)
+    def calculate_AEP(self, x_i, y_i, h_i=None, type_i=None, wd=None, ws=None):
+        self._run_wake_model(x_i=x_i, y_i=y_i, h_i=h_i, type_i=type_i, wd=wd, ws=ws)
         AEP_GWh_ilk = self.power_ilk * self.P_lk[na, :, :] * 24 * 365 * 1e-9
         return AEP_GWh_ilk
 

py_wake/tests/test_wake_models/test_iea37simple_bastankhah.py → py_wake/tests/test_wake_models/test_BastankhahGaussian.py

@@ -6,7 +6,65 @@ from py_wake.examples.data.iea37.iea37_reader import read_iea37_windrose,\
     read_iea37_windfarm
 from py_wake.site._site import UniformSite
 from py_wake.tests import npt
-from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian
+from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian,\
+    BastankhahGaussian
+
+
+def test_BastankhahGaussian_ex16():
+    _, _, freq = read_iea37_windrose(iea37_path + "iea37-windrose.yaml")
+    n_wt = 16
+    x, y, aep_ref = read_iea37_windfarm(iea37_path + 'iea37-ex%d.yaml' % n_wt)
+    if 0:
+        import matplotlib.pyplot as plt
+        plt.plot(x, y, '2k')
+        for i, (x_, y_) in enumerate(zip(x, y)):
+            plt.annotate(i, (x_, y_))
+        plt.axis('equal')
+        plt.show()
+    site = UniformSite(freq, ti=0.75)
+    windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
+    wake_model = BastankhahGaussian(windTurbines)
+
+    aep = AEP(site, windTurbines, wake_model, np.arange(0, 360, 22.5), [9.8])
+    aep_ilk = aep.calculate_AEP(x, y)
+    aep_MW_l = aep_ilk.sum((0, 2)) * 1000
+    # test that the result is equal to last run (no evidens that  these number are correct)
+    aep_ref = (355971.9717035484,
+               [9143.74048, 8156.71681, 11311.92915, 13955.06316, 19807.65346,
+                25196.64182, 39006.65223, 41463.31044, 23042.22602, 12978.30551,
+                14899.26913, 32320.21637, 67039.04091, 17912.40907, 12225.04134,
+                7513.75582])
+    npt.assert_almost_equal(aep_MW_l.sum(), aep_ref[0], 5)
+
+    npt.assert_array_almost_equal(aep_MW_l, aep_ref[1], 5)
+
+
+def test_BastankhahGaussian_wake_map():
+    _, _, freq = read_iea37_windrose(iea37_path + "iea37-windrose.yaml")
+    n_wt = 16
+    x, y, _ = read_iea37_windfarm(iea37_path + 'iea37-ex%d.yaml' % n_wt)
+    site = UniformSite(freq, ti=0.75)
+    windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
+
+    wake_model = BastankhahGaussian(windTurbines)
+    aep = AEP(site, windTurbines, wake_model, np.arange(0, 360, 22.5), [9.8])
+    x_j = np.linspace(-1500, 1500, 200)
+    y_j = np.linspace(-1500, 1500, 100)
+    X, Y, Z = aep.wake_map(x_j, y_j, 110, x, y, wd=[0], ws=[9])
+
+    # test that the result is equal to last run (no evidens that  these number are correct)
+    ref = [0.18, 3.6, 7.27, 8.32, 7.61, 6.64, 5.96, 6.04, 6.8, 7.69, 8.08, 7.87, 7.59, 7.46, 7.55, 7.84, 8.19]
+    if 0:
+        import matplotlib.pyplot as plt
+        c = plt.contourf(X, Y, Z, np.arange(-.25, 9.1, .01))
+        plt.colorbar(c)
+        plt.plot(x, y, '2k')
+        for i, (x_, y_) in enumerate(zip(x, y)):
+            plt.annotate(i, (x_, y_))
+        plt.plot(X[49, 100:133:2], Y[49, 100:133:2], '-.')
+        plt.axis('equal')
+        plt.show()
+    npt.assert_array_almost_equal(Z[49, 100:133:2], ref, 2)
 
 
 def test_IEA37SimpleBastankhahGaussian_ex16():
@@ -27,11 +85,12 @@ def test_IEA37SimpleBastankhahGaussian_ex16():
     aep = AEP(site, windTurbines, wake_model, np.arange(0, 360, 22.5), [9.8])
     aep_ilk = aep.calculate_AEP(x, y)
     aep_MW_l = aep_ilk.sum((0, 2)) * 1000
+    # test that the result is equal to results provided for IEA task 37
     npt.assert_almost_equal(aep_ref[0], aep_MW_l.sum(), 5)
     npt.assert_array_almost_equal(aep_ref[1], aep_MW_l, 5)
 
 
-def test_wake_map():
+def test_IEA37SimpleBastankhahGaussian_wake_map():
     _, _, freq = read_iea37_windrose(iea37_path + "iea37-windrose.yaml")
     n_wt = 16
     x, y, _ = read_iea37_windfarm(iea37_path + 'iea37-ex%d.yaml' % n_wt)
@@ -43,6 +102,7 @@ def test_wake_map():
     y_j = np.linspace(-1500, 1500, 100)
     X, Y, Z = aep.wake_map(x_j, y_j, 110, x, y, wd=[0], ws=[9])
 
+    # test that the result is equal to last run (no evidens that  these number are correct)
     ref = [3.32, 4.86, 7.0, 8.1, 7.8, 7.23, 6.86, 6.9, 7.3, 7.82, 8.11, 8.04, 7.87, 7.79, 7.85, 8.04, 8.28]
     npt.assert_array_almost_equal(Z[49, 100:133:2], ref, 2)
     if 0:
@@ -52,9 +112,6 @@ def test_wake_map():
         plt.plot(x, y, '2k')
         for i, (x_, y_) in enumerate(zip(x, y)):
             plt.annotate(i, (x_, y_))
+        plt.plot(X[49, 100:133:2], Y[49, 100:133:2], '-.')
         plt.axis('equal')
         plt.show()
-
-
-if __name__ == '__main__':
-    test_IEA37SimpleBastankhahGaussian_ex8()

py_wake/wake_models/gaussian.py

@@ -3,7 +3,12 @@ import numpy as np
 from numpy import newaxis as na
 
 
-class IEA37SimpleBastankhahGaussian(WakeModel, SquaredSum):
+class BastankhahGaussian(WakeModel, SquaredSum):
+    """Implemented according to
+    Bastankhah M and Porté-Agel F.
+    A new analytical model for wind-turbine wakes.
+    J. Renew. Energy. 2014;70:116-23.
+    """
     args4deficit = ['WS_lk', 'D_src_l', 'dw_jl', 'cw_jl', 'ct_lk']
 
     def __init__(self, windTurbines):
@@ -11,9 +16,29 @@ class IEA37SimpleBastankhahGaussian(WakeModel, SquaredSum):
         self.k = 0.0324555
 
     def calc_deficit(self, WS_lk, D_src_l, dw_jl, cw_jl, ct_lk):
+        sqrt1ct = np.sqrt(1 - ct_lk)
+        beta_lk = 1 / 2 * (1 + sqrt1ct) / sqrt1ct
+        sigma_sqr_jlk = (self.k * dw_jl[:, :, na] / D_src_l[na, :, na] + .2 * np.sqrt(beta_lk)[na])**2
+        exponent_jlk = -1 / (2 * sigma_sqr_jlk) * (cw_jl**2 / D_src_l[na, :]**2)[:, :, na]
+        # maximum added to avoid sqrt of negative number
+        radical_jlk = np.maximum(0, (1. - ct_lk[na, :, :] / (8. * sigma_sqr_jlk)))
+        deficit_jlk = (WS_lk[na, :, :] * (1. - np.sqrt(radical_jlk)) * np.exp(exponent_jlk))
+        return deficit_jlk
+
+
+class IEA37SimpleBastankhahGaussian(WakeModel, SquaredSum):
+    """Implemented according to
+    https://github.com/byuflowlab/iea37-wflo-casestudies/blob/master/iea37-wakemodel.pdf
+
+    Equivalent to BastankhahGaussian for beta=1/sqrt(8) ~ ct=0.9637188
+    """
+    args4deficit = ['WS_lk', 'D_src_l', 'dw_jl', 'cw_jl', 'ct_lk']
+
+    def __init__(self, windTurbines):
+        WakeModel.__init__(self, windTurbines)
+        self.k = 0.0324555
 
-        # Calculate the wake loss using
-        # simplified Bastankhah Gaussian wake model
+    def calc_deficit(self, WS_lk, D_src_l, dw_jl, cw_jl, ct_lk):
         sigma_jl = self.k * dw_jl + D_src_l[na, :] / np.sqrt(8.)
 
         exponent_jl = -0.5 * (cw_jl / sigma_jl)**2