- Steen Frandsen added turbulence implemented

- ti_map added
- default ti changed from 0.75 to 0.075

py_wake/aep_calculator.py

@@ -108,7 +108,7 @@ class AEPCalculator():
         AEP_GWh_ilk = self.power_ilk * self.P_ilk * 24 * 365 * 1e-9
         return AEP_GWh_ilk
 
-    def _WS_eff_map(self, x_j, y_j, h, x_i, y_i, type_i, h_i, wd=None, ws=None):
+    def _eff_map(self, type, x_j, y_j, h, x_i, y_i, type_i, h_i, wd=None, ws=None):
         h_i, type_i, wd, ws = self._get_defaults(x_i, h_i, type_i, wd, ws)
 
         def f(x, N=500, ext=.2):
@@ -131,12 +131,17 @@ class AEPCalculator():
         self._run_wake_model(x_i, y_i, h_i, type_i, wd, ws)
 
         h_j = np.zeros_like(x_j) + h
-        _, WS_jlk, _, P_ilk = self.site.local_wind(x_i=x_j, y_i=y_j, wd=wd, ws=ws)
+        _, WS_jlk, TI_jlk, P_ilk = self.site.local_wind(x_i=x_j, y_i=y_j, wd=wd, ws=ws)
         dw_ijl, cw_ijl, dh_ijl, _ = self.site.distances(x_i, y_i, h_i, x_j, y_j, h_j, self.WD_ilk.mean(2))
-        WS_eff_jlk = self.wake_model.wake_map(self.WS_ilk, self.WS_eff_ilk, dw_ijl,
+
+        if type == 'WS':
+            _eff_jlk = self.wake_model.ws_map(self.WS_ilk, self.WS_eff_ilk, self.TI_ilk, self.TI_ilk, dw_ijl,
                                               cw_ijl, dh_ijl, self.ct_ilk, type_i, WS_jlk)
+        elif type == 'TI':
+            _eff_jlk = self.wake_model.ti_map(self.WS_ilk, self.WS_eff_ilk, self.TI_ilk, self.TI_ilk, dw_ijl,
+                                              cw_ijl, dh_ijl, self.ct_ilk, type_i, TI_jlk)
 
-        return X_j, Y_j, WS_eff_jlk, P_ilk
+        return X_j, Y_j, _eff_jlk, P_ilk
 
     def wake_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None):
         """Calculate wake(effective wind speed) map
@@ -183,12 +188,62 @@ class AEPCalculator():
         --------
         plot_wake_map
         """
-        X_j, Y_j, WS_eff_jlk, P_ilk = self._WS_eff_map(x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
+        X_j, Y_j, WS_eff_jlk, P_ilk = self._eff_map(
+            'WS', x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
         if P_ilk.sum() > 0:
             WS_eff_jlk = WS_eff_jlk * (P_ilk / P_ilk.sum((1, 2)))
         return X_j, Y_j, WS_eff_jlk.sum((1, 2)).reshape(X_j.shape)
 
-    def plot_wake_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None, ax=None, levels=100):
+    def ti_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None, wd=None, ws=None):
+        """Calculate turbulence intensity map
+
+        Parameters
+        ----------
+        x_j : array_like or None, optional
+            X position map points
+        y_j : array_like
+            Y position of map points
+        h : int, float or None, optional
+            Height of wake map\n
+            If None, default, the mean hub height is used
+        wt_x : array_like, optional
+            X position of wind turbines
+        wt_y : array_like, optional
+            Y position of wind turbines
+        wt_type : array_like or None, optional
+            Type of the wind turbines
+        wt_height : array_like or None, optional
+            Hub height of the wind turbines\n
+            If None, default, the standard hub height is used
+        wd : int, float, array_like or None
+            Wind directions(s)\n
+            If None, default, the wake is calculated for site.default_wd
+        ws : int, float, array_like or None
+            Wind speed(s)\n
+            If None, default, the wake is calculated for site.default_ws
+
+        Returns
+        -------
+        X_j : array_like
+            2d array of map x positions
+        Y_j : array_like
+            2d array of map y positions
+        WS_eff_avg : array_like
+            2d array of average effective local wind speed taking into account
+            the probability of wind direction and speed
+
+        See Also
+        --------
+        plot_wake_map
+        """
+        X_j, Y_j, TI_eff_jlk, P_ilk = self._eff_map(
+            'TI', x_j, y_j, h, wt_x, wt_y, wt_type, wt_height, wd, ws)
+        if P_ilk.sum() > 0:
+            TI_eff_jlk = TI_eff_jlk * (P_ilk / P_ilk.sum((1, 2)))
+        return X_j, Y_j, TI_eff_jlk.sum((1, 2)).reshape(X_j.shape)
+
+    def plot_wake_map(self, x_j=None, y_j=None, h=None, wt_x=[], wt_y=[], wt_type=None, wt_height=None,
+                      wd=None, ws=None, ax=None, levels=100):
         """Plot wake(effective wind speed) map
 
         Parameters
@@ -268,7 +323,7 @@ class AEPCalculator():
             the probability of wind direction and speed
         """
         h_j = self.windTurbines.hub_height(type_j)
