Major revision of the code specifically address tip loss correction implementation

MITRotor/Aerodynamics.py

@@ -8,6 +8,7 @@
 
 from .RotorDefinition import RotorDefinition
 from .Geometry import BEMGeometry
+from .TipLoss import TipLossModel
 
 __all__ = [
     "AerodynamicModel",
@@ -35,7 +36,7 @@ class AerodynamicProperties:
         aoa (ArrayLike): Blade element angle of attack.
         Cl (ArrayLike): Blade element lift coefficient.
         Cd (ArrayLike): Blade element drag coefficient.
-        F (Optional[ArrayLike]): Blade element tip loss (optional).
+        F (ArrayLike): Blade element tip loss.
 
     Properties:
         W: Blade element inflow magnitude.
@@ -56,7 +57,7 @@ class AerodynamicProperties:
     aoa: ArrayLike
     Cl: ArrayLike
     Cd: ArrayLike
-    F: Optional[ArrayLike] = None
+    F: ArrayLike = None
 
     def __post_init__(self):
         pass
@@ -74,20 +75,50 @@ def phi(self):
         Blade element inflow direction.
         """
         return np.arctan2(self.Vax, self.Vtan)
+
+    @cached_property
+    def C_n(self):
+        """
+        Blade element axial blade force coefficient.
+        """
+        return self.Cl * np.cos(self.phi) + self.Cd * np.sin(self.phi)
 
     @cached_property
-    def Ctan(self):
+    def C_tan(self):
         """
-        Blade element tangengial force coefficient.
+        Blade element tangential blade force coefficient.
         """
         return self.Cl * np.sin(self.phi) - self.Cd * np.cos(self.phi)
+
+    @cached_property
+    def C_x(self):
+        """
+        Blade element axial area force coefficient.
+        """
+        return self.solidity * self.W**2 * self.C_n
 
     @cached_property
-    def Cax(self):
+    def C_tau(self):
         """
-        Blade element axial force coefficient.
+        Blade element tangential area force coefficient.
         """
-        return self.Cl * np.cos(self.phi) + self.Cd * np.sin(self.phi)
+        return self.solidity * self.W**2 * self.C_tan
+
+    @cached_property
+    def C_x_corr(self):
+        """
+        Corrected blade element area axial force coefficient.
+        """
+        return self.C_x / self.F
+
+    @cached_property
+    def C_tau_corr(self):
+        """
+        Corrected blade element area tangential force coefficient.
+        """
+        return self.C_tau / self.F
+
+
 
 
 class AerodynamicModel(ABC):
@@ -103,11 +134,10 @@ def __call__(
         geom: BEMGeometry,
         U: ArrayLike,
         wdir: ArrayLike,
+        # tiploss_model: TipLossModel
     ) -> AerodynamicProperties:
         """
-        Performs the aerodynamic calculations in a blade-element code using the
-        method outlined in Howland et al. 2020. (Influence of atmospheric conditions
-        on the power production of utility-scale wind turbines in yaw misalignment)
+        Performs the aerodynamic calculations in a blade-element code.
 
         Args:
             an (ArrayLike): Axial induction radial profile.
@@ -127,7 +157,7 @@ def __call__(
         ...
 
 
-class KraghAerodynamics(AerodynamicModel):
+class DefaultAerodynamics(AerodynamicModel):
     def __call__(
         self,
         an: ArrayLike,
@@ -139,6 +169,7 @@ def __call__(
         geom: BEMGeometry,
         U: ArrayLike,
         wdir: ArrayLike,
+        # tiploss_model: TipLossModel
     ) -> AerodynamicProperties:
         """
         Performs the aerodynamic calculations in a blade-element code using the
@@ -160,78 +191,40 @@ def __call__(
             AerodynamicProperties: Calculated aerodynamic properties stored in AerodynamicProperties object.
 
         """
-        local_yaw = wdir - yaw
-
-        Vax = U * ((1 - an) * np.cos(local_yaw * np.cos(geom.theta_mesh)) * np.cos(local_yaw * np.sin(geom.theta_mesh)))
-        Vtan = (1 + aprime) * tsr * geom.mu_mesh - U * (1 - an) * np.cos(local_yaw * np.sin(geom.theta_mesh)) * np.sin(
-            local_yaw * np.cos(geom.theta_mesh)
-        )
-
-        phi = np.arctan2(Vax, Vtan)
-        aoa = phi - rotor.twist(geom.mu_mesh) - pitch
-        aoa = np.clip(aoa, -np.pi / 2, np.pi / 2)
+        yaw_z = wdir - yaw
 
-        Cl, Cd = rotor.clcd(geom.mu_mesh, aoa)
-
-        solidity = rotor.solidity(geom.mu_mesh)
-
-        aero_props = AerodynamicProperties(
-            an, aprime, solidity, U * np.ones(geom.shape), wdir * np.ones(geom.shape), Vax, Vtan, aoa, Cl, Cd
+        Vax = (
+            U 
+            * (1 - an)
+            * np.cos(yaw_z * np.cos(geom.theta_mesh)) 
+            * np.cos(yaw_z * np.sin(geom.theta_mesh))
         )
 
-        return aero_props
-
-
-class DefaultAerodynamics(AerodynamicModel):
-    def __call__(
-        self,
-        an: ArrayLike,
-        aprime: ArrayLike,
-        pitch: float,
-        tsr: float,
-        yaw: float,
-        rotor: RotorDefinition,
-        geom: BEMGeometry,
-        U: ArrayLike,
-        wdir: ArrayLike,
-    ) -> AerodynamicProperties:
-        """
-        Performs the aerodynamic calculations in a blade-element code using the
-        method outlined in Howland et al. 2020. (Influence of atmospheric conditions
-        on the power production of utility-scale wind turbines in yaw misalignment)
-
-        Args:
-            an (ArrayLike): Axial induction radial profile.
-            aprime (ArrayLike): tangengial induction radial profile.
-            pitch (float): blade pitch angle [rad].
-            tsr (float): Rotor tip-speed ratio.
-            yaw (float): Rotor yaw angle [rad].
-            rotor (RotorDefinition): Turbine rotor definition object.
-            geom (BEMGeometry): Blade element geometry object.
-            U (ArrayLike): Inflow velocity on polar grid.
-            wdir (ArrayLike): Inflow direction on polar grid.
-
-        Returns:
-            AerodynamicProperties: Calculated aerodynamic properties stored in AerodynamicProperties object.
-
-        """
-        local_yaw = -yaw
-
-        Vax = U * ((1 - an) * np.cos(local_yaw))
-        Vtan = (1 + aprime) * tsr * geom.mu_mesh - U * np.cos(geom.theta_mesh) * np.sin(
-            local_yaw
+        Vtan = (
+            (1 + aprime) * tsr * geom.mu_mesh
+            - U * (1 - an) 
+            * np.cos(yaw_z * np.sin(geom.theta_mesh)) 
+            * np.sin(yaw_z * np.cos(geom.theta_mesh))
         )
 
         phi = np.arctan2(Vax, Vtan)
         aoa = phi - rotor.twist(geom.mu_mesh) - pitch
-        aoa = np.clip(aoa, -np.pi / 2, np.pi / 2)
 
         Cl, Cd = rotor.clcd(geom.mu_mesh, aoa)
 
         solidity = rotor.solidity(geom.mu_mesh)
 
         aero_props = AerodynamicProperties(
-            an, aprime, solidity, U * np.ones(geom.shape), wdir * np.ones(geom.shape), Vax, Vtan, aoa, Cl, Cd
+            an = an, 
+            aprime = aprime, 
+            solidity = solidity, 
+            U = U, 
+            wdir = wdir, 
+            Vax = Vax, 
+            Vtan = Vtan, 
+            aoa = aoa, 
+            Cl = Cl, 
+            Cd = Cd,
         )
 
-        return aero_props
+        return aero_props

MITRotor/BEMSolver.py

@@ -81,12 +81,21 @@ def Cl(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
 
     def Cd(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.Cd, grid)
-
-    def Cax(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        return average(self.geom, self.aero_props.Cax, grid)
-
-    def Ctan(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
-        return average(self.geom, self.aero_props.Ctan, grid)
+
+    def Cn(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_n, grid)
+
+    def Ctan(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_tan, grid)
+
+    def Cx(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_x_corr, grid)
+
+    def Ctau(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_tau_corr, grid)
+
+    def Ctau_uncorr(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_tau, grid)
 
     def F(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.F, grid)
@@ -96,28 +105,35 @@ def u4(self):
 
     def v4(self):
         return self._v4
-
+    
     def Cp(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         dCp = (
             self.tsr
-            * self.solidity(grid="sector")
             * self.geom.mu_mesh
-            * self.W(grid="sector") ** 2
-            * self.Ctan(grid="sector")
+            * self.Ctau_uncorr(grid="sector")
         )
         return average(self.geom, dCp, grid=grid)
-
+
+    def Cp_corr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        dCp = (
+            self.tsr
+            * self.geom.mu_mesh
+            * self.Ctau(grid="sector")
+        )
+        return average(self.geom, dCp, grid=grid)
+
     def Ct(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        _Ct = self.solidity(grid="sector") * self.W(grid="sector") ** 2 * self.Cax(grid="sector")
+        _Ct = self.aero_props.C_x
+        return average(self.geom, _Ct, grid=grid)
+
+    def Ct_corr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        _Ct = self.aero_props.C_x_corr
         return average(self.geom, _Ct, grid=grid)
 
     def Ctprime(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         Ctprime = self.Ct(grid="sector") / ((1 - self.a(grid="sector")) ** 2 * np.cos(self.yaw) ** 2)
         return average(self.geom, Ctprime, grid=grid)
 
-    def Cq(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        return average(self.geom, self.Cp(grid="sector") / self.tsr, grid=grid)
-
 
 @adaptivefixedpointiteration(max_iter=500, relaxations=[0.25, 0.5, 0.96])
 class BEM:
@@ -161,7 +177,7 @@ def sample_points(self, yaw: float = 0.0) -> tuple[ArrayLike, ArrayLike, ArrayLi
     def initial_guess(
         self, pitch: float, tsr: float, yaw: float = 0.0, U: ArrayLike = 1.0, wdir: ArrayLike = 0.0
     ) -> Tuple[ArrayLike, ...]:
-        a = 1 / 3 * np.ones(self.geometry.shape)
+        a = (1 / 3) * np.ones(self.geometry.shape)
         aprime = np.zeros(self.geometry.shape)
 
         return a, aprime
@@ -172,19 +188,33 @@ def residual(
         pitch: ArrayLike,
         tsr: ArrayLike,
         yaw: ArrayLike = 0.0,
-        U: ArrayLike = 1.0,
-        wdir: ArrayLike = 0.0,
+        U: ArrayLike = None,
+        wdir: ArrayLike = None,
     ) -> Tuple[ArrayLike, ...]:
         an, aprime = x
-
-        aero_props = self.aerodynamic_model(an, aprime, pitch, tsr, yaw, self.rotor, self.geometry, U, wdir)
+        U = U or np.ones(self.geometry.shape)
+        wdir = wdir or np.zeros(self.geometry.shape)
+
+        aero_props = self.aerodynamic_model(
+            an = an, 
+            aprime=aprime, 
+            pitch=pitch, 
+            tsr=tsr, 
+            yaw=yaw, 
+            rotor=self.rotor, 
+            geom=self.geometry, 
+            U=U, 
+            wdir=wdir)
+
         aero_props.F = self.tiploss_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry)
         e_an = self.momentum_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry) - an
         e_aprime = self.tangential_induction_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry) - aprime
 
         return e_an, e_aprime
 
-    def post_process(self, result: FixedPointIterationResult, pitch, tsr, yaw, U=1.0, wdir=0.0) -> BEMSolution:
+    def post_process(self, result: FixedPointIterationResult, pitch, tsr, yaw, U=None, wdir=None) -> BEMSolution:
+        U = U or np.ones(self.geometry.shape)
+        wdir = wdir or np.zeros(self.geometry.shape)
         an, aprime = result.x
         aero_props = self.aerodynamic_model(an, aprime, pitch, tsr, yaw, self.rotor, self.geometry, U, wdir)
         aero_props.F = self.tiploss_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry)