Removed some redefinitions of the axial induction factor

py_wake/deficit_models/rankinehalfbody.py

@@ -4,6 +4,7 @@ from py_wake.deficit_models import DeficitModel
 from py_wake.deficit_models import BlockageDeficitModel
 from py_wake.ground_models.ground_models import NoGround
 from py_wake.utils.gradients import hypot
+from py_wake.deficit_models.utils import a0
 
 
 class RankineHalfBody(BlockageDeficitModel):
@@ -22,23 +23,12 @@ class RankineHalfBody(BlockageDeficitModel):
                  upstream_only=False):
         DeficitModel.__init__(self, groundModel=groundModel)
         BlockageDeficitModel.__init__(self, upstream_only=upstream_only, superpositionModel=superpositionModel)
-        # coefficients for BEM approximation by Madsen (1997)
-        self.a0p = np.array([0.2460, 0.0586, 0.0883])
         # limiter to avoid singularities
         self.limiter = limiter
         # if used in a wind farm simulation, set deficit in wake region to
         # zero, as here the wake model is active
         self.exclude_wake = exclude_wake
 
-    def a0(self, ct_ilk):
-        """
-        BEM axial induction approximation by Madsen (1997).
-        """
-        # Evaluate with Horner's rule.
-        # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk
-        a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2]))
-        return a0_ilk
-
     def outside_body(self, WS_ilk, a0_ilk, R_il, dw_ijlk, cw_ijlk, r_ijlk):
         """
         Find all points lying outside Rankine Half Body, stagnation line given on p.3
@@ -54,7 +44,7 @@ class RankineHalfBody(BlockageDeficitModel):
 
     def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_):
         # source strength as given on p.7
-        a0_ilk = self.a0(ct_ilk)
+        a0_ilk = a0(ct_ilk)
         R_il = D_src_il / 2.
         m_ilk = 2. * WS_ilk * a0_ilk * np.pi * R_il[:, :, na]**2
         # radial distance
@@ -71,7 +61,7 @@ class RankineHalfBody(BlockageDeficitModel):
             deficit_ijlk = self.remove_wake(deficit_ijlk, dw_ijlk, cw_ijlk, D_src_il)
             # Close to the rotor the induced velocities become unphysical and are
             # limited to the induction in the rotor plane estimated by BEM.
-            induc = (WS_ilk * self.a0(ct_ilk))[:, na]
+            induc = (WS_ilk * a0(ct_ilk))[:, na]
             deficit_ijlk = np.where(deficit_ijlk > induc, induc * np.sign(deficit_ijlk), deficit_ijlk)
 
         return deficit_ijlk

py_wake/deficit_models/rathmann.py

@@ -4,6 +4,7 @@ from py_wake.deficit_models import DeficitModel
 from py_wake.deficit_models import BlockageDeficitModel
 from py_wake.ground_models.ground_models import NoGround
 from py_wake.utils.gradients import hypot
+from py_wake.deficit_models.utils import a0
 
 
 class Rathmann(BlockageDeficitModel):
@@ -29,8 +30,6 @@ class Rathmann(BlockageDeficitModel):
                  upstream_only=False):
         DeficitModel.__init__(self, groundModel=groundModel)
         BlockageDeficitModel.__init__(self, upstream_only=upstream_only, superpositionModel=superpositionModel)
-        # coefficients for BEM approximation by Madsen (1997)
-        self.a0p = np.array([0.2460, 0.0586, 0.0883])
         # limiter to avoid singularities
         self.limiter = limiter
         # coefficient for scaling the effective forcing
@@ -53,15 +52,6 @@ class Rathmann(BlockageDeficitModel):
         # layout term
         self.dmu_G_ijlk = dmu_ijlk * G_ijlk
 
-    def a0(self, ct_ilk):
-        """
-        BEM axial induction approximation by Madsen (1997).
-        """
-        # Evaluate with Horner's rule.
-        # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk
-        a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2]))
-        return a0_ilk
-
     def dmu(self, xi_ijlk):
         """
         Centreline deficit shape function. Same as for the vortex cylinder model.
@@ -96,7 +86,7 @@ class Rathmann(BlockageDeficitModel):
             self._calc_layout_terms(D_src_il, dw_ijlk, cw_ijlk)
 
         # circulation/strength of vortex dipole Eq. (1) in [1]
-        gammat_ilk = WS_ilk * 2. * self.a0(ct_ilk * self.sct)
+        gammat_ilk = WS_ilk * 2. * a0(ct_ilk * self.sct)
 
         deficit_ijlk = gammat_ilk[:, na] / 2. * self.dmu_G_ijlk
         # turn deficit into speed-up downstream

py_wake/deficit_models/utils.py