-        X_j, Y_j, WS_eff_jlk, P_jlk = self._WS_eff_map(x_j, y_j, h_j, wt_x, wt_y, wt_type, wt_height, wd, ws)
+        X_j, Y_j, WS_eff_jlk, P_jlk = self._eff_map('WS', x_j, y_j, h_j, wt_x, wt_y, wt_type, wt_height, wd, ws)
         P_jlk /= P_jlk.sum((1, 2))  # AEP if wind only comes from specified wd and ws
 
         # power_jlk = self.windTurbines.power_func(WS_eff_jlk, type_j)

py_wake/examples/data/iea37/_iea37.py

@@ -15,7 +15,7 @@ class IEA37_WindTurbines(OneTypeWindTurbines):
 
 
 class IEA37Site(UniformSite):
-    def __init__(self, n_wt, ti=.75):
+    def __init__(self, n_wt, ti=.075):
         assert n_wt in [9, 16, 36, 64]
 
         from py_wake.examples.data.iea37.iea37_reader import \
@@ -31,6 +31,7 @@ class IEA37Site(UniformSite):
 
 class IEA37AEPCalc():
     """Run the AEP calculator provided by IEA Task 37"""
+
     def __init__(self, n_wt):
         assert n_wt in [9, 16, 36, 64]
         self.n_wt = n_wt

py_wake/site/_site.py

@@ -238,8 +238,10 @@ class Site(ABC):
         theta = wd / 180 * np.pi
         if not ax.__class__.__name__ == 'PolarAxesSubplot':
             if hasattr(ax, 'subplot'):
+                ax.clf()
                 ax = ax.subplot(111, projection='polar')
             else:
+                ax.figure.clf()
                 ax = ax.figure.add_subplot(111, projection='polar')
         ax.set_theta_direction(-1)
         ax.set_theta_offset(np.pi / 2.0)
@@ -333,7 +335,8 @@ class UniformSite(Site):
         dh_ijl = np.zeros_like(dw_ijl)
         dh_ijl[:, :, :] = dh_ij[:, :, na]
 
-        dw_order_indices_l = np.argsort(-dw_il, 0).astype(np.int).T
+        # dw_order_indices_l = np.argsort(-dw_il, 0).astype(np.int).T
+        dw_order_indices_l = np.argsort((dw_ijl > 0).sum(0), 0).T
 
         return dw_ijl, hcw_ijl, dh_ijl, dw_order_indices_l
 

py_wake/tests/test_aep_calculator.py

@@ -9,6 +9,7 @@ from py_wake.examples.data import hornsrev1
 from py_wake.wake_models.gaussian import IEA37SimpleBastankhahGaussian
 import numpy as np
 from py_wake.aep_calculator import AEPCalculator
+from py_wake.turbulence_models.stf import NOJ_STF
 
 
 def test_aep_no_wake():
@@ -27,16 +28,18 @@ def test_wake_map():
     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])
-
+    m = 49, slice(100, 133, 2)
+    # print(np.round(Z[m], 2).tolist()) # ref
     if 0:
         import matplotlib.pyplot as plt
         c = plt.contourf(X, Y, Z)  # , np.arange(2, 10, .01))
         plt.colorbar(c)
         windTurbines.plot(x, y)
+        plt.plot(X[m], Y[m], '.-r')
         plt.show()
 
     ref = [3.27, 3.27, 9.0, 7.46, 7.46, 7.46, 7.46, 7.31, 7.31, 7.31, 7.31, 8.3, 8.3, 8.3, 8.3, 8.3, 8.3]
-    npt.assert_array_almost_equal(Z[49, 100:133:2], ref, 2)
+    npt.assert_array_almost_equal(Z[m], ref, 2)
 
 
 def test_aep_map():
@@ -50,16 +53,41 @@ def test_aep_map():
     x_j = np.arange(-150, 150, 20)
     y_j = np.arange(-250, 250, 20)
     X, Y, Z = aep.aep_map(x_j, y_j, 0, x, y, wd=[0], ws=np.arange(3, 25))
-    # print(y_j)
+    m = 17
     if 0:
         import matplotlib.pyplot as plt
         c = plt.contourf(X, Y, Z, 100)  # , np.arange(2, 10, .01))
         plt.colorbar(c)
         windTurbines.plot(x, y)
+        plt.plot(X[m], Y[m], '.-r')
         plt.show()
