Merge pull request #70 from ptrbortolotti/dev

fcl non negative, try except around 3d corrections

Author: ebranlard

pyFAST/airfoils/Polar.py

@@ -394,8 +394,11 @@ def correction3D(
                     / 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=1.0
+                fcl=0.0
         elif lift_method == "Snel":
             # Snel correction
             fcl = 3.0 * chord_over_r ** 2.0
@@ -730,7 +733,7 @@ def unsteadyParams(self, window_offset=None, nMin=720):
         else:
             window = [-20, 20]
         try:
-            alpha0cn = _find_alpha0(alpha, cn, window, direction='up')
+            alpha0cn = _find_alpha0(alpha, cn, window, direction='up', value_if_constant = 0.)
         except NoCrossingException:
             print("[WARN] Polar: Cn unsteady, cannot find zero crossing with up direction, trying down direction")
             alpha0cn = _find_alpha0(alpha, cn, window, direction='down')
@@ -1265,17 +1268,17 @@ def _alpha_window_in_bounds(alpha, window):
     return window
 
 
-def _find_alpha0(alpha, coeff, window, direction='up'):
+def _find_alpha0(alpha, coeff, window, direction='up', value_if_constant = np.nan):
     """Finds the point where coeff(alpha)==0 using interpolation.
     The search is narrowed to a window that can be specified by the user. The default window is yet enough for cases that make physical sense.
     The angle alpha0 is found by looking at a zero up crossing in this window, and interpolation is used to find the exact location.
     """
     # Constant case or only one value
-    if np.all(coeff == coeff[0]) or len(coeff) == 1:
+    if np.all(abs((coeff - coeff[0])<1e-8)) or len(coeff) == 1:
         if coeff[0] == 0:
             return 0
         else:
-            return np.nan
+            return value_if_constant
     # Ensuring window is within our alpha values
     window = _alpha_window_in_bounds(alpha, window)
 

pyFAST/converters/hawc2ToOpenfast.py

@@ -192,12 +192,15 @@ def AE_PC_C2_toAD(ae_filename, pc_filename, blade_c2def, r_AD=None, ADbldFilenam
         P = Polar(Re=Re, alpha=vAlpha, cl=M[:,1], cd=M[:,2], cm=M[:,3], radians=False)
         # Apply 3D correction
         if r>0 and correction3D:
-            P = P.correction3D(
-                r_over_R     = r/r_max,
-                chord_over_r = c/r,
-                tsr=tsr, 
-                lift_method="DuSelig", drag_method="None", blending_method="linear_25_45",
-                max_cl_corr=0.25)
+            try:
+                P = P.correction3D(
+                    r_over_R     = r/r_max,
+                    chord_over_r = c/r,
+                    tsr=tsr, 
+                    lift_method="DuSelig", drag_method="None", blending_method="linear_25_45",
+                    max_cl_corr=0.25)
+            except:
+                print('3D correction not applied at station ', ir)
 
         # Store polars
         M = np.column_stack((P.alpha, P.cl, P.cd, P.cm))