NACA 0012 Airfoil - Lift and Surface Pressure

This notebook validates the lift coefficient \(C_L(\alpha)\) and the chordwise surface-pressure distribution \(C_p(x/c)\) of a wall-modelled LES over a NACA 0012 section against the NASA Turbulence Modeling Resource (TMR) SA reference and the Ladson (NASA TM 4074) experiment.

Note: the simulation runs at a cost-constrained \(\mathrm{Re}_c = 5\times10^5\) while the vendored reference is \(\mathrm{Re}_c = 6\times10^6\). The attached-flow lift-curve slope is nearly Reynolds-independent (the primary check); the maximum lift and stall angle both drop at the lower Reynolds number (\(C_{L,\max}\approx1.0\) near \(\alpha\approx11^\circ\) at \(5\times10^5\) vs \(\approx1.5\) near \(16^\circ\) at \(6\times10^6\)) and, having no matched-Reynolds reference here, are read only qualitatively; the zero-lift drag and the suction-peak \(C_p\) magnitude likewise differ. Lift and drag are recovered from the immersed-boundary spread forces; the surface pressure is read from the body probe. See the case README.md for the physics and references.

Configuration

The angle-of-attack sweep and the grid-convergence pair are parsed straight from the case YAML files, so the notebook always reflects what was run.

[ ]:

import math
import pathlib

import numpy as np
import pandas as pd

from nassu.cfg.model import ConfigScheme

CASE_DIR = pathlib.Path("validation/external_aero/04_airfoil_naca0012")
REF_DIR = CASE_DIR / "reference"

# Angle-of-attack sweep (sim_id 0..5 -> alpha = 0,4,8,10,12,15 deg).
sweep_cfgs = ConfigScheme.sim_cfgs_from_file_dct(str(CASE_DIR / "04_airfoil_naca0012.nassu.yaml"))
# Grid-convergence pair at alpha = 10 deg.
grid_cfgs = ConfigScheme.sim_cfgs_from_file_dct(
    str(CASE_DIR / "04.1_airfoil_naca0012_grid_convergence.nassu.yaml")
)


def alpha_deg_of(cfg) -> float:
    """Angle of attack (deg) recovered from the body rotation about y."""
    rot_y = cfg.domain.bodies["naca0012"].transformation.rotation[1]
    return math.degrees(rot_y)


def inlet_speed(cfg) -> float:
    """Free-stream speed U_inf (lattice units) from the UniformFlow inlet BC."""
    return max(bc.params["ux"] for bc in cfg.models.BC.BC_map if "ux" in bc.params)


def chord_lu_of(cfg) -> float:
    """Chord in lattice units = STL native chord (48) x the config scale.

    The case is parametric: the ``chord_lu`` resolution knob drives the STL
    ``scale``, so reading it back keeps the normalisation correct at any
    resolution instead of assuming a fixed 48.
    """
    return 48.0 * float(cfg.domain.bodies["naca0012"].transformation.scale[0])


# Case constants (lattice units), read back from the (parametric) config so they
# track the `chord_lu` resolution knob rather than assuming a fixed grid.
_cfg0 = sweep_cfgs[("airfoilNACA0012", 0)]
CHORD = chord_lu_of(_cfg0)
RHO_INF = 1.0
U_INF = inlet_speed(_cfg0)
# Statistics window: discard the ~50-CTU spin-up, average to the run end
# (1 CTU = chord / U_inf steps; matches the config's derived spin-up).
STATS_START = round(50 * CHORD / U_INF)

sweep = sorted(sweep_cfgs.values(), key=alpha_deg_of)
for c in sweep:
    print(f"alpha = {alpha_deg_of(c):5.1f} deg   tau = {c.models.LBM.tau}   n_steps = {c.n_steps}")
print(f"\nU_inf = {U_INF}  chord = {CHORD}  Ma_LBM = {math.sqrt(3) * U_INF:.3f}")
print(f"stats start = {STATS_START} steps")

Reference data

The TMR files are Tecplot ASCII with zone, t="..." blocks; the helper below splits a file into its zones. The CFD lift/drag bracket is a small CSV.

[ ]:

def read_tecplot_zones(path: pathlib.Path) -> dict[str, pd.DataFrame]:
    """Parse a TMR Tecplot-ASCII .dat file into {zone_title: DataFrame}."""
    cols, zones, title, rows = None, {}, None, []
    for line in pathlib.Path(path).read_text().splitlines():
        s = line.strip()
        if not s or s.startswith("#"):
            continue
        if s.lower().startswith("variables"):
            cols = [c.strip().strip('"') for c in s.split("=", 1)[1].split('","')]
            cols = [c.strip().strip('"') for c in cols]
            continue
        if s.lower().startswith("zone"):
            if title is not None and rows:
                zones[title] = pd.DataFrame(rows, columns=cols)
            title = s.split("t=", 1)[1].strip().strip('"') if "t=" in s else f"zone{len(zones)}"
            rows = []
            continue
        parts = s.split()
        try:
            rows.append([float(p) for p in parts])
        except ValueError:
            continue
    if title is not None and rows:
        zones[title] = pd.DataFrame(rows, columns=cols)
    return zones


clcd_cfd = pd.read_csv(REF_DIR / "clcd_cfd_tmr_sa.csv", comment="#")
clcd_exp = read_tecplot_zones(REF_DIR / "CLCD_Ladson_expdata.dat")
cp_cfd = read_tecplot_zones(REF_DIR / "n0012cp_cfl3d_sa.dat")
print("CFD bracket:\n", clcd_cfd)
print("\nExperimental zones:", list(clcd_exp))
print("CFD Cp zones:", list(cp_cfd))

Lift and drag from the IBM forces

The total hydrodynamic force on the body is the spreading force the solver applied to enforce no-slip, summed over the Lagrangian surface nodes and negated. nassu.viz.read_body_ibm_force reads the per-node export_IBM_nodes force and rescales each node from its level-local lattice units to global units by \((2^{-l})^2\) (the \(\mathrm{d}x^2\) area ratio) before summing, so the near-wall multi-level refinement does not double-count. The streamwise component \(F_x\) gives drag and the vertical component \(F_z\) gives lift; both are normalised by the planform reference area \(A = c\,b\) (chord times span) and time-averaged over the statistics window (step \(\geq 24000\)).

[ ]:

import nassu.viz as common

common.use_style()


def planform_area(cfg) -> float:
    """Reference area A = chord * span in lattice units (alpha = 0 projection)."""
    # CHORD tracks the chord_lu resolution knob; span = the periodic y-extent.
    span = float(cfg.domain.domain_size.y)
    return CHORD * span


def lift_drag_coeffs(cfg, u_inf: float = U_INF, rho_inf: float = RHO_INF, stats_start=STATS_START):
    """Time-averaged (Cl, Cd) from the direct IBM spreading forces.

    Reads the per-node ``export_IBM_nodes`` force via ``nassu.viz.read_body_ibm_force``,
    which rescales each Lagrangian node's level-local force to global units
    (``(2**-lvl)**2``) and sums it into the net hydrodynamic force ON the body -
    so a multi-level near-wall refinement does NOT double-count. The force is
    time-averaged over the statistics window (step >= ``stats_start``); drag is the
    streamwise (x) component and lift the vertical (z, signed) component, each
    normalised by ``q A = 0.5 rho_inf u_inf**2 (c b)``.
    """
    bif = common.read_body_ibm_force(cfg, "naca0012", start_step=stats_start)
    f_body = common.time_mean_ibm_force(bif)  # net force ON the body, (fx, fy, fz)
    q_area = 0.5 * rho_inf * u_inf**2 * planform_area(cfg)
    cd = float(f_body[0] / q_area)  # drag along +x
    cl = float(f_body[2] / q_area)  # lift along +z (signed)
    return cl, cd


# Cl(alpha), Cd(alpha) over the sweep.
sweep_alpha = np.array([alpha_deg_of(c) for c in sweep])
sim_cl, sim_cd = [], []
for c in sweep:
    cl, cd = lift_drag_coeffs(c)
    sim_cl.append(cl)
    sim_cd.append(cd)
sim_cl, sim_cd = np.array(sim_cl), np.array(sim_cd)

import matplotlib.pyplot as plt

# Drag polar as a check: Nassu Cd(alpha) against the CFL3D / FUN3D drag bracket
# (Re=6e6, loaded in the reference cell). The absolute level differs (higher skin
# friction at the run's Re=5e5), so this is a consistency plot, not a matched-Re
# target.
fig, ax = common.fig_single()
ax.plot(sweep_alpha, sim_cd, **common.markers.sim(shape="o"), label="Nassu, Re=5e5")
ax.plot(
    clcd_cfd["alpha_deg"],
    clcd_cfd["cd_cfl3d"],
    **common.markers.exp_line(linestyle="-"),
    label="CFL3D (SA), Re=6e6",
)
ax.plot(
    clcd_cfd["alpha_deg"],
    clcd_cfd["cd_fun3d"],
    **common.markers.exp_line(linestyle="--"),
    label="FUN3D (SA), Re=6e6",
)
ax.set_xlabel(r"$\alpha$ (deg)")
ax.set_ylabel(r"$C_D$")
ax.legend()
plt.tight_layout()
plt.show(fig)

