Airfoil Script

# %%
# Imports
import matplotlib.pyplot as plt
from enum import Enum

import neuralfoil
import numpy as np

# %%
# INPUTS: Modyify these to change the script

# PARSEC airfoil parametrisation: http://pubs.sciepub.com/ajme/2/4/1.
# These values are the NACA0012 airfoil.
leading_edge_radius = 0.0155
upper_crest_x_location = 0.2966
upper_crest_y_location = 0.06002
upper_crest_curvature = -0.4515
lower_crest_x_location = 0.2966
lower_crest_y_location = -0.06002
lower_crest_curvature = 0.4515
trailing_edge_vertical_offset = 0.0
trailing_edge_thickness = 0.0
trailing_edge_angle_degrees = 0.0
trailing_edge_wedge_angle_degrees = 12.89

# Flow conditions
angle_of_attack = 0.0
reynolds_number = 1_000_000
mach = 0.4

# Constants
model_size = "large"
n_points_per_side = 100

# Toggle plotting on and off
plot = False


# %%
# Helper functions
def degrees_to_radians(value):
    return value * np.pi / 180.0


class Surface(Enum):
    UPPER = 0
    LOWER = 1


def _parsec_airfoil_surface_coords(surface, n_points):
    """Produce coordinates of a surface running from LE to TE. All angles in radians."""

    # Use cosine spacing to group points at the curvy bits of the airfoil
    x_coords = 0.5 * (1 - np.cos(np.linspace(0, np.pi, n_points)))

    x_matrix = np.zeros([n_points, 6], dtype=np.float64)
    for i in range(n_points):
        for j in range(6):
            x_matrix[i, j] = x_coords[i] ** (j + 0.5)

    if surface == Surface.UPPER:
        crest_x_loc = upper_crest_x_location
        crest_y_loc = upper_crest_y_location
        crest_curv = upper_crest_curvature
        te_y_coord = trailing_edge_vertical_offset + 0.5 * trailing_edge_thickness
        te_angle = degrees_to_radians(
            trailing_edge_angle_degrees - 0.5 * trailing_edge_wedge_angle_degrees
        )
        le_fac = 1.0
    else:
        crest_x_loc = lower_crest_x_location
        crest_y_loc = lower_crest_y_location
        crest_curv = lower_crest_curvature
        te_y_coord = trailing_edge_vertical_offset - 0.5 * trailing_edge_thickness
        te_angle = degrees_to_radians(
            trailing_edge_angle_degrees + 0.5 * trailing_edge_wedge_angle_degrees
        )
        le_fac = -1.0

    c = np.ones([6, 6], dtype=np.float64)
    for i in range(6):
        c[1, i] = crest_x_loc ** (0.5 + i)
        c[2, i] = 0.5 + i
        c[3, i] = crest_x_loc ** (-0.5 + i)
        c[4, i] = crest_x_loc ** (-1.5 + i)
        if i != 0:
            c[5, i] = 0.0
    c[3, :] = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5] * c[3, :]
    c[4, :] = [-0.25, 0.75, 3.75, 3.75, 15.75, 24.75] * c[4, :]

    le_term = le_fac * (2 * leading_edge_radius) ** 0.5

    b = np.array(
        [te_y_coord, crest_y_loc, np.tan(te_angle), 0.0, crest_curv, le_term], dtype=np.float64
    )

    a = np.linalg.solve(c, b)
    y_coords = np.dot(x_matrix, a)

    return list(zip(x_coords, y_coords, strict=True))


# %%
# Generate the airfoil coordinates
# Produce airfoil from PARSEC parameters (http://pubs.sciepub.com/ajme/2/4/1)
upper = _parsec_airfoil_surface_coords(Surface.UPPER, n_points=n_points_per_side)
lower = _parsec_airfoil_surface_coords(Surface.LOWER, n_points=n_points_per_side)
# Build coords anti-clockwise from upper trailing edge
airfoil_coords = list(reversed(upper)) + lower[1:]

# %%
# Plot the airfoil

# Close the trailing edge before plotting
x_vals = [p[0] for p in airfoil_coords] + [airfoil_coords[0][0]]
y_vals = [p[1] for p in airfoil_coords] + [airfoil_coords[0][1]]

fig = plt.figure()
ax = fig.subplots()
ax.plot(x_vals, y_vals, linewidth=2.5, solid_capstyle="round", color="#ECEEED")
ax.set_aspect("equal", "box")
ax.axis("off")
if plot:
    plt.show()


# %%
# Perform analysis using neuralfoil.
aero = neuralfoil.get_aero_from_coordinates(
    coordinates=np.array(airfoil_coords),
    alpha=angle_of_attack,
    Re=reynolds_number,
    model_size=model_size,
)
print(f"Lift Coefficient (Cl): {aero.get('CL')}")
print(f"Drag Coefficient (Cd): {aero.get('CD')}")
print(f"Moment Coefficient (Cm): {aero.get('CM')}")
print(f"Lift-to-Drag Ratio (Cl/Cd): {aero.get('CL') / aero.get('CD')}")
print(f"Analysis Confidence: {aero.get('analysis_confidence')}")