From a6eb895ff1a514d4d967dfeb7edd5d4b2845593d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrique=20Mill=C3=A1n=20Valbuena?= Date: Sun, 20 Sep 2020 18:09:45 +0200 Subject: [PATCH] Refactor BEM class (#30) * TST, CLN: Remove loss functions from BEM class - Easier to test. - Functions like this should be outside a class, to make it light and prevent technical debt. * BUG, BLD: Add __init__.py file to `bem` module * BLD: Add `bem` imports * BUG: Fix import * ENH: Include `model` module with the BEM implementation * TST: Use fixture for airfoil polars * ENH, TST: `Propeller` class updated with a new `Section` object * TST: Improve test coverage * BLD: Include `pytest-cov` in dev requirements * DOC: Fill docstrings * ENH: Include radii distribution of the sections * DOC, CLN: General changes * ENH: Refactor BEM class * TST: Remove versioneer files from test coverage * ENH: Make homogeneous the usage of degrees and radians The best solution for the moment that I have found is to only convert from deg to rad just before calling a trigonometric method (or approximation, like in the lift coefficient). Like that, it becomes easy to debug (we are used to read degrees, not radians) and the code is "robust'. Ideally there should be a class that defines `Angle` and computes correctly the trigonometric functions. * TST: Fix tests, they were not using correctly the new degrees/radians standard * ENH: Prevent division by zero * TST: Introduce unitary tests for model * Run tests: they do not assert against anything (is a nonlinear equation, that is why I have built a solver in the first place). * Reproduce: these tests check if we are getting the same results as the last time the code was editted. It is a guide to detect if existing bugs alter the result in a meaningful way. --- .coveragerc | 2 + requirements-dev.txt | 1 + src/pybem/__init__.py | 9 +- src/pybem/bem.py | 382 --------------------------------- src/pybem/bem/__init__.py | 2 + src/pybem/bem/loss.py | 81 +++++++ src/pybem/bem/model.py | 361 +++++++++++++++++++++++++++++++ src/pybem/models/__init__.py | 4 +- src/pybem/models/airfoils.py | 47 ++-- src/pybem/models/propeller.py | 272 ++++++++++++++++++++--- tests/bem/test_loss.py | 51 +++++ tests/bem/test_model.py | 138 ++++++++++++ tests/models/test_airfoils.py | 62 +++--- tests/models/test_propeller.py | 160 +++++++++++++- 14 files changed, 1100 insertions(+), 472 deletions(-) create mode 100644 .coveragerc delete mode 100644 src/pybem/bem.py create mode 100644 src/pybem/bem/__init__.py create mode 100644 src/pybem/bem/loss.py create mode 100644 src/pybem/bem/model.py create mode 100644 tests/bem/test_loss.py create mode 100644 tests/bem/test_model.py diff --git a/.coveragerc b/.coveragerc new file mode 100644 index 0000000..1e6a79a --- /dev/null +++ b/.coveragerc @@ -0,0 +1,2 @@ +[run] +omit = *version* \ No newline at end of file diff --git a/requirements-dev.txt b/requirements-dev.txt index 1a62f4d..9b81b93 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,4 +1,5 @@ pytest==6.0.2 +pytest-cov==2.10.1 black==20.8b1 pip-chill==1.0.0 pdbpp==0.10.2 diff --git a/src/pybem/__init__.py b/src/pybem/__init__.py index 64985c8..d692297 100644 --- a/src/pybem/__init__.py +++ b/src/pybem/__init__.py @@ -3,11 +3,10 @@ __version__ = get_versions()["version"] del get_versions +# Blade Element Method +from pybem.bem.model import BladeElementMethod + # Core -from pybem.models.propeller import Propeller from pybem.models.airfoils import Airfoil from pybem.models.flight_conditions import FlightConditions - - -# Blade Element Method -from pybem.bem import BladeElementMethod +from pybem.models.propeller import Propeller, Section diff --git a/src/pybem/bem.py b/src/pybem/bem.py deleted file mode 100644 index e6f5f2d..0000000 --- a/src/pybem/bem.py +++ /dev/null @@ -1,382 +0,0 @@ -from math import pi - -import numpy as np -from scipy.optimize import newton_krylov - -from pybem import Airfoil, FlightConditions, Propeller - - -class BladeElementMethod: - def __init__(self): - - # Flight conditions - self.flight = None - - # Airfoil - self.airfoil = None - - # Propeller - self.propeller = None - - self.F = 1 - - self.beta = None - self.r = None - self.dr = None - self.D = None - - self.J = None - - self.C_T = None - self.C_Q = None - - self.include_tip_loss = None - - def load_airfoil(self, alpha, polar_cl, polar_cd): - """ - Load airfoil polar. - - Parameters - ---------- - alpha: np.array - In degrees! - - polar_cl: np.array - - polar_cd: np.array - """ - - self.airfoil = Airfoil(alpha, polar_cl, polar_cd) - - def load_similarity(self, J): - """ - Load advance ratio similarity coefficient. - - Parameters - ---------- - J: float - """ - self.J = J - - def load_flight_conditions(self, airspeed, omega, altitude): - - self.flight = FlightConditions(airspeed, omega, altitude) - - def load_propeller(self, dist_r, dist_beta, dist_chord, n_blades=2): - """ - Create nondimensional propeller. - - Parameters - ---------- - dist_r: np.array - - dist_beta: np.array - - dist_chord: np.array - - n_blades: int - - Returns - ------- - Propeller object - """ - # Get geometry - r_hub = dist_r[0] - r_tip = dist_r[-1] - - _D = 2.0 * r_tip - - self.D = _D - - # Non-dimensional radius distribution - _dist_r = dist_r / _D - - # Create propeller - self.propeller = Propeller(r_hub, r_tip, _dist_r, dist_beta, dist_chord) - - # Load number of blades - self.B = n_blades - - return self.propeller - - def set_tip_loss(self, flag=True): - """ - Activate Prandtl tip loss model. - - Parameters - ---------- - flag: bool - """ - self.include_tip_loss = flag - - def induction_axial(self, r, phi, beta): - - _F = self.compute_tip_loss(r, phi) - _solidity = self.solidity(r) - - # Convert to radians - _phi, _beta = np.deg2rad([phi, beta]) - - # Compute angle of attack - alpha = beta - phi - - # Polar - _cl = self.airfoil.cl(alpha) - _cd = self.airfoil.cd(_cl) - - NUM = 4.0 * _F * (np.sin(_phi)) ** 2.0 - DEN = _solidity * (_cl * np.cos(_phi) - _cd * np.sin(_phi)) - - frac = NUM / DEN - 1 - - return 1.0 / frac - - def induction_tangential(self, r, phi, beta): - - # Convert to radians - _beta, _phi = np.deg2rad([beta, phi]) - - # Compute angle of attack - alpha = beta - phi - - F = self.compute_tip_loss(r, phi) - solidity = self.solidity(r) - - # Polar - _cl = self.airfoil.cl(alpha) - _cd = self.airfoil.cd(_cl) - - NUM = 4.0 * F * np.sin(_phi) * np.cos(_phi) - DEN = solidity * (_cl * np.sin(_phi) - _cd * np.cos(_phi)) - - frac = NUM / DEN + 1 - - return 1.0 / frac - - def solidity(self, r): - """ - Local element solidity. - - Parameters - ---------- - r: float - - Returns - ------- - float - """ - _D = self.D - _B = self.B - - _chord = self.propeller.chord(r) - - NUM = _B * _chord - DEN = 2.0 * pi * (r * _D) - - return NUM / DEN - - def compute_inflow_angle(self, r, x0=50): - """ - Solve nonlinear inflow angle equation. - - Parameters - ---------- - r: float - r / D parameter - - x0: float - Initial condition - - Returns - ------- - float - If NO convergence, returns np.nan - """ - self.beta = self.propeller.beta(r) - self.r = r - - try: - sol = newton_krylov(self._residual, x0) - except Exception as ex: - print(ex) - sol = np.array(np.nan) - - return sol.item() - - def compute_loads(self, dr=0.01): - """ - Compute blade loads. - - Parameters - ---------- - dr: float - Spacing between stations. - - Returns - ------- - tuple: thrust, torque - - """ - _J = self.J - _D = self.D - - _r_hub = self.propeller.r_hub / _D - _r_tip = self.propeller.r_tip / _D - - # Number of steps - N = np.floor((_r_tip - _r_hub) / dr) - - # Create nondimensional radius distribution - r_space = np.linspace(start=_r_hub, stop=_r_tip, num=N) - - # Initial condition - C_T = [0.0] - C_Q = [0.0] - - # Initial condition is 20% larger than the twist angle - phi0 = np.arctan(_J / _r_hub) - phi0 = np.rad2deg(phi0) - - phi_space = [phi0] - F_space = [1.0] - - idx = 0 # Index to control numpy arrays - for r in r_space[:-1]: - - # Compute induction angle - phi = self.compute_inflow_angle(r, phi_space[idx]) - - # Compute induction coefficients - axi = self.induction_axial(r, phi, self.beta) - tng = self.induction_tangential(r, phi, self.beta) - - # Tip loss - _F = self.compute_tip_loss(r, phi) - - # Compute Euler **implicit** derivative - F_T = 4.0 * pi * _J ** 2.0 * (r + dr) ** 1.0 * (1.0 + axi) * axi * _F - F_Q = 4.0 * pi * _J ** 1.0 * (r + dr) ** 3.0 * (1.0 + axi) * tng * _F - - # Integrate - C_T.append(C_T[idx] + dr * F_T) - C_Q.append(C_Q[idx] + dr * F_Q) - - # Save state - F_space.append(_F) - phi_space.append(phi) - - idx += 1 - - C_T = np.array(C_T) - C_Q = np.array(C_Q) - - self.C_T = C_T - self.C_Q = C_Q - - self.dr = dr - self.N = N - - # Pack up results - result = dict() - - result["r"] = r_space - result["C_T"] = C_T - result["C_Q"] = C_Q - result["F"] = F_space - result["phi"] = phi_space - - return result - - def _residual(self, phi): - """ - Nonlinear inflow angle equation residual. - - Parameters - ---------- - phi: float - Inflow angle, in deg - - Returns - ------- - float - Equation residual - """ - _phi = np.deg2rad(phi) - - _J = self.J - _r = self.r - - # Compute incidence coefficients - ind_ax = self.induction_axial(_r, phi, self.beta) - ind_tan = self.induction_tangential(_r, phi, self.beta) - - # Compute residual - return np.tan(_phi) - (_J / _r) * ((1.0 + ind_ax) / (1.0 - ind_tan)) - - def compute_tip_loss(self, r, phi): - """ - Prandtl tip loss coefficient. - - Parameters - ---------- - r: float - - phi: float - In degrees. - - Returns - ------- - F: float - """ - # Check correct use - if self.include_tip_loss is None: - raise ValueError("You need to invoke the set_tip_loss method first!") - - if self.include_tip_loss == False: - - return 1.0 - - else: - - _phi = np.deg2rad(phi) - - # Recover geometry - _B = self.B - _D = self.D - _R = self.propeller.r_tip - - _r = r * _D - - NUM = _B * (_R - _r) - DEN = 2.0 * _r * np.sin(_phi) - - f_tip = NUM / DEN - - return 2.0 * np.arccos(np.exp(-f_tip)) / pi - - def compute_hub_loss(self, r, phi): - """ - Hub loss coefficient. - - Parameters - ---------- - r: float - - phi: float - In degrees. - - Returns - ------- - F: float - """ - _phi = np.deg2rad(phi) - - _B = self.B - _R = self.propeller.r_hub - - N = _B * (r - _R) - D = 2 * r * np.sin(_phi) - - f_hub = N / D - - return 2.0 * np.arccos(np.exp(-f_hub)) / pi diff --git a/src/pybem/bem/__init__.py b/src/pybem/bem/__init__.py new file mode 100644 index 0000000..9654f5b --- /dev/null +++ b/src/pybem/bem/__init__.py @@ -0,0 +1,2 @@ +from .model import BladeElementMethod +from .loss import compute_hub_loss, compute_tip_loss diff --git a/src/pybem/bem/loss.py b/src/pybem/bem/loss.py new file mode 100644 index 0000000..1cf98d6 --- /dev/null +++ b/src/pybem/bem/loss.py @@ -0,0 +1,81 @@ +"""Prandtl loss models to account for hub and tip vortices. +""" + +from math import pi + +import numpy as np + + +def func(phi, x, a): + """Auxiliary function for the generalization of the tip losses. + + Parameters + ---------- + phi : float + Incidence angle in degrees. + x : float + Dimensionless radius r / R + a : float + Constant to distinguish between tip and hub loss. + + Returns + ------- + float + """ + + phi = np.deg2rad(phi) + + NUM = a - x + DEN = 2.0 * np.sin(phi) * x + + return NUM / DEN + + +def compute_tip_loss(phi, r, B): + """Compute Prandtl tip loss coefficient. + + Parameters + ---------- + phi : float + Incidence angle in degrees. + r : float + Dimensionless radius, r / R + B : int + Number of blades + + Returns + ------- + F_tip: float + """ + f_tip = func(phi=phi, x=r, a=1.0) + f_tip *= B + + F_tip = 2.0 * np.arccos(np.exp(-f_tip)) / pi + + return F_tip + + +def compute_hub_loss(phi, r, B, r_hubR): + """Compute Prandtl hub loss coefficient. + + Parameters + ---------- + phi : float + Incidence angle in degrees. + r : float + Dimensionless radius, r / R + B : int + Number of blades + r_hubR : float + Dimensionless location of the hub radius, R_hub / R + + Returns + ------- + F_hub: float + """ + f_hub = func(phi=phi, x=r, a=r_hubR) + f_hub *= -B + + F_hub = 2.0 * np.arccos(np.exp(-f_hub)) / pi + + return F_hub \ No newline at end of file diff --git a/src/pybem/bem/model.py b/src/pybem/bem/model.py new file mode 100644 index 0000000..a267620 --- /dev/null +++ b/src/pybem/bem/model.py @@ -0,0 +1,361 @@ +from functools import partial +from math import pi + +import numpy as np +from scipy.optimize import newton_krylov + +from .loss import compute_hub_loss, compute_tip_loss + + +class BladeElementMethod: + """Blade Element Method implementation. + + This class solves the BEM equations to compute the performance and + load distribution along a propeller. + + Parameters + ---------- + J : float + Advance ratio. + propeller + flight + tip_loss : bool + hub_loss : bool + + Attributes + ---------- + J : float + Advance ratio. + propeller + flight + tip_loss : bool + hub_loss : bool + r_dist : list of floats + phi : list of floats + C_T : np.array + C_Q : np.array + """ + + N_SECTIONS = 10 + EPSILON = 1e-6 + + def __init__(self, J, propeller, flight=None, tip_loss=False, hub_loss=False): + + self.J = J + + self.flight = flight + self.propeller = propeller + + self.tip_loss = tip_loss + self.hub_loss = hub_loss + + # Dimensionless thrust and torque distribution + self.phi = [] + self.C_T = None + self.C_Q = None + + def solve(self): + """Solve the BEM""" + + r_min = self.propeller.radii[0] + + # Create dimensionless radius distribution + start = r_min * (1.0 + self.EPSILON) + stop = 1.0 - self.EPSILON + r_dist = np.linspace(start=start, stop=stop, num=self.N_SECTIONS) + + phi0 = self.propeller.compute_beta(r_min) + phi = [] + for r in r_dist: + + # Solve inflow angle. + _phi = self.compute_inflow_angle(r=r, phi0=phi0) + + phi.append(_phi) + + phi0 = _phi + + self.phi = phi + self.r_dist = r_dist + + return r_dist, phi + + def _residual(self, phi, r): + """Nonlinear inflow angle equation residual. + + Parameters + ---------- + phi : float + Incidence angle, in degrees. + r : float + Dimensionless radius, r / R. + beta : float + Twist angle. + + Returns + ------- + residual : float + """ + + propeller = self.propeller + + # Compute angle of attack + beta = propeller.compute_beta(r) + alpha = beta - phi + + # Compute lift and drag coefficients + cl = propeller.compute_cl(r=r, alpha=alpha) + cd = propeller.compute_cd(r=r, alpha=alpha) + + # Compute incidence coefficients + sigma = propeller.compute_solidity(r) + + a = self.compute_axial_coefficient(r=r, phi=phi, cl=cl, cd=cd, sigma=sigma) + b = self.compute_angular_coefficient(r=r, phi=phi, cl=cl, cd=cd, sigma=sigma) + + # Compute residual + J = self.J + + phi = np.deg2rad(phi) + + SUM_1 = np.sin(phi) / (1.0 + a) + SUM_2 = J * np.cos(phi) / (r * (1.0 - b)) + + res = SUM_1 - SUM_2 + + return res + + def compute_prandtl_loss(self, r, phi): + """Compute tip and hub losses according to the Prandtl model. + + Parameters + ---------- + r : float + phi : float + Incidence angle in degrees + + Returns + ------- + F : float + Correction factor + """ + + propeller = self.propeller + + F = 1.0 # Assume no loss + + if self.tip_loss == True: + + F_tip = compute_tip_loss(B=propeller.B, r=r, phi=phi) + + F *= F_tip + + if self.hub_loss == True: + + r_hubR = propeller.R_hub / propeller.R + F_hub = compute_hub_loss(B=propeller.B, r=r, phi=phi, r_hubR=r_hubR) + + F *= F_hub + + return F + + @staticmethod + def compute_ct(cl, cd, phi): + """Compute section tangential force coefficient. + + Parameters + ---------- + cl : float + cd : float + phi : float + Incidence angle in degrees. + + Returns + ------- + ct : float + """ + + phi = np.deg2rad(phi) + + ct = cl * np.cos(phi) - cd * np.sin(phi) + + return ct + + @staticmethod + def compute_cn(cl, cd, phi): + """Compute section angular force coefficient. + + Parameters + ---------- + cl : float + cd : float + phi : float + Incidence angle in degrees. + + Returns + ------- + cn : float + """ + + phi = np.deg2rad(phi) + + cn = cl * np.sin(phi) + cd * np.cos(phi) + + return cn + + def compute_axial_coefficient(self, r, phi, cl, cd, sigma): + + ct = self.compute_ct(cl=cl, cd=cd, phi=phi) + + # Compute loss coefficients (if necessary) + F = self.compute_prandtl_loss(r=r, phi=phi) + + # Convert to radians + phi = np.deg2rad(phi) + + NUM = 4.0 * F * r * (np.sin(phi)) ** 2.0 + DEN = sigma * ct + + # Prevent division by zero warning + if np.isclose(DEN, 0.0): + frac = np.inf + else: + frac = NUM / DEN - 1.0 + + a = 1.0 / frac + + return a + + def compute_angular_coefficient(self, r, phi, cl, cd, sigma): + + cn = self.compute_cn(cl=cl, cd=cd, phi=phi) + + # Compute loss coefficients (if necessary) + F = self.compute_prandtl_loss(r=r, phi=phi) + + # Convert to radians + phi = np.deg2rad(phi) + + NUM = 4.0 * F * r * np.sin(phi) * np.cos(phi) + DEN = sigma * cn + + # Prevent division by zero warning + if np.isclose(DEN, 0.0): + frac = np.inf + else: + frac = NUM / DEN + 1.0 + + b = 1.0 / frac + + return b + + def compute_inflow_angle(self, r, phi0=0.0): + """Solve nonlinear inflow angle equation. + + Parameters + ---------- + r: float + r / D parameter + + x0: float + Initial condition + + Returns + ------- + float + If NO convergence, returns np.nan + """ + + # Fix section + func = partial(self._residual, r=r) + + # try: + phi = newton_krylov(func, phi0) + # except Exception as ex: + # print(ex) + # phi = np.array(np.nan) + + return phi.item() + + def compute_loads(self, dr=0.01): + """Compute blade loads. + + Parameters + ---------- + dr: float + Spacing between stations. + + Returns + ------- + tuple: thrust, torque + + """ + _J = self.J + _D = self.D + + _r_hub = self.propeller.r_hub / _D + _r_tip = self.propeller.r_tip / _D + + # Number of steps + N = np.floor((_r_tip - _r_hub) / dr) + + # Create nondimensional radius distribution + r_space = np.linspace(start=_r_hub, stop=_r_tip, num=N) + + # Initial condition + C_T = [0.0] + C_Q = [0.0] + + # Initial condition is 20% larger than the twist angle + phi0 = np.arctan(_J / _r_hub) + phi0 = np.rad2deg(phi0) + + phi_space = [phi0] + F_space = [1.0] + + idx = 0 # Index to control numpy arrays + for r in r_space[:-1]: + + # Compute induction angle + phi = self.compute_inflow_angle(r, phi_space[idx]) + + # Compute induction coefficients + axi = self.compute_axial_coefficient(r, phi, self.beta) + tng = self.compute_cn(r, phi, self.beta) + + # Tip loss + _F = self.compute_tip_loss(r, phi) + + # Compute Euler **implicit** derivative + F_T = 4.0 * pi * _J ** 2.0 * (r + dr) ** 1.0 * (1.0 + axi) * axi * _F + F_Q = 4.0 * pi * _J ** 1.0 * (r + dr) ** 3.0 * (1.0 + axi) * tng * _F + + # Integrate + C_T.append(C_T[idx] + dr * F_T) + C_Q.append(C_Q[idx] + dr * F_Q) + + # Save state + F_space.append(_F) + phi_space.append(phi) + + idx += 1 + + C_T = np.array(C_T) + C_Q = np.array(C_Q) + + self.C_T = C_T + self.C_Q = C_Q + + self.dr = dr + self.N = N + + # Pack up results + result = dict() + + result["r"] = r_space + result["C_T"] = C_T + result["C_Q"] = C_Q + result["F"] = F_space + result["phi"] = phi_space + + return result diff --git a/src/pybem/models/__init__.py b/src/pybem/models/__init__.py index 9d80376..9405c52 100644 --- a/src/pybem/models/__init__.py +++ b/src/pybem/models/__init__.py @@ -1,2 +1,2 @@ -from .airfoils import Airfoil -from .propeller import Propeller \ No newline at end of file +from .airfoils import BaseAirfoil, Airfoil +from .propeller import Propeller, Section \ No newline at end of file diff --git a/src/pybem/models/airfoils.py b/src/pybem/models/airfoils.py index a87dea5..f759bb3 100644 --- a/src/pybem/models/airfoils.py +++ b/src/pybem/models/airfoils.py @@ -1,41 +1,54 @@ from math import pi +import numpy as np from scipy.interpolate import interp1d -class AnalyticalAirfoil: - """Analytical implementation +class BaseAirfoil: + """Analytical airfoil implementation. + + Parameters + ---------- + cl_coeff : float + Constant to multiply the analytical lift coefficient. + cd_coeff : float + Constant to multiply the analytical drag coefficient. Attributes ---------- alpha : float - Angle of attack in radians. + Angle of attack in degrees. cl : float Lift coefficient. cd : float Drag coefficient. """ - def __init__(self): + def __init__(self, cl_coeff=1.0, cd_coeff=1.0): self.alpha = None self.cl = None self.cd = None + self.cl_coeff = cl_coeff + self.cd_coeff = cd_coeff + def compute_cl(self, alpha): """Compute lift coefficient. Parameters ---------- alpha : float - Angle of attack in radians. + Angle of attack in degrees. Returns ------- cl : float """ - cl = 2.0 * pi * alpha + alpha = np.deg2rad(alpha) + + cl = 2.0 * pi * alpha * self.cl_coeff self.alpha = alpha self.cl = cl @@ -48,7 +61,7 @@ def compute_cd(self, alpha): Parameters ---------- alpha : float - Angle of attack in radians. + Angle of attack in degrees. Returns ------- @@ -57,7 +70,7 @@ def compute_cd(self, alpha): cl = self.compute_cl(alpha=alpha) - cd = cl ** 2.0 + cd = self.cd_coeff * (cl ** 2.0) self.alpha = alpha self.cl = cl @@ -66,17 +79,24 @@ def compute_cd(self, alpha): return cd -class Airfoil(AnalyticalAirfoil): +class Airfoil(BaseAirfoil): """Airfoil section aerodynamic coefficients. Parameters ---------- alpha : array-like of floats - Angles of attack for `polar_cl` in radians. + Angles of attack for `polar_cl` in degrees. polar_cl : array-like of floats - Lift coefficients. + Lift coefficients for each angle present in alpha. polar_cd : array-like of floats - Drag coefficients. + Drag coefficients for each angle present in alpha. + + Attributes + ---------- + polar_cl : array-like of floats + polar_cd : array-like of floats + interpolant_cl : scipy.interpolate.interp1-like object + interpolant_cd : scipy.interpolate.interp1-like object """ def __init__(self, alpha, polar_cl, polar_cd): @@ -118,7 +138,10 @@ def compute_cd(self, cl=None, alpha=None): Parameters ---------- + alpha : float + Angle of attack. cl : float + Lift coefficient. Returns ------- diff --git a/src/pybem/models/propeller.py b/src/pybem/models/propeller.py index 26aa629..2a72c91 100644 --- a/src/pybem/models/propeller.py +++ b/src/pybem/models/propeller.py @@ -1,72 +1,274 @@ +from bisect import bisect_left +from math import pi + from scipy.interpolate import interp1d +class Section: + """Propeller blade section. + + This class is used by the `Propeller` class, to build the distribution of + chord, twist and solidity across the blade. + + Attributes + ---------- + r : float + Dimensionless radius. + beta : float + Geometrical twist. + chord : float + Airfoil section chord. + solidity : float + airfoil : AirfoilBase-like object + """ + + def __init__(self, r, beta, chord, airfoil, name=None): + + self.r = r + self.beta = beta + self.chord = chord + self.airfoil = airfoil + self.name = name + + def __repr__(self): + return ( + f"Section {self.name} (r = {self.r}, c = {self.chord}, beta = {self.beta})" + ) + + def cl(self, alpha): + """Section lift coefficient. + + Parameters + ---------- + alpha : float + + Returns + ------- + cl : float + """ + + _cl = self.airfoil.compute_cl(alpha=alpha) + + return _cl + + def cd(self, alpha): + """Section drag coefficient. + + Parameters + ---------- + alpha : float + + Returns + ------- + cd : float + """ + + _cd = self.airfoil.compute_cd(alpha=alpha) + + return _cd + + def solidity(self, B, R): + """Section solidity. + + Parameters + ---------- + B : int + Number of blades in the propeller. + R : float + Propeller radius. + + Returns + ------- + sigma : float + """ + + NUM = B * self.chord + DEN = 2.0 * pi * R + + sigma = NUM / DEN + + return sigma + + class Propeller: """Propeller definition. Parameters ---------- - r_hub : float - r_tip : float - r_dist : array-like of floats - beta_dist : array-like of floats - chord_dist : array-like of floats + B : int + Number of blades. + sections : list of Sections-like objects + + Attributes + ---------- + R : float + R_hub : float + sections : sorted list of Sections-like objects + radii : list-like of floats + Distribution of dimensionless radius distribution """ - def __init__(self, r_hub, r_tip, r_dist, beta_dist, chord_dist): + def __init__(self, B, sections): + + self.B = B # Get blade bounds - self.r_hub = r_hub - self.r_tip = r_tip + self.R = max(section.r for section in sections) + self.R_hub = min(section.r for section in sections) + + # Sort sections by ascending order of radius + _sections = sorted(sections, key=lambda section: section.r) + + self.sections = [] + dist_r = [] + dist_beta = [] + dist_chord = [] + for section in _sections: + + # Make radius dimensionless + section.r = section.r / self.R + + dist_r.append(section.r) + dist_beta.append(section.beta) + dist_chord.append(section.chord) - # Get airfoil properties - self.dist_twist = beta_dist - self.dist_chord = chord_dist - self.dist_r = r_dist + # Save + self.sections.append(section) - self.interpolant_twist = interp1d(r_dist, beta_dist) - self.interpolant_chord = interp1d(r_dist, chord_dist) + self._interpolant_chord = interp1d(dist_r, dist_chord) + self._interpolant_beta = interp1d(dist_r, dist_beta) - self.r = None - self._beta = None - self._chord = None + self.radii = dist_r - def get_beta(self, r): - """Airfoil twist at location r. + def _find_bracket(self, r): + """For a given dimensionless location `r`, find the + two surrounding airfoil Sections. Parameters ---------- - r: float - Location based on r_loc. + r : float + Dimensionless radius. Returns ------- - beta: float + left : Section-like object. + right : Section-like object. """ - self.r = r - _beta = self.interpolant_twist(r).item() + _radius = [section.r for section in self.sections] + + idx = bisect_left(_radius, r) - self._beta = _beta + left = self.sections[idx - 1] + right = self.sections[idx] - return _beta + return left, right - def get_chord(self, r): - """Airfoil chord at location r. + def compute_chord(self, r): + """Compute section chord via interpolation between sections. Parameters ---------- - r: float - Location based on r_loc. + r : float Returns ------- - chord: float + chord : float """ - self.r = r + chord = self._interpolant_chord(r).item() + + return chord + + def compute_beta(self, r): + """Compute section beta via interpolation between sections. + + Parameters + ---------- + r : float + + Returns + ------- + beta : float + """ + beta = self._interpolant_beta(r).item() + + return beta + + def compute_solidity(self, r): + """Compute section solidity via interpolation between sections. + + Parameters + ---------- + r : float + + Returns + ------- + sgima : float + """ + + left, right = self._find_bracket(r) + + radii = [left.r, right.r] + solidities = [ + left.solidity(B=self.B, R=self.R), + right.solidity(B=self.B, R=self.R), + ] + f = interp1d(radii, solidities) + + sigma = f(r).item() + + return sigma + + def compute_cl(self, r, alpha): + """Compute section lift coefficient via interpolation between sections. + + Parameters + ---------- + r : float + alpha : float + + Returns + ------- + cl : float + """ + + left, right = self._find_bracket(r) + + radii = [left.r, right.r] + coeffs = [ + left.cl(alpha=alpha), + right.cl(alpha=alpha), + ] + + f = interp1d(radii, coeffs) + + cl = f(r).item() + + return cl + + def compute_cd(self, r, alpha): + """Compute section drag coefficient via interpolation between sections. + + Parameters + ---------- + r : float + alpha : float + + Returns + ------- + cd : float + """ + + left, right = self._find_bracket(r) + + radii = [left.r, right.r] + coeffs = [ + left.cd(alpha=alpha), + right.cd(alpha=alpha), + ] - _chord = self.interpolant_chord(r).item() + f = interp1d(radii, coeffs) - self._chord = _chord + cd = f(r).item() - return _chord + return cd \ No newline at end of file diff --git a/tests/bem/test_loss.py b/tests/bem/test_loss.py new file mode 100644 index 0000000..49b2074 --- /dev/null +++ b/tests/bem/test_loss.py @@ -0,0 +1,51 @@ +from math import pi + +import numpy as np +import pytest +from numpy.testing import assert_allclose + +from pybem.bem.loss import compute_hub_loss, compute_tip_loss, func + + +def test_auxiliary_function(): + + phi = 30.0 # deg + x = 1.0 + a = 1.0 / 2.0 + + result = func(phi=phi, x=x, a=a) + + expected = -1.0 / 2.0 + + assert_allclose(expected, result) + + +@pytest.mark.parametrize( + argnames="B", argvalues=[(1), (2)], ids=["One blade", "Multiple blades"] +) +def test_tip_loss_factor(B): + + phi = 30.0 # deg + x = 1.0 / 2.0 + + result = compute_tip_loss(phi=phi, r=x, B=B) + + expected = 2.0 / pi * np.arccos(np.exp(-B)) + + assert_allclose(expected, result) + + +@pytest.mark.parametrize( + argnames="B", argvalues=[(1), (2)], ids=["One blade", "Multiple blades"] +) +def test_hub_loss_factor(B): + + phi = 30.0 # deg + x = 1.0 / 2.0 + r_hubR = 1.0 / 5.0 + + result = compute_hub_loss(phi=phi, r=x, B=B, r_hubR=r_hubR) + + expected = 2.0 / pi * np.arccos(np.exp(-B * 3.0 / 5.0)) + + assert_allclose(expected, result) diff --git a/tests/bem/test_model.py b/tests/bem/test_model.py new file mode 100644 index 0000000..ec034cf --- /dev/null +++ b/tests/bem/test_model.py @@ -0,0 +1,138 @@ +import pytest + +from pybem.models import Propeller, Section, BaseAirfoil +from pybem.bem import BladeElementMethod + +from numpy.testing import assert_allclose + + +@pytest.fixture +def propeller(): + + sections = [ + Section( + name="Hub", + r=0.3, + beta=60, + chord=0.4, + airfoil=BaseAirfoil(cl_coeff=1.0, cd_coeff=1e-2), + ), + Section( + name="Middle", + r=0.6, + beta=45, + chord=0.35, + airfoil=BaseAirfoil(cl_coeff=0.85, cd_coeff=1e-3), + ), + Section( + name="Tip", + r=1.2, + beta=30, + chord=0.2, + airfoil=BaseAirfoil(cl_coeff=0.5, cd_coeff=1e-3), + ), + ] + + B = 6 + + propeller = Propeller(B=B, sections=sections) + + return propeller + + +def test_run_bem_no_losses(propeller): + + bem = BladeElementMethod(J=1, propeller=propeller) + + bem.solve() + + pass # What could I assert? + + +def test_run_bem_with_losses(propeller): + + bem = BladeElementMethod(J=1, propeller=propeller, tip_loss=True, hub_loss=True) + + bem.solve() + + pass # What could I assert? + + +def test_reproduce_bem_no_losses(propeller): + """This is not a proper test, but it gives you some guarantees. + + With the actual state of the code, we are getting a result that seems + to be correct for the type of phenomena we are solving. + + Perhaps if we linearise the equations we could validate the solver, but + I haven't had time for the moment. So I check at least that for the same inputs, + I am getting the same outputs from here on. + + If a bug is found in the future, this test would fail, but that is not + necessarilly wrong. If the bug is proved to be correct, then the expected values + need to be updated. + """ + + bem = BladeElementMethod(J=1, propeller=propeller) + + bem.N_SECTIONS = 10 + + bem.solve() + + result_phi = bem.phi + + expected_phi = [ + 65.02457613402427, + 61.53158664709961, + 58.06287232558151, + 54.680737977050924, + 52.70649811115369, + 50.648517749244384, + 48.57169404768725, + 46.52659452619687, + 44.549823308619374, + 42.66534592922986, + ] + + assert_allclose(expected_phi, result_phi) + + +def test_reproduce_bem_with_losses(propeller): + """This is not a proper test, but it gives you some guarantees. + + With the actual state of the code, we are getting a result that seems + to be correct for the type of phenomena we are solving. + + Perhaps if we linearise the equations we could validate the solver, but + I haven't had time for the moment. So I check at least that for the same inputs, + I am getting the same outputs from here on. + + If a bug is found in the future, this test would fail, but that is not + necessarilly wrong. If the bug is proved to be correct, then the expected values + need to be updated. + """ + + bem = BladeElementMethod(J=1, propeller=propeller, tip_loss=True, hub_loss=True) + + bem.solve() + + bem.N_SECTIONS = 10 + + bem.solve() + + result_phi = bem.phi + + expected_phi = [ + 60.011868623160844, + 60.2642555093087, + 57.29588901630371, + 54.09221438071569, + 52.18007160086735, + 50.04596910839645, + 47.76712196831369, + 45.351112313287835, + 42.5933251947877, + 30.13349622416542, + ] + + assert_allclose(expected_phi, result_phi) \ No newline at end of file diff --git a/tests/models/test_airfoils.py b/tests/models/test_airfoils.py index 57a69d0..c8fc195 100644 --- a/tests/models/test_airfoils.py +++ b/tests/models/test_airfoils.py @@ -4,7 +4,27 @@ from numpy.testing import assert_allclose, assert_array_almost_equal import numpy as np -from pybem.models.airfoils import Airfoil, AnalyticalAirfoil +from pybem.models.airfoils import Airfoil, BaseAirfoil + + +@pytest.fixture +def polars(): + + analytical = BaseAirfoil() + + alphas = np.linspace(0, 1, num=1000) + polar_cl = [] + polar_cd = [] + + for alpha in alphas: + + _cl = analytical.compute_cl(alpha) + _cd = analytical.compute_cd(alpha) + + polar_cl.append(_cl) + polar_cd.append(_cd) + + return alphas, polar_cl, polar_cd @pytest.mark.parametrize( @@ -14,28 +34,20 @@ ) def test_analytical_airfoil_cl(alpha, expected_cl): - airfoil = AnalyticalAirfoil() + airfoil = BaseAirfoil() result_cl = airfoil.compute_cl(alpha) - assert_allclose(expected_cl, result_cl) - + deg2rad = pi / 180.0 -def test_airfoil_analytical_interpolation(): + expected_cl = 2.0 * pi * (alpha * deg2rad) - analytical = AnalyticalAirfoil() - - alphas = [-1.0, 0.0, 1.0] - expected_cl = [] - expected_cd = [] + assert_allclose(expected_cl, result_cl) - for alpha in alphas: - _cl = analytical.compute_cl(alpha) - _cd = analytical.compute_cd(alpha) +def test_airfoil_analytical_interpolation(polars): - expected_cl.append(_cl) - expected_cd.append(_cd) + alphas, expected_cl, expected_cd = polars empirical = Airfoil(alpha=alphas, polar_cl=expected_cl, polar_cd=expected_cd) @@ -54,21 +66,9 @@ def test_airfoil_analytical_interpolation(): assert_array_almost_equal(expected_cd, result_cd) -def test_airfoil_unseen_interpolation(): - - analytical = AnalyticalAirfoil() - - alphas = np.linspace(0, 1, num=1000) - polar_cl = [] - polar_cd = [] - - for alpha in alphas: - - _cl = analytical.compute_cl(alpha) - _cd = analytical.compute_cd(alpha) +def test_airfoil_unseen_interpolation(polars): - polar_cl.append(_cl) - polar_cd.append(_cd) + alphas, polar_cl, polar_cd = polars empirical = Airfoil(alpha=alphas, polar_cl=polar_cl, polar_cd=polar_cd) @@ -76,7 +76,9 @@ def test_airfoil_unseen_interpolation(): result_cl = empirical.compute_cl(alpha0) result_cd = empirical.compute_cd(alpha=alpha0) - expected_cl = 2.0 * pi * alpha0 + deg2rad = pi / 180.0 + + expected_cl = 2.0 * pi * (alpha0 * deg2rad) expected_cd = expected_cl ** 2.0 assert_allclose(expected_cl, result_cl) diff --git a/tests/models/test_propeller.py b/tests/models/test_propeller.py index 8d9e9c6..5dde9ef 100644 --- a/tests/models/test_propeller.py +++ b/tests/models/test_propeller.py @@ -1,10 +1,158 @@ +from math import pi + import pytest +from numpy.testing import assert_allclose + +from pybem.models.propeller import Propeller, Section +from pybem.models.airfoils import BaseAirfoil + + +@pytest.fixture +def propeller(): + + sections = [ + Section( + name="Hub", + r=0.3, + beta=60, + chord=0.2, + airfoil=BaseAirfoil(cl_coeff=1.0, cd_coeff=1.0), + ), + Section( + name="Middle", + r=0.6, + beta=45, + chord=0.35, + airfoil=BaseAirfoil(cl_coeff=2.0, cd_coeff=2.0), + ), + Section( + name="Tip", + r=1.2, + beta=30, + chord=0.5, + airfoil=BaseAirfoil(cl_coeff=3.0, cd_coeff=3.0), + ), + ] + + B = 6 + + propeller = Propeller(B=B, sections=sections) + + return propeller + + +@pytest.fixture +def section(): + + return Section(r=1.0, beta=0, chord=1.0, airfoil=BaseAirfoil()) + + +def test_section_solidity(section): + + result = section.solidity(1, 1.0) + + expected = 1.0 / (2.0 * pi) + + assert_allclose(expected, result) + + +def test_section_cl(section): + + alpha = 1.0 + deg2rad = pi / 180.0 + + result = section.cl(alpha=alpha) + + expected = 2.0 * pi * alpha * deg2rad + + assert_allclose(expected, result) + + +def test_section_cd(section): + + alpha = 1.0 + deg2rad = pi / 180.0 + + result = section.cd(alpha=1.0) + expected = (2.0 * pi * alpha * deg2rad) ** 2.0 + + assert_allclose(expected, result) + + +def test_propeller_init_radius(propeller): + + result = [propeller.R_hub, propeller.R] + expected = [0.3, 1.2] + + assert_allclose(expected, result) + + +def test_propeller_compute_chord(propeller): + + result = propeller.compute_chord(0.9) + expected = 0.47 + + assert_allclose(expected, result) + + +def test_propeller_compute_beta(propeller): + + result = propeller.compute_beta(0.9) + expected = 33.0 + + assert_allclose(expected, result) + + +def test_propeller_find_bracket_left(propeller): + + result_sections = propeller._find_bracket(0.3) + + expected_sections = (propeller.sections[0], propeller.sections[1]) + + assert expected_sections == result_sections + + +def test_propeller_find_bracket_right(propeller): + + result_sections = propeller._find_bracket(0.9) + + expected_sections = (propeller.sections[1], propeller.sections[2]) + + assert expected_sections == result_sections + + +def test_propeller_solidity(propeller): + + result_sigma = propeller.compute_solidity(0.9) + + expected_sigma = 0.37401411626595404 + + assert_allclose(expected_sigma, result_sigma) + + +def test_propeller_compute_cl(propeller): + + alpha = 1.0 + deg2rad = pi / 180.0 + + result_cl = propeller.compute_cl(r=0.75, alpha=alpha) + expected_cl = 2.0 * pi * (alpha * deg2rad) * (3.0 + 2.0) / 2.0 + + assert_allclose(expected_cl, result_cl) + + +def test_propeller_compute_cd(propeller): + + alpha = 1.0 + deg2rad = pi / 180.0 -from pybem import Propeller + result_cd = propeller.compute_cd(r=0.75, alpha=alpha) + expected_cd = ( + 2.0 + * pi ** 2.0 + * (alpha * deg2rad) ** 2.0 + * (3.0 * 3.0 ** 2.0 + 2.0 * 2.0 ** 2.0) + ) -@pytest.mark.xfail( - raises=NotImplementedError, reason="Propeller definition needs to be updated." -) -def test_propeller(): - raise NotImplementedError + assert_allclose(expected_cd, result_cd)