NACA 0012 Airfoil: Symmetric Airfoil Validation¶
This case evaluates Flow360 on the NACA 0012 airfoil from the NASA Turbulence Modeling Resource at Mach 0.15 and Reynolds number 6.0e6. It is a fundamental steady, subsonic, external-aerodynamics validation of the flow solver and its turbulence models.
The study is a sweep of 3 turbulence models (SA, SA-RC, SST) × 3 angles of attack (0°, 10°, 15°) = 9 cases. For each case we compute the lift and drag coefficients and the surface-pressure distribution, and compare them against experimental data and results from other CFD codes (CFL3D, FUN3D, OVERFLOW) run on the same workshop mesh.
The notebook runs top-to-bottom in a single kernel: load the mesh, submit the sweep, wait for completion, then post-process the in-kernel case objects into the three published figures.
Requirements¶
This notebook runs against the Flow360 Python API. Install it and configure your API key by following the installation & setup guide.
This notebook targets the latest Flow360 Python client. It fetches the case's
root asset and reference data with flow360.examples.download_benchmark_assets,
added in 25.10.4. Install the client (25.10.4 or newer) before running:
pip install "flow360>=25.10.4"
Beyond flow360, the notebook needs a few plotting/data packages:
pip install numpy pandas matplotlib pyvista
Run the cells top to bottom in a single kernel.
Imports¶
import json
import os
import re
import tarfile
from pathlib import Path
from tempfile import TemporaryDirectory
import flow360 as fl
import numpy as np
import pandas as pd
import pyvista as pv
import matplotlib.pyplot as plt
from flow360.examples import download_benchmark_assets
Input data¶
The post-processing compares Flow360 against experimental measurements and other
CFD codes. Fetch those reference files from the public benchmark bucket; this
recreates a local ./ref_data/ directory that the plotting cells read from.
download_benchmark_assets("NACA_0012", "ref_data")
Load project¶
The root asset for this case is a pre-generated volume mesh. Download it from the public benchmark bucket and start a fresh project from it, no private project id is needed.
root_asset_files = download_benchmark_assets("NACA_0012", "root_assets")
# The snapshot bundles the mesh with solver metadata and logs; select the
# volume mesh file (.cgns or .ugrid, possibly .zst-compressed) to load.
volume_mesh_file = next(
f for f in root_asset_files if f.removesuffix(".zst").endswith((".cgns", ".ugrid"))
)
project = fl.Project.from_volume_mesh(volume_mesh_file)
Physics setup¶
Every case in the sweep shares the same operating condition, reference geometry,
boundary conditions, and time stepping, only the turbulence model and angle of
attack change. We capture the shared setup in build_params, and the
model-specific pieces in three small turbulence builders.
Each builder returns the turbulence-model solver and its matching inlet turbulence quantities:
- SA: Spalart-Allmaras.
- SA-RC: Spalart-Allmaras with rotation/curvature correction.
- SST: Menter k-ω SST.
def make_sa():
return (
fl.SpalartAllmaras(
absolute_tolerance=1e-10,
equation_evaluation_frequency=1,
),
fl.TurbulenceQuantities(viscosity_ratio=0.210438),
)
def make_sa_rc():
return (
fl.SpalartAllmaras(
absolute_tolerance=1e-10,
rotation_correction=True,
equation_evaluation_frequency=1,
),
fl.TurbulenceQuantities(viscosity_ratio=0.210438),
)
def make_sst():
return (
fl.KOmegaSST(
absolute_tolerance=1e-10,
equation_evaluation_frequency=1,
),
fl.TurbulenceQuantities(
viscosity_ratio=0.009,
turbulent_intensity=0.052 * 0.01,
),
)
TURBULENCE_BUILDERS = {
"sa": make_sa,
"sa-rc": make_sa_rc,
"sst": make_sst,
}
build_params assembles the full SimulationParams for one case. The mesh
boundaries are addressed by name pattern: the airfoil surface (*5) is a viscous
wall, the symmetry planes (*1, *2) are slip walls, and the outer farfield
(*3, *4, *6) is a freestream carrying the turbulence quantities. The flow
solver uses a tight Navier-Stokes tolerance with the low-Mach preconditioner (the
freestream is nearly incompressible at M = 0.15), and a steady ramp-CFL schedule.
Surface Cp and Cf are written for the post-processing.
def build_params(project, alpha_deg, turbulence_model, turbulence_quantities):
with fl.SI_unit_system:
volume_mesh = project.volume_mesh
return fl.SimulationParams(
operating_condition=fl.AerospaceCondition.from_mach_reynolds(
mach=0.15,
reynolds_mesh_unit=6e6,
project_length_unit=1 * fl.u.m,
temperature=300 * fl.u.K,
alpha=alpha_deg * fl.u.deg,
),
reference_geometry=fl.ReferenceGeometry(
moment_center=(0, 0, 0),
moment_length=(1, 1, 1) * fl.u.m,
area=1 * fl.u.m**2,
),
models=[
fl.Wall(surfaces=volume_mesh["*5"]),
fl.SlipWall(entities=[volume_mesh["*1"], volume_mesh["*2"]]),
fl.Freestream(
surfaces=[volume_mesh["*3"], volume_mesh["*4"], volume_mesh["*6"]],
turbulence_quantities=turbulence_quantities,
),
fl.Fluid(
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-12,
low_mach_preconditioner=True,
),
turbulence_model_solver=turbulence_model,
),
],
time_stepping=fl.Steady(
max_steps=20000,
CFL=fl.RampCFL(
initial=1,
final=100,
ramp_steps=10000,
),
),
outputs=[
fl.SurfaceOutput(
entities=volume_mesh["*5"],
output_format="paraview",
output_fields=["Cp", "Cf"],
)
],
)
Submit cases¶
Sweep the three turbulence models over the three angles of attack and submit all
nine cases to Flow360. Each case is named
naca0012_<model>_a<alpha>_<solver_version>. Alongside each submitted case we
keep its turbulence label and angle so the post-processing can use the in-kernel
objects directly.
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")
ALPHAS_DEG = [0, 10, 15]
cases = []
for turbulence_name, builder in TURBULENCE_BUILDERS.items():
for alpha_deg in ALPHAS_DEG:
turbulence_model, turbulence_quantities = builder()
case = project.run_case(
params=build_params(
project=project,
alpha_deg=alpha_deg,
turbulence_model=turbulence_model,
turbulence_quantities=turbulence_quantities,
),
name=f"naca0012_{turbulence_name}_a{alpha_deg}_{SOLVER_VERSION}",
solver_version=SOLVER_VERSION,
)
cases.append({"turbulence": turbulence_name, "alpha_deg": alpha_deg, "case": case})
Wait for completion¶
Block until every submitted case has finished. This can take a long time (minutes to hours), since each case is a full steady solve of up to 20000 steps.
for record in cases:
record["case"].wait()
Postprocessing¶
The published results are three figures:
- Lift coefficient: CL vs α, one panel per turbulence model.
- Drag polar: CD vs CL, one panel per turbulence model.
- Surface pressure: Cp vs x/c, a grid over turbulence model × angle of attack.
Each is compared against experimental data and other CFD codes on the same mesh. The cells below read the results straight from the in-kernel case objects (no project is re-opened) and rebuild exactly these three figures.
Setup¶
Styling colors, the turbulence-model labels used across all plots, the reference
data file locations, and two small numeric helpers: mean_last_percent averages
the last 1% of the convergence history to get a settled coefficient, and
read_xy_csv reads a two-column reference curve.
os.makedirs("results", exist_ok=True)
FLOW360_COLOR = "#00643c" # Flow360 results
REFERENCE_COLOR = "black" # experimental / CFD-code references
TURBULENCE_LABELS = {"sa": "SA", "sa-rc": "SA-RC", "sst": "SST"}
TURBULENCE_MODELS = ("SA", "SA-RC", "SST")
CP_ALPHAS_DEG = (0.0, 10.0, 15.0)
# --- surface Cp extraction ---
CP_FIELD = "Cp"
CP_SLICE_ORIGIN = (0.0, -0.5, 0.0)
CP_SLICE_NORMAL = (0.0, 1.0, 0.0)
CP_SPLIT_VARIABLE = "Z" # coordinate that separates upper/lower surface
CP_SPLIT_THRESHOLD = 0.0
COORDINATE_INDEX = {"X": 0, "Y": 1, "Z": 2}
VTK_SUFFIXES = (".pvtu", ".pvtp", ".vtp", ".vtu", ".vtk")
# --- reference data ---
CFL3D_ZONE_RE = re.compile(r'^zone,\s*t="alpha=(?P<alpha>-?\d+(?:\.\d+)?)"$', re.IGNORECASE)
CP_CFL3D_REFERENCE_FILES = {
"SA": "ref_data/n0012cp_cfl3d_sa.dat",
"SA-RC": "ref_data/n0012cp_cfl3d_sarc.dat",
"SST": "ref_data/n0012cp_cfl3d_sst.dat",
}
CP_REFERENCE_FILES = {
0: "ref_data/exp_data_CP_alpha/exp_CP_alpha0.csv",
10: "ref_data/exp_data_CP_alpha/exp_CP_alpha10.csv",
15: "ref_data/exp_data_CP_alpha/exp_CP_alpha15.csv",
}
# --- plot series styling ---
CP_SERIES_STYLES = (
("Flow360 Upper Surface", FLOW360_COLOR, "-", 2),
("Flow360 Lower Surface", FLOW360_COLOR, "--", 2),
("CFL3D Upper Surface", REFERENCE_COLOR, ":", 1.75),
("CFL3D Lower Surface", REFERENCE_COLOR, "-.", 1.75),
)
COEFFICIENT_CODE_SERIES = ("CFL3D", "FUN3D", "OVERFLOW")
COEFFICIENT_SERIES_STYLES = (
("Flow360", FLOW360_COLOR, "d", 45, 1.2),
("CFL3D", REFERENCE_COLOR, "P", 55, 1.1),
("FUN3D", "#CC79A7", "s", 45, 1.1),
("OVERFLOW", "#D55E00", "^", 45, 1.1),
("Experiment", REFERENCE_COLOR, "x", 25, 1.0),
)
def mean_last_percent(series):
series = series.dropna()
if series.empty:
return None
start = max(int(len(series) * 0.99), 0)
return float(series.iloc[start:].mean())
def read_xy_csv(path, x_name, y_name):
df = pd.read_csv(path, header=None, usecols=[0, 1], names=[x_name, y_name], skipinitialspace=True)
df[x_name] = pd.to_numeric(df[x_name], errors="coerce")
df[y_name] = pd.to_numeric(df[y_name], errors="coerce")
return df.dropna(subset=[x_name, y_name])
Load the swept-case results¶
For each finished case, read the total-forces convergence history and take the
settled CL/CD (mean of the last 1% of steps). We collect these into case_df and
also key the case objects by (model, α) for the surface-pressure plot. Finally we
load the other-code reference coefficients from code_results.csv.
case_rows = []
cases_by_key = {}
for record in cases:
case = record["case"]
model_label = TURBULENCE_LABELS[record["turbulence"]]
alpha_deg = float(record["alpha_deg"])
cases_by_key[(model_label, alpha_deg)] = case
force_data = case.results.total_forces.as_dataframe()
if "physical_step" in force_data.columns:
force_data = force_data.sort_values("physical_step")
elif "pseudo_step" in force_data.columns:
force_data = force_data.sort_values("pseudo_step")
case_rows.append(
{
"turbulence_model": model_label,
"alpha_deg": alpha_deg,
"CL": mean_last_percent(force_data["CL"]),
"CD": mean_last_percent(force_data["CD"]),
}
)
case_df = pd.DataFrame(case_rows).sort_values(["turbulence_model", "alpha_deg"]).reset_index(drop=True)
code_df = pd.read_csv("ref_data/code_results.csv")
code_df["alpha_deg"] = pd.to_numeric(code_df["alpha_deg"])
code_df["CL"] = pd.to_numeric(code_df["CL"])
code_df["CD"] = pd.to_numeric(code_df["CD"])
Lift coefficient and drag polar¶
Assemble a tidy table of (turbulence_model, series, α, CL, CD) from the Flow360 results, the other CFD codes, and experiment. The experimental drag is interpolated onto the experimental lift curve so the polar is directly comparable. The same data drives both side-by-side figures, CL vs α and the CD vs CL drag polar, one panel per turbulence model.
def coefficient_plot_data(case_df, code_df):
flow_df = case_df[["turbulence_model", "alpha_deg", "CL", "CD"]].copy()
flow_df["series"] = "Flow360"
code_plot_df = (
code_df[code_df["Code"].isin(COEFFICIENT_CODE_SERIES)]
.rename(columns={"Turb": "turbulence_model", "Code": "series"})[
["turbulence_model", "series", "alpha_deg", "CL", "CD"]
]
.copy()
)
exp_cl_df = read_xy_csv("ref_data/exp_data_CD_CL/exp_CL_ALPHA2.csv", "alpha_deg", "CL")
exp_cd_df = (
read_xy_csv("ref_data/exp_data_CD_CL/exp_CD_CL.csv", "CL", "CD")
.sort_values("CL")
.drop_duplicates("CL")
)
exp_cl_df["CD"] = np.interp(exp_cl_df["CL"], exp_cd_df["CL"], exp_cd_df["CD"])
exp_cl_df["series"] = "Experiment"
exp_df = pd.concat(
[exp_cl_df.assign(turbulence_model=model_label) for model_label in TURBULENCE_MODELS],
ignore_index=True,
)
return pd.concat([flow_df, code_plot_df, exp_df], ignore_index=True)[
["turbulence_model", "series", "alpha_deg", "CL", "CD"]
]
def plot_coefficient_points(ax, panel_df, x_name, y_name):
handles = {}
for series_name, color, marker, size, line_width in COEFFICIENT_SERIES_STYLES:
series_df = panel_df[panel_df["series"] == series_name]
if series_df.empty:
continue
alpha = 0.7 if series_name in COEFFICIENT_CODE_SERIES else 1.0
if marker in ("x", "+", "*"):
scatter = ax.scatter(
series_df[x_name], series_df[y_name],
color=color, marker=marker, s=size, linewidths=line_width,
alpha=alpha, label=series_name,
)
else:
scatter = ax.scatter(
series_df[x_name], series_df[y_name],
facecolors=color, edgecolors=color, marker=marker, s=size,
linewidths=line_width, alpha=alpha, label=series_name,
)
handles[series_name] = scatter
ax.grid(True, alpha=0.3)
return handles
def write_side_by_side_plot(plot_df, path, x_name, y_name, x_label, y_label):
fig, axes = plt.subplots(1, len(TURBULENCE_MODELS), figsize=(18, 3.6), squeeze=False)
legend_handles = {}
for col_idx, model_label in enumerate(TURBULENCE_MODELS):
ax = axes[0][col_idx]
panel_df = plot_df[plot_df["turbulence_model"] == model_label]
ax.set_title(model_label)
ax.set_xlabel(x_label)
if col_idx == 0:
ax.set_ylabel(y_label)
legend_handles.update(plot_coefficient_points(ax, panel_df, x_name, y_name))
handles = [legend_handles[name] for name, *_ in COEFFICIENT_SERIES_STYLES if name in legend_handles]
fig.legend(
handles, [h.get_label() for h in handles],
loc="lower center", ncol=min(len(handles), 6), frameon=False,
)
fig.tight_layout(rect=(0, 0.12, 1, 1))
fig.savefig(path, dpi=300)
return fig
coefficients = coefficient_plot_data(case_df, code_df)
coefficients = coefficients.sort_values(["turbulence_model", "series", "alpha_deg"]).reset_index(drop=True)
write_side_by_side_plot(
coefficients, "results/CL_alpha_side_by_side.png",
"alpha_deg", "CL", r"$\alpha$ (deg)", r"$C_L$",
)
write_side_by_side_plot(
coefficients, "results/CD_vs_CL_side_by_side.png",
"CL", "CD", r"$C_L$", r"$C_D$",
)
plt.show()
Surface pressure¶
For each case, download the surface output, slice it along the airfoil mid-plane, and split it into upper and lower surfaces by the spanwise coordinate. Flow360's Cp is plotted against the CFL3D reference distributions and the experimental measurements, in a grid of turbulence model (rows) × angle of attack (columns). The y-axis is inverted, as is conventional for Cp.
def find_surface_files(root_dir):
root_dir = Path(root_dir)
for suffix in VTK_SUFFIXES:
matches = sorted(root_dir.rglob(f"*{suffix}"))
if matches:
return matches
return []
def load_surface_mesh(vtk_files):
meshes = []
for path in vtk_files:
mesh = pv.read(str(path))
if CP_FIELD in mesh.cell_data and CP_FIELD not in mesh.point_data:
mesh = mesh.cell_data_to_point_data(pass_cell_data=True)
meshes.append(mesh)
return meshes[0] if len(meshes) == 1 else pv.merge(meshes, merge_points=False)
def field_values(mesh, name):
coord_index = COORDINATE_INDEX.get(name.upper())
if coord_index is not None:
return mesh.points[:, coord_index]
if name in mesh.point_data:
values = mesh.point_data[name]
elif name in mesh.cell_data:
values = mesh.cell_data_to_point_data(pass_cell_data=True).point_data[name]
else:
raise KeyError(name)
if values.ndim == 2 and values.shape[1] == 1:
return values[:, 0]
if values.ndim != 1:
raise ValueError(f"Variable '{name}' is not scalar.")
return values
def extract_cp_slice(case):
with TemporaryDirectory() as tmp_dir:
archive_path = Path(tmp_dir) / "surfaces.tar.gz"
extract_dir = Path(tmp_dir) / "surfaces"
extract_dir.mkdir()
case.results.surfaces.to_file(str(archive_path), overwrite=True)
with tarfile.open(archive_path, "r:gz") as tar:
tar.extractall(extract_dir, filter="tar")
vtk_files = find_surface_files(extract_dir)
if not vtk_files:
raise RuntimeError(f"No surface files found for {case.name}")
surface_mesh = load_surface_mesh(vtk_files)
slice_mesh = surface_mesh.slice(normal=CP_SLICE_NORMAL, origin=CP_SLICE_ORIGIN)
if slice_mesh.n_points == 0:
raise RuntimeError(f"Slice plane did not intersect surface output for {case.name}")
slice_df = pd.DataFrame(
{
"x_over_c": field_values(slice_mesh, "X"),
"surface_coord": field_values(slice_mesh, CP_SPLIT_VARIABLE),
"Cp": field_values(slice_mesh, CP_FIELD),
}
)
return slice_df.dropna().sort_values("x_over_c").reset_index(drop=True)
def read_exp_cp_reference(alpha_deg):
path = CP_REFERENCE_FILES.get(int(round(alpha_deg)))
if path is None:
return pd.DataFrame(columns=["alpha_deg", "series", "x_over_c", "Cp"])
exp_df = pd.read_csv(path, header=None, usecols=[0, 1], names=["x_over_c", "Cp"], skipinitialspace=True)
exp_df["x_over_c"] = pd.to_numeric(exp_df["x_over_c"], errors="coerce")
exp_df["Cp"] = pd.to_numeric(exp_df["Cp"], errors="coerce")
exp_df = exp_df.dropna(subset=["x_over_c", "Cp"])
exp_df["series"] = "Experiment"
exp_df["alpha_deg"] = float(alpha_deg)
return exp_df[["alpha_deg", "series", "x_over_c", "Cp"]]
def read_cfl3d_cp_reference(alpha_deg, path):
alpha_int = int(round(alpha_deg))
current_alpha = None
rows = []
with open(path, encoding="utf-8") as f:
for line in f:
stripped = line.strip()
if not stripped or stripped.startswith("#") or stripped.startswith("variables="):
continue
zone_match = CFL3D_ZONE_RE.match(stripped)
if zone_match:
current_alpha = int(round(float(zone_match.group("alpha"))))
continue
if current_alpha != alpha_int:
continue
parts = stripped.split()
if len(parts) < 2:
continue
try:
x_val = float(parts[0])
cp_val = float(parts[1])
except ValueError:
continue
rows.append({"x_over_c": x_val, "Cp": cp_val})
if not rows:
return pd.DataFrame(columns=["alpha_deg", "series", "x_over_c", "Cp"])
cfl3d_df = pd.DataFrame(rows)
leading_edge_idx = int(cfl3d_df["x_over_c"].idxmin())
lower_df = cfl3d_df.iloc[: leading_edge_idx + 1].sort_values("x_over_c").copy()
upper_df = cfl3d_df.iloc[leading_edge_idx:].sort_values("x_over_c").copy()
lower_df["series"] = "CFL3D Lower Surface"
upper_df["series"] = "CFL3D Upper Surface"
lower_df["alpha_deg"] = float(alpha_deg)
upper_df["alpha_deg"] = float(alpha_deg)
return pd.concat(
[
upper_df[["alpha_deg", "series", "x_over_c", "Cp"]],
lower_df[["alpha_deg", "series", "x_over_c", "Cp"]],
],
ignore_index=True,
)
def cp_plot_data(model_label, alpha_deg, case, cfl3d_cp_reference_file):
slice_df = extract_cp_slice(case)
upper_df = slice_df[slice_df["surface_coord"] >= CP_SPLIT_THRESHOLD][["x_over_c", "Cp"]].copy()
lower_df = slice_df[slice_df["surface_coord"] <= CP_SPLIT_THRESHOLD][["x_over_c", "Cp"]].copy()
upper_df["series"] = "Flow360 Upper Surface"
lower_df["series"] = "Flow360 Lower Surface"
upper_df["alpha_deg"] = float(alpha_deg)
lower_df["alpha_deg"] = float(alpha_deg)
frames = [
upper_df[["alpha_deg", "series", "x_over_c", "Cp"]],
lower_df[["alpha_deg", "series", "x_over_c", "Cp"]],
read_exp_cp_reference(alpha_deg),
]
if cfl3d_cp_reference_file is not None:
frames.append(read_cfl3d_cp_reference(alpha_deg, cfl3d_cp_reference_file))
plot_df = pd.concat(frames, ignore_index=True)
plot_df["turbulence_model"] = model_label
return plot_df[["turbulence_model", "alpha_deg", "series", "x_over_c", "Cp"]]
def plot_cp_panel(ax, plot_df):
handles = {}
for series_name, color, line_style, line_width in CP_SERIES_STYLES:
series_df = plot_df[plot_df["series"] == series_name]
if not series_df.empty:
(line,) = ax.plot(
series_df["x_over_c"], series_df["Cp"],
color=color, linestyle=line_style, linewidth=line_width, label=series_name,
)
handles[series_name] = line
exp_df = plot_df[plot_df["series"] == "Experiment"]
if not exp_df.empty:
scatter = ax.scatter(
exp_df["x_over_c"], exp_df["Cp"],
facecolors="none", edgecolors=REFERENCE_COLOR, label="Experiment",
)
handles["Experiment"] = scatter
ax.invert_yaxis()
ax.grid(True, alpha=0.3)
return handles
plot_frames = []
for model_label in TURBULENCE_MODELS:
for alpha_deg in CP_ALPHAS_DEG:
case = cases_by_key.get((model_label, alpha_deg))
if case is None:
continue
plot_frames.append(
cp_plot_data(model_label, alpha_deg, case, CP_CFL3D_REFERENCE_FILES.get(model_label))
)
cp_df = pd.concat(plot_frames, ignore_index=True)
cp_df = cp_df.sort_values(["turbulence_model", "alpha_deg", "series", "x_over_c"]).reset_index(drop=True)
fig, axes = plt.subplots(len(TURBULENCE_MODELS), len(CP_ALPHAS_DEG), figsize=(21, 13.5), squeeze=False)
legend_handles = {}
for row_idx, model_label in enumerate(TURBULENCE_MODELS):
for col_idx, alpha_deg in enumerate(CP_ALPHAS_DEG):
ax = axes[row_idx][col_idx]
panel_df = cp_df[(cp_df["turbulence_model"] == model_label) & (cp_df["alpha_deg"] == alpha_deg)]
legend_handles.update(plot_cp_panel(ax, panel_df))
if row_idx == 0:
ax.set_title(rf"$\alpha = {alpha_deg:g}^\circ$")
if col_idx == 0:
ax.set_ylabel(f"{model_label}\n" + r"$C_p$")
if row_idx == len(TURBULENCE_MODELS) - 1:
ax.set_xlabel("x/c")
handles = [legend_handles[name] for name, *_ in CP_SERIES_STYLES if name in legend_handles]
if "Experiment" in legend_handles:
handles.append(legend_handles["Experiment"])
fig.legend(handles, [h.get_label() for h in handles], loc="lower center", ncol=len(handles), frameon=False)
fig.tight_layout(rect=(0, 0.04, 1, 1))
fig.savefig("results/Cp_vs_x.png", dpi=300)
plt.show()