-    # print(np.round(Z[17], 2).tolist())
+    # print(np.round(Z[m], 2).tolist()) # ref
     ref = [21.5, 21.4, 21.02, 20.34, 18.95, 16.54, 13.17, 10.17, 10.17, 13.17, 16.54, 18.95, 20.34, 21.02, 21.4]
-    npt.assert_array_almost_equal(Z[17], ref, 2)
+    npt.assert_array_almost_equal(Z[m], ref, 2)
+
+
+def test_ti_map():
+    site = IEA37Site(16)
+    x, y = site.initial_position.T
+    windTurbines = IEA37_WindTurbines()
+    wake_model = NOJ_STF(windTurbines)
+    aep = AEPCalculator(site, windTurbines, wake_model)
+    x_j = np.linspace(-1500, 1500, 200)
+    y_j = np.linspace(-1500, 1500, 100)
+    X, Y, Z = aep.ti_map(x_j, y_j, 110, x, y, wd=[0], ws=[9])
+    m = 49, slice(100, 133, 2)
+
+    # print(np.round(Z[m], 2).tolist())  # ref
+    if 0:
+        import matplotlib.pyplot as plt
+        c = plt.contourf(X, Y, Z, np.arange(.075, .50, .001))
+        plt.colorbar(c)
+        windTurbines.plot(x, y)
+        plt.plot(X[m], Y[m], '.-r')
+        plt.show()
+
+    ref = [0.43, 0.08, 0.08, 0.12, 0.15, 0.16, 0.18, 0.17, 0.16, 0.14, 0.12, 0.11, 0.12, 0.12, 0.12, 0.11, 0.11]
+    npt.assert_array_almost_equal(Z[m], ref, 2)
 
 
 def test_aep_map_no_turbines():
@@ -68,7 +96,7 @@ def test_aep_map_no_turbines():
     # n_wt = 16
     # x, y, _ = read_iea37_windfarm(iea37_path + 'iea37-ex%d.yaml' % n_wt)
 
-    site = UniformSite(freq, ti=0.75)
+    site = UniformSite(freq, ti=0.075)
     windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
     wake_model = IEA37SimpleBastankhahGaussian(windTurbines)
     aep = AEPCalculator(site, windTurbines, wake_model)

py_wake/tests/test_wake_model.py

@@ -16,7 +16,7 @@ from py_wake.tests.test_files.fuga import LUT_path_2MW_z0_0_03
 
 def test_wake_model():
     _, _, freq = read_iea37_windrose(iea37_path + "iea37-windrose.yaml")
-    site = UniformSite(freq, ti=0.75)
+    site = UniformSite(freq, ti=0.075)
     windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
     wake_model = NOJ(windTurbines)
     aep = AEPCalculator(site, windTurbines, wake_model)
@@ -30,7 +30,7 @@ def test_wec():
     wt_y = [433, 300, 0, 0, 0, -433, -433]
     wts = HornsrevV80()
 
-    site = UniformSite([1, 0, 0, 0], ti=0.75)
+    site = UniformSite([1, 0, 0, 0], ti=0.075)
 
     wake_model = Fuga(LUT_path_2MW_z0_0_03, wts)
     aep = AEPCalculator(site, wts, wake_model)

py_wake/tests/test_wake_models/test_fuga.py

@@ -14,7 +14,7 @@ def test_fuga():
     wts = HornsrevV80()
 
     path = tfp + 'fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+0/'
-    site = UniformSite([1, 0, 0, 0], ti=0.75)
+    site = UniformSite([1, 0, 0, 0], ti=0.075)
 
     wake_model = Fuga(path, wts)
     aep = AEPCalculator(site, wts, wake_model)
@@ -67,7 +67,7 @@ def cmp_fuga_with_colonel():
 
     path = tfp + 'fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+0/'
 
-    site = UniformSite([1, 0, 0, 0], ti=0.75)
+    site = UniformSite([1, 0, 0, 0], ti=0.075)
 
     wake_model = Fuga(path, wts)
     aep = AEPCalculator(site, wts, wake_model)

py_wake/turbulence_model.py

+from py_wake.wake_model import WakeModel
+import numpy as np
+
+
+class TurbulenceModel(WakeModel):
+    def calc_added_turbulence(self):
+        """Calculate added turbulence intensity caused by the x'th most upstream wind turbines
+        for all wind directions(l) and wind speeds(k) on a set of points(j)
+
+        This method must be overridden by subclass
+
+        Arguments required by this method must be added to the class list
+        args4addturb
+
+        See class documentation for examples and available arguments
+
+        Returns
+        -------
+        add_turb_jlk : array_like
+        """
+        pass
+
+    def calc_effective_TI(self, TI_lk, add_turb_jlk):
+        """Calculate effective turbulence intensity
+
+        Parameters
+        ----------
+        TI_lk : array_like
+            Local turbulence intensity at x'th most upstream turbines for all wind
+            directions(l) and wind speeds(k)
+        add_turb_jlk : array_like
+            deficit caused by upstream turbines(j) for all wind directions(l)
+            and wind speeds(k)
+
+        Returns
+        -------
+        TI_eff_lk : array_like
+            Effective wind speed at the x'th most upstream turbines for all wind
+            directions(l) and wind speeds(k)
+
+        """
+        pass
+
+    def ti_map(self, WS_ilk, WS_eff_ilk, TI_ilk, TI_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, TI_jlk):
+        """Calculate a wake (effecitve wind speed) map
+
+        Parameters
+        ----------
+        WS_ilk : array_like
+            Local wind speed [m/s] for each turbine(i), wind direction(l) and
+            wind speed(k)
+        WS_eff_ilk : array_like
+            Local effective wind speed [m/s] for each turbine(i), wind
+            direction(l) and wind speed(k)
+        dw_ijl : array_like
+            Down wind distance matrix between turbines(i) and map points (j) for
+            all wind directions(l) [m]
+        hcw_ijl : array_like
+            Horizontal cross wind distance matrix between turbines(i) and map
+            points(j) for all wind directions(l) [m]
+        dh_ijl : array_like
+            Vertical hub height distance matrix between turbines(i,i) for all
+            wind directions(l) [m]
+        ct_ilk : array_like
+            Thrust coefficient for all turbine(i), wind direction(l) and
+            wind speed(k)
+        types_i : array_like
+            Wind turbine type indexes
+        WS_jlk : array_like
+            Local wind speed [m/s] for map point(j), wind direction(l) and
+            wind speed(k)
+
+        Returns
+        -------
+        WS_eff_jlk : array_like
+            Local effective wind speed [m/s] for all map points(j),
+            wind direction(l) and wind speed(k)
+        """
+        return self._map(self.args4addturb, self.calc_added_turbulence, self.calc_effective_TI, WS_ilk, WS_eff_ilk, TI_ilk, TI_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, TI_jlk)
+
+
+class MaxSum():
+    def calc_effective_TI(self, TI_jlk, add_turb_ijlk):
+        return TI_jlk + np.max(add_turb_ijlk, 0)

py_wake/turbulence_models/stf.py

+from py_wake.turbulence_model import TurbulenceModel, MaxSum
+import numpy as np
+from numpy import newaxis as na
+from py_wake.wake_models.noj import NOJ
+
+
+class STFTurbulenceModel(MaxSum, TurbulenceModel):
+    args4addturb = ['dw_jl', 'cw_jl', 'D_src_l', 'ct_lk', 'TI_lk']
+
+    def calc_added_turbulence(self, dw_jl, cw_jl, D_src_l, ct_lk, TI_lk):
+        """ Calculate the added turbulence intensity at locations specified by
+        downstream distances (dw_jl) and crosswind distances (cw_jl)
+        caused by the wake of a turbine (diameter: D_src_l, thrust coefficient: Ct_lk).
+
+        Returns
+        -------
+        TI_eff_jlk: array:float
+            Effective turbulence intensity [-]
+        """
+        s_jl = dw_jl / D_src_l
+        # In the standard (see page 78), the maximal added TI is calculated as
+        # TI_add = 0.9/(1.5 + 0.3*d*sqrt(V_hub/c))
+        # where d is the downwind distance normalised by rotor diameter, c=1.0m/s
+        # Here it is assumed Ct = 7/V_hub (see Eq. (3.12) of ST Frandsen's thesis)
+        # thus, when using the acutal Ct, the function will be tranformed into:
+        # TI_add = 0.9/(1.5 + 0.8*d/sqrt(Ct))
+
+        TI_add_jlk = 0.9 / (1.5 + 0.8 * (dw_jl / D_src_l[na])[:, :, na] / np.sqrt(ct_lk)[na])
+
+        with np.warnings.catch_warnings():
+            np.warnings.filterwarnings('ignore', r'divide by zero encountered in true_divide')
+
+            # Theta_w is the characteristic view angle defined in Eq. (3.18) of
+            # ST Frandsen's thesis
+            theta_w = (180.0 / np.pi * np.arctan(1 / s_jl) + 10) / 2
+
+            # thetq denotes the acutally view angles
+            theta = np.arctan(cw_jl / dw_jl) * 180.0 / np.pi
+
+        # weights_jl = np.where(theta < 3 * theta_w, np.exp(-(theta / theta_w)**2), 0)
+        weights_jl = np.where(theta < theta_w, np.exp(-(theta / theta_w)**2), 0)
+
+        # the way effective added TI is calculated is derived from Eqs. (3.16-18)
+        # in ST Frandsen's thesis
+        TI_add_jlk = weights_jl[:, :, na] * (np.sqrt(TI_add_jlk**2 + TI_lk[na]**2) - TI_lk[na])
+        return TI_add_jlk
+
+
+class NOJ_STF(NOJ, STFTurbulenceModel):
+    pass
+
+
+def main():
+    if __name__ == '__main__':
+
+        from py_wake.aep_calculator import AEPCalculator
+        from py_wake.examples.data.iea37 import iea37_path
+        from py_wake.examples.data.iea37._iea37 import IEA37Site
+        from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
+        from py_wake.wake_models.noj import NOJ
+        # setup site, turbines and wakemodel
+        site = IEA37Site(16)
+        x, y = site.initial_position.T
+        windTurbines = IEA37_WindTurbines()
+
+        class NOJ_STF(NOJ, STFTurbulenceModel):
+            pass
+
+        wake_model = NOJ_STF(windTurbines)
+
+        # calculate AEP
+        aep_calculator = AEPCalculator(site, windTurbines, wake_model)
+        aep = aep_calculator.calculate_AEP(x, y)[0].sum()
+
+        # plot wake mape
+        import matplotlib.pyplot as plt
+
+        X, Y, Z = aep_calculator.ti_map(wt_x=x, wt_y=y, wd=[0], ws=[9])
+        c = plt.contourf(X, Y, Z, levels=100, cmap='Blues')
+        plt.colorbar(c, label='turbulence intensity [m/s]')
+        plt.title('AEP: %.2f GWh' % aep)
+        windTurbines.plot(x, y)
+        plt.show()
+
+
+main()

