zacheryasc/airfRANS-model-exploration
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
airfRANS-model-exploration/notebooks/airfrans_simulation_deep_dive.ipynb

1218 lines
53 KiB
Text
Raw Normal View History

init
2026-07-21 08:32:30 +00:00
{
"cells": [
{
"cell_type": "markdown",
"id": "f63c5c8a",
"metadata": {},
"source": [
"# AirfRANS simulation deep dive\n",
"\n",
"This notebook is a concrete map from the AirfRANS paper to the files in this repo's local raw OpenFOAM subset.\n",
"\n",
"Sources used while writing it:\n",
"\n",
"- Paper: Bonnet et al., **AirfRANS: High Fidelity Computational Fluid Dynamics Dataset for Approximating Reynolds-Averaged Navier-Stokes Solutions**, arXiv:2212.07564, v3.\n",
"- AirfRANS docs: introduction, simulation, dataset pages.\n",
"- `Extrality/NACA_simulation`: `README.md`, `params.yaml`, `main.py`, `dataset_generator.py`, `simulation_generator.py`.\n",
"\n",
"The paper's simulation is not a toy equation on a grid. It is a generated CFD case:\n",
"\n",
"```text\n",
"sample airfoil shape + Reynolds + angle of attack\n",
" ↓\n",
"convert Reynolds to freestream speed U_inf because chord = 1 m and ν is fixed\n",
" ↓\n",
"generate NACA 4/5-digit airfoil coordinates\n",
" ↓\n",
"generate an OpenFOAM C-grid blockMesh around the airfoil\n",
" ↓\n",
"run steady incompressible RANS with k-ω SST turbulence model using simpleFoam/SIMPLEC\n",
" ↓\n",
"postprocess fields, wall stresses, lift/drag/moment histories, VTK files\n",
" ↓\n",
"crop/slice/select fields for ML surrogate-model training\n",
"```\n",
"\n",
"You can run the code cells top-to-bottom. Change `SIM_INDEX` later to inspect another case.\n"
]
},
{
"cell_type": "markdown",
"id": "7caaa053",
"metadata": {},
"source": [
"## 1. Domain translation: what problem is being simulated?\n",
"\n",
"AirfRANS = airfoil data built from **Reynolds-Averaged Navier-Stokes** simulations.\n",
"\n",
"The simulated object is a 2D airfoil cross-section with chord length $c = 1\\ \\mathrm{m}$, not a full 3D wing. The surrounding fluid is air at sea-level-like conditions and $T = 298.15\\ \\mathrm{K}$. The flow is **subsonic** and treated as **incompressible** because the paper limits the regime to roughly $\\mathrm{Ma} < 0.3$.\n",
"\n",
"Each simulation solves for a steady **mean** flow around one airfoil. \"Mean\" matters because the real turbulent flow contains rapidly fluctuating eddies; RANS replaces those instant-by-instant fluctuations with time-averaged fields that are useful for engineering design.\n",
"\n",
"- **Mean velocity field** $\\bar u(x,y) = (\\bar u_x, \\bar u_y)$. This gives the average flow speed and direction at every point in the 2D domain. It shows acceleration over the suction side, stagnation near the leading edge, wake deficit behind the airfoil, and any separated/recirculating region. Velocity is the primary state variable: its gradients set viscous shear, and its deflection is directly related to lift.\n",
"- **Reduced pressure** $\\bar p/\\rho$. OpenFOAM's incompressible pressure variable is pressure divided by density, so its units are $\\mathrm{m^2/s^2}$ rather than Pa. This is convenient because the incompressible momentum equation uses pressure as a velocity-squared quantity. Pressure differences over the airfoil are crucial: low pressure on the upper/suction side and higher pressure on the lower/pressure side produce most of the lift, while pressure imbalance in the streamwise direction contributes to pressure drag.\n",
"- **Turbulent viscosity** $\\nu_t(x,y)$ from the turbulence model. This is not a material property like air's molecular viscosity $\\nu$; it is the $k$-$\\omega$ SST model's estimate of extra momentum mixing caused by unresolved turbulent eddies. Large $\\nu_t$ usually marks shear layers, wakes, and boundary-layer regions where turbulence transports momentum much faster than molecular diffusion alone. For surrogate models, $\\nu_t$ is important because it encodes where the closure model thinks turbulence is controlling the mean flow.\n",
"- **Surface pressure and viscous/shear stresses** on the airfoil wall. The volume fields describe the flow around the body, but forces come from stresses acting on the surface. Pressure acts normal to the wall; viscous shear acts tangentially through the near-wall velocity gradient. Integrating these wall quantities gives the aerodynamic loads.\n",
"- **Integrated force coefficients**: drag $C_D$, lift $C_L$, and moment coefficient. These nondimensionalize the total forces and torque by dynamic pressure and chord-based reference scales, so cases with different freestream speeds and airfoil shapes can be compared directly. $C_L$ measures useful vertical force, $C_D$ measures aerodynamic resistance, and the moment coefficient measures pitching tendency/stability.\n",
"\n",
"The original full dataset has 1000 simulations. This repo currently has a 50-case raw OpenFOAM subset under `data/raw/OF_dataset/`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1346d8f4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.234873Z",
"iopub.status.busy": "2026-07-20T19:11:12.234506Z",
"iopub.status.idle": "2026-07-20T19:11:12.518603Z",
"shell.execute_reply": "2026-07-20T19:11:12.518263Z"
}
},
"outputs": [],
"source": [
"from pathlib import Path\n",
"import gzip\n",
"import json\n",
"import math\n",
"import re\n",
"import sys\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"\n",
"def find_repo_root(start: Path) -> Path:\n",
" for candidate in (start, *start.parents):\n",
" if (candidate / \"pyproject.toml\").exists() and (candidate / \"src\" / \"airfrans_frontier\").exists():\n",
" return candidate\n",
" raise RuntimeError(f\"Could not find repository root from {start}\")\n",
"\n",
"\n",
"REPO_ROOT = find_repo_root(Path.cwd())\n",
"SRC = REPO_ROOT / \"src\"\n",
"if str(SRC) not in sys.path:\n",
" sys.path.insert(0, str(SRC))\n",
"\n",
"from airfrans_frontier.paths import DEFAULT_RAW_DATA_DIR, DEFAULT_RAW_MANIFEST_PATH\n",
"from airfrans_frontier.raw.inspect import format_raw_inspection, inspect_raw_subset\n",
"from airfrans_frontier.raw.manifest import load_raw_subset_manifest\n",
"\n",
"DATA_DIR = REPO_ROOT / DEFAULT_RAW_DATA_DIR\n",
"MANIFEST_PATH = REPO_ROOT / DEFAULT_RAW_MANIFEST_PATH\n",
"print(f\"repo: {REPO_ROOT}\")\n",
"print(f\"raw data: {DATA_DIR}\")\n"
]
},
{
"cell_type": "markdown",
"id": "5bb5d36c",
"metadata": {},
"source": [
"## 2. Constants and the key nondimensional numbers\n",
"\n",
"For incompressible flow, the key dimensionless control parameter is Reynolds number\n",
"\n",
"$$\n",
"\\mathrm{Re}=\\frac{U L}{\\nu}.\n",
"$$\n",
"\n",
"Here:\n",
"\n",
"- $U$ is the freestream/inlet velocity magnitude,\n",
"- $L$ is the characteristic length; AirfRANS uses the chord, $L = 1\\ \\mathrm{m}$, for the airfoil Reynolds number,\n",
"- $\\nu$ is kinematic viscosity. The dataset generation used about $1.56\\times 10^{-5}\\ \\mathrm{m^2/s}$.\n",
"\n",
"So, because $L=1$, AirfRANS effectively maps\n",
"\n",
"$$\n",
"U_\\infty \\approx \\mathrm{Re}\\,\\nu.\n",
"$$\n",
"\n",
"The paper samples $\\mathrm{Re}\\in[2,6]\\times10^6$, giving speeds around 31.2 to 93.6 m/s. The speed of sound at 298.15 K is about 346.1 m/s, so the upper speed gives $\\mathrm{Ma}\\approx0.27$, below the usual $0.3$ incompressible cutoff.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eb38682d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.520848Z",
"iopub.status.busy": "2026-07-20T19:11:12.520703Z",
"iopub.status.idle": "2026-07-20T19:11:12.523976Z",
"shell.execute_reply": "2026-07-20T19:11:12.523621Z"
}
},
"outputs": [],
"source": [
"NU_DATASET = 1.56e-5 # m^2/s, paper Table 6 / docs note\n",
"RHO = 1.184 # kg/m^3, paper Table 6\n",
"SOUND_SPEED = 346.1 # m/s, paper Table 6\n",
"CHORD = 1.0 # m\n",
"\n",
"for Re in [2e6, 3e6, 4e6, 5e6, 6e6]:\n",
" U = Re * NU_DATASET / CHORD\n",
" Ma = U / SOUND_SPEED\n",
" q_reduced = 0.5 * U**2 # reduced dynamic pressure q/rho for OpenFOAM incompressible pressure convention\n",
" q_physical = 0.5 * RHO * U**2\n",
" print(f\"Re={Re/1e6:.0f}e6 U_inf={U:6.2f} m/s Mach={Ma:.3f} q/rho={q_reduced:8.1f} m^2/s^2 q={q_physical:8.1f} Pa\")\n"
]
},
{
"cell_type": "markdown",
"id": "2631a6d7",
"metadata": {},
"source": [
"## 3. What equations are actually solved?\n",
"\n",
"The paper writes steady incompressible RANS in index notation. In plain terms:\n",
"\n",
"### Mass conservation\n",
"\n",
"$$\n",
"\\nabla\\cdot \\bar u = 0.\n",
"$$\n",
"\n",
"For incompressible flow, the mean velocity field has zero divergence: no local volume expansion/compression.\n",
"\n",
"### Momentum balance\n",
"\n",
"The paper gives, schematically,\n",
"\n",
"$$\n",
"\\partial_j(\\bar u_i\\bar u_j)\n",
"= -\\partial_i \\bar p + (\\nu + \\nu_t)\\partial_{jj}^2 \\bar u_i,\n",
"\\quad i\\in\\{1,2\\},\n",
"$$\n",
"\n",
"where pressure is the **reduced** pressure $p/\\rho$ in the incompressible convention. The indices are just compact notation for the two spatial directions: $i=1,2$ selects the velocity component being solved ($\\bar u_x$ or $\\bar u_y$), and a repeated $j$ means \"sum over both $x$ and $y$ directions.\"\n",
"\n",
"Interpretation:\n",
"\n",
"- $\\partial_j(\\bar u_i\\bar u_j)$: convection of momentum. It means the moving mean flow is carrying component-$i$ momentum through the $x$ and $y$ directions. For example, for $i=x$ it expands to $\\partial_x(\\bar u_x\\bar u_x)+\\partial_y(\\bar u_x\\bar u_y)$.\n",
"- $-\\partial_i \\bar p$: pressure-gradient force in the selected direction; pressure differences push and turn the flow.\n",
"- $(\\nu+\\nu_t)\\partial_{jj}^2\\bar u_i$: viscous/turbulent smoothing of the selected velocity component. Here $\\partial_{jj}^2$ means $\\partial_x^2+\\partial_y^2$, so this term responds to curvature/sharp spatial changes in the velocity field.\n",
"\n",
"The closure is the **$k$-$\\omega$ SST** turbulence model. That means OpenFOAM also solves extra turbulence equations for quantities such as turbulent kinetic energy $k$ and specific dissipation rate $\\omega$; those determine $\\nu_t$. The dataset exposes $\\nu_t$ as a target field because it indicates local turbulence intensity, but force computation only needs pressure and near-wall velocity gradients; on the airfoil wall $\\nu_t=0$ for this setup.\n",
"\n",
"OpenFOAM solves these equations with a finite-volume method: the domain is split into cells, integral conservation laws are enforced over each cell, and fluxes across cell faces couple neighboring cells.\n"
]
},
{
"cell_type": "markdown",
"id": "e41bd5ea",
"metadata": {},
"source": [
"## 4. What kinds of objects are involved?\n",
"\n",
"A raw AirfRANS case contains several distinct object types. Keeping them separate avoids confusion.\n",
"\n",
"| Object | Meaning in CFD | Where it appears in raw OpenFOAM |\n",
"|---|---|---|\n",
"| airfoil geometry | the solid obstacle boundary, chord 1 m | spline points inside `system/blockMeshDict`; wall patch named `aerofoil` |\n",
"| farfield/freestream boundary | outer domain boundary 200 chords away | patch named `freestream` |\n",
"| mesh vertices/points | coordinates used to define cell corners/faces | `constant/polyMesh/points.gz` |\n",
"| faces | polygon faces connecting vertices | `constant/polyMesh/faces.gz` |\n",
"| cells | finite volumes where equations are solved | OpenFOAM topology under `constant/polyMesh/`; volume fields have one row per cell |\n",
"| boundary patches | named sets of boundary faces | `constant/polyMesh/boundary` |\n",
"| volume fields | solved fields over cells | final time folder, e.g. `40000/U.gz`, `40000/p.gz`, turbulence fields |\n",
"| surface fields | wall quantities on airfoil faces | `40000/wallShearStress.gz`, `40000/forceCoeff.gz`, `40000/yPlus.gz` |\n",
"| force histories | integrated engineering outputs versus solver iteration | `postProcessing/forceCoeffs1/0/coefficient.dat` |\n",
"\n",
"The mesh is a C-grid/multiblock hexahedral mesh. It is deliberately very fine near the airfoil wall because drag depends strongly on wall shear stress, which depends on the velocity gradient immediately next to the wall.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a7370288",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.525617Z",
"iopub.status.busy": "2026-07-20T19:11:12.525521Z",
"iopub.status.idle": "2026-07-20T19:11:12.587426Z",
"shell.execute_reply": "2026-07-20T19:11:12.586915Z"
}
},
"outputs": [],
"source": [
"report = inspect_raw_subset(DATA_DIR, MANIFEST_PATH)\n",
"print(format_raw_inspection(report, sample_limit=8))\n",
"manifest = load_raw_subset_manifest(MANIFEST_PATH)\n",
"sim_names = list(manifest.simulation_names)\n"
]
},
{
"cell_type": "markdown",
"id": "8960a178",
"metadata": {},
"source": [
"## 5. Dataset design space: what is varied?\n",
"\n",
"Each simulation is defined by:\n",
"\n",
"1. turbulence model: AirfRANS dataset uses `SST` = $k$-$\\omega$ SST;\n",
"2. freestream speed $U_\\infty$;\n",
"3. angle of attack $\\alpha$ in degrees;\n",
"4. NACA airfoil parameters.\n",
"\n",
"The name stores these parameters. Example:\n",
"\n",
"```text\n",
"airFoil2D_SST_43.597_5.932_3.551_3.1_1.0_18.252\n",
" │ │ │ └──── NACA parameters; 4 numbers => 5-digit family\n",
" │ │ └────────── angle of attack α [deg]\n",
" │ └───────────────── U_inf [m/s]\n",
" └───────────────────── turbulence model\n",
"```\n",
"\n",
"NACA 4-digit family uses 3 parameters `(M, P, XX)`:\n",
"\n",
"- `M`: maximum camber in percent of chord,\n",
"- `P`: location of maximum camber in tenths of chord,\n",
"- `XX`: maximum thickness in percent of chord.\n",
"\n",
"NACA 5-digit family uses 4 parameters `(L, P, Q, XX)`:\n",
"\n",
"- `L`: controls design lift coefficient through $C_L = 0.15L$ in the classical definition,\n",
"- `P`: camber-location parameter in twentieths of chord,\n",
"- `Q`: 0 for standard camber, 1 for reflex camber,\n",
"- `XX`: maximum thickness in percent of chord.\n",
"\n",
"The AirfRANS paper samples 500 airfoils from the 4-digit family and 500 from the 5-digit family, for 1000 total simulations. The docs/code sample Reynolds uniformly in `[2e6, 6e6]` and angle of attack uniformly in `[-5°, 15°]`; high-|AoA| cases may run more iterations.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "473b211e",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.589141Z",
"iopub.status.busy": "2026-07-20T19:11:12.589030Z",
"iopub.status.idle": "2026-07-20T19:11:12.596837Z",
"shell.execute_reply": "2026-07-20T19:11:12.596532Z"
}
},
"outputs": [],
"source": [
"SIM_RE = re.compile(\n",
" r\"^airFoil2D_(?P<turbulence>[^_]+)_\"\n",
" r\"(?P<u_inf>-?\\d+(?:\\.\\d+)?)_\"\n",
" r\"(?P<alpha>-?\\d+(?:\\.\\d+)?)_\"\n",
" r\"(?P<params>.+)$\"\n",
")\n",
"\n",
"\n",
"def parse_sim_name(name: str) -> dict:\n",
" m = SIM_RE.match(name)\n",
" if not m:\n",
" raise ValueError(f\"unexpected simulation name: {name}\")\n",
" params = [float(x) for x in m.group(\"params\").split(\"_\")]\n",
" family = \"NACA-4-like\" if len(params) == 3 else \"NACA-5-like\" if len(params) == 4 else \"unknown\"\n",
" u_inf = float(m.group(\"u_inf\"))\n",
" return {\n",
" \"name\": name,\n",
" \"turbulence\": m.group(\"turbulence\"),\n",
" \"u_inf\": u_inf,\n",
" \"alpha_deg\": float(m.group(\"alpha\")),\n",
" \"naca_params\": params,\n",
" \"family\": family,\n",
" \"reynolds_from_name\": u_inf * CHORD / NU_DATASET,\n",
" \"mach_from_name\": u_inf / SOUND_SPEED,\n",
" }\n",
"\n",
"rows = [parse_sim_name(n) for n in sim_names]\n",
"u = np.array([r[\"u_inf\"] for r in rows])\n",
"a = np.array([r[\"alpha_deg\"] for r in rows])\n",
"re_vals = np.array([r[\"reynolds_from_name\"] for r in rows])\n",
"\n",
"print(f\"local simulations: {len(rows)}\")\n",
"print(f\"U_inf range: {u.min():.3f} .. {u.max():.3f} m/s\")\n",
"print(f\"Re range from U/nu: {re_vals.min():.3e} .. {re_vals.max():.3e}\")\n",
"print(f\"angle range: {a.min():.3f} .. {a.max():.3f} deg\")\n",
"print(f\"families: { {family: sum(r['family'] == family for r in rows) for family in sorted(set(r['family'] for r in rows))} }\")\n",
"rows[:3]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2f2e29bb",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.598192Z",
"iopub.status.busy": "2026-07-20T19:11:12.598103Z",
"iopub.status.idle": "2026-07-20T19:11:12.853202Z",
"shell.execute_reply": "2026-07-20T19:11:12.852448Z"
}
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(1, 3, figsize=(13, 3.5))\n",
"axes[0].hist(re_vals / 1e6, bins=12)\n",
"axes[0].set_title(\"Local subset Reynolds numbers\")\n",
"axes[0].set_xlabel(\"Re [millions]\")\n",
"axes[0].set_ylabel(\"cases\")\n",
"\n",
"axes[1].hist(a, bins=12)\n",
"axes[1].set_title(\"Angles of attack\")\n",
"axes[1].set_xlabel(\"α [deg]\")\n",
"\n",
"axes[2].scatter(re_vals / 1e6, a, s=28)\n",
"axes[2].set_title(\"Local subset coverage\")\n",
"axes[2].set_xlabel(\"Re [millions]\")\n",
"axes[2].set_ylabel(\"α [deg]\")\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "1ef32b82",
"metadata": {},
"source": [
"## 6. Pick one simulation and read its OpenFOAM dictionaries\n",
"\n",
"OpenFOAM cases are self-describing. Important files:\n",
"\n",
"- `system/controlDict`: solver application, time/iteration controls, force-coefficient postprocessor setup;\n",
"- `system/fvSchemes`: discretization schemes;\n",
"- `system/fvSolution`: linear solver and SIMPLE/SIMPLEC settings;\n",
"- `constant/transportProperties`: viscosity;\n",
"- `constant/turbulenceProperties`: RAS/turbulence model;\n",
"- `constant/polyMesh/boundary`: patch names and face ranges.\n",
"\n",
"Change `SIM_INDEX` to inspect a different case.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b1fa24ec",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.854493Z",
"iopub.status.busy": "2026-07-20T19:11:12.854338Z",
"iopub.status.idle": "2026-07-20T19:11:12.857336Z",
"shell.execute_reply": "2026-07-20T19:11:12.856721Z"
}
},
"outputs": [],
"source": [
"SIM_INDEX = 0\n",
"SIM_NAME = sim_names[SIM_INDEX]\n",
"SIM_DIR = DATA_DIR / SIM_NAME\n",
"meta = parse_sim_name(SIM_NAME)\n",
"print(SIM_NAME)\n",
"print(json.dumps(meta, indent=2))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e35d99c0",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.858618Z",
"iopub.status.busy": "2026-07-20T19:11:12.858481Z",
"iopub.status.idle": "2026-07-20T19:11:12.869906Z",
"shell.execute_reply": "2026-07-20T19:11:12.869229Z"
}
},
"outputs": [],
"source": [
"FLOAT_RE = re.compile(r\"[-+]?(?:\\d+(?:\\.\\d*)?|\\.\\d+)(?:[eE][-+]?\\d+)?\")\n",
"INT_RE = re.compile(r\"\\d+\")\n",
"\n",
"\n",
"def open_text(path: Path):\n",
" if path.suffix == \".gz\":\n",
" return gzip.open(path, \"rt\", errors=\"replace\")\n",
" return path.open(\"rt\", errors=\"replace\")\n",
"\n",
"\n",
"def read_text(path: Path, max_bytes: int = 200_000) -> str:\n",
" if path.suffix == \".gz\":\n",
" with gzip.open(path, \"rb\") as stream:\n",
" data = stream.read(max_bytes)\n",
" else:\n",
" data = path.read_bytes()[:max_bytes]\n",
" return data.decode(\"utf-8\", errors=\"replace\")\n",
"\n",
"\n",
"def assignment(text: str, key: str) -> str | None:\n",
" match = re.search(rf\"^\\s*{re.escape(key)}\\s+([^;]+);\", text, flags=re.MULTILINE)\n",
" return match.group(1).strip() if match else None\n",
"\n",
"\n",
"def vector_assignment(text: str, key: str) -> np.ndarray | None:\n",
" value = assignment(text, key)\n",
" if value is None:\n",
" return None\n",
" return np.array([float(item) for item in FLOAT_RE.findall(value)], dtype=float)\n",
"\n",
"\n",
"def parse_boundary(path: Path) -> dict[str, dict[str, int | str]]:\n",
" text = read_text(path, max_bytes=500_000)\n",
" patches: dict[str, dict[str, int | str]] = {}\n",
" for name, body in re.findall(r\"\\n\\s*([A-Za-z][A-Za-z0-9_]*)\\s*\\n\\s*\\{(.*?)\\n\\s*\\}\", text, flags=re.DOTALL):\n",
" patch_type = assignment(body, \"type\") or \"\"\n",
" n_faces = assignment(body, \"nFaces\")\n",
" start_face = assignment(body, \"startFace\")\n",
" if n_faces is not None and start_face is not None:\n",
" patches[name] = {\"type\": patch_type, \"nFaces\": int(n_faces), \"startFace\": int(start_face)}\n",
" return patches\n",
"\n",
"\n",
"def parse_foam_list(path: Path, columns: int) -> np.ndarray:\n",
" expected_count: int | None = None\n",
" in_values = False\n",
" rows = []\n",
" with open_text(path) as stream:\n",
" for line in stream:\n",
" stripped = line.strip()\n",
" if not in_values:\n",
" if expected_count is None and stripped.isdigit():\n",
" expected_count = int(stripped)\n",
" continue\n",
" if expected_count is not None and stripped == \"(\":\n",
" in_values = True\n",
" continue\n",
" continue\n",
" if stripped == \")\":\n",
" break\n",
" nums = [float(item) for item in FLOAT_RE.findall(stripped)]\n",
" if len(nums) >= columns:\n",
" rows.append(nums[:columns])\n",
" arr = np.array(rows, dtype=float)\n",
" if expected_count is not None and len(arr) != expected_count:\n",
" raise ValueError(f\"{path}: parsed {len(arr)} rows, expected {expected_count}\")\n",
" return arr.reshape(-1) if columns == 1 else arr\n",
"\n",
"\n",
"def parse_faces(path: Path, start_face: int, n_faces: int) -> list[list[int]]:\n",
" expected_count = None\n",
" in_values = False\n",
" face_index = -1\n",
" selected = []\n",
" stop_face = start_face + n_faces\n",
" with open_text(path) as stream:\n",
" for line in stream:\n",
" stripped = line.strip()\n",
" if not in_values:\n",
" if expected_count is None and stripped.isdigit():\n",
" expected_count = int(stripped)\n",
" continue\n",
" if expected_count is not None and stripped == \"(\":\n",
" in_values = True\n",
" continue\n",
" continue\n",
" if stripped == \")\":\n",
" break\n",
" face_index += 1\n",
" if face_index < start_face:\n",
" continue\n",
" if face_index >= stop_face:\n",
" break\n",
" values = [int(item) for item in INT_RE.findall(stripped)]\n",
" if values:\n",
" selected.append(values[1:])\n",
" if len(selected) != n_faces:\n",
" raise ValueError(f\"{path}: parsed {len(selected)} selected faces, expected {n_faces}\")\n",
" return selected\n",
"\n",
"\n",
"def load_force_coefficients(path: Path) -> tuple[list[str], np.ndarray]:\n",
" columns = []\n",
" with path.open(\"rt\", errors=\"replace\") as stream:\n",
" for line in stream:\n",
" if line.startswith(\"# Time\"):\n",
" columns = line[1:].split()\n",
" break\n",
" data = np.loadtxt(path, comments=\"#\")\n",
" return columns, data\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9f350525",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.871234Z",
"iopub.status.busy": "2026-07-20T19:11:12.871073Z",
"iopub.status.idle": "2026-07-20T19:11:12.876931Z",
"shell.execute_reply": "2026-07-20T19:11:12.876263Z"
}
},
"outputs": [],
"source": [
"control_text = read_text(SIM_DIR / \"system\" / \"controlDict\")\n",
"transport_text = read_text(SIM_DIR / \"constant\" / \"transportProperties\")\n",
"turbulence_text = read_text(SIM_DIR / \"constant\" / \"turbulenceProperties\")\n",
"\n",
"application = assignment(control_text, \"application\")\n",
"end_time = assignment(control_text, \"endTime\")\n",
"write_interval = assignment(control_text, \"writeInterval\")\n",
"u_inf_config = float(assignment(control_text, \"Uinf\"))\n",
"nu_case = float(assignment(transport_text, \"nu\"))\n",
"turbulence_model = assignment(turbulence_text, \"RASModel\")\n",
"drag_dir = vector_assignment(control_text, \"dragDir\")\n",
"lift_dir = vector_assignment(control_text, \"liftDir\")\n",
"alpha_from_drag = math.degrees(math.atan2(drag_dir[1], drag_dir[0])) if drag_dir is not None else float(\"nan\")\n",
"\n",
"print(f\"solver application: {application}\")\n",
"print(f\"final solver iteration/endTime: {end_time}\")\n",
"print(f\"writeInterval: {write_interval}\")\n",
"print(f\"RAS turbulence model: {turbulence_model}\")\n",
"print(f\"Uinf from controlDict: {u_inf_config:.3f} m/s\")\n",
"print(f\"nu from transportProperties: {nu_case:.3e} m^2/s\")\n",
"print(f\"Re = U*chord/nu: {u_inf_config / nu_case:.3e}\")\n",
"print(f\"Mach = U/c: {u_inf_config / SOUND_SPEED:.3f}\")\n",
"print(f\"dragDir: {drag_dir}\")\n",
"print(f\"liftDir: {lift_dir}\")\n",
"print(f\"angle inferred from dragDir: {alpha_from_drag:.3f} deg\")\n"
]
},
{
"cell_type": "markdown",
"id": "e1380985",
"metadata": {},
"source": [
"## 7. Boundary conditions in concrete terms\n",
"\n",
"For the incompressible SST cases in AirfRANS, the paper's Table 7 gives these patch meanings:\n",
"\n",
"| field | internal/initial | aerofoil wall | freestream boundary |\n",
"|---|---:|---|---|\n",
"| `U` | $U_\\infty$ | `noSlip`: wall velocity is zero | `freestreamVelocity` |\n",
"| `p` | 0 in reduced-pressure convention | `zeroGradient` | `freestreamPressure` |\n",
"| `nut` | $\\nu$ initially | low-Re wall function | `freestream` |\n",
"| `k` | about $0.001 U_\\infty^2 / Re_L$ | fixed value | `freestream` |\n",
"| `omega` | about $5U_\\infty/L$ | wall function | `freestream` |\n",
"\n",
"The important non-domain intuition:\n",
"\n",
"- **No slip wall** means the fluid sticks to the airfoil surface: velocity exactly zero at the wall. The rapid change from 0 at the wall to freestream away from it creates large velocity gradients.\n",
"- **Pressure zeroGradient at wall** means the solver does not impose a pressure value at the wall; pressure is solved consistently with momentum and continuity.\n",
"- **Freestream boundary** lets the outer C-grid behave like a far-away free flow boundary.\n",
"- The outer boundary is placed 200 chords away so the arbitrary boundary condition does not pollute the near-airfoil field too much.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1e0fc222",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.878190Z",
"iopub.status.busy": "2026-07-20T19:11:12.878049Z",
"iopub.status.idle": "2026-07-20T19:11:12.884898Z",
"shell.execute_reply": "2026-07-20T19:11:12.884089Z"
}
},
"outputs": [],
"source": [
"for field in [\"U\", \"p\", \"nut\", \"k\", \"omega\"]:\n",
" candidate = SIM_DIR / \"0\" / field\n",
" if not candidate.exists():\n",
" candidate = SIM_DIR / \"0\" / f\"turbulenceProperties:{field}\"\n",
" if candidate.exists():\n",
" text = read_text(candidate, max_bytes=2_500)\n",
" print(f\"\\n--- 0/{candidate.name} ---\")\n",
" # Show header plus boundaryField region start; enough to connect concepts to files.\n",
" lines = text.splitlines()\n",
" for line in lines[:80]:\n",
" print(line)\n",
" else:\n",
" print(f\"missing initial field for {field}\")\n"
]
},
{
"cell_type": "markdown",
"id": "f9513904",
"metadata": {},
"source": [
"## 8. Mesh: why so many cells for a 2D problem?\n",
"\n",
"The paper reports around 250k-300k cells per full simulation. That is large for a 2D teaching example, but small compared with industrial 3D CFD.\n",
"\n",
"The key reason is the **boundary layer** near the airfoil. Drag needs wall shear stress:\n",
"\n",
"$$\n",
"\\tau_w \\propto \\mu \\frac{\\partial u_\\text{tangent}}{\\partial n}\\bigg|_\\text{wall}.\n",
"$$\n",
"\n",
"That derivative is only trustworthy if the first cell center is extremely close to the wall. The AirfRANS mesh uses first-cell height $2\\ \\mu\\mathrm{m}$ near the airfoil to get $y^+\\approx 1$, i.e. low-Re wall resolution rather than a coarse wall-function-only treatment.\n",
"\n",
"`points.gz` are vertices. Volume field arrays such as `U.gz` are cell-centered. Boundary patch arrays are face-centered. Counts will not match one-to-one.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dc969024",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:12.886185Z",
"iopub.status.busy": "2026-07-20T19:11:12.886085Z",
"iopub.status.idle": "2026-07-20T19:11:14.072342Z",
"shell.execute_reply": "2026-07-20T19:11:14.071808Z"
}
},
"outputs": [],
"source": [
"points = parse_foam_list(SIM_DIR / \"constant\" / \"polyMesh\" / \"points.gz\", columns=3)\n",
"patches = parse_boundary(SIM_DIR / \"constant\" / \"polyMesh\" / \"boundary\")\n",
"print(f\"mesh vertices/points: {points.shape}\")\n",
"print(\"patches:\")\n",
"for name, info in patches.items():\n",
" print(f\" {name:12} type={info['type']:10} nFaces={info['nFaces']:8} startFace={info['startFace']}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "04c89f1f",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:14.074124Z",
"iopub.status.busy": "2026-07-20T19:11:14.074009Z",
"iopub.status.idle": "2026-07-20T19:11:14.347016Z",
"shell.execute_reply": "2026-07-20T19:11:14.346335Z"
}
},
"outputs": [],
"source": [
"rng = np.random.default_rng(20260720)\n",
"sample_count = min(80_000, len(points))\n",
"point_sample = points[rng.choice(len(points), size=sample_count, replace=False)]\n",
"\n",
"fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
"axes[0].scatter(point_sample[:, 0], point_sample[:, 1], s=0.15, alpha=0.25)\n",
"axes[0].set_title(\"Full C-grid vertex sample\")\n",
"axes[0].set_aspect(\"equal\", adjustable=\"box\")\n",
"axes[0].set_xlabel(\"x [m]\")\n",
"axes[0].set_ylabel(\"y [m]\")\n",
"\n",
"near = point_sample[(point_sample[:, 0] > -0.4) & (point_sample[:, 0] < 1.4) & (np.abs(point_sample[:, 1]) < 0.45)]\n",
"axes[1].scatter(near[:, 0], near[:, 1], s=0.4, alpha=0.45)\n",
"axes[1].set_title(\"Near-airfoil mesh concentration\")\n",
"axes[1].set_aspect(\"equal\", adjustable=\"box\")\n",
"axes[1].set_xlabel(\"x [m]\")\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "8adf2d5d",
"metadata": {},
"source": [
"## 9. Integrated forces: what $C_D$ and $C_L$ mean\n",
"\n",
"Surface stress force combines pressure and viscous stress. The paper writes a face force as\n",
"\n",
"$$\n",
"df = -p n\\,dS + 2\\mu S\\cdot n\\,dS.\n",
"$$\n",
"\n",
"OpenFOAM integrates over the airfoil wall. Then:\n",
"\n",
"- **drag** $D$ is the component parallel to freestream direction,\n",
"- **lift** $L$ is the component perpendicular to freestream direction,\n",
"- coefficients divide by dynamic pressure times reference area:\n",
"\n",
"$$\n",
"C_D = \\frac{D}{\\frac{1}{2}\\rho U_\\infty^2 A},\\qquad\n",
"C_L = \\frac{L}{\\frac{1}{2}\\rho U_\\infty^2 A}.\n",
"$$\n",
"\n",
"For this 2D dataset, the reference area is effectively chord $\\times$ unit span, so $A=1\\ \\mathrm{m^2}$.\n",
"\n",
"Why this matters for ML: a model can have visually plausible pressure/velocity fields but still rank airfoils incorrectly by $C_D$ or $C_L$. AirfRANS evaluates both fields and force coefficients.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ba58a7d4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:14.349038Z",
"iopub.status.busy": "2026-07-20T19:11:14.348886Z",
"iopub.status.idle": "2026-07-20T19:11:14.385211Z",
"shell.execute_reply": "2026-07-20T19:11:14.383778Z"
}
},
"outputs": [],
"source": [
"coeff_path = SIM_DIR / \"postProcessing\" / \"forceCoeffs1\" / \"0\" / \"coefficient.dat\"\n",
"coeff_columns, coeff_data = load_force_coefficients(coeff_path)\n",
"lookup = {name: idx for idx, name in enumerate(coeff_columns)}\n",
"time = coeff_data[:, lookup[\"Time\"]]\n",
"cd = coeff_data[:, lookup[\"Cd\"]]\n",
"cl = coeff_data[:, lookup[\"Cl\"]]\n",
"cm = coeff_data[:, lookup[\"CmPitch\"]]\n",
"\n",
"print(coeff_columns)\n",
"print(f\"rows: {len(coeff_data)}\")\n",
"print(f\"final iteration: {time[-1]:.0f}\")\n",
"print(f\"final Cd: {cd[-1]:.6f}\")\n",
"print(f\"final Cl: {cl[-1]:.6f}\")\n",
"print(f\"final CmPitch: {cm[-1]:.6f}\")\n",
"print(f\"last 500 Cd std: {np.std(cd[-min(500, len(cd)):]):.3e}\")\n",
"print(f\"last 500 Cl std: {np.std(cl[-min(500, len(cl)):]):.3e}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51cb325f",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:14.386355Z",
"iopub.status.busy": "2026-07-20T19:11:14.386250Z",
"iopub.status.idle": "2026-07-20T19:11:14.682440Z",
"shell.execute_reply": "2026-07-20T19:11:14.681608Z"
}
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
"axes[0].plot(time, cd, label=\"Cd\")\n",
"axes[0].plot(time, cl, label=\"Cl\")\n",
"axes[0].plot(time, cm, label=\"CmPitch\")\n",
"axes[0].set_title(\"Force/moment coefficient convergence history\")\n",
"axes[0].set_xlabel(\"simpleFoam iteration\")\n",
"axes[0].legend()\n",
"\n",
"window = min(750, len(time))\n",
"axes[1].plot(time[-window:], cd[-window:], label=\"Cd\")\n",
"axes[1].plot(time[-window:], cl[-window:], label=\"Cl\")\n",
"axes[1].set_title(f\"Final {window} iterations\")\n",
"axes[1].set_xlabel(\"simpleFoam iteration\")\n",
"axes[1].legend()\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "bc431803",
"metadata": {},
"source": [
"## 10. Airfoil-surface quantities: yPlus, wall shear, force coefficient density\n",
"\n",
"The `aerofoil` patch is a sequence of wall faces. For each face we can approximate a face center from its vertices, then plot surface fields:\n",
"\n",
"- `wallShearStress`: vector shear stress-like wall quantity from the OpenFOAM postprocessing field;\n",
"- `forceCoeff`: local force coefficient contribution on faces;\n",
"- `yPlus`: nondimensional first-cell wall distance. Values near or below 1 are what the paper wanted for resolved boundary layers.\n",
"\n",
"These surface quantities are central to drag. Lift is often dominated by pressure distribution; drag is more sensitive to velocity-gradient errors near the wall.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6d3ea796",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:14.683807Z",
"iopub.status.busy": "2026-07-20T19:11:14.683663Z",
"iopub.status.idle": "2026-07-20T19:11:14.777574Z",
"shell.execute_reply": "2026-07-20T19:11:14.776955Z"
}
},
"outputs": [],
"source": [
"aerofoil_patch = patches[\"aerofoil\"]\n",
"aerofoil_faces = parse_faces(\n",
" SIM_DIR / \"constant\" / \"polyMesh\" / \"faces.gz\",\n",
" start_face=int(aerofoil_patch[\"startFace\"]),\n",
" n_faces=int(aerofoil_patch[\"nFaces\"]),\n",
")\n",
"aerofoil_centers = np.array([points[face, :2].mean(axis=0) for face in aerofoil_faces])\n",
"\n",
"surface_force = parse_foam_list(SIM_DIR / \"40000\" / \"forceCoeff.gz\", columns=3)\n",
"wall_shear = parse_foam_list(SIM_DIR / \"40000\" / \"wallShearStress.gz\", columns=3)\n",
"y_plus = parse_foam_list(SIM_DIR / \"40000\" / \"yPlus.gz\", columns=1)\n",
"\n",
"force_norm = np.linalg.norm(surface_force[:, :2], axis=1)\n",
"shear_norm = np.linalg.norm(wall_shear[:, :2], axis=1)\n",
"\n",
"print(f\"aerofoil faces: {len(aerofoil_faces)}\")\n",
"print(f\"surface_force rows: {surface_force.shape}\")\n",
"print(f\"wall_shear rows: {wall_shear.shape}\")\n",
"print(f\"yPlus rows: {y_plus.shape}\")\n",
"print(f\"yPlus min/median/max: {y_plus.min():.4f} / {np.median(y_plus):.4f} / {y_plus.max():.4f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "29982f97",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:14.778789Z",
"iopub.status.busy": "2026-07-20T19:11:14.778638Z",
"iopub.status.idle": "2026-07-20T19:11:15.101946Z",
"shell.execute_reply": "2026-07-20T19:11:15.101306Z"
}
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(1, 3, figsize=(14, 4))\n",
"for ax, values, title in [\n",
" (axes[0], force_norm, \"|forceCoeff| on wall faces\"),\n",
" (axes[1], shear_norm, \"|wallShearStress|\"),\n",
" (axes[2], y_plus, \"yPlus\"),\n",
"]:\n",
" sc = ax.scatter(aerofoil_centers[:, 0], aerofoil_centers[:, 1], c=values, s=13, cmap=\"viridis\")\n",
" ax.set_title(title)\n",
" ax.set_aspect(\"equal\", adjustable=\"box\")\n",
" ax.set_xlabel(\"x [m]\")\n",
" fig.colorbar(sc, ax=ax, shrink=0.8)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a9582f0a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:15.103453Z",
"iopub.status.busy": "2026-07-20T19:11:15.103294Z",
"iopub.status.idle": "2026-07-20T19:11:15.282512Z",
"shell.execute_reply": "2026-07-20T19:11:15.282063Z"
}
},
"outputs": [],
"source": [
"idx = np.arange(len(y_plus))\n",
"fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)\n",
"axes[0].plot(idx, force_norm)\n",
"axes[0].set_ylabel(\"|forceCoeff|\")\n",
"axes[1].plot(idx, shear_norm)\n",
"axes[1].set_ylabel(\"|wall shear|\")\n",
"axes[2].plot(idx, y_plus)\n",
"axes[2].axhline(1.0, color=\"black\", linestyle=\"--\", linewidth=1, label=\"yPlus=1\")\n",
"axes[2].set_ylabel(\"yPlus\")\n",
"axes[2].set_xlabel(\"aerofoil boundary face index; not physical arc-length\")\n",
"axes[2].legend()\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "39e3601e",
"metadata": {},
"source": [
"## 11. Volume fields: what the surrogate model sees as targets\n",
"\n",
"At the final iteration folder (`40000/` here), the raw OpenFOAM case contains cell-centered fields. AirfRANS's preprocessed ML dataset keeps, after cropping/slicing:\n",
"\n",
"- velocity $\\bar u_x, \\bar u_y$,\n",
"- reduced pressure $\\bar p/\\rho$,\n",
"- turbulent kinematic viscosity $\\nu_t$,\n",
"- signed/implicit distance to the airfoil as input geometry signal,\n",
"- surface flag / normals as geometry features.\n",
"\n",
"The paper's supervised learning target per point is roughly\n",
"\n",
"$$\n",
"y_i = (\\bar u_x,\\bar u_y,\\bar p,\\nu_t)_i.\n",
"$$\n",
"\n",
"The input per point is roughly\n",
"\n",
"$$\n",
"x_i = (x,y,U_{\\infty,x},U_{\\infty,y},\\mathrm{distance},n_x,n_y)_i,\n",
"$$\n",
"\n",
"with normals zero away from the airfoil. Exact library attributes include `position`, `input_velocity`, `sdf`, `surface`, and `normals`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "187f7fea",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:15.284075Z",
"iopub.status.busy": "2026-07-20T19:11:15.283994Z",
"iopub.status.idle": "2026-07-20T19:11:17.039996Z",
"shell.execute_reply": "2026-07-20T19:11:17.039369Z"
}
},
"outputs": [],
"source": [
"U = parse_foam_list(SIM_DIR / \"40000\" / \"U.gz\", columns=3)\n",
"p = parse_foam_list(SIM_DIR / \"40000\" / \"p.gz\", columns=1)\n",
"nut_path = SIM_DIR / \"40000\" / \"turbulenceProperties:nut.gz\"\n",
"k_path = SIM_DIR / \"40000\" / \"turbulenceProperties:k.gz\"\n",
"omega_path = SIM_DIR / \"40000\" / \"turbulenceProperties:omega.gz\"\n",
"nut = parse_foam_list(nut_path, columns=1) if nut_path.exists() else np.array([])\n",
"k = parse_foam_list(k_path, columns=1) if k_path.exists() else np.array([])\n",
"omega = parse_foam_list(omega_path, columns=1) if omega_path.exists() else np.array([])\n",
"speed = np.linalg.norm(U[:, :2], axis=1)\n",
"\n",
"def summary(name, values):\n",
" finite = values[np.isfinite(values)]\n",
" return {\n",
" \"name\": name,\n",
" \"count\": int(values.size),\n",
" \"min\": float(np.min(finite)),\n",
" \"p01\": float(np.percentile(finite, 1)),\n",
" \"mean\": float(np.mean(finite)),\n",
" \"p99\": float(np.percentile(finite, 99)),\n",
" \"max\": float(np.max(finite)),\n",
" }\n",
"\n",
"summaries = [summary(\"|U_xy|\", speed), summary(\"Ux\", U[:,0]), summary(\"Uy\", U[:,1]), summary(\"p\", p)]\n",
"if nut.size: summaries.append(summary(\"nut\", nut))\n",
"if k.size: summaries.append(summary(\"k\", k))\n",
"if omega.size: summaries.append(summary(\"omega\", omega))\n",
"summaries\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2eccca27",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:17.041292Z",
"iopub.status.busy": "2026-07-20T19:11:17.041183Z",
"iopub.status.idle": "2026-07-20T19:11:17.671734Z",
"shell.execute_reply": "2026-07-20T19:11:17.671231Z"
}
},
"outputs": [],
"source": [
"plots = [(speed, \"speed |U_xy|\"), (p, \"reduced pressure p\"), (nut, \"turbulent viscosity nut\"), (k, \"turbulent kinetic energy k\")]\n",
"fig, axes = plt.subplots(2, 2, figsize=(12, 8))\n",
"for ax, (values, title) in zip(axes.flat, plots):\n",
" if values.size:\n",
" ax.hist(values[np.isfinite(values)], bins=80)\n",
" ax.set_yscale(\"log\")\n",
" ax.set_title(title)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "99baa471",
"metadata": {},
"source": [
"## 12. Airfoil geometry: generate a simple NACA 4-digit profile\n",
"\n",
"AirfRANS includes both NACA 4- and 5-digit families. This code implements only the 4-digit equations because they are compact and useful for intuition.\n",
"\n",
"For NACA `MPXX`:\n",
"\n",
"- $m=0.01M$ is maximum camber fraction,\n",
"- $p=0.1P$ is max-camber position as chord fraction,\n",
"- $t=XX/100$ is thickness fraction.\n",
"\n",
"The thickness envelope is\n",
"\n",
"$$\n",
"y_t(x)=\\frac{t}{0.2}\\left(0.2969\\sqrt{x}-0.1260x-0.3516x^2+0.2843x^3-0.1015x^4\\right).\n",
"$$\n",
"\n",
"The camber line bends the upper/lower surfaces by the local angle $\\theta=\\arctan(y_c')$.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fdb418b",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:17.673116Z",
"iopub.status.busy": "2026-07-20T19:11:17.672936Z",
"iopub.status.idle": "2026-07-20T19:11:17.777014Z",
"shell.execute_reply": "2026-07-20T19:11:17.776578Z"
}
},
"outputs": [],
"source": [
"def naca4(M=2.0, P=4.0, XX=12.0, n=250):\n",
" x = np.linspace(0, 1, n)\n",
" m = 0.01 * M\n",
" p = 0.1 * P\n",
" t = XX / 100\n",
" yt = (t / 0.2) * (0.2969*np.sqrt(x) - 0.1260*x - 0.3516*x**2 + 0.2843*x**3 - 0.1015*x**4)\n",
" if m == 0 or p == 0:\n",
" yc = np.zeros_like(x)\n",
" dyc = np.zeros_like(x)\n",
" else:\n",
" yc = np.where(x <= p, m / p**2 * (2*p*x - x**2), m / (1-p)**2 * ((1 - 2*p) + 2*p*x - x**2))\n",
" dyc = np.where(x <= p, 2*m / p**2 * (p - x), 2*m / (1-p)**2 * (p - x))\n",
" theta = np.arctan(dyc)\n",
" xu = x - yt * np.sin(theta)\n",
" yu = yc + yt * np.cos(theta)\n",
" xl = x + yt * np.sin(theta)\n",
" yl = yc - yt * np.cos(theta)\n",
" return x, yc, xu, yu, xl, yl\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 3))\n",
"for M, P, XX, label in [(0, 0, 12, \"0012 symmetric\"), (2, 4, 12, \"2412 cambered\"), (4, 4, 15, \"4415 thicker/cambered\")]:\n",
" x, yc, xu, yu, xl, yl = naca4(M, P, XX)\n",
" ax.plot(xu, yu, label=label)\n",
" ax.plot(xl, yl, color=ax.lines[-1].get_color())\n",
"ax.set_aspect(\"equal\", adjustable=\"box\")\n",
"ax.set_xlabel(\"x / chord\")\n",
"ax.set_ylabel(\"y / chord\")\n",
"ax.set_title(\"NACA 4-digit geometry intuition\")\n",
"ax.legend()\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "d7668533",
"metadata": {},
"source": [
"## 13. What the data is used for\n",
"\n",
"The dataset is built for **surrogate modeling**: replace a slow CFD solve with a learned function that maps geometry + boundary conditions to flow fields and force coefficients.\n",
"\n",
"Concrete ML tasks in the paper/docs:\n",
"\n",
"1. **Full data interpolation**: train on 800 simulations, test on 200 from same distribution.\n",
"2. **Scarce data interpolation**: train on 200, same 200-case test set as full regime.\n",
"3. **Reynolds extrapolation**: train only on $\\mathrm{Re}\\in[3,5]\\times10^6$, test below 3e6 and above 5e6.\n",
"4. **Angle-of-attack extrapolation**: train on $\\alpha\\in[-2.5^\\circ,12.5^\\circ]$, test on $[-5,-2.5]$ and $[12.5,15]$ degrees.\n",
"\n",
"Targets and evaluation:\n",
"\n",
"- field MSE for $\\bar u_x,\\bar u_y,\\bar p,\\nu_t$ in the volume;\n",
"- surface pressure MSE on airfoil;\n",
"- relative error for $C_D$ and $C_L$;\n",
"- Spearman rank correlation for $C_D$ and $C_L$, important for optimization because ranking airfoils correctly can matter more than exact coefficient calibration.\n",
"\n",
"The paper found drag harder than lift for the baselines. Reason: lift is often dominated by pressure distribution, which models predicted better; drag depends strongly on wall shear stress, i.e. near-wall velocity gradients, which are hard to recover from sampled point-cloud training.\n"
]
},
{
"cell_type": "markdown",
"id": "17c606a1",
"metadata": {},
"source": [
"## 14. Preprocessing: why raw OpenFOAM differs from the ML dataset\n",
"\n",
"Raw OpenFOAM cases contain the full farfield and many files. The preprocessed AirfRANS dataset:\n",
"\n",
"- clips the internal field to $x\\in[-2,4]$, $y\\in[-1.5,1.5]$, $z\\in[0,1]$;\n",
"- slices at $z=0.5$ to get a 2D representation from the thin 3D OpenFOAM mesh;\n",
"- computes an implicit distance / signed-distance-like field to the airfoil;\n",
"- keeps only pressure, velocity, turbulent viscosity, distance for internal nodes;\n",
"- for the airfoil patch, computes inward-pointing normals, slices to 1D, and keeps pressure, velocity, turbulent viscosity, normals.\n",
"\n",
"The raw case remains useful because it exposes solver dictionaries, gradients/postprocessing fields, and force histories. The preprocessed case is convenient for ML.\n"
]
},
{
"cell_type": "markdown",
"id": "2111d2fa",
"metadata": {},
"source": [
"## 15. Practical checklist when you inspect or model these simulations\n",
"\n",
"Use this checklist to avoid domain mistakes:\n",
"\n",
"1. **Do not treat the mesh as an image.** It is unstructured and highly nonuniform; near-wall density is intentional.\n",
"2. **Distinguish cells, vertices, and faces.** Raw OpenFOAM volume fields here are cell-centered; patch fields are face-centered; preprocessed VTK data may be node-based.\n",
"3. **Pressure is reduced pressure in incompressible OpenFOAM.** It is effectively $p/\\rho$, not absolute atmospheric pressure.\n",
"4. **Reynolds is encoded through velocity.** With fixed $\\nu$ and chord 1 m, $U_\\infty$ determines $\\mathrm{Re}$.\n",
"5. **Angle of attack rotates the freestream/force directions.** Check `dragDir` and `liftDir`, not only the filename.\n",
"6. **Drag needs near-wall gradients.** Good volume MSE does not guarantee good $C_D$.\n",
"7. **Lift/drag coefficients are integrated outputs.** They are computed from pressure and wall shear over the airfoil surface, normalized by dynamic pressure.\n",
"8. **Extrapolation tasks are materially harder.** A model good inside the sampled Re/AoA distribution may fail outside it.\n",
"\n",
"Next cells print raw dictionary snippets so you can connect the explanations above to actual files.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "336f680c",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:17.778411Z",
"iopub.status.busy": "2026-07-20T19:11:17.778320Z",
"iopub.status.idle": "2026-07-20T19:11:17.780856Z",
"shell.execute_reply": "2026-07-20T19:11:17.780315Z"
}
},
"outputs": [],
"source": [
"print(\"--- system/controlDict preview ---\")\n",
"print(read_text(SIM_DIR / \"system\" / \"controlDict\", max_bytes=6_000))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "940c0923",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:17.782129Z",
"iopub.status.busy": "2026-07-20T19:11:17.782013Z",
"iopub.status.idle": "2026-07-20T19:11:17.784645Z",
"shell.execute_reply": "2026-07-20T19:11:17.784234Z"
}
},
"outputs": [],
"source": [
"print(\"--- constant/turbulenceProperties preview ---\")\n",
"print(read_text(SIM_DIR / \"constant\" / \"turbulenceProperties\", max_bytes=3_000))\n",
"print(\"\\n--- constant/transportProperties preview ---\")\n",
"print(read_text(SIM_DIR / \"constant\" / \"transportProperties\", max_bytes=3_000))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3992406a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-20T19:11:17.785699Z",
"iopub.status.busy": "2026-07-20T19:11:17.785619Z",
"iopub.status.idle": "2026-07-20T19:11:17.787920Z",
"shell.execute_reply": "2026-07-20T19:11:17.787486Z"
}
},
"outputs": [],
"source": [
"print(\"--- final U field header and first values ---\")\n",
"print(read_text(SIM_DIR / \"40000\" / \"U.gz\", max_bytes=4_000))\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in a new issue Copy permalink