Imperial Front Wing: RANS¶
The Imperial Front Wing is a publicly documented F1-style multi-element front-wing geometry tested over a moving ground, with strongly three-dimensional flow dominated by the endplate and underwing vortices in the near wake. This case is a single steady RANS solution with the k-omega SST turbulence model on a Flow360 GAI (Geometry AI) automated mesh, at a freestream velocity of 25 m/s over a moving ground.
Flow360's surface flow and downstream wake velocity fields are compared against experimental flow-vis and PIV data and against an OpenFOAM reference. The notebook runs top-to-bottom in a single kernel: load the geometry, submit the case, wait for completion, then post-process the in-kernel case object into the eight published comparison figures (wing/endplate surface flow-vis and the Vy/Vz wake-plane slices) plus the Flow360-vs-OpenFOAM integrated wing-load table.
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, and rendering
packages:
pip install numpy pandas matplotlib vtk pillow
The surface flow-vis renders off-screen with VTK's surface-LIC mapper
(vtkRenderingLICOpenGL2, shipped with the standard VTK 9 wheels); on a
headless machine it needs a working off-screen OpenGL backend. Run the cells top
to bottom in a single kernel.
Imports¶
In [ ]:
import json
import os
import re
import shutil
import tarfile
from dataclasses import dataclass
from pathlib import Path
from typing import Optional, Tuple
from xml.etree import ElementTree as ET
import flow360 as fl
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import vtk
from matplotlib.cm import ScalarMappable
from matplotlib.colors import ListedColormap, Normalize, hsv_to_rgb
from matplotlib.ticker import MaxNLocator
from PIL import Image
from vtk.util import numpy_support as vtk_np
from flow360.examples import download_benchmark_assets
Input data¶
The post-processing compares Flow360 against experimental flow-vis and PIV
images and an OpenFOAM reference. Fetch those reference files from the public
benchmark bucket; this recreates a local ./ref_data/ directory (including the
openfoam/ sub-tree) that the plotting cells read from.
In [ ]:
download_benchmark_assets("Imperial_Front_Wing", "ref_data")
Case constants¶
The solver version, the freestream test condition, the case name, and the
downstream wake-plane X positions (metres) matching the experimental PIV planes
1, 3 and 5. The slice index is what the post-processing looks up in
ref_data/extents.json and the experimental images.
In [ ]:
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.10")
VELOCITY = 25.0 # freestream velocity [m/s], moving-ground front-wing condition
CASE_NAME = "ImperialFrontWing_SST_GAI"
WAKE_SLICE_X_M = {1: -0.294, 3: -0.150, 5: 0.150}
Load project¶
The root asset for this case is the uploaded geometry (wing + domain). Download
it from the public benchmark bucket and start a fresh project from it, no
private project id is needed. The faces are grouped by their faceId tag and the
bodies by file, so the setup can address surfaces by name pattern (e.g.
fw-mp*, wall-z-min-nearCar).
In [ ]:
root_asset_files = download_benchmark_assets("Imperial_Front_Wing", "root_assets")
# The snapshot includes processed copies under results/; select the source
# CAD file(s) to rebuild the geometry project from.
geometry_files = [
f for f in root_asset_files
if "/results/" not in f and "/logs/" not in f
and f.lower().endswith((".csm", ".egads", ".stp", ".step", ".iges", ".igs", ".stl"))
]
project = fl.Project.from_geometry(geometry_files, name="Imperial Front Wing")
geo = project.geometry
geo.group_faces_by_tag("faceId")
geo.group_bodies_by_tag("groupByFile")
Meshing setup¶
The mesh is built by the GAI (Geometry-AI) automated mesher, translated verbatim
from the exported GAI 4.json setup: global defaults, a user-defined farfield
plus a seeded custom near-body zone, per-element geometry/surface refinements, a
near-car refinement box, and boundary layers on the ground and wing surfaces.
In [ ]:
def build_meshing(geo):
with fl.SI_unit_system:
return fl.MeshingParams(
defaults=fl.MeshingDefaults(
surface_max_edge_length=500 * fl.u.mm,
curvature_resolution_angle=12 * fl.u.deg,
target_surface_node_count=1300000,
boundary_layer_first_layer_thickness=0.5 * fl.u.mm,
octree_spacing=fl.OctreeSpacing(base_spacing=1 * fl.u.mm),
),
volume_zones=[
fl.UserDefinedFarfield(),
fl.CustomZones(
name="Custom zone",
element_type="mixed",
entities=[
fl.SeedpointVolume(
name="Seed point volume",
point_in_mesh=[[0, 0, 2] * fl.u.m],
)
],
),
],
refinements=[
# Wing geometry accuracy (all FW elements except support-cv1).
fl.GeometryRefinement(
name="wing_ref",
geometry_accuracy=0.5 * fl.u.mm,
entities=[
geo["fw-ep*"],
geo["fw-flaps*"],
geo["fw-gf"],
geo["fw-mp*"],
geo["fw-nose*"],
geo["fw-support"],
],
),
# Near-car volume box.
fl.UniformRefinement(
name="nearcar",
spacing=4 * fl.u.mm,
entities=[
fl.Box(
name="nearcar",
center=[-0.15, 0, -0.2] * fl.u.m,
size=[1.75, 1.6, 1.2] * fl.u.m,
),
],
),
# Moving-ground surface refinement.
fl.SurfaceRefinement(
name="Surface refinement",
max_edge_length=28 * fl.u.mm,
entities=[geo["wall-z-min-nearCar"]],
),
# Front-wing surface refinement (all FW elements).
fl.SurfaceRefinement(
name="Surface refinement",
max_edge_length=7 * fl.u.mm,
curvature_resolution_angle=9 * fl.u.deg,
entities=[
geo["fw-nose*"],
geo["fw-support*"],
geo["fw-ep*"],
geo["fw-mp*"],
geo["fw-flaps*"],
geo["fw-gf"],
],
),
# Support-cv1 geometry accuracy.
fl.GeometryRefinement(
name="Geometry refinement",
geometry_accuracy=0.2 * fl.u.mm,
entities=[geo["fw-support-cv1"]],
),
# Passive spacing on the domain walls (unchanged).
fl.PassiveSpacing(
name="Passive spacing",
type="unchanged",
entities=[
geo["wall-z-min"],
geo["wall-z-max"],
geo["wall-y-min"],
geo["wall-y-max"],
geo["wall-x-min"],
geo["wall-x-max"],
],
),
# Boundary layer: moving ground.
fl.BoundaryLayer(
name="Boundary layer refinement",
first_layer_thickness=0.8 * fl.u.mm,
growth_rate=1.6,
entities=[geo["wall-z-min-nearCar"]],
),
# Boundary layer: main wing surfaces.
fl.BoundaryLayer(
name="Boundary layer refinement",
first_layer_thickness=0.05 * fl.u.mm,
growth_rate=1.2,
entities=[
geo["fw-mp*"],
geo["fw-gf"],
geo["fw-flaps*"],
geo["fw-ep*"],
],
),
# Boundary layer: nose + support.
fl.BoundaryLayer(
name="Boundary layer refinement",
first_layer_thickness=0.25 * fl.u.mm,
growth_rate=1.2,
entities=[geo["fw-support*"], geo["fw-nose*"]],
),
],
)
Physics setup¶
Boundary conditions and the fluid model. The six per-element front-wing walls
(passed in as walls, shared with the force outputs and convergence monitors)
use wall functions; the near-car ground is a moving wall at the freestream
speed; the side/top/bottom domain walls are slip walls and the streamwise walls
are freestream. The fluid is a steady k-omega SST solve with a low-Mach-
preconditioned Roe flux and rotation correction.
In [ ]:
def build_models(walls, geo):
fw_mp, fw_ep, fw_gf, fw_flaps, fw_support, fw_nose = walls
with fl.SI_unit_system:
return [
# Front-wing wall BCs.
fw_mp,
fw_ep,
fw_gf,
fw_flaps,
fw_support,
fw_nose,
# Moving ground.
fl.Wall(
surfaces=[geo["wall-z-min-nearCar"]],
use_wall_function=fl.WallFunction(),
velocity=[VELOCITY, 0, 0] * fl.u.m / fl.u.s,
),
# Domain BCs.
fl.SlipWall(
surfaces=[
geo["wall-y-min"],
geo["wall-y-max"],
geo["wall-z-max"],
geo["wall-z-min"],
],
),
fl.Freestream(surfaces=[geo["wall-x-min"], geo["wall-x-max"]]),
# Fluid model: steady k-omega SST.
fl.Fluid(
navier_stokes_solver=fl.NavierStokesSolver(
absolute_tolerance=1e-4,
riemann_solver=fl.RoeFlux(low_mach_preconditioner=True),
),
turbulence_model_solver=fl.KOmegaSST(
absolute_tolerance=1e-4,
rotation_correction=True,
),
),
]
Outputs setup¶
The solver outputs feeding the post-processing and convergence monitoring: a single-file surface output (used for the surface-LIC flow-vis), the X-normal wake slices carrying a free-stream-normalized velocity variable and total pressure, a symmetry-plane slice, a vorticity-tracking slice sweep, a Q-criterion isosurface, per-element and combined wing force outputs, and the CL/CD convergence monitors that also drive the stopping criteria.
In [ ]:
def build_outputs(geo, walls, all_fw_models, cd_convergence, cl_convergence):
fw_mp, fw_ep, fw_gf, fw_flaps, fw_support, fw_nose = walls
with fl.SI_unit_system:
# Custom variable: velocity normalized by the freestream speed.
vel_adim = fl.UserVariable(
name="velocity_adimensionalized",
value=[
fl.solution.velocity[0] / (VELOCITY * fl.u.m / fl.u.s),
fl.solution.velocity[1] / (VELOCITY * fl.u.m / fl.u.s),
fl.solution.velocity[2] / (VELOCITY * fl.u.m / fl.u.s),
],
)
# All patches except the domain walls -- the front-wing surfaces only.
wing_surfaces = [s for s in geo["*"] if not s.name.startswith("wall")]
return [
# Surface output (surface-LIC flow-vis).
fl.SurfaceOutput(
surfaces=wing_surfaces,
output_format=["paraview"],
output_fields=["primitiveVars", "Cp", "yPlus", "CfVec"],
write_single_file=True,
),
# Wake plane slices (X-normal) matching the experimental planes.
fl.SliceOutput(
name="wake_slices",
entities=[
fl.Slice(name=f"slice_x_{n}", normal=(1, 0, 0), origin=(x, 0, 0))
for n, x in WAKE_SLICE_X_M.items()
],
output_format=["paraview"],
output_fields=["Cp", vel_adim, fl.solution.Cpt_auto],
),
# Symmetry-plane slice.
fl.SliceOutput(
name="symmetry",
entities=[
fl.Slice(name="Y_symmetry", normal=(0, 1, 0), origin=(0, 0, 0)),
],
output_format=["paraview"],
output_fields=["Cp", vel_adim, fl.solution.Cpt_auto],
),
# Vorticity-tracking slice sweep.
fl.SliceOutput(
name="vorticity_tracking_slices",
entities=[
fl.Slice(
name=f"slice_vor_x_{x}", normal=(1, 0, 0), origin=(x, 0, 0)
)
for x in np.linspace(-0.9, 0.3, 9)
],
output_format=["paraview"],
output_fields=[
"vorticityMagnitude",
"vorticity_x",
"vorticity_y",
"vorticity_z",
],
),
# Q-criterion isosurface (vortex cores).
fl.IsosurfaceOutput(
name="qcriterion_isosurface",
isosurfaces=[
fl.Isosurface(
name="qcriterion_better",
field="qcriterion",
iso_value=1e-6,
),
],
output_format=["paraview"],
output_fields=["Cp", vel_adim],
),
# Per-component force outputs.
fl.ForceOutput(
name="wing_main_plane_forces",
models=[fw_mp],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
fl.ForceOutput(
name="wing_endplate_forces",
models=[fw_ep],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
fl.ForceOutput(
name="wing_gurney_flap_forces",
models=[fw_gf],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
fl.ForceOutput(
name="wing_flaps_forces",
models=[fw_flaps],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
fl.ForceOutput(
name="wing_support_forces",
models=[fw_support],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
fl.ForceOutput(
name="wing_nose_forces",
models=[fw_nose],
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
# Combined wing forces.
fl.ForceOutput(
name="wing_all",
models=all_fw_models,
output_fields=["CFx", "CFy", "CFz", "CMx", "CMy", "CMz"],
),
# Convergence monitors.
cd_convergence,
cl_convergence,
]
Simulation Params¶
Assemble the full SimulationParams. The per-element front-wing walls and the
CL/CD convergence monitors are created here (they are shared across the models,
outputs and stopping criteria), then the meshing, models and outputs builders are
called. The run stops early once the integrated CD and CL have both settled
within tolerance over a 50-step window.
In [ ]:
def build_params(geo):
"""Assemble the SimulationParams for the SST GAI-mesher run."""
with fl.SI_unit_system:
operating_condition = fl.AerospaceCondition(
velocity_magnitude=VELOCITY * fl.u.m / fl.u.s,
alpha=0 * fl.u.deg,
beta=0 * fl.u.deg,
)
# Per-element front-wing walls, reused across models, force outputs and
# convergence monitors.
fw_mp = fl.Wall(surfaces=[geo["fw-mp*"]], use_wall_function=fl.WallFunction())
fw_ep = fl.Wall(surfaces=[geo["fw-ep*"]], use_wall_function=fl.WallFunction())
fw_gf = fl.Wall(surfaces=[geo["fw-gf"]], use_wall_function=fl.WallFunction())
fw_flaps = fl.Wall(
surfaces=[geo["fw-flaps*"]], use_wall_function=fl.WallFunction()
)
fw_support = fl.Wall(
surfaces=[geo["fw-support*"]], use_wall_function=fl.WallFunction()
)
fw_nose = fl.Wall(
surfaces=[geo["fw-nose*"]], use_wall_function=fl.WallFunction()
)
walls = (fw_mp, fw_ep, fw_gf, fw_flaps, fw_support, fw_nose)
all_fw_models = [fw_gf, fw_mp, fw_ep, fw_flaps, fw_nose, fw_support]
cd_convergence = fl.ForceOutput(
name="CD_convergence", models=all_fw_models, output_fields=["CD"]
)
cl_convergence = fl.ForceOutput(
name="CL_convergence", models=all_fw_models, output_fields=["CL"]
)
return fl.SimulationParams(
meshing=build_meshing(geo),
operating_condition=operating_condition,
reference_geometry=fl.ReferenceGeometry(
moment_center=(1, 0, 0),
moment_length=(1, 1, 1),
area=1,
),
time_stepping=fl.Steady(max_steps=7000),
models=build_models(walls, geo),
outputs=build_outputs(
geo, walls, all_fw_models, cd_convergence, cl_convergence
),
run_control=fl.RunControl(
stopping_criteria=[
fl.StoppingCriterion(
monitor_field="CD",
monitor_output=cd_convergence,
tolerance=0.001,
tolerance_window_size=50,
),
fl.StoppingCriterion(
monitor_field="CL",
monitor_output=cl_convergence,
tolerance=0.01,
tolerance_window_size=50,
),
],
),
)
Submit cases¶
Submit the single steady SST case to Flow360 through the GAI beta mesher. The solver version is embedded in the case name so reruns across releases stay distinguishable. The length-scale-mismatch warning is bypassed because the source geometry length unit (mm) differs from the SI solve.
In [ ]:
case_name = f"{CASE_NAME}_{SOLVER_VERSION}"
with fl.warning_bypass("potential_length_scale_mismatch"):
case = project.run_case(
params=build_params(geo),
name=case_name,
use_beta_mesher=True,
use_geometry_AI=True,
solver_version=SOLVER_VERSION,
)
print(f"Submitted: {case_name} -> {case.id}")
Wait for completion¶
Block until the case has finished. This can take a long time (minutes to hours): it is a full steady solve of up to 7000 pseudo-steps on a ~52 M-node GAI mesh, ending early once the loads converge.
In [ ]:
case.wait()
Postprocessing¶
The published results are eight comparison figures and one table, all rebuilt
from the in-kernel case object (no project is re-opened):
- Surface flow-vis: off-screen VTK surface-LIC of the skin-friction field
coloured by
Cp, for the wing underside and the endplate, shown Experiment vs Flow360 vs OpenFOAM. - Wake-plane slices: free-stream-normalized
VyandVzon the X-normal wake planes 1, 3 and 5, shown Experiment vs Flow360 vs OpenFOAM. - Integrated wing loads: Flow360 vs OpenFOAM downforce and drag coefficients on the shared 1 m² reference area.
The OpenFOAM panels are drawn from ref_data/openfoam/; if that reference data
is absent, the OpenFOAM panel is simply skipped.
Setup¶
Render backend, off-screen surface-LIC mapper, reference paths, the per-slice crop windows and value ranges, the OpenFOAM reference loads, the experimental image pixel calibration, and the classic VTK rainbow colormap used for the slice contours.
In [ ]:
matplotlib.use("Agg")
try:
from vtkmodules.vtkRenderingLICOpenGL2 import vtkSurfaceLICMapper
except ImportError as exc:
raise ImportError(
"vtkSurfaceLICMapper is required. Install a VTK build that ships the "
"vtkRenderingLICOpenGL2 module (standard VTK 9 wheels include it)."
) from exc
REF_DIR = Path("ref_data")
EXTENTS_FILE = REF_DIR / "extents.json"
RESULTS_DIR = Path("results")
RESULTS_DIR.mkdir(parents=True, exist_ok=True)
SLICE_NUMBERS = (1, 3, 5)
COMPONENTS = ("Vy", "Vz")
COMPONENT_INDEX = {"Vy": 1, "Vz": 2}
# OpenFOAM reference wing loads (full wing), reported as force / (0.5 rho U^2) on
# the same 1 m^2 reference area as the Flow360 setup: CFz is downforce
# (negative-down), CFx is drag.
OPENFOAM_WING_LOADS = {"CFz": -1.955, "CFx": 0.208}
# Experimental images cover the LEFT half (y < 0). Crop windows per slice, in mm,
# as ((y_min, y_max), (z_min, z_max)).
SLICE_CROPS = {
1: ((-700.0, -400.0), (0.0, 250.0)),
3: ((-700.0, -400.0), (0.0, 250.0)),
5: ((-700.0, -250.0), (0.0, 250.0)),
}
# Lower Z cut applied to every slice panel (mm): trims the near-ground strip so
# all panels start at the same height.
Z_TRIM_MM = 5.0
OPENFOAM_DIR = REF_DIR / "openfoam"
_OF_SLICE_RE = re.compile(r"_slice_(\d+)__")
# Human-readable label for the field driving each render's colouring.
FIELD_LABELS = {"Cp": "$C_p$"}
# Pixel calibration of the experimental PIV images, per (component, slice):
# (pixels per 100 mm, x-pixel at Y=-500 mm, y-pixel at Z=0 mm). Used to place the
# image and its ticks in mm coordinates.
EXP_CALIB = {
("Vy", 1): (175, 325, 438),
("Vy", 3): (178, 331, 446),
("Vy", 5): (176, 327, 441),
("Vz", 1): (176, 327, 441),
("Vz", 3): (179, 332, 448),
("Vz", 5): (175, 326, 438),
}
def vtk_rainbow_colormap(n=256, name="vtk_rainbow"):
"""Classic VTK rainbow LUT (blue -> red) as a matplotlib colormap."""
hues = np.linspace(0.6667, 0.0, n)
hsv = np.column_stack([hues, np.ones(n), np.ones(n)])
return ListedColormap(hsv_to_rgb(hsv), name=name)
CMAP = vtk_rainbow_colormap()
Rendering cameras¶
Off-screen parallel-projection camera configs for the two surface-LIC views: the
wing underside and the endplate, both coloured by Cp over fixed value ranges.
In [ ]:
@dataclass
class CameraConfig:
position: Tuple[float, float, float]
focal_point: Tuple[float, float, float]
view_up: Tuple[float, float, float]
parallel_scale: float
clip_y_min: Optional[float] = None
clip_y_max: Optional[float] = None
rotate_degrees: int = 0
image_size: Tuple[int, int] = (1200, 1200)
color_array: str = "Cp"
scalar_min: Optional[float] = None
scalar_max: Optional[float] = None
# Geometry length unit is mm.
WING_UNDERSIDE_CAMERA = CameraConfig(
position=(-630.0, -350.0, -5000.0),
focal_point=(-630.0, -350.0, 230.0),
view_up=(-1.0, 0.0, 0.0),
parallel_scale=450.0,
clip_y_max=0.0,
rotate_degrees=-90,
scalar_min=-4.0,
scalar_max=0.5,
)
ENDPLATE_UNDERSIDE_CAMERA = CameraConfig(
position=(-590.0, 700.0, 5000.0),
focal_point=(-590.0, 700.0, 230.0),
view_up=(-1.0, 0.0, 0.0),
parallel_scale=400.0 / 1.2,
image_size=(750, 1200),
scalar_min=-0.5,
scalar_max=1.2,
)
Result readers¶
Download and extract the case's slice and surface archives, plus the VTK/EnSight readers that load the Flow360 slices and surfaces and the matching OpenFOAM reference data.
In [ ]:
def fetch_results(case, download_dir):
"""Download and extract the slices + surfaces tarballs."""
download_dir.mkdir(parents=True, exist_ok=True)
slice_dir = download_dir / "slices"
surface_dir = download_dir / "surfaces"
shutil.rmtree(slice_dir, ignore_errors=True)
archive = download_dir / "slices.tar.gz"
case.results.slices.download(str(archive), overwrite=True)
with tarfile.open(archive) as tf:
tf.extractall(slice_dir)
shutil.rmtree(surface_dir, ignore_errors=True)
archive = download_dir / "surfaces.tar.gz"
case.results.surfaces.download(str(archive), overwrite=True)
with tarfile.open(archive) as tf:
tf.extractall(surface_dir)
return slice_dir, surface_dir
def _resolve_pvd_to_vtu(pvd_path):
root = ET.parse(pvd_path).getroot()
datasets = root.findall(".//DataSet")
if not datasets:
raise RuntimeError(f"No DataSet entries in {pvd_path}")
file_attr = datasets[-1].get("file")
if not file_attr:
raise RuntimeError(f"DataSet in {pvd_path} has no 'file' attribute")
return (pvd_path.parent / file_attr).resolve()
def _find_slice_file(slice_dir, n):
"""Locate the data file for wake-plane slice slice_x_{n}."""
stem = f"slice_x_{n}"
for ext in (".pvtu", ".pvd", ".vtu", ".vtp"):
hits = sorted(p for p in slice_dir.rglob(f"*{ext}") if stem in p.name)
if hits:
return hits[0]
raise FileNotFoundError(
f"Could not find slice file for n={n} (stem '{stem}') under {slice_dir}"
)
def _read_dataset(path):
suffix = path.suffix.lower()
if suffix == ".vtu":
reader = vtk.vtkXMLUnstructuredGridReader()
elif suffix == ".vtp":
reader = vtk.vtkXMLPolyDataReader()
elif suffix == ".pvtu":
reader = vtk.vtkXMLPUnstructuredGridReader()
elif suffix == ".pvd":
return _read_dataset(_resolve_pvd_to_vtu(path))
else:
raise RuntimeError(f"Unsupported file type: {path}")
reader.SetFileName(str(path))
reader.Update()
return reader.GetOutput()
def _ensure_point_data(ds, array_name):
if ds.GetPointData().GetArray(array_name) is not None:
return ds
if ds.GetCellData().GetArray(array_name) is not None:
c2p = vtk.vtkCellDataToPointData()
c2p.SetInputData(ds)
c2p.PassCellDataOff()
c2p.Update()
return c2p.GetOutput()
available = [
ds.GetPointData().GetArrayName(i)
for i in range(ds.GetPointData().GetNumberOfArrays())
] + [
ds.GetCellData().GetArrayName(i)
for i in range(ds.GetCellData().GetNumberOfArrays())
]
raise RuntimeError(f"Array '{array_name}' not found. Available: {available}")
def _find_openfoam_slice(n):
"""EnSight case path for the OpenFOAM wake-exp slice numbered n, or None."""
if not OPENFOAM_DIR.is_dir():
return None
for d in OPENFOAM_DIR.iterdir():
if not d.is_dir() or "Wake_exp" not in d.name:
continue
m = _OF_SLICE_RE.search(d.name)
if m and int(m.group(1)) == n:
return d / f"{d.name}.case"
return None
def read_openfoam_slice(n):
"""Read the OpenFOAM wake-exp slice n, or None if unavailable.
Returns (Y_mm, Z_mm, Unorm) with coordinates converted to the mm slice frame.
"""
case_path = _find_openfoam_slice(n)
if case_path is None or not case_path.is_file():
return None
reader = vtk.vtkGenericEnSightReader()
reader.SetCaseFileName(str(case_path))
reader.ReadAllVariablesOn()
reader.Update()
blk = reader.GetOutput().GetBlock(0)
while blk is not None and blk.IsA("vtkMultiBlockDataSet"):
blk = blk.GetBlock(0)
if blk is None or blk.GetPoints() is None:
return None
pts = vtk_np.vtk_to_numpy(blk.GetPoints().GetData())
unorm = vtk_np.vtk_to_numpy(blk.GetPointData().GetArray("Unorm"))
return pts[:, 1] * 1000.0, pts[:, 2] * 1000.0, unorm
def _find_surface_file(surface_dir):
pvds = list(surface_dir.rglob("*.pvd"))
if pvds:
return _resolve_pvd_to_vtu(pvds[0])
for ext in ("*.vtp", "*.pvtu", "*.vtu"):
hits = list(surface_dir.rglob(ext))
if hits:
return hits[0]
raise FileNotFoundError(f"No surface file found under {surface_dir}")
def read_wing_surfaces(surface_dir):
ds = _read_dataset(_find_surface_file(surface_dir))
ds = _ensure_point_data(ds, "Cp")
ds = _ensure_point_data(ds, "CfVec")
if ds.IsA("vtkUnstructuredGrid"):
gf = vtk.vtkGeometryFilter()
gf.SetInputData(ds)
gf.Update()
ds = gf.GetOutput()
return ds
def read_openfoam_wing_surface():
"""OpenFOAM wing surface as polydata in the mm frame, or None.
Reads the EnSight surfaces/surfaces.case export directly (scaled m -> mm). The
LIC vector is UwNorm (near-wall velocity, the OpenFOAM analogue of Flow360's
CfVec); colouring uses the shared Cp field so the colormap matches Flow360.
"""
case = OPENFOAM_DIR / "surfaces" / "surfaces.case"
if not case.is_file():
return None
reader = vtk.vtkGenericEnSightReader()
reader.SetCaseFileName(str(case))
reader.ReadAllVariablesOn()
reader.Update()
blk = reader.GetOutput().GetBlock(0)
while blk is not None and blk.IsA("vtkMultiBlockDataSet"):
blk = blk.GetBlock(0)
if blk is None or blk.GetPoints() is None:
return None
gf = vtk.vtkGeometryFilter()
gf.SetInputData(blk)
gf.Update()
poly = gf.GetOutput()
pts = poly.GetPoints()
coords = vtk_np.vtk_to_numpy(pts.GetData()) * 1000.0 # metres -> mm
pts.SetData(vtk_np.numpy_to_vtk(coords, deep=1))
poly.Modified()
return poly
def _mirror_surface_y(poly):
"""Copy of poly mirrored about the symmetry plane y=0.
The OpenFOAM wing export covers only the left half (y <= 0); mirroring places
it on the right half to match the Flow360 endplate view. The y component of
UwNorm is flipped to stay physical and the triangle winding reversed.
"""
out = vtk.vtkPolyData()
out.DeepCopy(poly)
coords = vtk_np.vtk_to_numpy(out.GetPoints().GetData()).copy()
coords[:, 1] *= -1.0
out.GetPoints().SetData(vtk_np.numpy_to_vtk(coords, deep=1))
uw = out.GetPointData().GetArray("UwNorm")
if uw is not None:
vec = vtk_np.vtk_to_numpy(uw).copy()
vec[:, 1] *= -1.0
new = vtk_np.numpy_to_vtk(vec, deep=1)
new.SetName("UwNorm")
out.GetPointData().RemoveArray("UwNorm")
out.GetPointData().AddArray(new)
rev = vtk.vtkReverseSense()
rev.SetInputData(out)
rev.ReverseCellsOn()
rev.ReverseNormalsOff()
rev.Update()
return rev.GetOutput()
Rendering and composition helpers¶
The off-screen surface-LIC renderer and the figure-composition helpers: a side-by-side comparison composer for the surface renders (with a matched colorbar) and the 2×2 wake-slice figure that overlays the experimental image, a slice locator, and the Flow360 and OpenFOAM contours.
In [ ]:
def _clip_polydata_y(polydata, *, y_min=None, y_max=None):
out = polydata
for origin_y, normal_y in ((y_min, 1.0), (y_max, -1.0)):
if origin_y is None:
continue
plane = vtk.vtkPlane()
plane.SetOrigin(0.0, origin_y, 0.0)
plane.SetNormal(0.0, normal_y, 0.0)
clipper = vtk.vtkClipPolyData()
clipper.SetInputData(out)
clipper.SetClipFunction(plane)
clipper.InsideOutOff()
clipper.Update()
out = clipper.GetOutput()
return out
def _build_lut(vmin, vmax):
"""Diverging blue-white-red LUT for the Cp surface renders."""
lut = vtk.vtkColorTransferFunction()
lut.AddRGBPoint(vmin, 0.23, 0.30, 0.75)
lut.AddRGBPoint(0.5 * (vmin + vmax), 0.95, 0.95, 0.95)
lut.AddRGBPoint(vmax, 0.71, 0.02, 0.15)
return lut
def _lut_to_cmap(vmin, vmax, n=256):
"""Sample the VTK render LUT into a matplotlib colormap for the colorbar."""
lut = _build_lut(vmin, vmax)
colors = [lut.GetColor(v) for v in np.linspace(vmin, vmax, n)]
return ListedColormap(colors)
def render_surface_lic(polydata, camera, out_png, vector_array="CfVec"):
poly = polydata
if camera.clip_y_min is not None or camera.clip_y_max is not None:
poly = _clip_polydata_y(poly, y_min=camera.clip_y_min, y_max=camera.clip_y_max)
poly.GetPointData().SetActiveVectors(vector_array)
poly.GetPointData().SetActiveScalars(camera.color_array)
scalar_data = vtk_np.vtk_to_numpy(poly.GetPointData().GetArray(camera.color_array))
vmin = (
float(camera.scalar_min)
if camera.scalar_min is not None
else float(np.percentile(scalar_data, 2))
)
vmax = (
float(camera.scalar_max)
if camera.scalar_max is not None
else float(np.percentile(scalar_data, 98))
)
lut = _build_lut(vmin, vmax)
mapper = vtkSurfaceLICMapper()
mapper.SetInputData(poly)
mapper.SetLookupTable(lut)
mapper.SetScalarModeToUsePointData()
mapper.SetColorModeToMapScalars()
mapper.ScalarVisibilityOn()
mapper.SelectColorArray(camera.color_array)
lic = mapper.GetLICInterface()
lic.SetNumberOfSteps(40)
lic.SetStepSize(0.4)
lic.SetEnhancedLIC(1)
lic.SetEnhanceContrast(1)
actor = vtk.vtkActor()
actor.SetMapper(mapper)
renderer = vtk.vtkRenderer()
renderer.AddActor(actor)
renderer.SetBackground(1.0, 1.0, 1.0)
cam = renderer.GetActiveCamera()
cam.SetParallelProjection(True)
cam.SetPosition(*camera.position)
cam.SetFocalPoint(*camera.focal_point)
cam.SetViewUp(*camera.view_up)
cam.SetParallelScale(camera.parallel_scale)
renderer.ResetCameraClippingRange()
render_window = vtk.vtkRenderWindow()
render_window.SetOffScreenRendering(1)
render_window.AddRenderer(renderer)
render_window.SetSize(*camera.image_size)
render_window.Render()
w2i = vtk.vtkWindowToImageFilter()
w2i.SetInput(render_window)
w2i.SetInputBufferTypeToRGB()
w2i.ReadFrontBufferOff()
w2i.Update()
writer = vtk.vtkPNGWriter()
writer.SetFileName(str(out_png))
writer.SetInputConnection(w2i.GetOutputPort())
writer.Write()
if camera.rotate_degrees:
img = Image.open(out_png).rotate(camera.rotate_degrees, expand=True)
img.save(out_png)
return camera.color_array, vmin, vmax
def compose_comparison(panels, out_png, title=None):
"""Compose framed, labelled renders side-by-side.
panels is a list of dicts {"label": str, "image": Path,
"scale": (color_array, vmin, vmax) | None}. Panels carrying a scale (the
numerical results) get a horizontal colorbar built from the same range.
"""
imgs = [Image.open(p["image"]).convert("RGB") for p in panels]
aspects = [im.width / im.height for im in imgs]
panel_h = 4.6
fig = plt.figure(figsize=(sum(aspects) * panel_h + 0.4, panel_h + 0.9))
gs = fig.add_gridspec(
2,
len(panels),
height_ratios=[1.0, 0.045],
width_ratios=aspects,
hspace=0.06,
wspace=0.12,
)
for i, (panel, img) in enumerate(zip(panels, imgs)):
ax = fig.add_subplot(gs[0, i])
ax.imshow(img, aspect="auto")
ax.set_xticks([])
ax.set_yticks([])
for spine in ax.spines.values():
spine.set_visible(True)
spine.set_linewidth(1.4)
spine.set_edgecolor("0.25")
ax.set_title(panel["label"], fontsize=13, pad=8)
if panel.get("scale") is not None:
color_array, vmin, vmax = panel["scale"]
sub = gs[1, i].subgridspec(1, 3, width_ratios=[1, 6, 1])
cax = fig.add_subplot(sub[0, 1])
sm = ScalarMappable(
norm=Normalize(vmin=vmin, vmax=vmax),
cmap=_lut_to_cmap(vmin, vmax),
)
cbar = fig.colorbar(sm, cax=cax, orientation="horizontal")
ticks = MaxNLocator(nbins=4, steps=[1, 2, 2.5, 5, 10]).tick_values(
vmin, vmax
)
ticks = ticks[(ticks >= vmin) & (ticks <= vmax)]
cbar.set_ticks(ticks)
cbar.ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda v, _: f"{v:g}"))
cbar.ax.tick_params(labelsize=7)
cbar.set_label(FIELD_LABELS.get(color_array, color_array), fontsize=10)
if title:
fig.suptitle(title, fontsize=15, y=1.02)
fig.savefig(out_png, dpi=200, bbox_inches="tight", facecolor="white")
plt.close(fig)
def _add_slice_grid(ax):
ax.grid(True, which="major", color="0.5", linewidth=0.6, alpha=0.55)
ax.set_axisbelow(False)
def _style_slice_axes(ax, y_min, y_max, z_min, z_max, xticks, yticks):
ax.set_xlabel("$Y$ Position [mm]")
ax.set_ylabel("$Z$ Position [mm]")
ax.set_xticks(xticks)
ax.set_yticks(yticks)
ax.set_xlim(y_min, y_max)
ax.set_ylim(z_min, z_max)
ax.set_aspect("equal")
_add_slice_grid(ax)
def _ticks_100(lo, hi):
"""100 mm tick locations spanning [lo, hi]."""
return np.arange(np.ceil(lo / 100.0) * 100.0, hi + 1.0, 100.0)
def _exp_extent(comp, n, exp_png):
"""imshow extent (y_min, y_max, z_min, z_max) in mm for an experimental image
from its pixel calibration, or None if no calibration is registered."""
calib = EXP_CALIB.get((comp, n))
if calib is None:
return None
px_per_100mm, x_at_ym500, y_at_z0 = calib
width, height = Image.open(exp_png).size
ppm = px_per_100mm / 100.0
y_min = -500.0 - x_at_ym500 / ppm
y_max = -500.0 + (width - x_at_ym500) / ppm
z_max = y_at_z0 / ppm
z_min = (y_at_z0 - height) / ppm
return (y_min, y_max, z_min, z_max)
def compose_slice_figure(slice_dir, n, comp, vmin, vmax, exp_png, locator_png, out_png):
"""2x2 figure: experimental contour, slice locator + colorbar, and the
Flow360 and OpenFOAM contours (OpenFOAM omitted if reference data absent)."""
ds = _read_dataset(_find_slice_file(slice_dir, n))
ds = _ensure_point_data(ds, "velocity_adimensionalized")
pts = vtk_np.vtk_to_numpy(ds.GetPoints().GetData())
vec = vtk_np.vtk_to_numpy(ds.GetPointData().GetArray("velocity_adimensionalized"))
if vec.ndim == 1:
raise RuntimeError(
"velocity_adimensionalized is scalar; expected a 3-component vector"
)
values = vec[:, COMPONENT_INDEX[comp]]
(y_min, y_max), (z_min, z_max) = SLICE_CROPS[n]
z_min = max(z_min, Z_TRIM_MM)
y, z = pts[:, 1], pts[:, 2]
mask = (y >= y_min - 50) & (y <= y_max + 50) & (z >= z_min - 50) & (z <= z_max + 50)
y, z, values = y[mask], z[mask], values[mask]
fig = plt.figure(figsize=(10, 8))
gs = fig.add_gridspec(2, 2, hspace=0.30, wspace=0.20)
levels = np.linspace(vmin, vmax, 41)
xticks = _ticks_100(y_min, y_max)
yticks = _ticks_100(z_min, z_max)
comp_label = f"$V_{comp[1].lower()}$ normalized"
# Top-left: experimental contour, placed/ticked from its pixel calibration.
ax_exp = fig.add_subplot(gs[0, 0])
ext = _exp_extent(comp, n, exp_png) or (y_min, y_max, z_min, z_max)
ey_min, ey_max, ez_min, ez_max = ext
ax_exp.imshow(
plt.imread(exp_png),
extent=[ey_min, ey_max, ez_min, ez_max],
aspect="equal",
interpolation="bilinear",
origin="upper",
)
ax_exp.set_title("Experimental Result")
ax_exp.set_xlabel("$Y$ Position [mm]")
ax_exp.set_ylabel("$Z$ Position [mm]")
ez_min = max(ez_min, Z_TRIM_MM)
ax_exp.set_xticks(_ticks_100(ey_min, ey_max))
ax_exp.set_yticks(_ticks_100(ez_min, ez_max))
ax_exp.set_xlim(ey_min, ey_max)
ax_exp.set_ylim(ez_min, ez_max)
_add_slice_grid(ax_exp)
# Top-right: slice locator + horizontal colorbar.
gs_tr = gs[0, 1].subgridspec(2, 1, height_ratios=[3, 1], hspace=0.05)
ax_loc = fig.add_subplot(gs_tr[0])
ax_loc.imshow(plt.imread(locator_png))
ax_loc.set_axis_off()
ax_cbar = fig.add_subplot(gs_tr[1])
sm = ScalarMappable(norm=Normalize(vmin=vmin, vmax=vmax), cmap=CMAP)
cbar = fig.colorbar(sm, cax=ax_cbar, orientation="horizontal")
cbar.set_ticks(np.linspace(vmin, vmax, 5))
cbar.set_label(comp_label, fontsize=11)
# Bottom-left: Flow360 contour.
ax_sim = fig.add_subplot(gs[1, 0])
ax_sim.tricontourf(
y, z, values, levels=levels, cmap=CMAP, vmin=vmin, vmax=vmax, extend="both"
)
ax_sim.set_title("Simulation (Flow360)")
_style_slice_axes(ax_sim, y_min, y_max, z_min, z_max, xticks, yticks)
# Bottom-right: OpenFOAM contour from the matching wake-exp slice.
of = read_openfoam_slice(n)
if of is not None:
of_y, of_z, of_unorm = of
of_values = of_unorm[:, COMPONENT_INDEX[comp]]
of_mask = (
(of_y >= y_min - 50)
& (of_y <= y_max + 50)
& (of_z >= z_min - 50)
& (of_z <= z_max + 50)
)
ax_of = fig.add_subplot(gs[1, 1])
ax_of.tricontourf(
of_y[of_mask],
of_z[of_mask],
of_values[of_mask],
levels=levels,
cmap=CMAP,
vmin=vmin,
vmax=vmax,
extend="both",
)
ax_of.set_title("Simulation (OpenFOAM)")
_style_slice_axes(ax_of, y_min, y_max, z_min, z_max, xticks, yticks)
fig.savefig(out_png, dpi=200, bbox_inches="tight", facecolor="white")
plt.close(fig)
def _final_force_coeffs(custom_forces, force_name):
"""Final converged force coefficients for a ForceOutput, or None if absent."""
if force_name not in set(custom_forces.custom_force_names):
return None
force = custom_forces.get_custom_force_by_name(force_name)
force.reload_data()
final = force.as_dataframe().iloc[-1].to_dict()
def _coeff(name):
for key in (f"total{name}", name, f"{force_name}_{name}"):
if key in final:
return final[key]
for key in final:
if key.lower().endswith(name.lower()):
return final[key]
return None
return {c: _coeff(c) for c in ("CFx", "CFy", "CFz", "CMx", "CMy", "CMz")}
def write_loads_comparison(case):
"""Flow360 (wing_all surface forces) vs OpenFOAM reference wing loads, on the
same 1 m^2 reference area. Saves the published results/loads_vs_openfoam.csv
and returns the table."""
rows = []
f360 = _final_force_coeffs(case.results.custom_forces, "wing_all")
if f360 is not None:
rows.append({"source": "Flow360", "CFz": f360["CFz"], "CFx": f360["CFx"]})
rows.append(
{
"source": "OpenFOAM",
"CFz": OPENFOAM_WING_LOADS["CFz"],
"CFx": OPENFOAM_WING_LOADS["CFx"],
}
)
df = pd.DataFrame(rows)
df.to_csv(RESULTS_DIR / "loads_vs_openfoam.csv", index=False)
return df
Download and extract the case results¶
Pull the slice and surface archives from the finished in-kernel case and load
the per-slice value ranges used by the wake-plane contours.
In [ ]:
download_dir = RESULTS_DIR / "_downloads"
render_dir = RESULTS_DIR / "_renders"
render_dir.mkdir(parents=True, exist_ok=True)
slice_dir, surface_dir = fetch_results(case, download_dir)
with open(EXTENTS_FILE) as f:
extents = json.load(f)
Integrated wing loads: Flow360 vs OpenFOAM¶
Table of total (pressure + viscous) wing-load coefficients on the shared 1 m²
reference area: Flow360's wing_all surface forces against the OpenFOAM
reference. Produces results/loads_vs_openfoam.csv.
In [ ]:
loads_df = write_loads_comparison(case)
loads_df
Surface flow-vis: wing underside and endplate¶
Off-screen surface-LIC of the skin-friction field coloured by Cp, composed
Experiment | Flow360 | OpenFOAM. Flow360 uses CfVec for the LIC streaks;
OpenFOAM uses UwNorm on the same Cp colormap (the endplate view mirrors the
left-half OpenFOAM export onto the right). Produces
results/compare_wing_underside.png and results/compare_endplate_underside.png.
In [ ]:
surfaces = read_wing_surfaces(surface_dir)
flow_wing = render_dir / "flow360_wing_underside.png"
flow_endplate = render_dir / "flow360_endplate_underside.png"
wing_scale = render_surface_lic(surfaces, WING_UNDERSIDE_CAMERA, flow_wing)
endplate_scale = render_surface_lic(surfaces, ENDPLATE_UNDERSIDE_CAMERA, flow_endplate)
# OpenFOAM surface flow-vis: same cameras and Cp colormap, LIC driven by UwNorm.
of_surface = read_openfoam_wing_surface()
of_wing = of_endplate = None
of_wing_scale = of_endplate_scale = None
if of_surface is not None:
of_wing = render_dir / "openfoam_wing_underside.png"
of_endplate = render_dir / "openfoam_endplate_underside.png"
of_wing_scale = render_surface_lic(
of_surface, WING_UNDERSIDE_CAMERA, of_wing, vector_array="UwNorm"
)
of_endplate_scale = render_surface_lic(
_mirror_surface_y(of_surface),
ENDPLATE_UNDERSIDE_CAMERA,
of_endplate,
vector_array="UwNorm",
)
wing_panels = [
{"label": "Experiment", "image": REF_DIR / "wing_underside.png", "scale": None},
{"label": "Flow360", "image": flow_wing, "scale": wing_scale},
]
endplate_panels = [
{"label": "Experiment", "image": REF_DIR / "endplate_underside.png", "scale": None},
{"label": "Flow360", "image": flow_endplate, "scale": endplate_scale},
]
if of_wing is not None:
wing_panels.append({"label": "OpenFOAM", "image": of_wing, "scale": of_wing_scale})
endplate_panels.append(
{"label": "OpenFOAM", "image": of_endplate, "scale": of_endplate_scale}
)
compose_comparison(
wing_panels,
RESULTS_DIR / "compare_wing_underside.png",
title="Wing underside -- surface skin-friction flow-vis",
)
compose_comparison(
endplate_panels,
RESULTS_DIR / "compare_endplate_underside.png",
title="Endplate -- surface skin-friction flow-vis",
)
Wake-plane slices: Vy and Vz at planes 1, 3 and 5¶
For each of the three wake planes and both velocity components, build the 2×2
comparison of the experimental PIV image, the slice locator, and the Flow360 and
OpenFOAM free-stream-normalized contours over the per-slice value range from
extents.json. Produces the six results/compare_V{y,z}_slice{1,3,5}.png.
In [ ]:
for n in SLICE_NUMBERS:
locator = REF_DIR / f"slice{n}.png"
for comp in COMPONENTS:
vmin, vmax = extents[str(n)][comp]
exp = REF_DIR / f"{comp}-slice{n}-exp.png"
out = RESULTS_DIR / f"compare_{comp}_slice{n}.png"
compose_slice_figure(slice_dir, n, comp, vmin, vmax, exp, locator, out)
print(f"Done. See {RESULTS_DIR}")