NLF Airfoil: Natural Laminar Flow Validation¶
This case validates Flow360 on the NLF(1)-0416 natural-laminar-flow airfoil, a benchmark for boundary-layer transition modeling that drives laminar-flow extent and drag on modern wings. All cases run at Mach 0.1, Re = 4×10⁶ (mesh unit) on a structured C-grid volume mesh.
The study is an angle-of-attack sweep from −6° up to ~15°, run twice with two turbulence configurations so the effect of transition is isolated:
- SA-AFT: Spalart-Allmaras with the γ-Reθ transition model (N_crit = 7.2).
- SA: Spalart-Allmaras, fully turbulent (no transition model).
Each sweep forks sequentially outward from α = 0° to warm-start every case. Integrated loads and skin-friction distributions are compared against experiment (NASA TP-1861) and reference CFD codes (OVERFLOW, FUN3D). The notebook runs top-to-bottom in one kernel: load the mesh, submit both sweeps, wait, then post-process the in-kernel case objects into the 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 uses a few data/plotting/mesh packages:
pip install numpy pandas matplotlib pyvista
Run the cells top to bottom in a single kernel.
Imports¶
import os
import tarfile
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¶
Fetch the case's reference data (experiment, OVERFLOW, FUN3D) from the public
benchmark bucket; this recreates a local ./ref_data/ directory that both the sweep
definition and the post-processing read from. The experimental CL_alpha file also
sets the upper bound of the angle-of-attack sweep, and constants name the mesh
boundaries and solver settings.
download_benchmark_assets("NLF_airfoil", "ref_data")
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")
WALL = "blk-1/airfoil"
FARFIELD = "blk-1/farfield"
SLIPWALL = "blk-1/sides"
MAX_STEPS = 2e4
# Alpha sweep: 2° steps from -6° to 14° (incl. 5°), then 0.5° steps above 14°
# (upper limit taken from the experimental data).
alpha_max = pd.read_csv(
os.path.join("ref_data", "CL_alpha_Experiment.csv"), skipinitialspace=True
)["N"].max()
alphas_coarse = np.sort(np.unique(np.append(np.arange(-6.0, 16.0, 2.0), 5.0)))
alphas_fine = np.arange(14.5, alpha_max + 0.5, 0.5) if alpha_max > 14.0 else np.array([])
alphas = np.unique(np.concatenate([alphas_coarse, alphas_fine]))
Load project¶
The root asset for this case is a pre-generated structured C-grid volume mesh (built in Pointwise). 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("NLF_airfoil", "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="NLF Airfoil")
vm = project.volume_mesh
Physics setup¶
The airfoil surface is a viscous wall, the outer boundary is a freestream, and the
spanwise sides are slip walls. Two fluid models drive the two sweeps: fluid_transition
adds the γ-Reθ transition model (N_crit = 7.2) on top of Spalart-Allmaras, while
fluid_sa is fully turbulent Spalart-Allmaras. build_models assembles the boundary
conditions and the chosen fluid into the model list for a run.
with fl.SI_unit_system:
wall_surf = vm[WALL]
wall_model = fl.Wall(surfaces=wall_surf)
fluid_transition = fl.Fluid(
turbulence_model_solver=fl.SpalartAllmaras(
absolute_tolerance=1e-10,
update_jacobian_frequency=1,
equation_evaluation_frequency=1,
linear_solver=fl.LinearSolver(max_iterations=35),
),
transition_model_solver=fl.TransitionModelSolver(
absolute_tolerance=1e-10,
N_crit=7.2,
update_jacobian_frequency=1,
equation_evaluation_frequency=1,
linear_solver=fl.LinearSolver(max_iterations=35),
),
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-12,
update_jacobian_frequency=4,
equation_evaluation_frequency=1,
linear_solver=fl.LinearSolver(max_iterations=35),
),
)
fluid_sa = fl.Fluid(
turbulence_model_solver=fl.SpalartAllmaras(
absolute_tolerance=1e-10,
linear_solver=fl.LinearSolver(max_iterations=35),
),
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-12,
linear_solver=fl.LinearSolver(max_iterations=35),
),
)
def build_models(fluid):
with fl.SI_unit_system:
return [
wall_model,
fl.Freestream(surfaces=[vm[FARFIELD]]),
fl.SlipWall(surfaces=[vm[SLIPWALL]]),
fluid,
]
Outputs setup¶
Write surface Cp, Cf, yPlus, and CfVec on the airfoil so the post-processing
can extract the skin-friction distribution.
def build_outputs():
with fl.SI_unit_system:
return [
fl.SurfaceOutput(
surfaces=wall_surf,
output_format="both",
output_fields=["Cp", "Cf", "yPlus", "CfVec"],
),
]
Simulation Params¶
make_params assembles the full SimulationParams shared by every case: the
Mach 0.1 / Re 4×10⁶ operating condition, a steady ramp-CFL schedule, the boundary
conditions and fluid, the surface outputs, and the reference geometry. Only the angle
of attack changes between cases. run_alpha_sweep submits one sweep: it runs α = 0°
first, then marches outward to positive and negative angles, forking each case from
the previous one to warm-start it.
def make_params(fluid):
with fl.SI_unit_system:
return fl.SimulationParams(
operating_condition=fl.AerospaceCondition.from_mach_reynolds(
mach=0.1,
reynolds_mesh_unit=4e6,
project_length_unit=1 * fl.u.m,
temperature=540.0 * fl.u.R,
alpha=0.0 * fl.u.deg,
beta=0.0 * fl.u.deg,
reference_mach=0.1,
),
time_stepping=fl.Steady(
max_steps=MAX_STEPS,
CFL=fl.RampCFL(initial=1, final=100, ramp_steps=2000),
),
models=build_models(fluid),
outputs=build_outputs(),
reference_geometry=fl.ReferenceGeometry(moment_center=(0, 0.5, 0.5) * fl.u.m),
)
def run_alpha_sweep(params, suffix):
collected = []
pos = sorted(a for a in alphas if a > 0)
neg = sorted((a for a in alphas if a < 0), reverse=True)
params.operating_condition.alpha = 0.0 * fl.u.deg
base = project.run_case(
params, name=f"alpha_0.00{suffix}_{SOLVER_VERSION}", solver_version=SOLVER_VERSION
)
collected.append(base)
parent = base
for a in pos:
params.operating_condition.alpha = a * fl.u.deg
parent = project.run_case(
params,
name=f"alpha_{a:.2f}{suffix}_{SOLVER_VERSION}",
fork_from=parent,
solver_version=SOLVER_VERSION,
)
collected.append(parent)
parent = base
for a in neg:
params.operating_condition.alpha = a * fl.u.deg
parent = project.run_case(
params,
name=f"alpha_{a:.2f}{suffix}_{SOLVER_VERSION}",
fork_from=parent,
solver_version=SOLVER_VERSION,
)
collected.append(parent)
return collected
Submit cases¶
Run both sweeps, first the transition (SA-AFT) configuration, then the fully turbulent
SA configuration, and collect every submitted case into submitted for the
post-processing to reuse directly.
submitted = []
params_transition = make_params(fluid_transition)
params_sa = make_params(fluid_sa)
submitted.extend(run_alpha_sweep(params_transition, "_with_transition"))
submitted.extend(run_alpha_sweep(params_sa, "_no_transition"))
print(f"Launched {len(submitted)} cases")
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 case in submitted:
case.wait()
Postprocessing¶
The published results are three figures, all comparing Flow360 against experiment and reference CFD codes:
- CL vs α: lift curve over the sweep (SA-AFT and SA).
- CL vs CD: drag polar (SA-AFT and SA).
- Cf vs x/c at α = 5°: skin-friction distribution, showing the laminar-to-turbulent transition captured by the SA-AFT model.
The cells below read results straight from the in-kernel case objects (no project is re-opened) and rebuild exactly these three figures.
Setup¶
Plot colors (Flow360 brand green; black/red/blue for the references), the sweep range
used for the polar plots, and two small helpers: case_alpha reads a case's angle of
attack, and load_ref reads a reference CSV and normalizes its column names.
FLOW360_COLOR = "#00643c" # Flow360 results
REFERENCE_COLOR = "black" # experiment
OVERFLOW_COLOR = "#C00" # OVERFLOW (transition / SA-AFT)
OVERFLOW_SA_COLOR = "#1f77b4" # OVERFLOW (SA)
ref_data_dir = "ref_data"
os.makedirs("results", exist_ok=True)
ALPHA_MIN, ALPHA_MAX = -6.0, 15.0
def case_alpha(case):
return case.params.operating_condition.alpha.to("deg").to_value()
def load_ref(filename, col_map):
df = pd.read_csv(os.path.join(ref_data_dir, filename), skipinitialspace=True)
df.columns = df.columns.str.strip()
return df.rename(columns={k: v for k, v in col_map.items() if k in df.columns})
Load forces from the swept cases¶
For every finished case, take the settled lift and drag coefficients (average over the last 1/12 of the convergence history) and tag each with its angle of attack and whether it used the transition model. Then load the experiment and OVERFLOW reference coefficients for the two polar plots.
cases = []
for case in submitted:
if case.status.value != "completed":
print(f"Skipping {case.name} (status={case.status.value})")
continue
cases.append((case, case.name, "_with_transition" in case.name))
n_tr = sum(tr for _, _, tr in cases)
print(f"Cases: {n_tr} with transition, {len(cases) - n_tr} without")
rows = []
for case, case_name, has_transition in cases:
CF = case.results.total_forces.get_averages(1 / 12)
cl = float(CF["CL"]) if "CL" in CF else None
cd = float(CF["CD"]) if "CD" in CF else None
if cl is not None and cd is not None:
rows.append(
{
"alpha": case_alpha(case),
"CL": cl,
"CD": cd,
"transition": "yes" if has_transition else "no",
}
)
df_all = (
pd.DataFrame(rows).sort_values(["transition", "alpha"])
if rows
else pd.DataFrame(columns=["alpha", "CL", "CD", "transition"])
)
df_tr = df_all[df_all["transition"] == "yes"].copy()
df_no_tr = df_all[df_all["transition"] == "no"].copy()
# Reference coefficients: experiment, OVERFLOW (transition / SA-AFT), OVERFLOW (SA).
df_cl_alpha_exp = load_ref("CL_alpha_Experiment.csv", {"N": "alpha", "CL": "cl"})
df_cl_alpha_overflow = load_ref("CL_alpha_OVERFLOW_Coder.csv", {"N": "alpha", "CL": "cl"})
df_cl_alpha_overflow_sa = load_ref("CL_alpha_OVERFLOW_Coder_SA.csv", {"N": "alpha", "CL": "cl"})
df_cl_cd_exp = load_ref("CL_CD_Experiment.csv", {"CL": "cl", "CD": "cd"})
df_cl_cd_overflow = load_ref("CL_CD_OVERFLOW_Coder.csv", {"CL": "cl", "CD": "cd"})
df_cl_cd_overflow_sa = load_ref("CL_CD_OVERFLOW_Coder_SA.csv", {"CL": "cl", "CD": "cd"})
CL vs α¶
Lift curve over the sweep range (−6° to 15°). Flow360 SA-AFT and SA results are drawn as lines with markers; the experiment and both OVERFLOW references are scatter points.
df_alpha_exp_r = (
df_cl_alpha_exp[(df_cl_alpha_exp["alpha"] >= ALPHA_MIN) & (df_cl_alpha_exp["alpha"] <= ALPHA_MAX)]
if not df_cl_alpha_exp.empty
else pd.DataFrame()
)
df_alpha_of_r = (
df_cl_alpha_overflow[(df_cl_alpha_overflow["alpha"] >= ALPHA_MIN) & (df_cl_alpha_overflow["alpha"] <= ALPHA_MAX)]
if not df_cl_alpha_overflow.empty
else pd.DataFrame()
)
df_alpha_of_sa_r = (
df_cl_alpha_overflow_sa[(df_cl_alpha_overflow_sa["alpha"] >= ALPHA_MIN) & (df_cl_alpha_overflow_sa["alpha"] <= ALPHA_MAX)]
if not df_cl_alpha_overflow_sa.empty
else pd.DataFrame()
)
df_tr_alpha_r = (
df_tr[(df_tr["alpha"] >= ALPHA_MIN) & (df_tr["alpha"] <= ALPHA_MAX)].sort_values("alpha")
if not df_tr.empty
else pd.DataFrame()
)
df_no_tr_alpha_r = (
df_no_tr[(df_no_tr["alpha"] >= ALPHA_MIN) & (df_no_tr["alpha"] <= ALPHA_MAX)].sort_values("alpha")
if not df_no_tr.empty
else pd.DataFrame()
)
fig, ax = plt.subplots(figsize=(10, 6))
if not df_alpha_exp_r.empty:
ax.scatter(df_alpha_exp_r["alpha"], df_alpha_exp_r["cl"], color=REFERENCE_COLOR, s=35, marker="o", label="Experiment", alpha=0.9, zorder=2, edgecolors="none")
if not df_alpha_of_r.empty:
ax.scatter(df_alpha_of_r["alpha"], df_alpha_of_r["cl"], color=OVERFLOW_COLOR, s=30, marker="s", label="OVERFLOW (SA-AFT)", alpha=0.9, zorder=2, edgecolors="none")
if not df_alpha_of_sa_r.empty:
ax.scatter(df_alpha_of_sa_r["alpha"], df_alpha_of_sa_r["cl"], color=OVERFLOW_SA_COLOR, s=30, marker="^", label="OVERFLOW (SA)", alpha=0.9, zorder=2, edgecolors="none")
if not df_tr_alpha_r.empty:
ax.plot(df_tr_alpha_r["alpha"], df_tr_alpha_r["CL"], color=FLOW360_COLOR, linestyle="-", marker="o", markersize=5, label="Flow360 (SA-AFT)", linewidth=2, zorder=2)
if not df_no_tr_alpha_r.empty:
ax.plot(df_no_tr_alpha_r["alpha"], df_no_tr_alpha_r["CL"], color=FLOW360_COLOR, linestyle="--", marker="s", markersize=5, label="Flow360 (SA)", zorder=2)
ax.set(xlabel="Alpha (deg)", ylabel="CL", xlim=(ALPHA_MIN, ALPHA_MAX))
ax.legend(loc="best", fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("results/CL_alpha_comparison.png", dpi=150)
plt.show()
CL vs CD (drag polar)¶
Drag polar over the same range. Flow360 is drawn as a line ordered by angle of attack; the references are scatter points (their source files do not carry α).
df_tr_cd_plot = (
df_tr[(df_tr["alpha"] >= ALPHA_MIN) & (df_tr["alpha"] <= ALPHA_MAX)].sort_values("alpha")
if not df_tr.empty
else pd.DataFrame()
)
df_no_tr_cd_plot = (
df_no_tr[(df_no_tr["alpha"] >= ALPHA_MIN) & (df_no_tr["alpha"] <= ALPHA_MAX)].sort_values("alpha")
if not df_no_tr.empty
else pd.DataFrame()
)
fig, ax = plt.subplots(figsize=(10, 6))
if not df_cl_cd_exp.empty:
ax.scatter(df_cl_cd_exp["cd"], df_cl_cd_exp["cl"], color=REFERENCE_COLOR, s=35, marker="o", label="Experiment", alpha=0.9, zorder=2, edgecolors="none")
if not df_cl_cd_overflow.empty:
ax.scatter(df_cl_cd_overflow["cd"], df_cl_cd_overflow["cl"], color=OVERFLOW_COLOR, s=30, marker="s", label="OVERFLOW (SA-AFT)", alpha=0.9, zorder=2, edgecolors="none")
if not df_cl_cd_overflow_sa.empty:
ax.scatter(df_cl_cd_overflow_sa["cd"], df_cl_cd_overflow_sa["cl"], color=OVERFLOW_SA_COLOR, s=30, marker="^", label="OVERFLOW (SA)", alpha=0.9, zorder=2, edgecolors="none")
if not df_tr_cd_plot.empty:
ax.plot(df_tr_cd_plot["CD"], df_tr_cd_plot["CL"], color=FLOW360_COLOR, linestyle="-", marker="o", markersize=5, label="Flow360 (SA-AFT)", linewidth=2, zorder=2)
if not df_no_tr_cd_plot.empty:
ax.plot(df_no_tr_cd_plot["CD"], df_no_tr_cd_plot["CL"], color=FLOW360_COLOR, linestyle="--", marker="s", markersize=5, label="Flow360 (SA)", zorder=2)
ax.set(xlabel="CD", ylabel="CL", xlim=(0, 0.05))
ax.legend(loc="best", fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("results/CL_CD_comparison.png", dpi=150)
plt.show()
Cf vs x/c at α = 5°¶
Download the surface output of the transition (SA-AFT) case at α = 5°, extract the
airfoil-section skin-friction distribution at the mid-plane (y = 0), and compare it
against the OVERFLOW and FUN3D references. sort_by_arc_length walks the extracted
surface points into a continuous chordwise path (upper then lower surface) so the
distribution plots cleanly.
def sort_by_arc_length(x, y, z, field, start_from_le=True, x_direction=1):
if len(x) <= 1:
return x, field
start_idx = np.argmin(x) if start_from_le else np.argmax(x)
sample_size = min(100, len(x))
all_distances = [
np.sqrt((x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2 + (z[i] - z[j]) ** 2)
for i in range(sample_size)
for j in range(i + 1, min(i + 10, sample_size))
]
max_jump = np.percentile(all_distances, 50) * 5 if all_distances else np.inf
visited = np.zeros(len(x), dtype=bool)
path = [start_idx]
visited[start_idx] = True
cur = start_idx
while len(path) < len(x):
mask = ~visited
if not np.any(mask):
break
dx = x[mask] - x[cur]
dy = y[mask] - y[cur]
dz = z[mask] - z[cur]
dist = np.sqrt(dx**2 + dy**2 + dz**2)
idxs = np.where(mask)[0]
valid = dist <= max_jump
if np.any(valid):
vi = idxs[valid]
vd = dist[valid]
x_score = np.where((x[vi] - x[cur]) * x_direction > 0, 1.0, 0.1)
dist_score = 1.0 / (vd + 1e-10)
dist_score /= dist_score.max()
cur = vi[np.argmax(0.6 * x_score + 0.4 * dist_score)]
else:
cur = idxs[np.argmin(dist)]
path.append(cur)
visited[cur] = True
p = np.array(path)
return x[p], field[p]
def extract_surface_data(vtu_path, field_name="Cp", target_y=0.0, tolerance=1e-6):
mesh = pv.read(vtu_path)
x, y, z = mesh.points[:, 0], mesh.points[:, 1], mesh.points[:, 2]
mask = np.abs(y - target_y) < tolerance
xs, ys, zs = x[mask], y[mask], z[mask]
if field_name in mesh.point_data:
fv = mesh.point_data[field_name][mask]
elif field_name in mesh.cell_data:
fv = mesh.cell_data_to_point_data().point_data[field_name][mask]
else:
return None, None
if len(xs) == 0:
return None, None
up = zs > 1e-6
lo = zs < -1e-6
if np.sum(up) < len(zs) * 0.1 or np.sum(lo) < len(zs) * 0.1:
med = np.median(zs)
up, lo = zs > med, zs < med
xu, fu = sort_by_arc_length(xs[up], ys[up], zs[up], fv[up], True, 1) if np.sum(up) else (np.array([]), np.array([]))
xl, fl_ = sort_by_arc_length(xs[lo], ys[lo], zs[lo], fv[lo], False, -1) if np.sum(lo) else (np.array([]), np.array([]))
return np.concatenate([xu, xl]), np.concatenate([fu, fl_])
def find_vtu_file(case_dir):
for d in [case_dir, os.path.join(case_dir, "surfaces")]:
if not os.path.exists(d):
continue
canonical = os.path.join(d, "surface_Block_Aerofoil_proc0.vtu")
if os.path.exists(canonical):
return canonical
vtu = [f for f in os.listdir(d) if f.endswith(".vtu")]
if vtu:
return os.path.join(d, vtu[0])
return None
def extract_tar_if_needed(case_dir):
tar = os.path.join(case_dir, "surfaces.tar.gz")
out = os.path.join(case_dir, "surfaces")
if os.path.exists(tar) and not os.path.exists(out):
with tarfile.open(tar, "r:gz") as t:
t.extractall(case_dir, filter="tar")
# Download the transition (SA-AFT) case at alpha = 5° and extract its Cf distribution.
SURFACE_ALPHA = 5.0
surface_data_dir = os.path.join("results", "surface_data")
os.makedirs(surface_data_dir, exist_ok=True)
df_cf = None
for case, case_name, has_transition in cases:
if not has_transition or abs(case_alpha(case) - SURFACE_ALPHA) >= 0.01:
continue
dest = os.path.join(surface_data_dir, f"case_{case.id}")
print(f" {case_name} (alpha={SURFACE_ALPHA:.1f}deg)...")
case.results.download(surface=True, destination=dest)
extract_tar_if_needed(dest)
vtu = find_vtu_file(dest)
if vtu:
x_cf, cf_vals = extract_surface_data(vtu, "Cf")
if x_cf is not None and len(x_cf) > 0:
df_cf = pd.DataFrame({"x": x_cf, "cf": cf_vals}).sort_values("x").reset_index(drop=True)
break
deg = int(SURFACE_ALPHA)
fig, ax = plt.subplots(figsize=(10, 6))
if df_cf is not None:
ax.scatter(df_cf["x"], df_cf["cf"], s=4, color=FLOW360_COLOR, label="Flow360 (SA-AFT)", zorder=2)
for ref_name, ref_file, color, marker in [
("OVERFLOW (kw-SST γ-Reθ)", f"CF_{deg}deg_OVERFLOW.csv", OVERFLOW_COLOR, "s"),
("FUN3D (kw-SST-γ-Reθ)", f"CF_{deg}deg_FUN3D.csv", "gray", "^"),
]:
ref_path = os.path.join(ref_data_dir, ref_file)
if os.path.exists(ref_path):
df_ref = pd.read_csv(ref_path, skipinitialspace=True)
df_ref.columns = df_ref.columns.str.strip()
xc = next((c for c in df_ref.columns if c.lower() in ["x", "x/c"]), None)
fc = next((c for c in df_ref.columns if c.upper() == "CF"), None)
if xc and fc:
ax.scatter(df_ref[xc], df_ref[fc], s=6, marker=marker, color=color, label=ref_name, alpha=0.8, zorder=2)
ax.set(xlabel="x/c", ylabel="Cf", xlim=(0, 1))
ax.legend(loc="best", fontsize=10, markerscale=3)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f"results/Cf_vs_x_{deg}deg.png", dpi=150)
plt.show()