py_wake/wake_model.py

@@ -45,7 +45,9 @@ class WakeModel(ABC):
     Arguments available for calc_deficit (specifiy in args4deficit):
 
     - WS_lk: Local wind speed without wake effects
+    - TI_lk: local turbulence intensity without wake effects
     - WS_eff_lk: Local wind speed with wake effects
+    - WS_eff_lk: local turbulence intensity with wake effects
     - D_src_l: Diameter of source turbine
     - D_dst_jl: Diameter of destination turbine
     - dw_jl: Downwind distance from turbine i to point/turbine j
@@ -140,10 +142,14 @@ class WakeModel(ABC):
 
         K = WS_ilk.shape[2]
         deficit_nk = np.zeros((I * I * L, K))
+        if hasattr(self, 'calc_added_turbulence'):
+            add_turb_nk = np.zeros((I * I * L, K))
 
         indices = np.arange(I * I * L).reshape((I, I, L))
         WS_mk = WS_ilk.astype(np.float).reshape((I * L, K))
         WS_eff_mk = WS_mk.copy()
+        TI_mk = TI_ilk.astype(np.float).reshape((I * L, K))
+        TI_eff_mk = TI_mk.copy()
         dw_n = dw_iil.flatten()
         hcw_n = hcw_iil.flatten()
         if self.wec != 1:
@@ -172,6 +178,8 @@ class WakeModel(ABC):
             if j < I - 1:
                 arg_funcs = {'WS_lk': lambda: WS_mk[m],
                              'WS_eff_lk': lambda: WS_eff_mk[m],
+                             'TI_lk': lambda: TI_mk[m],
+                             'TI_eff_lk': lambda: TI_eff_mk[m],
                              'D_src_l': lambda: D_i[i_wt_l],
                              'D_dst_jl': lambda: D_i[dw_order_indices_dl[:, j + 1:]].T,
                              'dw_jl': lambda: dw_n[n_dw],
@@ -182,12 +190,50 @@ class WakeModel(ABC):
                              'ct_lk': lambda: ct_lk}
                 args = {k: arg_funcs[k]() for k in self.args4deficit}
                 deficit_nk[n_dw] = self.calc_deficit(**args)
+                if hasattr(self, 'calc_added_turbulence'):
+                    args = {k: arg_funcs[k]() for k in self.args4addturb}
+                    add_turb_nk[n_dw] = self.calc_added_turbulence(**args)
 
         WS_eff_ilk = WS_eff_mk.reshape((I, L, K))
         TI_eff_ilk = TI_ilk  # Effective TI is not yet implemented
         return WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk, WD_ilk, WS_ilk, TI_ilk, P_ilk
 
