ONERA M6 Wing: Transonic Grid Convergence Study¶

The ONERA M6 wing is a canonical transonic CFD validation case: a 30-degree swept, semi-span wing with a symmetric ONERA D airfoil section that develops a characteristic lambda shock on the upper surface. Conditions are Mach 0.84, Reynolds number 11.72 million (based on mean aerodynamic chord) and 3.06 degrees angle of attack, solved with the Spalart-Allmaras turbulence model.

This notebook runs a grid-convergence study: the same case is solved on four mesh-refinement levels (coarse, medium, fine, extra-fine) taken from the NASA Turbulence Modeling Resource ONERA M6 grid family. Each level is its own volume-mesh project. We then assess grid convergence of the integrated forces and compare surface pressure and skin-friction distributions against reference CFD (FUN3D) and experiment.

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 vtk

Run the cells top to bottom in a single kernel.

Imports¶

Standard scientific stack plus flow360. download_benchmark_assets fetches the public benchmark snapshots (mesh root assets and reference data). vtk reads the surface-slice output files.

In [ ]:
import os
import tarfile

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import vtk

import flow360 as fl

# Shared plot style (Flow360 = brand green, reference CFD = grey, experiment = black).
FLOW360_COLOR = "#00643c"
FUN3D_COLOR = "#666666"
REFERENCE_COLOR = "black"

# Reference-solver visual style: same colour, different markers / dashes.
_REF_MARKERS = ["o", "s", "^", "v", "D", "P"]
_REF_LINESTYLES = ["--", "-.", ":", (0, (3, 1, 1, 1)), "--", "-."]

from flow360.examples import download_benchmark_assets

Input data¶

Flow conditions, reference geometry and the per-level mesh projects. The MESH_LEVELS values are public benchmark snapshot ids (not private cloud project handles): each names a saved root-asset snapshot of one refinement level. Downloading the reference data recreates a local ./ref_data/ directory that the post-processing cells read from (grid-convergence CSV, ONERA experiment, FUN3D distributions).

In [ ]:
# ── Flow conditions ──────────────────────────────────────────────────────────
MACH = 0.84
RE = 11.72e6
ALPHA = 3.06
T_REF = 297.78
MAC = 0.80167
SEMI_SPAN = 1.47602
REF_AREA = 1.15315

# Boundary names in the volume meshes.
WING = "WING"
FARFIELD = "FARFIELD"
SYMMETRY = "SYMMETRY"

# Span stations (y/b) at which Cp / Cf are extracted.
Y_OVER_B = [0.20, 0.44, 0.65, 0.80, 0.90, 0.96, 0.99]

# Public benchmark snapshot ids, one per mesh-refinement level.
MESH_LEVELS = {
    "coarse": "prj-f7e7c5cb-ca38-43b6-9335-4ee9bcf414fa",
    "medium": "prj-8363174b-1f1c-4c42-bdae-705a5fe42803",
    "fine": "prj-5f9d0e44-17e3-4eb0-a3a4-17eda0a7f974",
    "extra_fine": "prj-10e1160f-7b2e-4b92-a103-70e1fe5b1894",
}

SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")

# Reference data -> ./ref_data/
download_benchmark_assets("Onera_M6", "ref_data")

GRID_CONV_CSV = "ref_data/combined_forces_pitchmom_maxmut.csv"
ONERA_EXP_CSV = "ref_data/ONERAb114.csv"
FUN3D_CSV = "ref_data/FUN3D_OM6_A3p06_CPCF.csv"

Load projects (per-level meshes)¶

Each refinement level has its own volume-mesh root asset. Download the snapshot for each level and build a Flow360 project from its volume mesh.

In [ ]:
projects = {}
for level_name, project_id in MESH_LEVELS.items():
    files = download_benchmark_assets("Onera_M6", "root_assets", project_id)
    # Select the volume mesh file (.ugrid here); its .mapbc sidecar is
    # downloaded alongside and picked up automatically.
    volume_mesh_file = next(
        f for f in files if f.removesuffix(".zst").endswith((".cgns", ".ugrid"))
    )
    projects[level_name] = fl.Project.from_volume_mesh(volume_mesh_file)
    print(f"Loaded {level_name} mesh project.")