# Tabular supplement (below the figure, never the sole output).
pd.DataFrame({"alpha_deg": sweep_alpha, "Cl": sim_cl, "Cd": sim_cd})

Lift curve \(C_L(\alpha)\)

[ ]:

fig, ax = plt.subplots()
# Experimental scatter (all grit zones), source named once in the legend.
for name, df in clcd_exp.items():
    ax.plot(
        df["alpha, deg"],
        df["cl"],
        **common.markers.exp(shape="o"),
        label="Ladson (NASA TM 4074), Re=6e6" if name == list(clcd_exp)[0] else None,
    )
# CFD bracket.
ax.plot(
    clcd_cfd["alpha_deg"],
    clcd_cfd["cl_cfl3d"],
    **common.markers.exp_line(linestyle="-"),
    label="CFL3D (SA), Re=6e6",
)
ax.plot(
    clcd_cfd["alpha_deg"],
    clcd_cfd["cl_fun3d"],
    **common.markers.exp_line(linestyle="--"),
    label="FUN3D (SA), Re=6e6",
)
# Thin-airfoil slope reference.
aa = np.linspace(0, 12, 50)
ax.plot(
    aa,
    2 * np.pi * np.deg2rad(aa),
    color=common.colors.refline,
    lw=1,
    ls=":",
    label=r"$2\pi\,\alpha$ (thin airfoil)",
)
# Nassu.
ax.plot(sweep_alpha, sim_cl, **common.markers.sim(), label="Nassu, Re=5e5")
ax.set_xlabel(r"$\alpha$ (deg)")
ax.set_ylabel(r"$C_L$")
ax.set_xlim(-1, 20)
ax.legend()
plt.tight_layout()
plt.show(fig)

Surface pressure \(C_p(x/c)\) at \(\alpha = 10^\circ\)

The body probe stores the near-wall density; \(C_p = (\rho - \rho_\infty)/(\tfrac{1}{2}\rho_\infty U_\infty^2 \cdot 3)\) since \(p = c_s^2\rho\) and \(c_s^2 = 1/3\). Points are split into upper/lower surface (by z, the lift direction) and plotted against \(x/c\).

[ ]:

def surface_cp(cfg, u_inf: float = U_INF, rho_inf: float = RHO_INF, t_start: int = STATS_START):
    """Mean Cp at the body-surface probe points, returned with chordwise x/c.

    Reads the ``body_pressure`` surface-density series
    (``cfg.output.exports["body_pressure"].series.bodies["naca0012"].inst``); with
    ``p = c_s^2 rho`` and ``c_s^2 = 1/3``, ``Cp = (rho - rho_inf) / (0.5 rho_inf
    u_inf^2 * 3)``. Points are split into upper/lower surface by z (the lift
    direction) and mapped to ``x/c`` about the leading edge.
    """
    hs = cfg.output.exports["body_pressure"].series.bodies["naca0012"].inst
    pts = pd.read_csv(hs.points_filename)
    rho = hs.read_full_data("rho")
    rho = rho[rho["time_step"] >= t_start].drop(columns="time_step")
    # read_full_data columns follow the points-file row order, so the mean over
    # time aligns element-wise with the pts rows.
    cp = ((rho - rho_inf) / (0.5 * rho_inf * u_inf**2 * 3.0)).mean().to_numpy()
    x_le = pts["x"].min()
    xc = (pts["x"].to_numpy() - x_le) / CHORD
    upper = pts["z"].to_numpy() >= pts["z"].mean()
    return xc, cp, upper


cfg10 = next(c for c in sweep if abs(alpha_deg_of(c) - 10.0) < 0.5)
xc, cp, upper = surface_cp(cfg10)

fig, ax = plt.subplots()
ax.plot(xc[upper], cp[upper], **common.markers.sim(shape="o"), label="Nassu (suction), Re=5e5")
ax.plot(xc[~upper], cp[~upper], **common.markers.sim(shape="s"), label="Nassu (pressure), Re=5e5")
cp_ref = cp_cfd.get("alpha=10")
if cp_ref is not None:
    ax.plot(cp_ref["x"], cp_ref["cp"], **common.markers.exp_line(), label="CFL3D (SA), Re=6e6")
ax.invert_yaxis()  # -Cp up, aerodynamics convention
ax.set_xlabel(r"$x/c$")
ax.set_ylabel(r"$C_p$")
ax.set_title(r"$\alpha = 10^\circ$")
ax.legend()
plt.tight_layout()
plt.show(fig)