-    def wake_map(self, WS_ilk, WS_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, WS_jlk):
+    def _map(self, args4func, func, sum_func, WS_ilk, WS_eff_ilk, TI_ilk, TI_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, TI_jlk):
+        D_i = self.windTurbines.diameter(types_i)
+        H_i = self.windTurbines.hub_height(types_i)
+        if self.wec != 1:
+            hcw_ijl = hcw_ijl / self.wec
+        I, J, L = dw_ijl.shape
+        K = WS_ilk.shape[2]
+
+        deficit_ijlk = []
+        for i in range(I):
+            deficit_jlk = np.zeros((J, L, K))
+
+            for l in range(L):
+                m = dw_ijl[i, :, l] > 0
+
+                arg_funcs = {'WS_lk': lambda: WS_ilk[i, l][na],
+                             'WS_eff_lk': lambda: WS_eff_ilk[i, l][na],
+                             'TI_lk': lambda: TI_ilk[i, l][na],
+                             'TI_eff_lk': lambda: TI_eff_ilk[i, l][na],
+                             'D_src_l': lambda: D_i[i][na],
+                             'D_dst_jl': lambda: None,
+                             'dw_jl': lambda: dw_ijl[i, :, l][m][:, na],
+                             'cw_jl': lambda: np.sqrt(hcw_ijl[i, :, l][m]**2 + dh_ijl[i, :, l][m]**2)[:, na],
+                             'hcw_jl': lambda: hcw_ijl[i, :, l][m][:, na],
+                             'dh_jl': lambda: dh_ijl[i, :, l][m][:, na],
+                             'h_l': lambda: H_i[i][na],
+                             'ct_lk': lambda: ct_ilk[i, l][na]}
+
+                args = {k: arg_funcs[k]() for k in args4func}
+                deficit_jlk[:, l][m] = func(**args)[:, 0]
+
+            deficit_ijlk.append(deficit_jlk)
+        deficit_ijlk = np.array(deficit_ijlk)
+        return sum_func(TI_jlk, deficit_ijlk)
+
+    def ws_map(self, WS_ilk, WS_eff_ilk, TI_ilk, TI_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, WS_jlk):
         """Calculate a wake (effecitve wind speed) map
 
         Parameters
@@ -198,6 +244,12 @@ class WakeModel(ABC):
         WS_eff_ilk : array_like
             Local effective wind speed [m/s] for each turbine(i), wind
             direction(l) and wind speed(k)
+        TI_ilk : array_like
+            Local turbulence intensity for each turbine(i), wind direction(l) and
+            wind speed(k)
+        TI_eff_ilk : array_like
+            Local effective turbulence intensity for each turbine(i), wind
+            direction(l) and wind speed(k)
         dw_ijl : array_like
             Down wind distance matrix between turbines(i) and map points (j) for
             all wind directions(l) [m]
@@ -222,38 +274,8 @@ class WakeModel(ABC):
             Local effective wind speed [m/s] for all map points(j),
             wind direction(l) and wind speed(k)
         """
-        D_i = self.windTurbines.diameter(types_i)
-        H_i = self.windTurbines.hub_height(types_i)
-        if self.wec != 1:
-            hcw_ijl = hcw_ijl / self.wec
-        I, J, L = dw_ijl.shape
-        K = WS_ilk.shape[2]
-
-        deficit_ijlk = []
-        for i in range(I):
-            deficit_jlk = np.zeros((J, L, K))
-
-            for l in range(L):
-                m = dw_ijl[i, :, l] > 0
-
-                arg_funcs = {'WS_lk': lambda: WS_ilk[i, l][na],
-                             'WS_eff_lk': lambda: WS_eff_ilk[i, l][na],
-                             'D_src_l': lambda: D_i[i][na],
-                             'D_dst_jl': lambda: None,
-                             'dw_jl': lambda: dw_ijl[i, :, l][m][:, na],
-                             'cw_jl': lambda: np.sqrt(hcw_ijl[i, :, l][m]**2 + dh_ijl[i, :, l][m]**2)[:, na],
-                             'hcw_jl': lambda: hcw_ijl[i, :, l][m][:, na],
-                             'dh_jl': lambda: dh_ijl[i, :, l][m][:, na],
-                             'h_l': lambda: H_i[i][na],
-                             'ct_lk': lambda: ct_ilk[i, l][na]}
-
-                args = {k: arg_funcs[k]() for k in self.args4deficit}
-
-                deficit_jlk[:, l][m] = self.calc_deficit(**args)[:, 0]
-
-            deficit_ijlk.append(deficit_jlk)
-        deficit_ijlk = np.array(deficit_ijlk)
-        return self.calc_effective_WS(WS_jlk, deficit_ijlk)
+        return self._map(self.args4deficit, self.calc_deficit, self.calc_effective_WS, WS_ilk, WS_eff_ilk, TI_ilk,
+                         TI_eff_ilk, dw_ijl, hcw_ijl, dh_ijl, ct_ilk, types_i, WS_jlk)
 
     @abstractmethod
     def calc_deficit(self):
@@ -274,7 +296,7 @@ class WakeModel(ABC):
         pass
 
     @abstractmethod
