Nassu validation #677: NACA 0012 wing geometry - parametric generator (analytic, sharp/TMR TE) + ready-to-use ASCII STL (chord=1, span=1, 4108 tris, watertight)

make_naca0012_stl.py

#!/usr/bin/env python3
"""Generate a watertight NACA 0012 wing STL (analytic, sharp trailing edge).

The thickness uses the closed-TE ("sharp") NACA 0012 coefficient (-0.1036),
matching the NASA Turbulence Modeling Resource (TMR) airfoil definition, so the
geometry is consistent with the Cl/Cp/Cf reference data used for validation.

The airfoil is generated at zero angle of attack with unit chord; rotate to the
target alpha and scale to the target chord (in lattice units) in the case config.
Cross-section is convex, so end caps are a centroid triangle fan (watertight).

Usage:
    python make_naca0012_stl.py --chord 1.0 --span 1.0 \
        --n-chord 120 --n-span 30 --out naca0012.stl [--ascii]
"""
import argparse
import struct

import numpy as np


def naca0012_thickness(x: np.ndarray) -> np.ndarray:
    # 5 * t with t = 0.12; last coeff -0.1036 closes the TE exactly (sharp TE).
    return 0.6 * (
        0.2969 * np.sqrt(x)
        - 0.1260 * x
        - 0.3516 * x**2
        + 0.2843 * x**3
        - 0.1036 * x**4
    )


def section_loop(n_chord: int) -> np.ndarray:
    """Closed, ordered 2D perimeter (x, y) of the section, unique points."""
    # Cosine spacing: dense at LE and TE.
    theta = np.linspace(0.0, np.pi, n_chord)
    x = 0.5 * (1.0 - np.cos(theta))  # 0 (LE) -> 1 (TE)
    yt = naca0012_thickness(x)
    # Upper TE->LE (x decreasing), then lower LE->TE (x increasing).
    upper = np.column_stack([x[::-1], yt[::-1]])  # (1,0) ... (0,0)
    lower = np.column_stack([x, -yt])  # (0,0) ... (1,0)
    # Drop shared LE (start of lower) and shared TE (end of lower) to keep unique.
    loop = np.vstack([upper, lower[1:-1]])
    return loop  # (M, 2), open loop (last connects back to first)


def build_mesh(chord, span, n_chord, n_span):
    loop = section_loop(n_chord) * chord  # (M, 2)
    m = loop.shape[0]
    zs = np.linspace(0.0, span, n_span + 1)
    center = np.array([0.25 * chord, 0.0, 0.5 * span])  # interior reference

    tris = []

    def add(a, b, c):
        a, b, c = np.asarray(a), np.asarray(b), np.asarray(c)
        n = np.cross(b - a, c - a)
        # Orient outward: normal should point away from the interior reference.
        if np.dot(n, (a + b + c) / 3.0 - center) < 0.0:
            b, c = c, b
        tris.append((a, b, c))

    # Side skin: quads between consecutive perimeter edges and span stations.
    for k in range(n_span):
        z0, z1 = zs[k], zs[k + 1]
        for i in range(m):
            j = (i + 1) % m
            p0 = np.array([loop[i, 0], loop[i, 1], z0])
            p1 = np.array([loop[j, 0], loop[j, 1], z0])
            p2 = np.array([loop[j, 0], loop[j, 1], z1])
            p3 = np.array([loop[i, 0], loop[i, 1], z1])
            add(p0, p1, p2)
            add(p0, p2, p3)

    # End caps: centroid triangle fan at z=0 and z=span (section is convex).
    cen2d = loop.mean(axis=0)
    for z in (0.0, span):
        c = np.array([cen2d[0], cen2d[1], z])
        for i in range(m):
            j = (i + 1) % m
            a = np.array([loop[i, 0], loop[i, 1], z])
            b = np.array([loop[j, 0], loop[j, 1], z])
            add(c, a, b)

    return tris


def write_binary_stl(tris, path):
    with open(path, "wb") as f:
        f.write(b"\0" * 80)
        f.write(struct.pack("<I", len(tris)))
        for a, b, c in tris:
            n = np.cross(b - a, c - a)
            ln = np.linalg.norm(n)
            n = n / ln if ln > 0 else n
            f.write(struct.pack("<3f", *n))
            for v in (a, b, c):
                f.write(struct.pack("<3f", *v))
            f.write(struct.pack("<H", 0))


def write_ascii_stl(tris, path, name="naca0012"):
    with open(path, "w") as f:
        f.write(f"solid {name}\n")
        for a, b, c in tris:
            n = np.cross(b - a, c - a)
            ln = np.linalg.norm(n)
            n = n / ln if ln > 0 else n
            f.write(f"  facet normal {n[0]:.6e} {n[1]:.6e} {n[2]:.6e}\n")
            f.write("    outer loop\n")
            for v in (a, b, c):
                f.write(f"      vertex {v[0]:.6e} {v[1]:.6e} {v[2]:.6e}\n")
            f.write("    endloop\n  endfacet\n")
        f.write(f"endsolid {name}\n")


def main():
    p = argparse.ArgumentParser()
    p.add_argument("--chord", type=float, default=1.0)
    p.add_argument("--span", type=float, default=1.0)
    p.add_argument("--n-chord", type=int, default=120, help="stations per surface")
    p.add_argument("--n-span", type=int, default=30)
    p.add_argument("--out", default="naca0012.stl")
    p.add_argument("--ascii", action="store_true")
    a = p.parse_args()
    tris = build_mesh(a.chord, a.span, a.n_chord, a.n_span)
    if a.ascii:
        write_ascii_stl(tris, a.out)
    else:
        write_binary_stl(tris, a.out)
    print(f"wrote {a.out}: {len(tris)} triangles "
          f"(chord={a.chord}, span={a.span}, n_chord={a.n_chord}, n_span={a.n_span})")


if __name__ == "__main__":
    main()

naca0012_coarse_ascii.stl