Grid convergence at \(\alpha = 10^\circ\)

The lift from the coarse (level 3) and fine (level 4 near-wall shell) runs brackets the resolution sensitivity at the headline angle.

[ ]:

rows = []
for c in grid_cfgs.values():
    cl, cd = lift_drag_coeffs(c)
    rows.append({"run": c.name, "Cl": cl, "Cd": cd})
grid_df = pd.DataFrame(rows).sort_values("run").reset_index(drop=True)

cl_ref = float(clcd_cfd.loc[clcd_cfd["alpha_deg"] == 10, "cl_cfl3d"].iloc[0])
cd_ref = float(clcd_cfd.loc[clcd_cfd["alpha_deg"] == 10, "cd_cfl3d"].iloc[0])

# Grid-convergence bracket at alpha=10 deg as bars: coarse vs fine near-wall
# resolution, with the CFL3D SA reference (Re=6e6) drawn as a line. The lift
# should change little between resolutions (the resolution-sensitivity bracket);
# the reference is a different-Re anchor, not a matched-Re target.
labels = grid_df["run"].tolist()
xb = np.arange(len(labels))
fig, ax = common.fig_double()
for axq, q, ref, ttl in ((ax[0], "Cl", cl_ref, r"$C_L$"), (ax[1], "Cd", cd_ref, r"$C_D$")):
    bars = axq.bar(xb, grid_df[q], color=common.colors.sim, width=0.6, label="Nassu, Re=5e5")
    axq.bar_label(bars, fmt="%.4f", padding=2, fontsize=9)
    axq.axhline(ref, color=common.colors.exp, ls="--", lw=1.5, label="CFL3D (SA), Re=6e6")
    axq.set_xticks(xb)
    axq.set_xticklabels(labels, rotation=20, ha="right")
    common.bar_axis(axq)
    axq.set_ylabel(ttl)
    axq.set_title(f"{ttl} at " + r"$\alpha = 10^\circ$")
    axq.legend()
plt.tight_layout()
plt.show(fig)

grid_df

Flow field

Instantaneous velocity magnitude on the mid-span plane (plane_series.mid_span, y-normal at the spanwise centre) at \(\alpha = 10^\circ\), framed from the airfoil geometry: chordwise flow over the section, the suction-side acceleration and the wake.

[ ]:

common.enable_offscreen()

# alpha = 10 deg section (the headline angle). The airfoil STL span overhangs the
# periodic y-domain, so the mid-span plane is pinned explicitly at y = 12.
body, geom = common.read_body(cfg10, "naca0012")
view = common.frame_body(geom, "y", slice_coord=12.0, downstream=2.0, half_extent=3.0)

panel = common.Panel(
    r"$\alpha = 10^\circ$",
    common.PlaneSource.from_cfg(cfg10, series="plane_series", plane="mid_span"),
    view,
)
steps = [panel.source.steps[-1]]
plotter = common.render_grid(
    [panel],
    steps=steps,
    scalar="u_mag",
    cmap="viridis",
    clim=(0.0, 0.15),
    bar_title="|u|",
    bodies=[body],
)
plotter.show()

Summary

Success criteria for this case:

  • \(C_L(\alpha)\) tracks the CFL3D/FUN3D bracket and the Ladson scatter in the linear (attached) regime, with slope close to \(2\pi\) per radian (\(\approx 0.11\) per deg) - the Reynolds-robust primary check.

  • The stall onset (\(C_{L,\max}\) and its angle) is read only qualitatively: both fall at the run’s \(\mathrm{Re}_c = 5\times10^5\) (expect \(C_{L,\max}\approx1.0\) near \(\alpha\approx11^\circ\)) and are not directly comparable to the \(6\times10^6\) reference.

  • Drag reported as a convergence check (the \(6\times10^6\) value \(C_{D,0}\approx 0.0082\) is not the target at \(5\times10^5\), where skin friction is higher).

  • \(C_p(x/c)\) at \(\alpha=10^\circ\) matches the CFD reference shape (suction peak and stagnation captured); the suction-peak magnitude sits below the \(6\times10^6\) reference at the lower Reynolds number.

  • Lift changes little between the coarse and fine near-wall resolutions.

GPU runs are required to populate the cells above.

Version

[ ]:

sim_info = sweep[0].output.read_info()
print("Version:", sim_info["version"])
print("Commit hash:", sim_info["commit"])

Configuration

[ ]:

from IPython.display import Code

Code(filename=str(CASE_DIR / "04_airfoil_naca0012.nassu.yaml"), language="yaml")