Add 3D stall delay correction, fix tangential induction model, add CachedLUT

MITRotor/Aerodynamics.py

@@ -122,6 +122,10 @@ def C_tau_corr(self):
 
 
 class AerodynamicModel(ABC):
+    def __init__(self, apply_3D_stall_correction: bool = False):
+        self.apply_3D_stall_correction = apply_3D_stall_correction
+
+
     @abstractmethod
     def __call__(
         self,
@@ -215,10 +219,10 @@ def __call__(
         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)
 
+        Cl, Cd = rotor.clcd(geom.mu_mesh, aoa)
+
         aero_props = AerodynamicProperties(
             an = an,
             aprime = aprime,
@@ -264,6 +268,8 @@ def __call__(
             geom (BEMGeometry): Blade element geometry object.
             U (ArrayLike): Inflow velocity on polar grid.
             wdir (ArrayLike): Inflow direction on polar grid.
+            tilt (float): Rotor tilt angle [rad].
+            apply_3D_stall_correction (bool): Whether to apply 3D stall correction to angle of attack.
 
         Returns:
             AerodynamicProperties: Calculated aerodynamic properties stored in AerodynamicProperties object.
@@ -283,10 +289,11 @@ def __call__(
         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)
 
+        Cl, Cd = rotor.clcd(geom.mu_mesh, aoa, tsr=tsr, apply_3D_stall_correction=self.apply_3D_stall_correction)
+
+
         aero_props = AerodynamicProperties(
             an = an,
             aprime = aprime,

MITRotor/BEMSolver.py

@@ -89,26 +89,21 @@ 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)
+        return average(self.geom, self.aero_props.C_x, 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 Cx_corr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_x_corr, grid)
+
+    def Ctau_corr(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+        return average(self.geom, self.aero_props.C_tau_corr, grid)
 
     def F(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.F, grid)
 
     def Cp(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        dCp = (
-            self.tsr
-            * self.geom.mu_mesh
-            * 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
@@ -119,18 +114,14 @@ def Cp_corr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
     def Ct(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         _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"):
         eff_yaw = calc_eff_yaw(self.yaw, self.tilt)
         Ctprime = self.Ct(grid="sector") / ((1 - self.a(grid="sector")) ** 2 * np.cos(eff_yaw) ** 2)
         return average(self.geom, Ctprime, grid=grid)
 
 
-@adaptivefixedpointiteration(max_iter=500, relaxations=[0.25, 0.5, 0.96])
+@adaptivefixedpointiteration(max_iter=1000, relaxations=[0.25, 0.5, 0.96])
 class BEM:
     """
     A generic BEM class which facilitates dependency injection for various models.

MITRotor/RotorDefinition.py

@@ -61,8 +61,68 @@ def Cl(self, x, inflow):
     def Cd(self, x, inflow):
         return self.cd_interp(x, inflow, grid=False)
 
-    def __call__(self, x, inflow):
-        return self.Cl(x, inflow), self.Cd(x, inflow)
+    def __call__(self, x, inflow, tsr=None, apply_3D_stall_correction=False, c_on_r=None):
+        '''
+        Inputs:
+            - x: radial position normalized by rotor radius (mu)
+            - inflow: angle of attack (radians)
+            - tsr: tip speed ratio (only needed if apply_3D_stall_correction=True)
+            - apply_3D_stall_correction: whether to apply 3D stall correction (Du and Selig 1998)
+            - c_on_r: chord-to-radius ratio (only needed if apply_3D_stall_correction=True)
+
+        Returns:
+            - Cl, Cd with optional 3D stall correction applied
+        '''
+        a = 1
+        b = 1
+        d = 1
+        Cl_2D = self.Cl(x, inflow)
+        Cd_2D = self.Cd(x, inflow)
+
+        if apply_3D_stall_correction:
+            # Apply 3D stall correction (Du and Selig 1998)
+
+            # Compute slope of Cl curve
+            CL1 = self.Cl(x, np.radians(0) * np.ones_like(x))
+            CL2 = self.Cl(x, np.radians(5) * np.ones_like(x))
+            m = (CL2 - CL1) / np.radians(5 - 0)
+
+            # Compute alpha_0 (angle of attack at which Cl=0)
+            alpha_0 = -CL1 / m
+
+            tsr_mod = 1 / np.sqrt(1 + tsr**2)
+
+            exp1 = np.minimum((d/(tsr_mod * x + 1e-3)), 10)
+            exp2 = np.minimum((d/(2*tsr_mod * x + 1e-3)), 10)
+
+            fl = (
+                (1 / (m)) * 
+                ((1.6 * c_on_r / 0.1267) * 
+                ((a - c_on_r**exp1) / 
+                 (b + c_on_r**exp1)) - 1)
+            )
+            fd = (
+                (1 / (m)) * 
+                ((1.6 * c_on_r / 0.1267) * 
+                ((a - c_on_r**exp2) / 
+                 (b + c_on_r**exp2)) - 1)
+            )
+
+            Clp = m * (inflow - alpha_0)
+
+            Cd0 = self.Cd(x, alpha_0)
+
+            del_Cl = np.maximum(0, fl * (Clp - Cl_2D))
+            del_Cd = np.maximum(0, fd * (Cd_2D - Cd0))
+
+            Cl_3D = Cl_2D + del_Cl
+            Cd_3D = Cd_2D + del_Cd
+
+            return Cl_3D, Cd_3D
+
+        else:
+            # No 3D stall correction
+            return Cl_2D, Cd_2D
 
 
 class RotorDefinition:
@@ -153,5 +213,7 @@ def twist(self, mu):
     def solidity(self, mu):
         return self.solidity_func(mu)
 
-    def clcd(self, mu, aoa):
-        return self.airfoil_func(mu, aoa)
+    def clcd(self, mu, aoa, tsr=None, apply_3D_stall_correction=False):
+        solidity = self.solidity(mu)
+        c_on_r = solidity * (2 * np.pi) / (self.N_blades)
+        return self.airfoil_func(mu, aoa, tsr, apply_3D_stall_correction, c_on_r)

MITRotor/TangentialInduction.py

@@ -56,8 +56,8 @@ def __call__(
     ) -> ArrayLike:
         eff_yaw = calc_eff_yaw(yaw, tilt)
         aprime = (
-            np.clip(aero_props.C_tau_corr, -2, 2)
-            / (4 * np.maximum(geom.mu_mesh, 0.1) ** 2 * tsr * (1 - aero_props.an) * np.cos(eff_yaw))
+            np.clip(aero_props.C_tau_corr, -10, 10)
+            / (4 * np.maximum(geom.mu_mesh, 0.1) * tsr * np.maximum((1 - aero_props.an), 0.001) * np.cos(eff_yaw))
         )
 
         return aprime

examples/example_02_rotor_distributions.py

@@ -3,21 +3,23 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from MITRotor import BEM, BEMGeometry, IEA10MW, BEMSolution
+from MITRotor import BEM, BEMGeometry, IEA10MW, BEMSolution, HeckMomentum
 
 figdir = Path("fig")
 figdir.mkdir(exist_ok=True, parents=True)
 
 
 def plot_radial_distributions(sol: BEMSolution, save_to: Path):
     fig, axes = plt.subplots(3, 1, sharex=True)
-    axes[0].plot(sol.geom.mu, sol.Cp(grid="annulus"), label="$C_P$")
-    axes[1].plot(sol.geom.mu, sol.Ct(grid="annulus"), label="$C_T$")
-    axes[2].plot(sol.geom.mu, sol.a(grid="annulus"), label="$a_n$")
+    axes[0].plot(sol.geom.mu, sol.Cp(grid="annulus"))
+    axes[1].plot(sol.geom.mu, sol.Ct(grid="annulus"))
+    axes[2].plot(sol.geom.mu, sol.a(grid="annulus"))
 
-    [ax.legend(loc="lower center") for ax in axes]
 
-    axes[-1].set_xlabel("Radial position, $\mu$ [-]")
+    axes[0].set_ylabel("$C_P$")
+    axes[1].set_ylabel("$C_T$")
+    axes[2].set_ylabel("$a_n$")
+    axes[-1].set_xlabel("$\mu$")
     plt.xlim(0, 1)
     plt.savefig(save_to, dpi=300, bbox_inches="tight")
 
@@ -45,10 +47,10 @@ def plot_azimuthal_variations(sol: BEMSolution, save_to: Path):
 if __name__ == "__main__":
     # Initialize rotor with increased radial resolution.
     rotor = IEA10MW()
-    bem = BEM(rotor=rotor, geometry=BEMGeometry(Nr=100, Ntheta=10))
+    bem = BEM(rotor=rotor, geometry=BEMGeometry(Nr=100, Ntheta=10), momentum_model=HeckMomentum(averaging='sector'))
 
     # solve BEM for a control set point.
-    pitch, tsr, yaw = np.deg2rad(0), 7.0, np.deg2rad(0.0)
+    pitch, tsr, yaw = np.deg2rad(0), 7.0, np.deg2rad(20.0)
     sol = bem(pitch, tsr, yaw)
 
     # Plot

examples/example_03_pitch_tsr_contour.py

@@ -72,13 +72,12 @@ def plot(df_surface: pl.DataFrame):
     cbar.set_label(label=r"$C_T~$(-)")
 
     # Save figure to file
-    plt.savefig(figdir / "example_3_pitch_tsr_contour.png", dpi=500, bbox_inches="tight")
+    plt.savefig(figdir / "example_03_pitch_tsr_contour.png", dpi=500, bbox_inches="tight")
 
 
 if __name__ == "__main__":
     print("Generate BEM data. This may take a couple of minutes.")
     df = generate()
     print(df)
     print("plotting contour.")
-    plot(df)
-    plt.show()
+    plot(df)

examples/example_05_3D_stall_delay_correction.py

@@ -0,0 +1,48 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from pathlib import Path
+from MITRotor import (
+    BEM,
+    IEA10MW,
+    DefaultAerodynamics
+)
+
+figdir = Path("fig")
+figdir.mkdir(exist_ok=True, parents=True)
+
+def test():
+    model = BEM(
+        aerodynamic_model=DefaultAerodynamics(apply_3D_stall_correction=False),
+        rotor=IEA10MW(),
+    )
+
+    aoa = np.radians(np.linspace(-5, 30, 100))
+    mu = 0.4 * np.ones_like(aoa)
+    tsr = 2.0
+    Cl_uncorr, Cd_uncorr = model.rotor.clcd(
+        mu, aoa, tsr=tsr, apply_3D_stall_correction=False
+    )
+    Cl_corr, Cd_corr = model.rotor.clcd(
+        mu, aoa, tsr=tsr, apply_3D_stall_correction=True
+    )
+
+
+    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
+    ax[0].plot(np.degrees(aoa), Cl_uncorr, label="2D Airfoil")
+    ax[0].plot(np.degrees(aoa), Cl_corr, label="3D Stall Corrected")
+    ax[0].set_xlabel("Angle of Attack (degrees)")
+    ax[0].set_ylabel("Cl")
+    ax[0].set_title("Lift Coefficient with 3D Stall Correction")
+    ax[0].legend()
+
+    ax[1].plot(np.degrees(aoa), Cd_uncorr, label="2D Airfoil")
+    ax[1].plot(np.degrees(aoa), Cd_corr, label="3D Stall Corrected", color='orange')
+    ax[1].set_xlabel("Angle of Attack (degrees)")
+    ax[1].set_ylabel("Cd")
+    ax[1].set_title("Drag Coefficient with 3D Stall Correction")
+    ax[1].legend()
+
+    plt.savefig(figdir / "3D_stall_correction_test.png")
+
+if __name__ == "__main__":
+    test()