Physics / Outputs setup¶

Builds the simulation parameters for one mesh level: the flow solver (Navier-Stokes + Spalart-Allmaras), boundary conditions, operating condition, and the requested outputs (surface Cp/Cf, span-station slices, and a CL force monitor that drives an automatic stopping criterion). The parameters are identical across all four levels.

In [ ]:
def build_case_params(vm):
    """Build the SimulationParams for one mesh level's volume mesh `vm`."""
    y_slices = [round(yb * SEMI_SPAN, 5) for yb in Y_OVER_B]

    with fl.SI_unit_system:
        wall = fl.Wall(surfaces=vm[WING], name="wing")
        wall_out = fl.ForceOutput(
            output_fields=["CL"],
            name="cl_monitor",
            models=[wall],
        )
        run_control = fl.RunControl(
            stopping_criteria=[
                fl.StoppingCriterion(
                    monitor_field="CL",
                    monitor_output=wall_out,
                    tolerance=0.001,
                    tolerance_window_size=500,
                )
            ]
        )
        operating_condition = fl.AerospaceCondition.from_mach_reynolds(
            mach=MACH,
            reynolds_mesh_unit=RE,
            temperature=T_REF * fl.u.K,
            alpha=ALPHA * fl.u.deg,
            beta=0.0 * fl.u.deg,
            project_length_unit=MAC * fl.u.m,
        )
        reference_geometry = fl.ReferenceGeometry(
            area=REF_AREA,
            moment_center=(0.0, 0.0, 0.0),
            moment_length=(MAC, MAC, MAC),
        )
        time_stepping = fl.Steady(max_steps=10000)

        models = [
            wall,
            fl.Fluid(
                navier_stokes_solver=fl.NavierStokesSolver(
                    limit_pressure_density=True,
                    absolute_tolerance=1e-4,
                ),
                turbulence_model_solver=fl.SpalartAllmaras(absolute_tolerance=1e-4),
            ),
            fl.Freestream(surfaces=vm[FARFIELD], name="farfield"),
            fl.SlipWall(surfaces=vm[SYMMETRY], name="symmetry"),
        ]

        outputs = [
            fl.SurfaceOutput(
                surfaces=vm[WING],
                output_fields=["Cp", "Cf"],
            ),
            fl.SurfaceSliceOutput(
                name="wing_sections",
                entities=[
                    fl.Slice(
                        name=f"section_y{yb:.2f}b",
                        origin=(0.0, y_coord, 0.0),
                        normal=(0, 1, 0),
                    )
                    for yb, y_coord in zip(Y_OVER_B, y_slices)
                ],
                target_surfaces=vm[WING],
                output_fields=["Cp", "Cf"],
            ),
            fl.SliceOutput(
                slices=[
                    fl.Slice(
                        origin=[0, 0, 0] * fl.u.m,
                        normal=[0, 1, 0],
                        name="Yslice",
                    )
                ],
                output_fields=["Cp"],
            ),
            wall_out,
        ]

        params = fl.SimulationParams(
            operating_condition=operating_condition,
            reference_geometry=reference_geometry,
            time_stepping=time_stepping,
            models=models,
            outputs=outputs,
            run_control=run_control,
        )
    return params

Submit cases (the level sweep)¶

Submit one case per mesh level. Every case uses the same parameters and solver version; only the underlying mesh changes.

In [ ]:
cases = {}
for level_name, project in projects.items():
    vm = project.volume_mesh
    params = build_case_params(vm)
    case = project.run_case(
        params=params,
        name=f"oneraM6_gridConvergence_{level_name}_{SOLVER_VERSION}",
        solver_version=SOLVER_VERSION,
    )
    cases[level_name] = case
    print(f"Submitted: {level_name:12s} -> {case.id}")

Wait for completion¶

Block until every case has finished before post-processing its results.

In [ ]:
for level_name, case in cases.items():
    print(f"Waiting for {level_name} ...")
    case.wait()
print("All cases complete.")

Postprocessing¶

Reproduce the four published figures: CL and CD grid convergence, and Cp / Cf distributions at seven span stations. Results are read from the in-kernel case objects submitted above, no projects are re-opened.

First, extract the converged force coefficients (and node count) for each level.

In [ ]:
rows = []
for level_name, case in cases.items():
    avg = case.results.total_forces.averages
    n = case.volume_mesh.stats.n_nodes
    rows.append({"level": level_name, "N": n, "CL": avg["CL"], "CD": avg["CD"]})
    print(f"  {level_name:12s}  N={n:>12,}  CL={avg['CL']:.5f}  CD={avg['CD']:.5f}")

df_fl = pd.DataFrame(rows)

Grid-convergence plot helper¶

Load the reference grid-convergence data and define a helper that plots a force coefficient versus node count, overlaying the reference CFD solvers and the Flow360 result.

In [ ]:
df_ref = pd.read_csv(GRID_CONV_CSV).rename(columns={"C_L": "CL", "C_D": "CD"})

# Expand the cryptic solver tokens in the reference CSV for a readable legend.
solver_label = {
    "USM3D SA-neg, FO turb, Prism_Hex": "USM3D SA-neg (1st-order turb.), Prism-Hex",
    "USM3D SA-neg, FO turb, Tetrahedral": "USM3D SA-neg (1st-order turb.), Tet",
}


def grid_convergence_plot(metric, ylabel, ylim, savepath):
    fig, ax = plt.subplots(figsize=(7, 5))
    ax.set_title(
        f"ONERA M6 Wing, Grid Convergence, {ylabel}\n"
        r"Mach 0.84, $\alpha$ = 3.06°",
        fontsize=11,
    )

    # Reference solvers.
    for i, (solver, grp) in enumerate(df_ref.groupby("solver")):
        g = grp.sort_values("N")
        ax.plot(
            g["N"], g[metric],
            linestyle=_REF_LINESTYLES[i % len(_REF_LINESTYLES)],
            marker=_REF_MARKERS[i % len(_REF_MARKERS)],
            markersize=7, linewidth=1.3, color=REFERENCE_COLOR, alpha=0.7,
            label=solver_label.get(solver, solver),
        )

    # Flow360.
    ds = df_fl.sort_values("N")
    ax.plot(
        ds["N"], ds[metric],
        linestyle="-", marker="D", markersize=9, linewidth=1.8,
        color=FLOW360_COLOR, label="Flow360", zorder=5,
    )

    ax.set_xscale("log")
    ax.set_xlim(20_000, 120_000_000)
    if ylim is not None:
        ax.set_ylim(*ylim)
    ax.set_xlabel("N (nodes)", fontsize=11)
    ax.set_ylabel(ylabel, fontsize=11)
    ax.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.5)
    ax.legend(fontsize=8, loc="best", framealpha=0.9)
    plt.tight_layout()

    # Write the published artifact to its exact path, then show it inline.
    os.makedirs("results", exist_ok=True)
    fig.savefig(savepath, dpi=150, bbox_inches="tight")
    plt.show()

CL Grid Convergence¶

In [ ]:
grid_convergence_plot("CL", "$C_L$", (0.22, 0.30), "results/grid_convergence_CL.png")

CD Grid Convergence¶

In [ ]:
grid_convergence_plot("CD", "$C_D$", (0.015, 0.025), "results/grid_convergence_CD.png")

Surface-slice extraction¶

Helpers to download and read the surface-slice outputs (.pvtu) of the finest mesh, plus the wing-geometry relations used to normalise the chordwise coordinate at each span station.

In [ ]:
# ── Wing geometry (for chordwise normalisation) ──────────────────────────────
TAPER_RATIO = 0.562
LE_SWEEP = np.radians(30.0)
TE_SWEEP = np.radians(15.8)
ROOT_CHORD = SEMI_SPAN * (np.tan(TE_SWEEP) - np.tan(LE_SWEEP)) / (TAPER_RATIO - 1)


def local_le_x(y_abs):
    return y_abs * np.tan(LE_SWEEP)


def local_chord(y_abs):
    return ROOT_CHORD + y_abs * (np.tan(TE_SWEEP) - np.tan(LE_SWEEP))


