Inconsistent output for correction3D() method when self._radians == True vs. False

[jbdv-no]

If I am using all functionality correctly, I believe there may be issues with the correction3D method of the openfast_toolbox.airfoils.Polar class, when the self._radians property of an instance is True.

Consider the following example:

import numpy as np
import matplotlib.pyplot as plt
from openfast_toolbox.airfoils.Polar import Polar

# aoa_rad = ...
# cl = ...
# cd = ...
# cm = ...

# may also create an object, with the same outcome, without passing `radians=True`  
polarRAD = Polar(alpha=aoa_rad, cl=cl, cd=cd, cm=cm, Re=1.0, radians=True)

# for reference
polarDEG = Polar(alpha=np.rad2deg(aoa_rad), cl=cl, cd=cd, cm=cm, Re=1.0, radians=False)

r_R = 0.197
c_r = 0.232
tsr = 9.153
polarRAD3D = polarRAD.correction3D(r_over_R=r_R, chord_over_r=c_r, tsr=tsr)

# for reference
polarDEG3D = polarDEG.correction3D(r_over_R=r_R, chord_over_r=c_r, tsr=tsr)

ax = plt.subplots(1, 1)[1]
ax.plot(np.rad2deg(aoa_rad), cl, label='original')
ax.plot(polarRAD3D.alpha, polarRAD3D.cl, ls='-', label='polarRAD3D')
ax.plot(polarDEG3D.alpha, polarDEG3D.cl, ls='--', label='polarDEG3D')
ax.legend()
ax.set_xlim([-50, 50])

plt.show()

Executing this code results in the following error:

Exception has occurred: UnboundLocalError
cannot access local variable 'alpha' where it is not associated with a value
File "...path...\openfast_toolbox\openfast_toolbox\airfoils\Polar.py", line 360, in correction3D

where the relevant line(s) of code in the openfast_toolbox is:

        # rename and convert units for convenience
        if self._radians:
            alpha = alpha  # <- problematic line
        else:
            alpha = np.radians(self.alpha)

I believe this line should read:

        if self._radians:
            alpha = self.alpha

Making this change the program execution continues, however, the output from the polars generated with input angle in degree/radian are not the same.

As far as I can tell, this is due to the following:

• In lines 346-249, a number of variables related to the linear region of the polar are calculated. If self._radians == True, the angular alpha_* variables have radians as unit, and cl_slope 1/radians -- as expected I believe:

              alpha_linear_region, _, cl_slope, alpha0 = self.linear_region()
              alpha_linear_min = alpha_linear_region[0]
              alpha_linear_max = alpha_linear_region[-1]
              _, alpha_max_corr = self.cl_max()

• In lines 365-367, these alpha_* variables are converted to radians again:

          alpha_max_corr = np.radians(alpha_max_corr)
          alpha_linear_min = np.radians(alpha_linear_min)
          alpha_linear_max = np.radians(alpha_linear_max)

• Also in lines 386/7, alpha0 and cl_slope are converted as if they have units degree and 1/degree, respectively:

              cl_slope = np.degrees(cl_slope)
              alpha0 = np.radians(alpha0)

