diff --git a/MITRotor/Aerodynamics.py b/MITRotor/Aerodynamics.py index 265872c..a2538c4 100644 --- a/MITRotor/Aerodynamics.py +++ b/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, diff --git a/MITRotor/BEMSolver.py b/MITRotor/BEMSolver.py index daea9bc..0d8be3f 100644 --- a/MITRotor/BEMSolver.py +++ b/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,10 +114,6 @@ 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) @@ -130,7 +121,7 @@ def Ctprime(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"): 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. diff --git a/MITRotor/RotorDefinition.py b/MITRotor/RotorDefinition.py index 4e8d487..a1065af 100644 --- a/MITRotor/RotorDefinition.py +++ b/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) diff --git a/MITRotor/TangentialInduction.py b/MITRotor/TangentialInduction.py index cd493ec..6d5dbcf 100644 --- a/MITRotor/TangentialInduction.py +++ b/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 diff --git a/examples/example_02_rotor_distributions.py b/examples/example_02_rotor_distributions.py index 517ea21..8df3084 100644 --- a/examples/example_02_rotor_distributions.py +++ b/examples/example_02_rotor_distributions.py @@ -3,7 +3,7 @@ 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) @@ -11,13 +11,15 @@ 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 diff --git a/examples/example_03_pitch_tsr_contour.py b/examples/example_03_pitch_tsr_contour.py index 0b6923e..32eb154 100644 --- a/examples/example_03_pitch_tsr_contour.py +++ b/examples/example_03_pitch_tsr_contour.py @@ -72,7 +72,7 @@ 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__": @@ -80,5 +80,4 @@ def plot(df_surface: pl.DataFrame): df = generate() print(df) print("plotting contour.") - plot(df) - plt.show() + plot(df) \ No newline at end of file diff --git a/examples/example_05_3D_stall_delay_correction.py b/examples/example_05_3D_stall_delay_correction.py new file mode 100644 index 0000000..ae65bf2 --- /dev/null +++ b/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() \ No newline at end of file