def read_pvtu(pvtu_path, field):
    """Read a .pvtu slice. Returns (points, values) ordered as a closed airfoil
    curve: upper surface by x ascending, then lower surface by x descending."""
    reader = vtk.vtkXMLPUnstructuredGridReader()
    reader.SetFileName(pvtu_path)
    reader.Update()
    output = reader.GetOutput()
    n = output.GetNumberOfPoints()
    pts_vtk = output.GetPoints()
    pts = np.array([pts_vtk.GetPoint(i) for i in range(n)], dtype=np.float64)
    arr = output.GetPointData().GetArray(field)
    if arr is None:
        raise ValueError(f"Field '{field}' not found in {pvtu_path}")
    if arr.GetNumberOfComponents() == 1:
        values = np.array([arr.GetValue(i) for i in range(n)], dtype=np.float64)
    else:
        values = np.array([arr.GetTuple(i) for i in range(n)], dtype=np.float64)

    z = pts[:, 2]
    upper = z >= 0
    idx_upper = np.where(upper)[0][np.argsort(pts[upper, 0])]
    idx_lower = np.where(~upper)[0][np.argsort(pts[~upper, 0])[::-1]]
    order = np.concatenate([idx_upper, idx_lower])
    return pts[order], values[order]


def download_and_extract_surfaces(case, level):
    surface_dir = os.path.join("results", f"surfaces_{level}")
    os.makedirs(surface_dir, exist_ok=True)
    tgz_path = os.path.join(surface_dir, "surfaces.tar.gz")
    if not os.path.exists(tgz_path):
        case.results.surfaces.download(to_file=tgz_path)
    marker = os.path.join(surface_dir, ".extracted")
    if not os.path.exists(marker):
        with tarfile.open(tgz_path) as tar:
            tar.extractall(surface_dir)
        open(marker, "w").close()
    return surface_dir


def read_surface_slices(surface_dir, field="Cp"):
    """Read slice PVTUs -> {section_idx: DataFrame(x_norm, <field>)}."""
    data = {}
    for idx, yb in enumerate(Y_OVER_B):
        fpath = os.path.join(surface_dir, f"surface_slice_section_y{yb:.2f}b.pvtu")
        if not os.path.exists(fpath):
            continue
        pts, values = read_pvtu(fpath, field)
        y_abs = yb * SEMI_SPAN
        x_norm = (pts[:, 0] - local_le_x(y_abs)) / local_chord(y_abs)
        mask = (x_norm >= -0.02) & (x_norm <= 1.02)
        data[idx] = pd.DataFrame({"x_norm": x_norm[mask], field: values[mask]})
    return data


# Use the finest mesh that was run for the slice comparisons.
LEVEL_ORDER = ["extra_fine", "fine", "medium", "coarse"]
slice_level = next(lvl for lvl in LEVEL_ORDER if lvl in cases)
surface_dir = download_and_extract_surfaces(cases[slice_level], slice_level)
cp_slices = read_surface_slices(surface_dir, field="Cp")
print(f"Using {slice_level} mesh for span-station slices.")

Cp at Span Stations¶

Chordwise surface-pressure distribution at each span station: Flow360 (finest mesh) versus the ONERA experiment and FUN3D.

In [ ]:
df_exp = pd.read_csv(ONERA_EXP_CSV)
df_fun3d = pd.read_csv(FUN3D_CSV)
df_fun3d_l1 = df_fun3d[df_fun3d["grid_level"] == "L1"]

ncols = 4
nrows = (len(Y_OVER_B) + ncols - 1) // ncols
fig, axes = plt.subplots(nrows, ncols, figsize=(16, nrows * 4.5))
axes = axes.flatten()

