{ "cells": [ { "cell_type": "markdown", "id": "5b164eeb", "metadata": {}, "source": "# NACA 0012 Airfoil - Lift and Surface Pressure\n\nThis notebook validates the lift coefficient $C_L(\\alpha)$ and the chordwise\nsurface-pressure distribution $C_p(x/c)$ of a wall-modelled LES over a NACA 0012\nsection against the NASA Turbulence Modeling Resource (TMR) SA reference and the\nLadson (NASA TM 4074) experiment.\n\n**Note:** the simulation runs at a cost-constrained $\\mathrm{Re}_c = 5\\times10^5$\nwhile the vendored reference is $\\mathrm{Re}_c = 6\\times10^6$. The attached-flow\nlift-curve slope is nearly Reynolds-independent (the primary check); the maximum\nlift and stall angle both drop at the lower Reynolds number ($C_{L,\\max}\\approx1.0$\nnear $\\alpha\\approx11^\\circ$ at $5\\times10^5$ vs $\\approx1.5$ near $16^\\circ$ at\n$6\\times10^6$) and, having no matched-Reynolds reference here, are read only\nqualitatively; the zero-lift drag and the suction-peak $C_p$ magnitude likewise\ndiffer. Lift and drag are recovered from the immersed-boundary spread forces; the\nsurface pressure is read from the body probe. See the case `README.md` for the\nphysics and references." }, { "cell_type": "markdown", "id": "4636564c", "metadata": {}, "source": [ "## Configuration\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "85dfcf52", "metadata": {}, "outputs": [], "source": [ "import math\n", "import pathlib\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from nassu.cfg.model import ConfigScheme\n", "\n", "CASE_DIR = pathlib.Path(\"validation/external_aero/04_airfoil_naca0012\")\n", "REF_DIR = CASE_DIR / \"reference\"\n", "\n", "# Angle-of-attack sweep (sim_id 0..5 -> alpha = 0,4,8,10,12,15 deg).\n", "sweep_cfgs = ConfigScheme.sim_cfgs_from_file_dct(str(CASE_DIR / \"04_airfoil_naca0012.nassu.yaml\"))\n", "# Grid-convergence pair at alpha = 10 deg.\n", "grid_cfgs = ConfigScheme.sim_cfgs_from_file_dct(\n", " str(CASE_DIR / \"04.1_airfoil_naca0012_grid_convergence.nassu.yaml\")\n", ")\n", "\n", "\n", "def alpha_deg_of(cfg) -> float:\n", " \"\"\"Angle of attack (deg) recovered from the body rotation about y.\"\"\"\n", " rot_y = cfg.domain.bodies[\"naca0012\"].transformation.rotation[1]\n", " return math.degrees(rot_y)\n", "\n", "\n", "def inlet_speed(cfg) -> float:\n", " \"\"\"Free-stream speed U_inf (lattice units) from the UniformFlow inlet BC.\"\"\"\n", " return max(bc.params[\"ux\"] for bc in cfg.models.BC.BC_map if \"ux\" in bc.params)\n", "\n", "\n", "def chord_lu_of(cfg) -> float:\n", " \"\"\"Chord in lattice units = STL native chord (48) x the config scale.\n", "\n", " The case is parametric: the ``chord_lu`` resolution knob drives the STL\n", " ``scale``, so reading it back keeps the normalisation correct at any\n", " resolution instead of assuming a fixed 48.\n", " \"\"\"\n", " return 48.0 * float(cfg.domain.bodies[\"naca0012\"].transformation.scale[0])\n", "\n", "\n", "# Case constants (lattice units), read back from the (parametric) config so they\n", "# track the `chord_lu` resolution knob rather than assuming a fixed grid.\n", "_cfg0 = sweep_cfgs[(\"airfoilNACA0012\", 0)]\n", "CHORD = chord_lu_of(_cfg0)\n", "RHO_INF = 1.0\n", "U_INF = inlet_speed(_cfg0)\n", "# Statistics window: discard the ~50-CTU spin-up, average to the run end\n", "# (1 CTU = chord / U_inf steps; matches the config's derived spin-up).\n", "STATS_START = round(50 * CHORD / U_INF)\n", "\n", "sweep = sorted(sweep_cfgs.values(), key=alpha_deg_of)\n", "for c in sweep:\n", " print(f\"alpha = {alpha_deg_of(c):5.1f} deg tau = {c.models.LBM.tau} n_steps = {c.n_steps}\")\n", "print(f\"\\nU_inf = {U_INF} chord = {CHORD} Ma_LBM = {math.sqrt(3) * U_INF:.3f}\")\n", "print(f\"stats start = {STATS_START} steps\")" ] }, { "cell_type": "markdown", "id": "e5cc6c4e", "metadata": {}, "source": [ "## Reference data\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "85e124ca", "metadata": {}, "outputs": [], "source": [ "def read_tecplot_zones(path: pathlib.Path) -> dict[str, pd.DataFrame]:\n", " \"\"\"Parse a TMR Tecplot-ASCII .dat file into {zone_title: DataFrame}.\"\"\"\n", " cols, zones, title, rows = None, {}, None, []\n", " for line in pathlib.Path(path).read_text().splitlines():\n", " s = line.strip()\n", " if not s or s.startswith(\"#\"):\n", " continue\n", " if s.lower().startswith(\"variables\"):\n", " cols = [c.strip().strip('\"') for c in s.split(\"=\", 1)[1].split('\",\"')]\n", " cols = [c.strip().strip('\"') for c in cols]\n", " continue\n", " if s.lower().startswith(\"zone\"):\n", " if title is not None and rows:\n", " zones[title] = pd.DataFrame(rows, columns=cols)\n", " title = s.split(\"t=\", 1)[1].strip().strip('\"') if \"t=\" in s else f\"zone{len(zones)}\"\n", " rows = []\n", " continue\n", " parts = s.split()\n", " try:\n", " rows.append([float(p) for p in parts])\n", " except ValueError:\n", " continue\n", " if title is not None and rows:\n", " zones[title] = pd.DataFrame(rows, columns=cols)\n", " return zones\n", "\n", "\n", "clcd_cfd = pd.read_csv(REF_DIR / \"clcd_cfd_tmr_sa.csv\", comment=\"#\")\n", "clcd_exp = read_tecplot_zones(REF_DIR / \"CLCD_Ladson_expdata.dat\")\n", "cp_cfd = read_tecplot_zones(REF_DIR / \"n0012cp_cfl3d_sa.dat\")\n", "print(\"CFD bracket:\\n\", clcd_cfd)\n", "print(\"\\nExperimental zones:\", list(clcd_exp))\n", "print(\"CFD Cp zones:\", list(cp_cfd))" ] }, { "cell_type": "markdown", "id": "0a85eaf3", "metadata": {}, "source": "## Lift and drag from the IBM forces\n\nThe total hydrodynamic force on the body is the spreading force the solver applied\nto enforce no-slip, summed over the Lagrangian surface nodes and negated.\n`nassu.viz.read_body_ibm_force` reads the per-node `export_IBM_nodes` force and\nrescales each node from its level-local lattice units to global units by\n$(2^{-l})^2$ (the $\\mathrm{d}x^2$ area ratio) before summing, so the near-wall\nmulti-level refinement does not double-count. The streamwise component $F_x$ gives\ndrag and the vertical component $F_z$ gives lift; both are normalised by the\nplanform reference area $A = c\\,b$ (chord times span) and time-averaged over the\nstatistics window (step $\\geq 24000$)." }, { "cell_type": "code", "execution_count": null, "id": "3b590190", "metadata": {}, "outputs": [], "source": [ "import nassu.viz as common\n", "\n", "common.use_style()\n", "\n", "\n", "def planform_area(cfg) -> float:\n", " \"\"\"Reference area A = chord * span in lattice units (alpha = 0 projection).\"\"\"\n", " # CHORD tracks the chord_lu resolution knob; span = the periodic y-extent.\n", " span = float(cfg.domain.domain_size.y)\n", " return CHORD * span\n", "\n", "\n", "def lift_drag_coeffs(cfg, u_inf: float = U_INF, rho_inf: float = RHO_INF, stats_start=STATS_START):\n", " \"\"\"Time-averaged (Cl, Cd) from the direct IBM spreading forces.\n", "\n", " Reads the per-node ``export_IBM_nodes`` force via ``nassu.viz.read_body_ibm_force``,\n", " which rescales each Lagrangian node's level-local force to global units\n", " (``(2**-lvl)**2``) and sums it into the net hydrodynamic force ON the body -\n", " so a multi-level near-wall refinement does NOT double-count. The force is\n", " time-averaged over the statistics window (step >= ``stats_start``); drag is the\n", " streamwise (x) component and lift the vertical (z, signed) component, each\n", " normalised by ``q A = 0.5 rho_inf u_inf**2 (c b)``.\n", " \"\"\"\n", " bif = common.read_body_ibm_force(cfg, \"naca0012\", start_step=stats_start)\n", " f_body = common.time_mean_ibm_force(bif) # net force ON the body, (fx, fy, fz)\n", " q_area = 0.5 * rho_inf * u_inf**2 * planform_area(cfg)\n", " cd = float(f_body[0] / q_area) # drag along +x\n", " cl = float(f_body[2] / q_area) # lift along +z (signed)\n", " return cl, cd\n", "\n", "\n", "# Cl(alpha), Cd(alpha) over the sweep.\n", "sweep_alpha = np.array([alpha_deg_of(c) for c in sweep])\n", "sim_cl, sim_cd = [], []\n", "for c in sweep:\n", " cl, cd = lift_drag_coeffs(c)\n", " sim_cl.append(cl)\n", " sim_cd.append(cd)\n", "sim_cl, sim_cd = np.array(sim_cl), np.array(sim_cd)\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "# Drag polar as a check: Nassu Cd(alpha) against the CFL3D / FUN3D drag bracket\n", "# (Re=6e6, loaded in the reference cell). The absolute level differs (higher skin\n", "# friction at the run's Re=5e5), so this is a consistency plot, not a matched-Re\n", "# target.\n", "fig, ax = common.fig_single()\n", "ax.plot(sweep_alpha, sim_cd, **common.markers.sim(shape=\"o\"), label=\"Nassu, Re=5e5\")\n", "ax.plot(\n", " clcd_cfd[\"alpha_deg\"],\n", " clcd_cfd[\"cd_cfl3d\"],\n", " **common.markers.exp_line(linestyle=\"-\"),\n", " label=\"CFL3D (SA), Re=6e6\",\n", ")\n", "ax.plot(\n", " clcd_cfd[\"alpha_deg\"],\n", " clcd_cfd[\"cd_fun3d\"],\n", " **common.markers.exp_line(linestyle=\"--\"),\n", " label=\"FUN3D (SA), Re=6e6\",\n", ")\n", "ax.set_xlabel(r\"$\\alpha$ (deg)\")\n", "ax.set_ylabel(r\"$C_D$\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show(fig)\n", "\n", "# Tabular supplement (below the figure, never the sole output).\n", "pd.DataFrame({\"alpha_deg\": sweep_alpha, \"Cl\": sim_cl, \"Cd\": sim_cd})" ] }, { "cell_type": "markdown", "id": "0d75c813", "metadata": {}, "source": [ "### Lift curve $C_L(\\alpha)$" ] }, { "cell_type": "code", "execution_count": null, "id": "c918c6f0", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "# Experimental scatter (all grit zones), source named once in the legend.\n", "for name, df in clcd_exp.items():\n", " ax.plot(\n", " df[\"alpha, deg\"],\n", " df[\"cl\"],\n", " **common.markers.exp(shape=\"o\"),\n", " label=\"Ladson (NASA TM 4074), Re=6e6\" if name == list(clcd_exp)[0] else None,\n", " )\n", "# CFD bracket.\n", "ax.plot(\n", " clcd_cfd[\"alpha_deg\"],\n", " clcd_cfd[\"cl_cfl3d\"],\n", " **common.markers.exp_line(linestyle=\"-\"),\n", " label=\"CFL3D (SA), Re=6e6\",\n", ")\n", "ax.plot(\n", " clcd_cfd[\"alpha_deg\"],\n", " clcd_cfd[\"cl_fun3d\"],\n", " **common.markers.exp_line(linestyle=\"--\"),\n", " label=\"FUN3D (SA), Re=6e6\",\n", ")\n", "# Thin-airfoil slope reference.\n", "aa = np.linspace(0, 12, 50)\n", "ax.plot(\n", " aa,\n", " 2 * np.pi * np.deg2rad(aa),\n", " color=common.colors.refline,\n", " lw=1,\n", " ls=\":\",\n", " label=r\"$2\\pi\\,\\alpha$ (thin airfoil)\",\n", ")\n", "# Nassu.\n", "ax.plot(sweep_alpha, sim_cl, **common.markers.sim(), label=\"Nassu, Re=5e5\")\n", "ax.set_xlabel(r\"$\\alpha$ (deg)\")\n", "ax.set_ylabel(r\"$C_L$\")\n", "ax.set_xlim(-1, 20)\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show(fig)" ] }, { "cell_type": "markdown", "id": "b710d960", "metadata": {}, "source": [ "### Surface pressure $C_p(x/c)$ at $\\alpha = 10^\\circ$\n", "\n", "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$." ] }, { "cell_type": "code", "execution_count": null, "id": "8e919ffb", "metadata": {}, "outputs": [], "source": [ "def surface_cp(cfg, u_inf: float = U_INF, rho_inf: float = RHO_INF, t_start: int = STATS_START):\n", " \"\"\"Mean Cp at the body-surface probe points, returned with chordwise x/c.\n", "\n", " Reads the ``body_pressure`` surface-density series\n", " (``cfg.output.exports[\"body_pressure\"].series.bodies[\"naca0012\"].inst``); with\n", " ``p = c_s^2 rho`` and ``c_s^2 = 1/3``, ``Cp = (rho - rho_inf) / (0.5 rho_inf\n", " u_inf^2 * 3)``. Points are split into upper/lower surface by z (the lift\n", " direction) and mapped to ``x/c`` about the leading edge.\n", " \"\"\"\n", " hs = cfg.output.exports[\"body_pressure\"].series.bodies[\"naca0012\"].inst\n", " pts = pd.read_csv(hs.points_filename)\n", " rho = hs.read_full_data(\"rho\")\n", " rho = rho[rho[\"time_step\"] >= t_start].drop(columns=\"time_step\")\n", " # read_full_data columns follow the points-file row order, so the mean over\n", " # time aligns element-wise with the pts rows.\n", " cp = ((rho - rho_inf) / (0.5 * rho_inf * u_inf**2 * 3.0)).mean().to_numpy()\n", " x_le = pts[\"x\"].min()\n", " xc = (pts[\"x\"].to_numpy() - x_le) / CHORD\n", " upper = pts[\"z\"].to_numpy() >= pts[\"z\"].mean()\n", " return xc, cp, upper\n", "\n", "\n", "cfg10 = next(c for c in sweep if abs(alpha_deg_of(c) - 10.0) < 0.5)\n", "xc, cp, upper = surface_cp(cfg10)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(xc[upper], cp[upper], **common.markers.sim(shape=\"o\"), label=\"Nassu (suction), Re=5e5\")\n", "ax.plot(xc[~upper], cp[~upper], **common.markers.sim(shape=\"s\"), label=\"Nassu (pressure), Re=5e5\")\n", "cp_ref = cp_cfd.get(\"alpha=10\")\n", "if cp_ref is not None:\n", " ax.plot(cp_ref[\"x\"], cp_ref[\"cp\"], **common.markers.exp_line(), label=\"CFL3D (SA), Re=6e6\")\n", "ax.invert_yaxis() # -Cp up, aerodynamics convention\n", "ax.set_xlabel(r\"$x/c$\")\n", "ax.set_ylabel(r\"$C_p$\")\n", "ax.set_title(r\"$\\alpha = 10^\\circ$\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show(fig)" ] }, { "cell_type": "markdown", "id": "61efd96a", "metadata": {}, "source": [ "## Grid convergence at $\\alpha = 10^\\circ$\n", "\n", "The lift from the coarse (level 3) and fine (level 4 near-wall shell) runs brackets the resolution sensitivity at the headline angle." ] }, { "cell_type": "code", "execution_count": null, "id": "0e877067", "metadata": {}, "outputs": [], "source": [ "rows = []\n", "for c in grid_cfgs.values():\n", " cl, cd = lift_drag_coeffs(c)\n", " rows.append({\"run\": c.name, \"Cl\": cl, \"Cd\": cd})\n", "grid_df = pd.DataFrame(rows).sort_values(\"run\").reset_index(drop=True)\n", "\n", "cl_ref = float(clcd_cfd.loc[clcd_cfd[\"alpha_deg\"] == 10, \"cl_cfl3d\"].iloc[0])\n", "cd_ref = float(clcd_cfd.loc[clcd_cfd[\"alpha_deg\"] == 10, \"cd_cfl3d\"].iloc[0])\n", "\n", "# Grid-convergence bracket at alpha=10 deg as bars: coarse vs fine near-wall\n", "# resolution, with the CFL3D SA reference (Re=6e6) drawn as a line. The lift\n", "# should change little between resolutions (the resolution-sensitivity bracket);\n", "# the reference is a different-Re anchor, not a matched-Re target.\n", "labels = grid_df[\"run\"].tolist()\n", "xb = np.arange(len(labels))\n", "fig, ax = common.fig_double()\n", "for axq, q, ref, ttl in ((ax[0], \"Cl\", cl_ref, r\"$C_L$\"), (ax[1], \"Cd\", cd_ref, r\"$C_D$\")):\n", " bars = axq.bar(xb, grid_df[q], color=common.colors.sim, width=0.6, label=\"Nassu, Re=5e5\")\n", " axq.bar_label(bars, fmt=\"%.4f\", padding=2, fontsize=9)\n", " axq.axhline(ref, color=common.colors.exp, ls=\"--\", lw=1.5, label=\"CFL3D (SA), Re=6e6\")\n", " axq.set_xticks(xb)\n", " axq.set_xticklabels(labels, rotation=20, ha=\"right\")\n", " common.bar_axis(axq)\n", " axq.set_ylabel(ttl)\n", " axq.set_title(f\"{ttl} at \" + r\"$\\alpha = 10^\\circ$\")\n", " axq.legend()\n", "plt.tight_layout()\n", "plt.show(fig)\n", "\n", "grid_df" ] }, { "cell_type": "markdown", "id": "995242f5", "metadata": {}, "source": "## Flow field\n\nInstantaneous velocity magnitude on the mid-span plane (`plane_series.mid_span`,\ny-normal at the spanwise centre) at $\\alpha = 10^\\circ$, framed from the airfoil\ngeometry: chordwise flow over the section, the suction-side acceleration and the\nwake." }, { "cell_type": "code", "execution_count": null, "id": "4ef94317", "metadata": {}, "outputs": [], "source": [ "common.enable_offscreen()\n", "\n", "# alpha = 10 deg section (the headline angle). The airfoil STL span overhangs the\n", "# periodic y-domain, so the mid-span plane is pinned explicitly at y = 12.\n", "body, geom = common.read_body(cfg10, \"naca0012\")\n", "view = common.frame_body(geom, \"y\", slice_coord=12.0, downstream=2.0, half_extent=3.0)\n", "\n", "panel = common.Panel(\n", " r\"$\\alpha = 10^\\circ$\",\n", " common.PlaneSource.from_cfg(cfg10, series=\"plane_series\", plane=\"mid_span\"),\n", " view,\n", ")\n", "steps = [panel.source.steps[-1]]\n", "plotter = common.render_grid(\n", " [panel],\n", " steps=steps,\n", " scalar=\"u_mag\",\n", " cmap=\"viridis\",\n", " clim=(0.0, 0.15),\n", " bar_title=\"|u|\",\n", " bodies=[body],\n", ")\n", "plotter.show()" ] }, { "cell_type": "markdown", "id": "3e16a257", "metadata": {}, "source": "## Summary\n\nSuccess criteria for this case:\n\n- $C_L(\\alpha)$ tracks the CFL3D/FUN3D bracket and the Ladson scatter in the\n linear (attached) regime, with slope close to $2\\pi$ per radian ($\\approx 0.11$\n per deg) - the Reynolds-robust primary check.\n- The stall onset ($C_{L,\\max}$ and its angle) is read only qualitatively: both\n fall at the run's $\\mathrm{Re}_c = 5\\times10^5$ (expect $C_{L,\\max}\\approx1.0$\n near $\\alpha\\approx11^\\circ$) and are not directly comparable to the\n $6\\times10^6$ reference.\n- Drag reported as a convergence check (the $6\\times10^6$ value\n $C_{D,0}\\approx 0.0082$ is not the target at $5\\times10^5$, where skin friction\n is higher).\n- $C_p(x/c)$ at $\\alpha=10^\\circ$ matches the CFD reference shape (suction peak and\n stagnation captured); the suction-peak magnitude sits below the $6\\times10^6$\n reference at the lower Reynolds number.\n- Lift changes little between the coarse and fine near-wall resolutions.\n\nGPU runs are required to populate the cells above." }, { "cell_type": "markdown", "id": "a4522e4d", "metadata": {}, "source": [ "## Version" ] }, { "cell_type": "code", "execution_count": null, "id": "32b7b518", "metadata": {}, "outputs": [], "source": [ "sim_info = sweep[0].output.read_info()\n", "print(\"Version:\", sim_info[\"version\"])\n", "print(\"Commit hash:\", sim_info[\"commit\"])" ] }, { "cell_type": "markdown", "id": "f55d501c", "metadata": {}, "source": [ "## Configuration" ] }, { "cell_type": "code", "execution_count": null, "id": "30c6776a", "metadata": {}, "outputs": [], "source": [ "from IPython.display import Code\n", "\n", "Code(filename=str(CASE_DIR / \"04_airfoil_naca0012.nassu.yaml\"), language=\"yaml\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }