From db9c4790906f61292398eb0dac4dca7cb10f2563 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Mon, 19 Feb 2024 12:36:24 -0500 Subject: [PATCH] add new tools for assessing airfoil quality --- .../compute_airfoil_quality.py | 107 +++++++++++++----- 1 file changed, 77 insertions(+), 30 deletions(-) diff --git a/studies/AirfoilDatabaseQuality/compute_airfoil_quality.py b/studies/AirfoilDatabaseQuality/compute_airfoil_quality.py index 5711bd06..885dcffb 100644 --- a/studies/AirfoilDatabaseQuality/compute_airfoil_quality.py +++ b/studies/AirfoilDatabaseQuality/compute_airfoil_quality.py @@ -7,9 +7,11 @@ airfoil_database_path = asb._asb_root / "geometry" / "airfoil" / "airfoil_database" current_airfoil_database_names = airfoil_database_path.glob("*.dat") + class QualityError(Exception): pass + def compute_airfoil_quality(af: asb.Airfoil): if af.name in current_airfoil_database_names: raise QualityError("Airfoil already exists in the database.") @@ -21,17 +23,17 @@ def compute_airfoil_quality(af: asb.Airfoil): raise QualityError("Airfoil has no coordinates (0).") # Check that the airfoil has enough coordinates - # if len(af.coordinates) < 40: - # raise QualityError("Airfoil has too few coordinates.") + if len(af.coordinates) <= 25: + raise QualityError(f"Airfoil has too few coordinates ({len(af.coordinates)}).") # Check that airfoil is roughly normalized - if np.any(af.x() <= -0.2): + if np.any(af.x() <= -0.05): raise QualityError("Airfoil has abnormally low x-coordinates.") - if np.any(af.x() >= 1.2): + if np.any(af.x() >= 1.05): raise QualityError("Airfoil has abnormally high x-coordinates.") - if np.any(af.y() <= -0.7): + if np.any(af.y() <= -0.5): raise QualityError("Airfoil has abnormally low y-coordinates.") - if np.any(af.y() >= 0.7): + if np.any(af.y() >= 0.5): raise QualityError("Airfoil has abnormally high y-coordinates.") # Check if the airfoil is self-intersecting @@ -47,27 +49,59 @@ def compute_airfoil_quality(af: asb.Airfoil): i = np.argwhere(ds == 0)[0][0] raise QualityError(f"Doubled point at ({af.x()[i]}, {af.y()[i]}).") - # if np.any(ds > 0.12): - # raise QualityError("Airfoil has abnormally long edge lengths.") + if np.any(ds > 0.18): + raise QualityError(f"Airfoil has abnormally long edge lengths ({np.max(ds)}).") angles = np.arctan2d(dy, dx) - if np.any(np.abs(np.diff(angles, period=360)) > 150): - raise QualityError("Airfoil has abnormally large changes in angle.") - + d_angles = np.abs(np.diff(angles, period=360)) + # if np.any(d_angles > 120): + # i = np.argmax(d_angles) + # raise QualityError(f"Airfoil has abnormally large changes in angle at (" + # f"{af.x()[i]}, " + # f"{af.y()[i]}" + # f"), {d_angles[i]:.3g} deg.") + + # is_mid_section = af.x() > 0.05 is_mid_section = np.logical_and( af.x() > 0.05, - af.x() < 0.95 + af.x() < 0.99 ) - if np.any(np.abs(np.diff(angles, period=360)[is_mid_section[1:-1]]) > 30): - # print(np.abs(np.diff(angles, period=360)])) - raise QualityError("Airfoil has abnormally large changes in angle in mid-section.") + if np.any(d_angles[is_mid_section[1:-1]] > 30): + i = np.argmax(d_angles[is_mid_section[1:-1]]) + raise QualityError(f"Airfoil has abnormally large changes in angle in mid-section at (" + f"{af.x()[is_mid_section][i - 1]}, " + f"{af.y()[is_mid_section][i - 1]}" + f"), {d_angles[is_mid_section[1:-1]][i]:.3g} deg.") # Normalize the airfoil af = af.normalize() # Check that the airfoil is representable by a Kulfan airfoil - ka = af.to_kulfan_airfoil() + ka: asb.KulfanAirfoil = af.to_kulfan_airfoil() + + indices_upper = np.logical_and( + af.upper_coordinates()[:, 0] >= 0, + af.upper_coordinates()[:, 0] <= 1 + ) + y_deviance_upper = af.upper_coordinates()[indices_upper, 1] - ka.upper_coordinates( + af.upper_coordinates()[indices_upper, 0])[:, 1] + indices_lower = np.logical_and( + af.lower_coordinates()[:, 0] >= 0, + af.lower_coordinates()[:, 0] <= 1 + ) + y_deviance_lower = af.lower_coordinates()[indices_lower, 1] - ka.lower_coordinates( + af.lower_coordinates()[indices_lower, 0])[:, 1] + + if np.max(np.abs(y_deviance_upper)) > 0.01: + i = np.argmax(np.abs(y_deviance_upper)) + raise QualityError( + f"Airfoil is not representable by a Kulfan airfoil (upper deviance of {y_deviance_upper[i]:.3g} at x={af.upper_coordinates()[i, 0]:.3g}).") + + if np.max(np.abs(y_deviance_lower)) > 0.01: + i = np.argmax(np.abs(y_deviance_lower)) + raise QualityError( + f"Airfoil is not representable by a Kulfan airfoil (lower deviance of {y_deviance_lower[i]:.3g} at x={af.lower_coordinates()[i, 0]:.3g}).") # if not ka.as_shapely_polygon().is_valid: # raise QualityError("Kulfan airfoil is self-intersecting.") @@ -76,9 +110,6 @@ def compute_airfoil_quality(af: asb.Airfoil): # raise QualityError("Airfoil is not representable by a Kulfan airfoil.") - - - if __name__ == '__main__': airfoil_database_path = asb._asb_root / "geometry" / "airfoil" / "airfoil_database" airfoil_database = [ @@ -93,16 +124,32 @@ def compute_airfoil_quality(af: asb.Airfoil): try: compute_airfoil_quality(af) except QualityError as e: - print(f"Airfoil {af.name} failed quality checks: {e}") + print(f"Airfoil {af.name.ljust(20)} failed quality checks: {e}") af.draw() - # all_kulfan_parameters = { - # k: np.stack([ - # af.kulfan_parameters[k] - # for af in UIUC_airfoils - # ], axis=0) - # for k in UIUC_airfoils[0].kulfan_parameters.keys() - # } - # for i in np.argsort(np.abs(kulfan_parameters["leading_edge_weight"]))[-10:]: - # UIUC_airfoils[i].draw() - # asb.Airfoil(UIUC_airfoils[i].name).draw() + airfoil_database_kulfan = [ + af.to_kulfan_airfoil() + for af in airfoil_database + ] + + all_kulfan_parameters = { + k: np.stack([ + af.kulfan_parameters[k] + for af in airfoil_database_kulfan + ], axis=0) + for k in airfoil_database_kulfan[0].kulfan_parameters.keys() + } + # for i in np.argsort(np.abs(all_kulfan_parameters["TE_thickness"]))[-10:]: + # fig, ax = plt.subplots() + # plt.plot( + # airfoil_database[i].normalize().x(), + # airfoil_database[i].normalize().y(), + # ".k", markersize=5, zorder=4, + # ) + # plt.plot( + # airfoil_database_kulfan[i].x(), + # airfoil_database_kulfan[i].y(), + # "-r", linewidth=2, + # ) + # plt.title("Airfoil: " + airfoil_database[i].name) + # plt.show()