+import numpy as np
+
+
+def a0(ct_ilk, a0p=np.array([0.2460, 0.0586, 0.0883])):
+    """
+    BEM axial induction approximation by
+    Madsen, H. A., Larsen, T. J., Pirrung, G. R., Li, A., and Zahle, F.: Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact, Wind Energ. Sci., 5, 1–27, https://doi.org/10.5194/wes-5-1-2020, 2020.
+    """
+    # Evaluate with Horner's rule.
+    # a0_ilk = a0p[2] * ct_ilk**3 + a0p[1] * ct_ilk**2 + a0p[0] * ct_ilk
+    a0_ilk = ct_ilk * (a0p[0] + ct_ilk * (a0p[1] + ct_ilk * a0p[2]))
+    return a0_ilk

py_wake/deficit_models/vortexcylinder.py

@@ -6,6 +6,7 @@ from py_wake.utils.elliptic import ellipticPiCarlson
 from py_wake.deficit_models import DeficitModel
 from py_wake.deficit_models import BlockageDeficitModel
 from py_wake.utils.gradients import hypot, cabs
+from py_wake.deficit_models.utils import a0
 
 
 class VortexCylinder(BlockageDeficitModel):
@@ -27,8 +28,6 @@ class VortexCylinder(BlockageDeficitModel):
                  upstream_only=False):
         DeficitModel.__init__(self, groundModel=groundModel)
         BlockageDeficitModel.__init__(self, upstream_only=upstream_only, superpositionModel=superpositionModel)
-        # coefficients for BEM approximation by Madsen (1997)
-        self.a0p = np.array([0.2460, 0.0586, 0.0883])
         # limiter to avoid singularities
         self.limiter = limiter
         # if used in a wind farm simulation, set deficit in wake region to
@@ -72,15 +71,6 @@ class VortexCylinder(BlockageDeficitModel):
         # deficit shape function
         self.dmu_ijlk = term1_ijlk + term2_ijlk
 
-    def a0(self, ct_ilk):
-        """
-        BEM axial induction approximation by Madsen (1997).
-        """
-        # Evaluate with Horner's rule.
-        # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk
-        a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2]))
-        return a0_ilk
-
     def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_):
         """
         The analytical relationships can be found in [1,2], in particular equations (7-8) from [1].
@@ -90,7 +80,7 @@ class VortexCylinder(BlockageDeficitModel):
             self._calc_layout_terms(D_src_il, dw_ijlk, cw_ijlk)
 
         # circulation/strength of vortex cylinder
-        gammat_ilk = WS_ilk * 2. * self.a0(ct_ilk)
+        gammat_ilk = WS_ilk * 2. * a0(ct_ilk)
 
         deficit_ijlk = gammat_ilk[:, na] / 2. * self.dmu_ijlk
 

py_wake/deficit_models/vortexdipole.py

@@ -4,6 +4,7 @@ from py_wake.deficit_models import DeficitModel
 from py_wake.deficit_models import BlockageDeficitModel
 from py_wake.ground_models.ground_models import NoGround
 from py_wake.utils.gradients import hypot, cabs
+from py_wake.deficit_models.utils import a0
 
 
 class VortexDipole(BlockageDeficitModel):
@@ -27,8 +28,6 @@ class VortexDipole(BlockageDeficitModel):
                  upstream_only=False):
         DeficitModel.__init__(self, groundModel=groundModel)
         BlockageDeficitModel.__init__(self, upstream_only=upstream_only, superpositionModel=superpositionModel)
-        # coefficients for BEM approximation by Madsen (1997)
-        self.a0p = np.array([0.2460, 0.0586, 0.0883])
         # limiter to avoid singularities
         self.limiter = limiter
         # coefficient for scaling the effective forcing
@@ -37,15 +36,6 @@ class VortexDipole(BlockageDeficitModel):
         # zero, as here the wake model is active
         self.exclude_wake = exclude_wake
 
-    def a0(self, ct_ilk):
-        """
-        BEM axial induction approximation by Madsen (1997).
-        """
-        # Evaluate with Horner's rule.
-        # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk
-        a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2]))
-        return a0_ilk
-
     def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_):
         """
         The analytical relationships can be found in [1,2].
@@ -54,7 +44,7 @@ class VortexDipole(BlockageDeficitModel):
         # radial distance
         r_ijlk = hypot(dw_ijlk, cw_ijlk)
         # circulation/strength of vortex dipole Eq. (1) in [1]
-        gammat_ilk = WS_ilk * 2. * self.a0(ct_ilk * self.sct)
+        gammat_ilk = WS_ilk * 2. * a0(ct_ilk * self.sct)
         # Eq. (2) in [1], induced velocities away from centreline, however
         # here it is simplified. Effectively the equations are the same as for
         # a Rankine Half Body.