NACA 4412: Trailing-Edge Separation¶
This case evaluates Flow360 on the NACA 4412 airfoil from the NASA Turbulence Modeling Resource, a classic trailing-edge separation benchmark. The flow is steady and subsonic at Mach 0.09, Reynolds number 1.52e6, and an angle of attack of 13.87°, conditions under which the boundary layer separates near the trailing edge on the suction side of the airfoil.
The study runs the same case with three turbulence models (SA, SA-RC, and k-ω SSTm), and compares Flow360 against CFL3D, FUN3D, and experimental data from the workshop. For each model we look at the surface-pressure distribution, the boundary-layer velocity and turbulent-shear profiles at six chordwise stations near the trailing edge, and the integrated aerodynamic loads.
The notebook runs top-to-bottom in a single kernel: load the mesh, submit the three cases, wait for completion, then post-process the in-kernel case objects into the published table and 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¶
In [ ]:
import os
import re
import shutil
import tarfile
from pathlib import Path
import flow360 as fl
import numpy as np
import pandas as pd
import pyvista as pv
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from flow360.examples import download_benchmark_assets
Input data¶
The post-processing compares Flow360 against CFL3D, FUN3D, and experimental
data. Fetch those reference files from the public benchmark bucket; this
recreates a local ./ref_data/ directory that the plotting cells read from.
In [ ]:
download_benchmark_assets("NACA_4412", "ref_data")
Load project¶
The root asset for this case is a pre-generated volume mesh of the NACA 4412 airfoil. Download it from the public benchmark bucket and start a fresh project from it, no private project id is needed.
In [ ]:
root_asset_files = download_benchmark_assets("NACA_4412", "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, name="NACA 4412")
Physics setup¶
Every case shares the same boundary conditions and flow solver; only the
turbulence model changes across the sweep. The airfoil surface (*Wall) is a
viscous wall, the two spanwise planes are slip walls, and the outer boundary
(*Farfield) is a freestream. The Navier-Stokes solver runs at a tight
tolerance with the low-Mach preconditioner disabled.
In [ ]:
def build_models(turbulence_model, volume_mesh):
with fl.SI_unit_system:
return [
fl.Fluid(
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-10,
linear_solver=fl.LinearSolver(max_iterations=35),
low_mach_preconditioner=False,
),
turbulence_model_solver=turbulence_model,
),
fl.SlipWall(surfaces=[volume_mesh["Block/Symm1"], volume_mesh["Block/Symm2"]]),
fl.Wall(surfaces=volume_mesh["*Wall"]),
fl.Freestream(surfaces=volume_mesh["*Farfield"]),
]
Outputs setup¶
The case writes two solver outputs. Six chordwise slices near the trailing edge
(defined by build_slices) carry the fields needed for the boundary-layer
velocity and turbulent-shear profiles, and a surface output on the airfoil wall
carries the pressure coefficient for the surface-pressure comparison.
In [ ]:
def build_slices():
return [
fl.Slice(
name="x=.6753",
normal=(-0.9867857316, 0, 0.1620306139),
origin=(0.6986141272, -0.5, 0.2125983647) * fl.u.m,
),
fl.Slice(
name="x=.7308",
normal=(-0.9849131727, 0, 0.1730492481),
origin=(0.7557469103, -0.5, 0.2033860993) * fl.u.m,
),
fl.Slice(
name="x=.7863",
normal=(-0.9821066699, 0, 0.1883254867),
origin=(0.8135267289, -0.5, 0.1930202262) * fl.u.m,
),
fl.Slice(
name="x=.8418",
normal=(-0.9794377168, 0, 0.2017467692),
origin=(0.8710465558, -0.5, 0.1814956895) * fl.u.m,
),
fl.Slice(
name="x=.8973",
normal=(-0.9747913323, 0, 0.2231184855),
origin=(0.9297989129, -0.5, 0.1687878009) * fl.u.m,
),
fl.Slice(
name="x=.9528",
normal=(-0.9706600354, 0, 0.2404560163),
origin=(0.9879697115, -0.5, 0.1548373658) * fl.u.m,
),
]
def build_outputs(volume_mesh):
with fl.SI_unit_system:
return [
fl.SliceOutput(
slices=build_slices(),
output_fields=["mut", "primitiveVars", "gradW"],
output_format="paraview",
),
fl.SurfaceOutput(
entities=volume_mesh["*Wall"],
output_fields=["Cp", "Cf", "CfVec"],
output_format="paraview",
),
]
Simulation Params¶
create_params assembles the full SimulationParams for one turbulence model:
the operating condition (Mach 0.09, Re = 1.52e6, α = 13.87°), the reference
geometry, the models above, a steady time-stepping schedule, and the outputs.
In [ ]:
def create_params(project, turbulence_model):
with fl.SI_unit_system:
volume_mesh = project.volume_mesh
return fl.SimulationParams(
operating_condition=fl.AerospaceCondition.from_mach_reynolds(
mach=0.09,
reynolds_mesh_unit=1.52e6,
project_length_unit=1 * fl.u.m,
temperature=298 * fl.u.K,
alpha=13.87 * 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=build_models(turbulence_model, volume_mesh),
time_stepping=fl.Steady(max_steps=20000),
outputs=build_outputs(volume_mesh),
)
Submit cases¶
Submit the three cases, SA, SA-RC, and k-ω SSTm, to Flow360. The SSTm case
is given a longer step budget (35000 vs 20000) to converge the separated
boundary layer. Each submitted case is kept in submitted_cases alongside its
model key so the post-processing can use the in-kernel objects directly.
In [ ]:
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")
turbulence_setups = {
"sa": fl.SpalartAllmaras(
absolute_tolerance=1e-10,
linear_solver=fl.LinearSolver(max_iterations=25),
),
"sa_rc": fl.SpalartAllmaras(
absolute_tolerance=1e-10,
rotation_correction=True,
linear_solver=fl.LinearSolver(max_iterations=25),
),
"sstm": fl.KOmegaSST(
absolute_tolerance=1e-10,
linear_solver=fl.LinearSolver(max_iterations=25),
update_jacobian_frequency=1,
equation_evaluation_frequency=1,
),
}
submitted_cases = []
for model_key, turbulence_model in turbulence_setups.items():
params = create_params(project, turbulence_model)
if model_key == "sstm":
params.time_stepping.max_steps = 35000
case = project.run_case(
params=params,
name=f"NACA_4412_{model_key}_{SOLVER_VERSION}",
solver_version=SOLVER_VERSION,
tags=[SOLVER_VERSION],
)
submitted_cases.append({"model_key": model_key, "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 35000 steps.
In [ ]:
for record in submitted_cases:
record["case"].wait()
Postprocessing¶
The published results are one table and two figures:
- Integrated loads: CL, CD, CDp, CDv for each turbulence model, tabulated next to the CFL3D, FUN3D, and OVERFLOW reference values.
- Surface pressure: Cp vs x/c, one panel per turbulence model, compared against CFL3D, FUN3D, and experiment.
- Velocity and turbulent-shear profiles: u, v, and u'v' vs wall-normal distance at six chordwise stations near the trailing edge, in a grid of turbulence model × velocity component.
The cells below read the results straight from the in-kernel case objects (no project is re-opened) and rebuild exactly these artifacts.
Setup¶
Styling colors (Flow360 brand green for the solver results, black/blue for the reference codes), the turbulence-model labels used across every artifact, the reference-data directory, and the slice-profile constants (normalization velocity, spanwise binning, and the axis limits used for each velocity component).
In [ ]:
os.makedirs("results", exist_ok=True)
FLOW360_COLOR = "#00643c" # Flow360 results
REFERENCE_COLOR = "black" # experimental / CFD-code references
REF_DIR = Path("ref_data")
VTK_SUFFIXES = (".pvtu", ".pvtp", ".vtp", ".vtu", ".vtk")
MODEL_LABELS = {"sa": "SA", "sa_rc": "SA-RC", "sstm": "SSTm"}
MODEL_ORDER = ["sa", "sa_rc", "sstm"]
FORCE_COLUMNS = {"CL": "CL", "CD": "CD", "CDp": "CDPressure", "CDv": "CDSkinFriction"}
# --- integrated-loads table ---
LOAD_TABLE_COLUMNS = ["Code", "Turb. model", "CL", "CD", "CDp", "CDv"]
LOAD_TABLE_MODEL_LABELS = {"sa": "SA", "sa_rc": "SA-RC", "sstm": "k-omega SSTm"}
LOAD_TABLE_DECIMALS = {"CL": 4, "CD": 5, "CDp": 5, "CDv": 6}
# --- surface-pressure reference styling ---
SURFACE_REFERENCE_COLORS = {"CFL3D": REFERENCE_COLOR, "FUN3D": "#4c72b0"}
SURFACE_REFERENCE_LINESTYLES = {"CFL3D": "-", "FUN3D": (0, (4, 3))}
# --- boundary-layer slice profiles ---
CFL3D_SLICE_SUFFIXES = {"sa": "sa", "sa_rc": "sarc", "sstm": "sst"}
SLICE_PLOT_COMPONENTS = ("vaxial", "vnormal", "turbulent_shear")
SLICE_COMPONENT_LABELS = {
"vaxial": r"u/U$_{ref}$",
"vnormal": r"v/U$_{ref}$",
"turbulent_shear": r"u'v'/U$_{ref}^2$",
}
SLICE_COMPONENT_XLIMS = {
"vaxial": (-0.2, 1.3),
"vnormal": (-0.2, 0.2),
"turbulent_shear": (-0.02, 0.02),
}
SLICE_COMPONENT_XTICKS = {
"vaxial": np.arange(-0.2, 1.2001, 0.2),
"vnormal": np.arange(-0.2, 0.2001, 0.05),
"turbulent_shear": np.arange(-0.02, 0.0201, 0.01),
}
SLICE_AX_BOX_ASPECT = 6.0 / 5.0
SLICE_Z_BIN_SIZE = 0.001
SLICE_Z_MAX = 0.12
SLICE_UREF = 0.09 * 0.93
Result-loading helpers¶
Small utilities for reading Flow360 result archives (_extract_vtk_files,
_prepare_mesh), for the settled integrated loads (_mean_last_percent averages
the last 1% of the convergence history), for extracting the mid-span surface Cp
curve (_extract_surface_curve), and for assembling and rendering the surface-
pressure figure and the integrated-loads table.
In [ ]:
def _prepare_mesh(mesh, fields):
if mesh is None:
return None
if isinstance(mesh, pv.MultiBlock):
blocks = [_prepare_mesh(block, fields) for block in mesh if block is not None]
blocks = [b for b in blocks if b is not None and getattr(b, "n_points", 0) > 0]
if not blocks:
return None
if len(blocks) == 1:
return blocks[0]
return pv.merge(blocks, merge_points=False)
if getattr(mesh, "n_points", 0) == 0:
return None
if any(name in mesh.cell_data and name not in mesh.point_data for name in fields):
mesh = mesh.cell_data_to_point_data(pass_cell_data=True)
elif not mesh.point_data and mesh.cell_data:
mesh = mesh.cell_data_to_point_data(pass_cell_data=True)
return mesh
def _extract_vtk_files(results_obj, stem):
archive_path = Path("results") / f"{stem}.tar.gz"
extract_dir = Path("results") / stem
results_obj.to_file(str(archive_path), overwrite=True)
extract_dir.mkdir(parents=True, exist_ok=True)
with tarfile.open(archive_path, mode="r:*") as archive:
archive.extractall(path=extract_dir)
vtk_files = []
for suffix in VTK_SUFFIXES:
vtk_files.extend(sorted(extract_dir.rglob(f"*{suffix}")))
return archive_path, extract_dir, vtk_files
def _cleanup_paths(*paths):
for path in paths:
path = Path(path)
if path.is_dir():
shutil.rmtree(path)
elif path.exists():
path.unlink()
def _mean_last_percent(series, fraction=0.01):
if series is None:
return None
values = pd.to_numeric(series, errors="coerce").dropna()
if values.empty:
return None
start = max(int(len(values) * (1.0 - fraction)), 0)
return float(values.iloc[start:].mean())
def _extract_surface_curve(case):
archive_path, extract_dir, vtk_files = _extract_vtk_files(
case.results.surfaces, f"{case.id}_surfaces"
)
try:
if not vtk_files:
raise FileNotFoundError(f"No VTK surfaces found for case {case.name!r}.")
meshes = []
for vtk_path in vtk_files:
mesh = _prepare_mesh(pv.read(str(vtk_path)), ("Cp",))
if mesh is not None:
meshes.append(mesh)
if not meshes:
raise RuntimeError(f"No readable surface mesh found for case {case.name!r}.")
merged = meshes[0] if len(meshes) == 1 else pv.merge(meshes, merge_points=False)
section = _prepare_mesh(merged.slice(normal=(0, 1, 0), origin=(0, -0.5, 0)), ("Cp",))
if section is None or section.n_points == 0:
raise RuntimeError(f"Mid-span surface slice is empty for case {case.name!r}.")
curve = pd.DataFrame(
{
"x": section.points[:, 0],
"z": section.points[:, 2],
"Cp": pd.to_numeric(section.point_data["Cp"], errors="coerce"),
}
).dropna(subset=["x", "z"])
upper = curve[curve["z"] >= 0].sort_values("x").copy()
upper["surface"] = "upper"
lower = curve[curve["z"] < 0].sort_values("x", ascending=False).copy()
lower["surface"] = "lower"
return pd.concat([upper, lower], ignore_index=True)
finally:
_cleanup_paths(archive_path, extract_dir)
def _ordered_flow_surface(model_df):
upper = model_df[model_df["surface"] == "upper"].sort_values("x")
lower = model_df[model_df["surface"] == "lower"].sort_values("x", ascending=False)
return pd.concat([upper, lower], ignore_index=True)
def _load_surface_reference(filename):
reference_df = pd.read_csv(REF_DIR / filename)
reference_df["model"] = reference_df["model_key"].map(MODEL_LABELS)
return reference_df
def _load_cp_reference():
return pd.read_csv(
REF_DIR / "exp_cp.csv", header=None, names=["x", "value"], skipinitialspace=True
)
def _plot_surface_quantity(surface_df, quantity, out_base, reference_df, experiment_df=None):
fig, axes = plt.subplots(1, len(MODEL_ORDER), figsize=(14.5, 4.4), sharex=True)
for idx, (ax, model_key) in enumerate(zip(axes, MODEL_ORDER)):
model_flow = _ordered_flow_surface(surface_df[surface_df["model_key"] == model_key])
model_ref = reference_df[reference_df["model_key"] == model_key]
if not model_flow.empty:
ax.plot(
model_flow["x"], model_flow[quantity],
color=FLOW360_COLOR, linewidth=2.0,
label="Flow360" if idx == 0 else None,
)
for source in ("CFL3D", "FUN3D"):
source_df = model_ref[model_ref["source"] == source].sort_values("point_index")
if source_df.empty:
continue
ax.plot(
source_df["x"], source_df["value"],
color=SURFACE_REFERENCE_COLORS[source],
linestyle=SURFACE_REFERENCE_LINESTYLES[source],
linewidth=1.3, label=source if idx == 0 else None,
)
if experiment_df is not None and not experiment_df.empty:
ax.scatter(
experiment_df["x"], experiment_df["value"],
color=REFERENCE_COLOR, marker="s", s=22, facecolors="white",
linewidths=0.8, label="Experiment" if idx == 0 else None, zorder=3,
)
ax.set_title(MODEL_LABELS[model_key])
ax.set_xlim(0.0, 1.0)
ax.grid(True, alpha=0.3)
ax.set_xlabel("x/c")
if quantity == "Cp":
ax.set_ylim(2.0, -8.0)
axes[0].set_ylabel(quantity)
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="upper center", ncol=len(labels), bbox_to_anchor=(0.5, 1.03), frameon=True)
fig.tight_layout(rect=(0, 0, 1, 0.93))
fig.savefig(f"results/{out_base}.png", dpi=220, bbox_inches="tight")
return fig
def _normalize_reference_model(name):
value = str(name).strip().upper()
if value == "SA":
return "sa"
if value == "SA-RC":
return "sa_rc"
if value in {"SST", "SSTM"}:
return "sstm"
return None
def _load_reference_load_table():
reference_df = pd.read_csv(REF_DIR / "reference_data.csv", dtype=str).fillna("")
reference_df["model_key"] = reference_df["Turb. model"].map(_normalize_reference_model)
return reference_df.dropna(subset=["model_key"]).copy()
def _format_reference_load_value(value):
value = str(value).strip()
return value if value else "\N{EN DASH}"
def _format_flow_load_value(quantity, value):
if pd.isna(value):
return "\N{EN DASH}"
return f"{float(value):.{LOAD_TABLE_DECIMALS[quantity]}f}"
def _make_integrated_loads_table(flow_df):
reference_df = _load_reference_load_table()
rows = []
for model_key in MODEL_ORDER:
model_label = LOAD_TABLE_MODEL_LABELS[model_key]
model_reference = reference_df[reference_df["model_key"] == model_key]
for _, ref_row in model_reference.iterrows():
row = {
column: _format_reference_load_value(ref_row.get(column, ""))
for column in LOAD_TABLE_COLUMNS
}
row["Turb. model"] = model_label
rows.append(row)
flow_row = flow_df[flow_df["model_key"] == model_key]
if not flow_row.empty:
flow_row = flow_row.iloc[0]
rows.append(
{
"Code": "Flow360",
"Turb. model": model_label,
**{
quantity: _format_flow_load_value(quantity, flow_row.get(quantity))
for quantity in ["CL", "CD", "CDp", "CDv"]
},
}
)
return pd.DataFrame(rows, columns=LOAD_TABLE_COLUMNS)
Boundary-layer slice-profile helpers¶
The velocity/turbulent-shear figure compares Flow360 profiles at each slice
against CFL3D and experimental data read from Tecplot-style zoned .dat files.
These helpers parse those zones, extract and spanwise-average the Flow360 slice
fields, match each Flow360 slice to its reference zone by chordwise station, and
draw the model × component grid.
In [ ]:
def _slug(text):
return re.sub(r"[^a-zA-Z0-9._-]+", "_", str(text)).strip("_")
def _slice_sort_value(text):
match = re.search(r"x\s*[=:_]?\s*(-?\d*\.?\d+(?:[eE][-+]?\d+)?)", str(text))
return float(match.group(1)) if match else float("inf")
def _slice_label(text):
value = _slice_sort_value(text)
if np.isfinite(value):
return f"x={value:.4f}".replace("0.", ".").replace("-0.", "-.")
return str(text)
def _load_zoned_data(dat_path, columns):
zones = {}
current_name = None
rows = []
with open(dat_path, encoding="utf-8") as handle:
for raw in handle:
line = raw.strip()
if not line or line.startswith(("#", "VARIABLES")):
continue
if line.startswith("ZONE"):
if current_name and rows:
zones[current_name] = pd.DataFrame(rows, columns=columns)
match = re.search(r'T\s*=\s*"([^"]+)"', line)
current_name = match.group(1) if match else line
rows = []
continue
parts = line.split()
if len(parts) >= len(columns):
rows.append([float(value) for value in parts[: len(columns)]])
if current_name and rows:
zones[current_name] = pd.DataFrame(rows, columns=columns)
return zones
def _load_experimental_zones():
path = REF_DIR / "experimental_data.dat"
if not path.is_file():
raise FileNotFoundError(f"Experimental slice data not found: {path}")
return _load_zoned_data(path, ["x", "z", "u", "v", "uv"])
def _load_cfl3d_slice_zones(model_key):
suffix = CFL3D_SLICE_SUFFIXES[model_key]
vel_zones = _load_zoned_data(REF_DIR / f"n4412_cfl3d_vel_{suffix}.dat", ["x", "z", "u", "v"])
uv_zones = _load_zoned_data(REF_DIR / f"n4412_cfl3d_uv_{suffix}.dat", ["x", "z", "uv"])
zones = {}
for name, vel_df in vel_zones.items():
uv_df = uv_zones.get(name)
if uv_df is None:
continue
df = vel_df.copy()
df["uv"] = uv_df["uv"].to_numpy(dtype=float)
zones[name] = df
if not zones:
raise RuntimeError(f"No paired CFL3D velocity/uv slice zones found for {MODEL_LABELS[model_key]}.")
return zones
def _slice_dataframe(mesh):
velocity = np.asarray(mesh.point_data["velocity"], dtype=float)
grad_u = np.asarray(mesh.point_data["grad_u"], dtype=float)
mut = np.asarray(mesh.point_data["mut"], dtype=float)
return pd.DataFrame(
{
"z": np.asarray(mesh.points[:, 2], dtype=float),
"vnormal": velocity[:, 2],
"vaxial": velocity[:, 0],
"turbulent_shear": -(mut * grad_u[:, 2]),
}
)
def _spanwise_average_by_z(df, bin_size=SLICE_Z_BIN_SIZE):
if df.empty:
return df.copy()
if bin_size <= 0.0:
raise ValueError(f"bin_size must be positive. Got {bin_size}.")
return (
df.assign(_z_key=(df["z"] / float(bin_size)).round().astype(int))
.groupby("_z_key", as_index=False)
.agg({"z": "mean", "vnormal": "mean", "vaxial": "mean", "turbulent_shear": "mean"})
.sort_values("z")
.reset_index(drop=True)
)
def _match_experimental_zone(slice_name, zones):
if slice_name in zones:
return zones[slice_name]
slug = _slug(slice_name)
if slug in zones:
return zones[slug]
x_match = re.search(r"x\s*[=:_]?\s*(-?\d*\.?\d+(?:[eE][-+]?\d+)?)", slice_name)
if x_match:
x_value = float(x_match.group(1))
best_key = None
best_delta = None
for key in zones:
key_match = re.search(r"x\s*[=:_]?\s*(-?\d*\.?\d+(?:[eE][-+]?\d+)?)", key)
if not key_match:
continue
delta = abs(float(key_match.group(1)) - x_value)
if best_delta is None or delta < best_delta:
best_delta = delta
best_key = key
if best_key is not None and best_delta is not None and best_delta < 1e-4:
return zones[best_key]
for key, zone in zones.items():
if key in slice_name or _slug(key) in slug:
return zone
return None
def _project_experimental(zone_df, z0):
return {
"z_rel": zone_df["z"].to_numpy(dtype=float) - z0,
"vaxial": zone_df["u"].to_numpy(dtype=float),
"vnormal": zone_df["v"].to_numpy(dtype=float),
"turbulent_shear": zone_df["uv"].to_numpy(dtype=float),
}
def _flow_slice_data(df, z0):
z_rel = (df["z"] - z0).to_numpy(dtype=float)
return {
"z_rel": z_rel,
"vaxial": (df["vaxial"] / SLICE_UREF).to_numpy(dtype=float),
"vnormal": (df["vnormal"] / SLICE_UREF).to_numpy(dtype=float),
"turbulent_shear": (df["turbulent_shear"] / (SLICE_UREF * SLICE_UREF)).to_numpy(dtype=float),
}
def _add_slice_rows(rows, model_key, slice_name, source, data):
for component in SLICE_PLOT_COMPONENTS:
z_values = np.asarray(data["z_rel"], dtype=float)
values = np.asarray(data[component], dtype=float)
mask = (z_values >= 0.0) & (z_values <= SLICE_Z_MAX + 1e-12)
for z_rel, value in zip(z_values[mask], values[mask]):
rows.append(
{
"model_key": model_key,
"model": MODEL_LABELS[model_key],
"slice": slice_name,
"source": source,
"component": component,
"z_rel": z_rel,
"value": value,
}
)
def _collect_slice_profile_rows(case, model_key, experimental_zones, cfl3d_zones):
archive_path, extract_dir, vtk_files = _extract_vtk_files(case.results.slices, f"{case.id}_slices")
try:
if not vtk_files:
return []
rows = []
for vtk_path in vtk_files:
mesh = _prepare_mesh(pv.read(str(vtk_path)), ("velocity", "mut", "grad_u"))
if mesh is None or mesh.n_points == 0:
continue
df = _spanwise_average_by_z(_slice_dataframe(mesh))
if df.empty:
continue
vtk_slice_name = Path(vtk_path).stem
exp_zone = _match_experimental_zone(vtk_slice_name, experimental_zones)
if exp_zone is None or exp_zone.empty:
raise RuntimeError(f"No experimental zone match for slice {vtk_slice_name!r}.")
z0 = float(exp_zone["z"].iloc[0])
cfl3d_zone = _match_experimental_zone(vtk_slice_name, cfl3d_zones)
if cfl3d_zone is None or cfl3d_zone.empty:
raise RuntimeError(f"No CFL3D zone match for {MODEL_LABELS[model_key]} slice {vtk_slice_name!r}.")
slice_name = _slice_label(vtk_slice_name)
_add_slice_rows(rows, model_key, slice_name, "Flow360", _flow_slice_data(df, z0))
_add_slice_rows(rows, model_key, slice_name, "Experimental data", _project_experimental(exp_zone, z0))
_add_slice_rows(rows, model_key, slice_name, "CFL3D", _project_experimental(cfl3d_zone, z0))
return rows
finally:
_cleanup_paths(archive_path, extract_dir)
def _plot_v_uv_field_comparison(profile_df):
if profile_df.empty:
return None
slices = sorted(profile_df["slice"].unique(), key=_slice_sort_value)
color_cycle = plt.rcParams["axes.prop_cycle"].by_key()["color"]
colors = {name: color_cycle[idx % len(color_cycle)] for idx, name in enumerate(slices)}
fig, axes = plt.subplots(len(MODEL_ORDER), len(SLICE_PLOT_COMPONENTS), figsize=(15.5, 18.0), sharey=True)
for row_idx, model_key in enumerate(MODEL_ORDER):
for col_idx, component in enumerate(SLICE_PLOT_COMPONENTS):
ax = axes[row_idx, col_idx]
subset = profile_df[
(profile_df["model_key"] == model_key) & (profile_df["component"] == component)
]
for slice_name in slices:
slice_subset = subset[subset["slice"] == slice_name]
if slice_subset.empty:
continue
color = colors[slice_name]
for source, style in (
("Flow360", {"linestyle": "-", "linewidth": 1.1}),
("CFL3D", {"linestyle": "--", "linewidth": 1.0}),
("Experimental data", {"linestyle": "none", "marker": "x", "markersize": 3}),
):
source_df = slice_subset[slice_subset["source"] == source].sort_values("z_rel")
if source_df.empty:
continue
ax.plot(source_df["value"], source_df["z_rel"], color=color, label="_nolegend_", **style)
if row_idx == 0:
ax.set_title(SLICE_COMPONENT_LABELS[component])
if col_idx == 0:
ax.set_ylabel(f"{MODEL_LABELS[model_key]}\nz - z$_0$")
if row_idx == len(MODEL_ORDER) - 1:
ax.set_xlabel(SLICE_COMPONENT_LABELS[component])
ax.set_xlim(*SLICE_COMPONENT_XLIMS[component])
ax.set_xticks(SLICE_COMPONENT_XTICKS[component])
ax.set_ylim(-0.005, SLICE_Z_MAX)
ax.set_yticks(np.arange(0.0, SLICE_Z_MAX + 0.0001, 0.02))
ax.set_box_aspect(SLICE_AX_BOX_ASPECT)
ax.grid(True, alpha=0.3)
source_handles = [
Line2D([], [], color="black", linewidth=1.1, label="Flow360"),
Line2D([], [], color="black", linestyle="--", linewidth=1.0, label="CFL3D"),
Line2D([], [], color="black", linestyle="none", marker="x", markersize=5, label="Experimental data"),
]
slice_handles = [Line2D([], [], color=colors[name], linewidth=1.2, label=name) for name in slices]
fig.legend(handles=source_handles + slice_handles, loc="upper center", ncol=5, bbox_to_anchor=(0.5, 1.0))
fig.tight_layout(rect=(0, 0, 1, 0.94))
fig.savefig("results/v_uv_field_comparison.png", dpi=250, bbox_inches="tight")
return fig
Load the case results¶
For each finished case, read the settled integrated loads from the total-forces
history, extract the mid-span surface Cp curve, and collect the boundary-layer
slice profiles (matched against the experimental and CFL3D reference zones). The
results are gathered into flow_df, surface_df, and slice_profile_df, which
drive the table and figures below.
In [ ]:
experimental_zones = _load_experimental_zones()
cfl3d_slice_zones = {model_key: _load_cfl3d_slice_zones(model_key) for model_key in MODEL_ORDER}
flow_rows = []
surface_rows = []
slice_profile_rows = []
for record in submitted_cases:
case = record["case"]
model_key = record["model_key"]
force_df = case.results.total_forces.as_dataframe()
if "physical_step" in force_df.columns:
force_df = force_df.sort_values("physical_step")
elif "pseudo_step" in force_df.columns:
force_df = force_df.sort_values("pseudo_step")
row = {"model_key": model_key, "model": MODEL_LABELS[model_key]}
for out_key, column_name in FORCE_COLUMNS.items():
row[out_key] = _mean_last_percent(force_df.get(column_name))
flow_rows.append(row)
curve_df = _extract_surface_curve(case)
curve_df["model_key"] = model_key
curve_df["model"] = MODEL_LABELS[model_key]
surface_rows.append(curve_df)
slice_profile_rows.extend(
_collect_slice_profile_rows(case, model_key, experimental_zones, cfl3d_slice_zones[model_key])
)
flow_df = pd.DataFrame(flow_rows).sort_values("model_key").reset_index(drop=True)
surface_df = pd.concat(surface_rows, ignore_index=True)
slice_profile_df = pd.DataFrame(slice_profile_rows)
Integrated loads by turbulence model¶
Tabulate the Flow360 lift and drag coefficients (total, pressure, and viscous) for each turbulence model next to the CFL3D, FUN3D, and OVERFLOW reference values.
In [ ]:
loads_table = _make_integrated_loads_table(flow_df)
loads_table.to_csv("results/integrated_loads_table.csv", index=False)
loads_table
Surface pressure coefficient by turbulence model¶
Plot the Flow360 mid-span Cp distribution against the CFL3D and FUN3D reference curves and the experimental measurements, one panel per turbulence model. The y-axis is inverted, as is conventional for Cp.
In [ ]:
cp_surface_reference_df = _load_surface_reference("cp_reference_curves.csv")
cp_experiment_df = _load_cp_reference()
_plot_surface_quantity(surface_df, "Cp", "cp_comparison", cp_surface_reference_df, cp_experiment_df)
plt.show()
Velocity and turbulent-shear profiles by turbulence model¶
For each turbulence model, plot the wall-normal profiles of axial velocity, normal velocity, and turbulent shear stress at the six chordwise stations near the trailing edge, overlaid with the CFL3D and experimental references. Colors distinguish the chordwise stations; line styles distinguish the data sources.
In [ ]:
_plot_v_uv_field_comparison(slice_profile_df)
plt.show()