If I modify the logic of the relevant routine to the following:

    def correction3D(
        self,
        r_over_R,
        chord_over_r,
        tsr,
        lift_method="DuSelig",
        drag_method="None",
        blending_method="linear_25_45",
        max_cl_corr=0.25,
        alpha_max_corr=None,
        alpha_linear_min=None,
        alpha_linear_max=None,
    ):
        """Applies 3-D corrections for rotating sections from the 2-D data.

        Parameters
        ----------
        r_over_R : float
            local radial position / rotor radius
        chord_over_r : float
            local chord length / local radial location
        tsr : float
            tip-speed ratio
        lift_method : string, optional
            flag switching between Du-Selig and Snel corrections
        drag_method : string, optional
            flag switching between Eggers correction and None
        blending_method: string:
             blending method used to blend from 3D to 2D polar. default 'linear_25_45'
        max_cl_corr: float, optional
             maximum correction allowed, default is 0.25.
        alpha_max_corr : float, optional (deg)
            maximum angle of attack to apply full correction
        alpha_linear_min : float, optional (deg)
            angle of attack where linear portion of lift curve slope begins
        alpha_linear_max : float, optional (deg)
            angle of attack where linear portion of lift curve slope ends

        Returns
        -------
        polar : Polar
            A new Polar object corrected for 3-D effects

        Notes
        -----
        The Du-Selig method :cite:`Du1998A-3-D-stall-del` is used to correct lift, and
        the Eggers method :cite:`Eggers-Jr2003An-assessment-o` is used to correct drag.

        """

        if alpha_max_corr == None and alpha_linear_min == None and alpha_linear_max == None:
            # ASSUMPTIONS:
            #   - No alpha_* are provided
            #   - All are calculated in self's native angular units (DEG or RAD)
            alpha_linear_region, _, cl_slope, alpha0 = self.linear_region()
            alpha_linear_min = alpha_linear_region[0]
            alpha_linear_max = alpha_linear_region[-1]
            _, alpha_max_corr = self.cl_max()
            find_linear_region = False
        elif alpha_max_corr * alpha_linear_min * alpha_linear_max == None:
            raise Exception(
                "Define all or none of the keyword arguments alpha_max_corr, alpha_linear_min, and alpha_linear_max"
            )
        else:
            # ASSUMPTIONS:
            #   - All alpha_* are provided
            #   - Given in DEG, since that is indicated in the docstring
            #   - May/may not be the same as self's native angular units (DEG or RAD)
            find_linear_region = True

        # rename and convert units for convenience
        cl_2d = self.cl
        cd_2d = self.cd
        if self._radians:
            # changed from `alpha = alpha` to `alpha = self.alpha`
            alpha = self.alpha
        else:
            alpha = np.radians(self.alpha)

        if find_linear_region or not self._radians:
            # ASSUMPTIONS:
            #   Here because:
            #   - alpha_* values were all provided in DEG (as per docstring), or
            #   - alpha_* values were calculated in self's native angluar units, which were set to DEG
            alpha_max_corr = np.radians(alpha_max_corr)
            alpha_linear_min = np.radians(alpha_linear_min)
            alpha_linear_max = np.radians(alpha_linear_max)

        # ASSUMPTIONS:
        #   - from this point forward, alpha* quantities (aside from alpha0, which may still be unknown) are in RADIANS

        # parameters in Du-Selig model
        a = 1
        b = 1
        d = 1
        lam = tsr / (1 + tsr ** 2) ** 0.5  # modified tip speed ratio
        if np.abs(r_over_R)>1e-4:
            expon = d / lam / r_over_R
        else:
            expon = d / lam / 1e-4

        # find linear region with numpy polyfit
        if find_linear_region:
            idx = np.logical_and(alpha >= alpha_linear_min, alpha <= alpha_linear_max)
            p = np.polyfit(alpha[idx], cl_2d[idx], 1)
            cl_slope = p[0]
            alpha0 = -p[1] / cl_slope
            # ASSUMPTIONS:
            # - alpha0 is in rad
            # - cl_slope is in 1/rad
        else:
            # added test for self._radians
            if not self._radians:
                # ASSUMPTION:
                #  - Here because alpha0 and cl_slope were calculated w.r.t. self's native angular unit, which was DEG
                cl_slope = np.degrees(cl_slope)
                alpha0 = np.radians(alpha0)


        # ASSUMPTIONS:
        #   - all alpha* quantities are in RAD
        #   - cl_slope is in 1/RAD

        if lift_method == "DuSelig":
            # Du-Selig correction factor
            if np.abs(cl_slope)>1e-4:
                fcl = (
                    1.0
                    / cl_slope
                    * (1.6 * chord_over_r / 0.1267 * (a - chord_over_r ** expon) / (b + chord_over_r ** expon) - 1)
                )
                # Force fcl to stay non-negative
                if fcl < 0.:
                    fcl = 0.
            else:
                fcl=0.0
        elif lift_method == "Snel":
            # Snel correction
            fcl = 3.0 * chord_over_r ** 2.0
        else:
            raise Exception("The keyword argument lift_method (3d correction for lift) can only be DuSelig or Snel.")

        # 3D correction for lift
        cl_linear = cl_slope * (alpha - alpha0)
        cl_corr = fcl * (cl_linear - cl_2d)
        # Bound correction +/- max_cl_corr
        cl_corr = np.clip(cl_corr, -max_cl_corr, max_cl_corr)
        # Blending
        if blending_method == "linear_25_45":
            # We adjust fully between +/- 25 deg, linearly to +/- 45
            adj_alpha = np.radians([-180, -45, -25, 25, 45, 180])
            adj_value = np.array([0, 0, 1, 1, 0, 0])
            adj = np.interp(alpha, adj_alpha, adj_value)
        elif blending_method == "heaviside":
            # Apply (arbitrary!) smoothing function to smoothen the 3D corrections and zero them out away from alpha_max_corr
            delta_corr = 10
            y1 = 1.0 - smooth_heaviside(alpha, k=1, rng=(alpha_max_corr, alpha_max_corr + np.deg2rad(delta_corr)))
            y2 = smooth_heaviside(alpha, k=1, rng=(0.0, np.deg2rad(delta_corr)))
            adj = y1 * y2
        else:
            raise NotImplementedError("blending :", blending_method)
        cl_3d = cl_2d + cl_corr * adj

        # Eggers 2003 correction for drag
        if drag_method == "Eggers":
            delta_cd = cl_corr * (np.sin(alpha) - 0.12 * np.cos(alpha)) / (np.cos(alpha) + 0.12 * np.sin(alpha)) * adj
        elif drag_method == "None":
            delta_cd = 0.0
        else:
            raise Exception("The keyword argument darg_method (3d correction for drag) can only be Eggers or None.")

        cd_3d = cd_2d + delta_cd

        return type(self)(Re=self.Re, alpha=np.degrees(alpha), cl=cl_3d, cd=cd_3d, cm=self.cm, radians=False)

I get the same result for a 3D corrected airfoil, regardless of whether the initial Polar is created with its angle of attack input in degree/radians (which I assume is the expected behaviour)

Am I doing something wrong, using the toolbox incorrectly, or is this possibly a bug?