for idx, yb in enumerate(Y_OVER_B):
    ax = axes[idx]

    exp_sec = df_exp[df_exp["Y/b"] == yb]
    if not exp_sec.empty:
        ax.scatter(
            exp_sec["X/L"], exp_sec["Cp"], s=20, color=REFERENCE_COLOR,
            marker="o", zorder=4, label="Experiment (ONERA)",
        )

    fun3d_sec = df_fun3d_l1[np.isclose(df_fun3d_l1["eta"].astype(float), yb, atol=0.015)]
    if not fun3d_sec.empty:
        ax.plot(
            fun3d_sec["X_over_C"].astype(float), fun3d_sec["CP"].astype(float),
            linewidth=1.0, color=FUN3D_COLOR, alpha=0.6, zorder=3, label="FUN3D",
        )

    if idx in cp_slices:
        d = cp_slices[idx]
        ax.plot(
            d["x_norm"], d["Cp"], linewidth=1.2, color=FLOW360_COLOR,
            alpha=0.8, zorder=5, label="Flow360",
        )

    ax.set_title(f"Section {idx + 1}, y/b = {yb}", fontsize=10)
    ax.set_xlabel("x/c", fontsize=9)
    ax.set_ylabel("$C_p$", fontsize=9)
    ax.invert_yaxis()
    ax.set_xlim(-0.02, 1.02)
    ax.grid(True, linestyle="--", alpha=0.3)
    ax.legend(fontsize=7)

for j in range(len(Y_OVER_B), len(axes)):
    axes[j].set_visible(False)

plt.tight_layout()

# Write the published artifact to its exact path, then show it inline.
os.makedirs("results", exist_ok=True)
fig.savefig("results/cp_slices_comparison.png", dpi=150, bbox_inches="tight")
plt.show()

Cf at Span Stations¶

Chordwise skin-friction magnitude at each span station: Flow360 (finest mesh) versus FUN3D (no experimental data available for Cf).

In [ ]:
df_fun3d = pd.read_csv(FUN3D_CSV)
df_fun3d_l1 = df_fun3d[df_fun3d["grid_level"] == "L1"]

ncols = 4
nrows = (len(Y_OVER_B) + ncols - 1) // ncols
fig, axes = plt.subplots(nrows, ncols, figsize=(16, nrows * 4.5))
axes = axes.flatten()

for idx, yb in enumerate(Y_OVER_B):
    ax = axes[idx]

    fun3d_sec = df_fun3d_l1[np.isclose(df_fun3d_l1["eta"].astype(float), yb, atol=0.015)]
    if not fun3d_sec.empty:
        cfx = fun3d_sec["CFX"].astype(float)
        cfy = fun3d_sec["CFY"].astype(float)
        cfz = fun3d_sec["CFZ"].astype(float)
        cf_mag = np.sqrt(cfx**2 + cfy**2 + cfz**2)
        ax.plot(
            fun3d_sec["X_over_C"].astype(float), cf_mag,
            linewidth=1.0, color=FUN3D_COLOR, alpha=0.6, zorder=3, label="FUN3D",
        )

    fpath = os.path.join(surface_dir, f"surface_slice_section_y{yb:.2f}b.pvtu")
    if os.path.exists(fpath):
        pts, cf_vec = read_pvtu(fpath, "Cf")
        y_abs = yb * SEMI_SPAN
        x_norm = (pts[:, 0] - local_le_x(y_abs)) / local_chord(y_abs)
        mask = (x_norm >= -0.02) & (x_norm <= 1.02)
        if cf_vec.ndim == 2:
            cf_magnitude = np.linalg.norm(cf_vec[mask], axis=1)
        else:
            cf_magnitude = np.abs(cf_vec[mask])
        ax.plot(
            x_norm[mask], cf_magnitude, linewidth=1.2, color=FLOW360_COLOR,
            alpha=0.8, zorder=5, label="Flow360",
        )

    ax.set_title(f"Section {idx + 1}, y/b = {yb}", fontsize=10)
    ax.set_xlabel("x/c", fontsize=9)
    ax.set_ylabel("$|C_f|$", fontsize=9)
    ax.set_xlim(-0.02, 1.02)
    ax.grid(True, linestyle="--", alpha=0.3)
    ax.legend(fontsize=7)

for j in range(len(Y_OVER_B), len(axes)):
    axes[j].set_visible(False)

plt.tight_layout()

# Write the published artifact to its exact path, then show it inline.
os.makedirs("results", exist_ok=True)
fig.savefig("results/cf_slices_comparison.png", dpi=150, bbox_inches="tight")
plt.show()