-    def calc_effective_WS(self, WS_lk, deficit_jlk):
+    def calc_effective_WS(self, WS_jlk, deficit_ijlk):
         """Calculate effective wind speed
 
         This method must be overridden by subclass or by adding SquaredSum or
@@ -282,11 +304,11 @@ class WakeModel(ABC):
 
         Parameters
         ----------
-        WS_lk : array_like
-            Local wind speed at x'th most upstream turbines for all wind
+        WS_jlk : array_like
+            Local wind speed at turbine/site(j) for all wind
             directions(l) and wind speeds(k)
-        deficit_jlk : array_like
-            deficit caused by upstream turbines(j) for all wind directions(l)
+        deficit_ijlk : array_like
+            deficit caused by upstream turbines(i) on all downstream turbines/points (j) for all wind directions(l)
             and wind speeds(k)
 
         Returns
@@ -304,13 +326,13 @@ class WakeModel(ABC):
 
 
 class SquaredSum():
-    def calc_effective_WS(self, WS_lk, deficit_jlk):
-        return WS_lk - np.sqrt(np.sum(deficit_jlk**2, 0))
+    def calc_effective_WS(self, WS_jlk, deficit_ijlk):
+        return WS_jlk - np.sqrt(np.sum(deficit_ijlk**2, 0))
 
 
 class LinearSum():
-    def calc_effective_WS(self, WS_lk, deficit_jlk):
-        return WS_lk - np.sum(deficit_jlk, 0)
+    def calc_effective_WS(self, WS_jlk, deficit_ijlk):
+        return WS_jlk - np.sum(deficit_ijlk, 0)
 
 
 def main():
@@ -333,7 +355,7 @@ def main():
         n_wt = 16
         x, y, _ = read_iea37_windfarm(iea37_path + 'iea37-ex%d.yaml' % n_wt)
 
-        site = UniformSite(freq, ti=0.75)
+        site = UniformSite(freq, ti=0.075)
         windTurbines = IEA37_WindTurbines(iea37_path + 'iea37-335mw.yaml')
 
         wake_model = MyWakeModel(windTurbines)

py_wake/wake_models/noj.py

@@ -3,35 +3,36 @@ import numpy as np
 from numpy import newaxis as na
 
 
-class NOJ(SquaredSum, WakeModel):
-    args4deficit = ['WS_lk', 'D_src_l', 'D_dst_jl', 'dw_jl', 'cw_jl', 'ct_lk']
-
-    def __init__(self, windTurbines, k=.1, **kwargs):
-        super().__init__(windTurbines, **kwargs)
+class AreaOverlappingFactor():
+    def __init__(self, k=.1):
         self.k = k
 
-    def calc_deficit(self, WS_lk, D_src_l, D_dst_jl, dw_jl, cw_jl, ct_lk):
+    def overlapping_area_factor(self, dw_jl, cw_jl, D_src_l, D_dst_jl):
+        """Calculate overlapping factor
 
-        # Calculate the wake loss using NOJ
-        # Jensen, Niels Otto. "A note on wind generator interaction." (1983)
-
-        # In NOJensen wake model:
-        #      V_def = v*(1-sqrt(1-Ct))/(1+k*dist_down/R)**2
-        R_l = D_src_l / 2
-        term_denominator_jl = (1 + self.k * dw_jl / R_l[na, :])**2
+        Parameters
+        ----------
+        dw_jl : array_like
+            down wind distance [m]
+        cw_jl : array_like
+            cross wind distance [m]
+        D_src_l : array_like
+            Diameter of source turbines [m]
+        D_dst_jl : array_like or None
+            Diameter of destination turbines [m]. If None destination is assumed to be a point
 
-        wake_radius_jl = (self.k * dw_jl + R_l[na, :])
+        Returns
+        -------
+        A_ol_factor_jl : array_like
+            area overlaping factor
+        """
+        wake_radius_jl = (self.k * dw_jl + D_src_l[na, :] / 2)
         if D_dst_jl is None:
-            A_ol_factor_jl = wake_radius_jl > cw_jl
+            return wake_radius_jl > cw_jl
         else:
-            A_ol_factor_jl = self.cal_overlapping_area_factor(wake_radius_jl,
-                                                              (D_dst_jl / 2),
-                                                              cw_jl)
-        term_numerator_lk = WS_lk * (1 - np.sqrt(1 - ct_lk))
-        with np.warnings.catch_warnings():
-            np.warnings.filterwarnings('ignore', r'invalid value encountered in true_divide')
-            deficit_jlk = term_numerator_lk[na] * (A_ol_factor_jl / term_denominator_jl)[:, :, na]
-        return deficit_jlk
+            return self.cal_overlapping_area_factor(wake_radius_jl,
+                                                    (D_dst_jl / 2),
+                                                    cw_jl)
 
     def cal_overlapping_area_factor(self, R1, R2, d):
         """ Calculate the overlapping area of two circles with radius R1 and
