2D CRM High-Lift: Multi-Element Airfoil Validation¶
This case evaluates Flow360 on the multi-element 2D CRM high-lift airfoil (slat / main / flap) from the NASA Turbulence Modeling Resource, a standard benchmark for high-lift RANS predictions. The flow is a steady, subsonic solve at Mach 0.2, Reynolds number 5.0e6, and angle of attack 16°, using the Spalart-Allmaras turbulence model on the HLPW-4 L7 workshop mesh with translational periodic boundaries.
It is a single case. Its results are validated against FUN3D reference data from the NASA TMR: the surface pressure and skin-friction distributions along the airfoil, and boundary-layer velocity/eddy-viscosity profiles at three chordwise stations (x/c = -0.03, 0.4, 0.95). The three FUN3D profile datasets use different grid families and turbulence discretizations.
The notebook runs top-to-bottom in a single kernel: load the mesh, build the simulation, submit the case, wait for completion, then post-process the in-kernel case object 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, including
pyvista to read the surface and slice output meshes:
pip install numpy pandas matplotlib pyvista
Run the cells top to bottom in a single kernel.
Imports¶
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 flow360.examples import download_benchmark_assets
Input data¶
The post-processing compares Flow360 against FUN3D reference data from the NASA
TMR. Fetch those reference files from the public benchmark bucket; this recreates
a local ./ref_data/ directory (surface Cp/Cf and velocity-profile datasets)
that the plotting cells read from.
download_benchmark_assets("2D_CRM", "ref_data")
Load project¶
The root asset for this case is the pre-generated HLPW-4 L7 workshop 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("2D_CRM", "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="2D CRM High-Lift")
Simulation parameters¶
All of the physics and outputs are assembled in build_params. The operating
condition is Mach 0.2 at Reynolds number 5.0e6 (per mesh unit), T = 272.1 K, and
α = 16°, β = 0°; the reference geometry places the moment center at quarter-chord.
The models tag the mesh boundaries by name: every *wall* surface is a viscous
wall, *farfield is the freestream, and the two y-symmetry planes form a
translational periodic pair. The fluid uses Spalart-Allmaras with a steady
adaptive-CFL ramp to 15000 steps.
Two outputs are written for post-processing: slice outputs (primitive variables
and eddy-viscosity ratio) on three x-normal planes at the profile stations, and a
surface output carrying Cp, Cf, and the friction vector on the walls.
SLICE_X_POSITIONS = (-0.03, 0.4, 0.95)
def build_params(project):
volume_mesh = project.volume_mesh
slices = [
fl.Slice(name=f"slice_{x}", normal=(1, 0, 0), origin=(x, 0, 0) * fl.u.m)
for x in SLICE_X_POSITIONS
]
with fl.SI_unit_system:
return fl.SimulationParams(
reference_geometry=fl.ReferenceGeometry(
moment_center=[0.25, 0, 0],
moment_length=[1, 1, 1],
area=1,
),
operating_condition=fl.AerospaceCondition.from_mach_reynolds(
mach=0.2,
reynolds_mesh_unit=5e6,
project_length_unit=1 * fl.u.m,
temperature=272.1,
alpha=16 * fl.u.deg,
beta=0 * fl.u.deg,
),
time_stepping=fl.Steady(
max_steps=15000,
CFL=fl.AdaptiveCFL(
max=100,
max_relative_change=50,
convergence_limiting_factor=0.9,
),
),
models=[
fl.Wall(surfaces=volume_mesh["*wall*"]),
fl.Freestream(surfaces=volume_mesh["*farfield"]),
fl.Periodic(
surface_pairs=[(volume_mesh["y-symmetry plane1"], volume_mesh["y-symmetry plane2"])],
spec=fl.Translational(),
),
fl.Fluid(
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-11,
linear_solver=fl.LinearSolver(max_iterations=35),
kappa_MUSCL=0.33,
),
turbulence_model_solver=fl.SpalartAllmaras(
absolute_tolerance=1e-10,
linear_solver=fl.LinearSolver(max_iterations=25),
equation_evaluation_frequency=1,
),
),
],
outputs=[
fl.SliceOutput(
output_fields=["primitiveVars", "mutRatio"],
output_format="paraview",
slices=slices,
),
fl.SurfaceOutput(
surfaces=volume_mesh["*wall*"],
output_fields=["Cp", "Cf", "CfVec"],
output_format="paraview",
),
],
)
Submit case¶
Submit the single case to Flow360. It is named L7-WMesh-15k_<solver_version>.
The submitted case object is kept in the kernel so the post-processing can read
its results directly.
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")
CASE_NAME = f"L7-WMesh-15k_{SOLVER_VERSION}"
case = project.run_case(
params=build_params(project),
name=CASE_NAME,
solver_version=SOLVER_VERSION,
)
Wait for completion¶
Block until the case has finished. This can take a long time (minutes to hours), since it is a full steady solve of up to 15000 steps.
case.wait()
Postprocessing¶
The published results are three figures:
- Surface pressure coefficient: Cp vs x/c along the airfoil, Flow360 vs FUN3D.
- Surface skin-friction coefficient: Cf vs x/c along the airfoil, Flow360 vs FUN3D.
- Velocity profiles: u, w, and eddy-viscosity ratio vs z/c at three x/c stations, Flow360 vs three FUN3D reference datasets.
All three read the results straight from the in-kernel case object (no project
is re-opened). The cells below build exactly these figures.
Setup and shared helpers¶
Styling colors, the reference-data file locations, and helpers used by both the
surface and profile plots: extract_archive downloads a result archive and
returns the VTK files inside, and prepare_mesh promotes cell data to point data
and merges multi-block meshes. read_tecplot_zones parses the FUN3D Tecplot
reference files into per-zone dataframes.
FLOW360_COLOR = "#00643c" # Flow360 results
REFERENCE_COLOR = "black" # FUN3D reference data
RESULTS_DIR = Path("results")
RESULTS_DIR.mkdir(exist_ok=True)
SURFACE_REFERENCE_PATH = Path("ref_data") / "fun3d_cp_cf.dat"
PROFILE_REFERENCE_PATH = Path("ref_data") / "fun3d_velocity_profiles.dat"
def extract_archive(results_obj, archive_path, extract_dir):
archive_path.parent.mkdir(parents=True, exist_ok=True)
extract_dir.mkdir(parents=True, exist_ok=True)
results_obj.to_file(str(archive_path), overwrite=True)
with tarfile.open(archive_path, "r:gz") as tar:
tar.extractall(extract_dir)
aggregate_files = []
for suffix in (".pvtu", ".pvtp"):
aggregate_files.extend(sorted(extract_dir.rglob(f"*{suffix}")))
if aggregate_files:
return aggregate_files
files = []
for suffix in (".vtu", ".vtp", ".vtk"):
files.extend(sorted(extract_dir.rglob(f"*{suffix}")))
return files
def prepare_mesh(mesh, required_point_arrays=()):
if isinstance(mesh, pv.MultiBlock):
blocks = [
prepare_mesh(block, required_point_arrays)
for block in mesh
if block is not None and getattr(block, "n_points", 0) > 0
]
return pv.merge(blocks, merge_points=False) if blocks else pv.PolyData()
if not mesh.point_data and mesh.cell_data:
return mesh.cell_data_to_point_data(pass_cell_data=True)
missing = [name for name in required_point_arrays if name in mesh.cell_data and name not in mesh.point_data]
return mesh.cell_data_to_point_data(pass_cell_data=True) if missing else mesh
def read_tecplot_zones(path):
lines = Path(path).read_text(encoding="utf-8").splitlines()
variables = []
zones = []
i = 0
while i < len(lines):
line = lines[i].strip()
if not line:
i += 1
continue
if line.upper().startswith("VARIABLES"):
variables = re.findall(r'"([^"]+)"', line)
i += 1
continue
if not line.upper().startswith("ZONE"):
i += 1
continue
match = re.search(r'T\s*=\s*"([^"]+)"', line, flags=re.IGNORECASE)
zone_name = match.group(1) if match else f"zone_{len(zones) + 1}"
nodes_match = re.search(r"Nodes\s*=\s*(\d+)", line, flags=re.IGNORECASE)
elements_match = re.search(r"Elements\s*=\s*(\d+)", line, flags=re.IGNORECASE)
nodes = int(nodes_match.group(1)) if nodes_match else None
elements = int(elements_match.group(1)) if elements_match else 0
i += 1
while i < len(lines):
header = lines[i].strip()
if not header:
i += 1
continue
if header.upper().startswith("ZONE") or re.match(r"^[+\-]?\d", header):
break
match = re.search(r"Nodes\s*=\s*(\d+)", header, flags=re.IGNORECASE)
if match:
nodes = int(match.group(1))
match = re.search(r"Elements\s*=\s*(\d+)", header, flags=re.IGNORECASE)
if match:
elements = int(match.group(1))
i += 1
if nodes is None or not variables:
continue
rows = []
read_rows = 0
while i < len(lines) and read_rows < nodes:
row = lines[i].strip()
if row:
values = row.split()
if len(values) >= len(variables):
rows.append({key: float(val) for key, val in zip(variables, values[: len(variables)])})
read_rows += 1
i += 1
i += elements
zones.append((zone_name, pd.DataFrame(rows)))
return zones
Extract the surface Cp/Cf slice¶
Download the surface output, read every wall VTK, group them into slat / main /
flap, and slice each along the airfoil mid-plane (a y-normal plane at y = -0.5).
The result is a list of (group, dataframe) with X, Z, Cp, and Cf on the
sliced surface, the raw data for both surface plots.
SLICE_NORMAL = (0.0, 1.0, 0.0)
SLICE_ORIGIN = (0.0, -0.5, 0.0)
def group_name(name):
lowered = str(name).lower()
if "flap" in lowered:
return "flap"
if "slat" in lowered:
return "slat"
if "wing" in lowered or "main" in lowered:
return "wing"
return None
def extract_surface_data(case):
tmp_dir = RESULTS_DIR / "_surface_tmp"
archive_path = tmp_dir / "surfaces.tar.gz"
extract_dir = tmp_dir / "extract"
files = extract_archive(case.results.surfaces, archive_path, extract_dir)
grouped = {"flap": [], "slat": [], "wing": []}
for path in files:
group = group_name(path.stem)
if group is None:
continue
mesh = prepare_mesh(pv.read(path), ("Cp", "Cf"))
sliced = mesh.slice(normal=SLICE_NORMAL, origin=SLICE_ORIGIN)
if sliced.n_points == 0:
continue
grouped[group].append(
pd.DataFrame(
{
"X": sliced.points[:, 0],
"Z": sliced.points[:, 2],
"Cp": np.asarray(sliced.point_data["Cp"]).reshape(-1),
"Cf": np.asarray(sliced.point_data["Cf"]).reshape(-1),
}
)
)
output = []
for group, parts in grouped.items():
if parts:
output.append((group, pd.concat(parts, ignore_index=True)))
shutil.rmtree(tmp_dir, ignore_errors=True)
return output
surface_flow_data = extract_surface_data(case)
Surface plotting helpers¶
order_points reconnects the sliced points into a continuous curve along the
airfoil surface (nearest-neighbor walk, oriented clockwise), and split_by_jump
breaks the curve at large gaps so distinct surfaces are not joined by spurious
lines. parse_surface_reference reads the FUN3D Cp/Cf reference for a variable.
def is_clockwise(data):
if len(data) < 3:
return True
area2 = np.sum(data["X"] * np.roll(data["Z"], -1) - np.roll(data["X"], -1) * data["Z"])
return area2 < 0.0
def order_points(data):
work = data.reset_index(drop=True)
if len(work) <= 2:
return work
points = work[["X", "Z"]].to_numpy()
order = [0]
unvisited = set(range(1, len(points)))
while unvisited:
last = order[-1]
idx = np.fromiter(unvisited, dtype=int)
deltas = points[idx] - points[last]
next_idx = int(idx[int(np.argmin(np.einsum("ij,ij->i", deltas, deltas)))])
order.append(next_idx)
unvisited.remove(next_idx)
ordered = work.iloc[order].reset_index(drop=True)
if not is_clockwise(ordered):
ordered = ordered.iloc[[0] + list(range(len(ordered) - 1, 0, -1))].reset_index(drop=True)
return ordered
def split_by_jump(data, jump_factor=12.0):
if len(data) <= 2:
return [data]
distances = np.sqrt(np.sum(np.diff(data[["X", "Z"]].to_numpy(), axis=0) ** 2, axis=1))
positive = distances[distances > 0]
if len(positive) == 0:
return [data]
threshold = float(np.median(positive)) * jump_factor
breaks = np.where(distances > threshold)[0]
if len(breaks) == 0:
return [data]
start = 0
segments = []
for idx in breaks:
segment = data.iloc[start : idx + 1].reset_index(drop=True)
if not segment.empty:
segments.append(segment)
start = idx + 1
tail = data.iloc[start:].reset_index(drop=True)
if not tail.empty:
segments.append(tail)
return segments
SURFACE_LIMITS = {"Cp": (-14, 2), "Cf": (-0.005, 0.06)}
SURFACE_LABELS = {"Cp": r"$C_p$", "Cf": r"$C_f$"}
def parse_surface_reference(variable):
column = "cp" if variable == "Cp" else "abs(cf)"
zones = []
for name, data in read_tecplot_zones(SURFACE_REFERENCE_PATH):
group = group_name(name)
if group is None or "x/c" not in data.columns or column not in data.columns:
continue
zone = data[["x/c", column]].rename(columns={"x/c": "X", column: "value"})
zone["group"] = group
zone["source"] = "FUN3D"
zones.append(zone)
return zones
def plot_surface_variable(flow_data, variable):
reference = parse_surface_reference(variable)
plt.figure(figsize=(8, 5))
flow_label_used = False
for group, data in flow_data:
ordered = order_points(data[["X", "Z", variable]].rename(columns={variable: "value"}))
for segment in split_by_jump(ordered):
plt.plot(
segment["X"],
segment["value"],
color=FLOW360_COLOR,
linewidth=1.8,
label="Flow360" if not flow_label_used else None,
)
flow_label_used = True
ref_label_used = False
for zone in reference:
plt.plot(
zone["X"],
zone["value"],
linestyle="None",
marker=".",
markersize=2.0,
color=REFERENCE_COLOR,
label="FUN3D" if not ref_label_used else None,
)
ref_label_used = True
plt.xlabel("x/c")
plt.ylabel(SURFACE_LABELS[variable])
plt.ylim(*SURFACE_LIMITS[variable])
if variable == "Cp":
plt.gca().invert_yaxis()
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
stem = "cp_x" if variable == "Cp" else "cf_x"
plt.savefig(f"results/{stem}.png", dpi=300)
plt.show()
Surface pressure coefficient¶
Cp vs x/c along the airfoil, Flow360 vs the FUN3D reference. The y-axis is
inverted, as is conventional for Cp. Saved to results/cp_x.png.
plot_surface_variable(surface_flow_data, "Cp")
Surface skin-friction coefficient¶
Cf vs x/c along the airfoil, Flow360 vs the FUN3D reference. Saved to
results/cf_x.png.
plot_surface_variable(surface_flow_data, "Cf")
Velocity profiles at three x/c stations¶
Download the slice outputs and, at each of the three stations (slat x/c = -0.03,
main x/c = 0.4, flap x/c = 0.95), extract non-dimensional u, w (velocity
components normalized by the freestream Mach 0.2), and eddy-viscosity ratio versus
z/c. These are binned in z to a single profile per station and compared against
three FUN3D reference datasets (two grid families, one with second-order
turbulence). The figure is a 3×3 grid of station × quantity, saved to
results/velocity_profiles.png.
PROFILE_ROWS = ("slat", "main", "flap")
PROFILE_QUANTITIES = ("u", "w", "mut")
PROFILE_TITLES = {"slat": "x/c = -0.03", "main": "x/c = 0.4", "flap": "x/c = 0.95"}
PROFILE_WINDOWS = {
"slat": (-0.0448, -0.0378),
"main": (0.06, 0.10),
"flap": (-0.025, 0.10),
}
PROFILE_X_LIMITS = {
("slat", "u"): (0.0, 2.0),
("slat", "w"): (0.0, 2.0),
("slat", "mut"): (0.0, 100.0),
("main", "u"): (0.0, 2.0),
("main", "w"): (0.0, 0.1),
("main", "mut"): (0.0, 600.0),
("flap", "u"): (0.0, 2.0),
("flap", "w"): (-1.0, 0.0),
("flap", "mut"): (0.0, 800.0),
}
PROFILE_X_TICKS = {
("slat", "u"): [0.0, 0.5, 1.0, 1.5, 2.0],
("slat", "w"): [0.0, 0.5, 1.0, 1.5, 2.0],
("slat", "mut"): [0.0, 50.0, 100.0],
("main", "u"): [0.0, 0.5, 1.0, 1.5, 2.0],
("main", "w"): [0.0, 0.02, 0.04, 0.06, 0.08, 0.10],
("main", "mut"): [0.0, 200.0, 400.0, 600.0],
("flap", "u"): [0.0, 0.5, 1.0, 1.5, 2.0],
("flap", "w"): [-1.0, -0.8, -0.6, -0.4, -0.2, 0.0],
("flap", "mut"): [0.0, 200.0, 400.0, 600.0, 800.0],
}
PROFILE_LABELS = {"u": r"u/U$_{ref}$", "w": r"w/U$_{ref}$", "mut": r"$\mu_t/\mu_{ref}$"}
REFERENCE_COLORS = {"family1": "#cc0000", "family1_2nd": "#111111", "family2": "#00b020"}
REFERENCE_LINESTYLES = {"family1": "-", "family1_2nd": "--", "family2": ":"}
def profile_family(name):
lowered = name.lower()
if "2nd order turb" in lowered:
return "family1_2nd"
if "family 2" in lowered:
return "family2"
return "family1"
def profile_row_from_x(value):
target = {"slat": -0.03, "main": 0.4, "flap": 0.95}
return min(target, key=lambda key: abs(target[key] - value))
def profile_row_from_name(name):
match = re.search(r"(-?\d+(?:\.\d+)?)", name)
return profile_row_from_x(float(match.group(1))) if match else "main"
def parse_profile_reference():
zones = {}
for name, data in read_tecplot_zones(PROFILE_REFERENCE_PATH):
if {"x/c", "z/c", "u/Uref", "w/Uref", "mu_t/mu_ref"} - set(data.columns):
continue
row = profile_row_from_x(float(data["x/c"].mean()))
family = profile_family(name)
zones[(row, family)] = data.rename(
columns={"z/c": "z", "u/Uref": "u", "w/Uref": "w", "mu_t/mu_ref": "mut"}
)[["z", "u", "w", "mut"]]
return zones
def extract_profiles(case):
tmp_dir = RESULTS_DIR / "_slice_tmp"
archive_path = tmp_dir / "slices.tar.gz"
extract_dir = tmp_dir / "extract"
files = extract_archive(case.results.slices, archive_path, extract_dir)
profiles = {}
for path in files:
mesh = prepare_mesh(pv.read(path), ("velocity", "mutRatio"))
if mesh is None or getattr(mesh, "n_points", 0) == 0:
continue
velocity = np.asarray(mesh.point_data["velocity"])
if velocity.ndim != 2 or velocity.shape[1] < 3:
continue
row = profile_row_from_name(path.stem)
data = pd.DataFrame(
{
"z": mesh.points[:, 2],
"u": velocity[:, 0] / 0.2,
"w": velocity[:, 2] / 0.2,
"mut": np.asarray(mesh.point_data["mutRatio"]).reshape(-1),
}
)
data["z_bin"] = (data["z"] / 0.0001).round().astype(int)
profiles[row] = (
data.groupby("z_bin", as_index=False)[["z", "u", "w", "mut"]]
.mean()
.sort_values("z")
.reset_index(drop=True)
)
shutil.rmtree(tmp_dir, ignore_errors=True)
return profiles
def profile_window(data, row, quantity):
zmin, zmax = PROFILE_WINDOWS[row]
subset = data[(data["z"] >= zmin) & (data["z"] <= zmax)][["z", quantity]].dropna()
return subset.sort_values("z")
reference_profiles = parse_profile_reference()
flow_profiles = extract_profiles(case)
fig, axes = plt.subplots(3, 3, figsize=(12, 11))
for i, row in enumerate(PROFILE_ROWS):
for j, quantity in enumerate(PROFILE_QUANTITIES):
ax = axes[i][j]
for family in ("family1", "family1_2nd", "family2"):
ref = reference_profiles.get((row, family))
if ref is None:
continue
window = profile_window(ref, row, quantity)
if window.empty:
continue
ax.plot(
window[quantity],
window["z"],
color=REFERENCE_COLORS[family],
linestyle=REFERENCE_LINESTYLES[family],
linewidth=1.5,
label=f"FUN3D {family.replace('_', ' ')}",
)
if row in flow_profiles:
window = profile_window(flow_profiles[row], row, quantity)
if not window.empty:
ax.plot(
window[quantity],
window["z"],
color=FLOW360_COLOR,
linewidth=1.8,
label="Flow360",
)
ax.set_xlim(*PROFILE_X_LIMITS[(row, quantity)])
ax.set_ylim(*PROFILE_WINDOWS[row])
ax.set_xticks(PROFILE_X_TICKS[(row, quantity)])
ax.grid(True, alpha=0.3)
ax.set_title(PROFILE_TITLES[row], fontsize=12)
ax.set_xlabel(PROFILE_LABELS[quantity], fontsize=11)
ax.set_ylabel("z/c", fontsize=11)
handles, labels = axes[0][0].get_legend_handles_labels()
dedup_handles = []
dedup_labels = []
for handle, label in zip(handles, labels):
if label not in dedup_labels:
dedup_handles.append(handle)
dedup_labels.append(label)
axes[0][0].legend(dedup_handles, dedup_labels, loc="upper left", fontsize=9)
fig.subplots_adjust(left=0.07, right=0.98, top=0.97, bottom=0.06, wspace=0.42, hspace=0.42)
fig.savefig("results/velocity_profiles.png", dpi=300)
plt.show()