@@ -91,6 +92,31 @@ class NOJ(SquaredSum, WakeModel):
         return A_ol_f
 
 
+class NOJ(SquaredSum, WakeModel, AreaOverlappingFactor):
+    args4deficit = ['WS_lk', 'D_src_l', 'D_dst_jl', 'dw_jl', 'cw_jl', 'ct_lk']
+
+    def __init__(self, windTurbines, k=.1, **kwargs):
+        WakeModel.__init__(self, windTurbines, **kwargs)
+        AreaOverlappingFactor.__init__(self, k)
+
+    def calc_deficit(self, WS_lk, D_src_l, D_dst_jl, dw_jl, cw_jl, ct_lk):
+
+        # Calculate the wake loss using NOJ
+        # Jensen, Niels Otto. "A note on wind generator interaction." (1983)
+        # V_def = v*(1-sqrt(1-Ct))/(1+k*dist_down/R)**2
+
+        R_src_l = D_src_l / 2
+
+        term_denominator_jl = (1 + self.k * dw_jl / R_src_l[na, :])**2
+
+        A_ol_factor_jl = self.overlapping_area_factor(dw_jl, cw_jl, D_src_l, D_dst_jl)
+        term_numerator_lk = WS_lk * (1 - np.sqrt(1 - ct_lk))
+        with np.warnings.catch_warnings():
+            np.warnings.filterwarnings('ignore', r'invalid value encountered in true_divide')
+            deficit_jlk = term_numerator_lk[na] * (A_ol_factor_jl / term_denominator_jl)[:, :, na]
+        return deficit_jlk
+
+
 def main():
     if __name__ == '__main__':
 

py_wake/wind_turbines.py

@@ -87,8 +87,9 @@ class WindTurbines():
 
     def _ct_power(self, ws_i, type_i=0):
         if np.any(type_i != 0):
-            CT = np.zeros_like(ws_i)
-            P = np.zeros_like(ws_i)
+            CT = np.zeros_like(ws_i, dtype=np.float)
+            P = np.zeros_like(ws_i, dtype=np.float)
+            type_i = np.zeros(np.asarray(ws_i).shape[0]) + type_i
             for t in np.unique(type_i).astype(np.int):
                 m = type_i == t
                 CT[m] = self.ct_funcs[t](ws_i[m])
@@ -170,7 +171,7 @@ def cube_power(ws_cut_in=3, ws_cut_out=25, ws_rated=12, power_rated=5000):
     return power_func
 
 
-def dummy_thrust(ws_cut_in=3, ws_cut_out=25, ws_rated=12, ct_rated=0.88):
+def dummy_thrust(ws_cut_in=3, ws_cut_out=25, ws_rated=12, ct_rated=8 / 9):
     # temporary thrust curve fix
     def ct_func(ws):
         ws = np.asarray(ws)
@@ -181,7 +182,8 @@ def dummy_thrust(ws_cut_in=3, ws_cut_out=25, ws_rated=12, ct_rated=0.88):
             ct[m] = ct_rated
             idx = (ws >= ws_rated) & (ws <= ws_cut_out)
             # second order polynomial fit for above rated
-            ct[idx] = np.polyval(np.polyfit([ws_rated, (ws_rated + ws_cut_out) / 2, ws_cut_out], [ct_rated, 0.4, 0.03], 2), ws[idx])
+            ct[idx] = np.polyval(np.polyfit([ws_rated, (ws_rated + ws_cut_out) / 2,
+                                             ws_cut_out], [ct_rated, 0.4, 0.03], 2), ws[idx])
         return ct
     return ct_func
 
@@ -192,16 +194,20 @@ def main():
         wts = WindTurbines(names=['tb1', 'tb2'],
                            diameters=[80, 120],
                            hub_heights=[70, 110],
-                           ct_funcs=[lambda _: 8 / 9,
-                                     lambda _: 8 / 9],
+                           ct_funcs=[lambda ws: ws * 0 + 8 / 9,
+                                     dummy_thrust()],
                            power_funcs=[cube_power(ws_cut_in=3, ws_cut_out=25, ws_rated=12, power_rated=2000),
                                         cube_power(ws_cut_in=3, ws_cut_out=25, ws_rated=12, power_rated=3000)],
                            power_unit='kW')
 
         ws = np.arange(25)
         import matplotlib.pyplot as plt
-        power = wts.power(ws)
-        plt.plot(ws, power, label=wts.name(0))
+        plt.plot(ws, wts.power(ws, 0), label=wts.name(0))
+        plt.plot(ws, wts.power(ws, 1), label=wts.name(1))
+        plt.legend()
+        plt.show()
+        plt.plot(ws, wts.ct(ws, 0), label=wts.name(0))
+        plt.plot(ws, wts.ct(ws, 1), label=wts.name(1))
         plt.legend()
         plt.show()