SCWS Common Research Model: Static & Dynamic Rolling Aerodynamics¶

This case reproduces the Second AIAA Stability & Control Prediction Workshop (S&CPW2) validation study on the NASA/Boeing Common Research Model, comparing Flow360 against experimental data from NASA Langley's 12-Foot Low-Speed Tunnel.

The study has two branches. The static branch is a steady RANS (Spalart-Allmaras) sweep over angle of attack from 0° to 20° at Mach 0.052 and a mean-aerodynamic-chord Reynolds number of 200,000; each case forks from the previous one to warm-start the solve. The dynamic branch is a single unsteady SA-DDES case at α = 3°, Mach 0.0358, and Re = 140,000, with the aircraft forced through a ±5° sinusoidal roll oscillation at 0.0264 Hz over 20 full cycles at 1° per time step. The notebook runs top-to-bottom in one kernel: load the geometry, submit both branches, wait for completion, then post-process the in-kernel case objects into the four published figures (static coefficient and lift/drag curves, and the dynamic roll time-histories and roll-phase hysteresis).

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 standard scientific-Python packages:

pip install numpy pandas matplotlib

Run the cells top to bottom in a single kernel.

Imports¶

In [ ]:
import json
import math
import os
import re
from pathlib import Path

import flow360 as fl
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from flow360.examples import download_benchmark_assets

Helper modules¶

The forced-roll case drives the aircraft with an angle expression evaluated by the solver. The two sibling helper modules that ship with this case are inlined here so the notebook is self-contained: the first derives the roll-angle expression from the physical oscillation parameters, and the second is the expression string that fl.AngleExpression consumes.

In [ ]:
# Inlined from the case's sibling helpers so the notebook needs no local imports
# (no `import rolling_expression` / `import rolling_function_definition`).

# --- rolling_function_definition.py: derive the roll-angle expression ---
# The aircraft rolls sinusoidally by +/-5 deg at 0.0264 Hz. The solver evaluates
# the angle expression in nondimensional time, so the physical angular frequency
# is scaled by (length_unit / reference_speed).
amplitude_deg = 5.0                  # roll amplitude [deg]
frequency_hz = 0.0264                # oscillation frequency [Hz]
length_unit = 0.0254                 # length scale [m]
reference_speed = 340.074144758659   # reference velocity [m/s]

amplitude_rad = amplitude_deg * math.pi / 180.0
omega = 2.0 * math.pi * frequency_hz          # angular frequency [rad/s]
time_scale = length_unit / reference_speed    # nondimensional time scaling
expression_string = f"{amplitude_rad}*sin({omega * time_scale}*t)"

# --- rolling_expression.py: the angle expression consumed by fl.AngleExpression ---
Expression = f"\n{expression_string}\n"

Input data¶

The post-processing compares Flow360 against wind-tunnel measurements and the workshop participants' submissions. Fetch those reference files from the public benchmark bucket; this recreates a local ./ref_data/ directory that the plotting cells read from.

In [ ]:
download_benchmark_assets("SCWS", "ref_data")

Load project¶

The root asset for this case is the published CRM geometry. Download it from the public benchmark bucket and start a fresh project from it, no private project id is needed.

In [ ]:
root_asset_files = download_benchmark_assets("SCWS", "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="SCWS Common Research Model")

Case constants¶

Reference geometry, operating points, and sweep definitions for both branches. SOLVER_VERSION and the SCWS_RUN_STATIC / SCWS_RUN_DYNAMIC toggles can be overridden through the environment; by default both branches run.

In [ ]:
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")

REFERENCE_AREA_IN2 = 594720
MOMENT_CENTER_IN = (1325.9, 0.0, 177.95)
MAC_IN = 275.8
SPAN_IN = 2313.5

STATIC_MACH = 0.052
STATIC_RE = 200000
STATIC_WARMUP_ALPHA = 0.0
STATIC_ALPHAS = (4.0, 8.0, 12.0, 16.0, 20.0)

DYNAMIC_MACH = 0.0358
DYNAMIC_RE = 140000
DYNAMIC_ALPHA = 3.0
DYNAMIC_FREQUENCY_HZ = 0.0264
DYNAMIC_ANGLE_PER_STEP_DEG = 1.0
DYNAMIC_TOTAL_CYCLES = 20

RUN_STATIC = os.environ.get("SCWS_RUN_STATIC", "1") != "0"
RUN_DYNAMIC = os.environ.get("SCWS_RUN_DYNAMIC", "1") != "0"

Meshing setup¶

Both branches share the same beta-mesher setup. The wall surfaces are discovered from the geometry (everything that is not the farfield or the rotor/static interface), and two custom volume zones are built: a rotating zone bounded by the walls and the interface, and a farfield zone. make_common_meshing sets the surface, boundary-layer, and curvature spacing and keeps the farfield/interface spacing unchanged.

In [ ]:
def get_wall_surfaces(geometry):
    entity_info = getattr(geometry, "entity_info", None)
    if entity_info is not None and hasattr(entity_info, "get_boundaries"):
        names = []
        for surface in entity_info.get_boundaries():
            surface_name = getattr(surface, "private_attribute_id", None) or getattr(
                surface, "name", None
            )
            if not surface_name:
                continue
            if "Farfield" in surface_name or "Interface" in surface_name:
                continue
            names.append(surface_name)
        if names:
            return [geometry[name] for name in names]
    return [geometry["Aircraft"]]


def make_volumes(geometry, wall_surfaces):
    rotating_volume = fl.CustomVolume(
        boundaries=[*wall_surfaces, geometry["Interface"]],
        name="RotatingVolume",
    )
    farfield_volume = fl.CustomVolume(
        boundaries=[geometry["Farfield"], geometry["Interface"]],
        name="FarfieldVolume",
    )
    return rotating_volume, farfield_volume


def make_common_meshing(project, rotating_volume, farfield_volume, wall_surfaces):
    farfield = fl.UserDefinedFarfield()
    return fl.MeshingParams(
        defaults=fl.MeshingDefaults(
            surface_max_edge_length=(125000 * 2 * 3.13 / 12.5) * fl.u.inch,
            boundary_layer_first_layer_thickness=0.125 * fl.u.mm,
            curvature_resolution_angle=10 * fl.u.deg,
        ),
        volume_zones=[
            farfield,
            fl.CustomZones(entities=[farfield_volume, rotating_volume]),
        ],
        refinements=[
            fl.SurfaceRefinement(
                faces=wall_surfaces,
                max_edge_length=8 * fl.u.inch,
                curvature_resolution_angle=5 * fl.u.deg,
            ),
            fl.SurfaceRefinement(
                faces=[project.geometry["Interface"]],
                max_edge_length=100 * fl.u.inch,
                curvature_resolution_angle=5 * fl.u.deg,
            ),
            fl.PassiveSpacing(
                faces=[project.geometry["Farfield"], project.geometry["Interface"]],
                type="unchanged",
            ),
        ],
    )

Physics setup¶

The reference geometry and the operating condition are shared builders. The operating condition is defined from Mach and a Reynolds number per mesh unit, using the model's inch length unit. The boundary-condition and solver stacks differ between the two branches, steady RANS Spalart-Allmaras for the static sweep, unsteady SA-DDES plus a rolling fl.Rotation model for the dynamic case, so they are assembled inside the two SimulationParams builders below.

In [ ]:
def make_reference_geometry():
    return fl.ReferenceGeometry(
        area=REFERENCE_AREA_IN2 * fl.u.inch**2,
        moment_center=MOMENT_CENTER_IN * fl.u.inch,
        moment_length=(SPAN_IN, MAC_IN, SPAN_IN) * fl.u.inch,
    )


def make_operating_condition(mach, reynolds_length_unit, alpha_deg):
    return fl.AerospaceCondition.from_mach_reynolds(
        mach=mach,
        reynolds_mesh_unit=reynolds_length_unit,
        temperature=518 * fl.u.R,
        project_length_unit=1 * fl.u.inch,
        alpha=alpha_deg * fl.u.deg,
    )

Simulation Params¶

make_static_params assembles a steady RANS case at a given angle of attack: adaptive-CFL steady time stepping, a tightly converged Navier-Stokes solver with the low-Mach preconditioner, and the Spalart-Allmaras turbulence model.

make_dynamic_params assembles the unsteady rolling case: the rotating zone is given its roll axis and center, the time step and step count are derived from the oscillation frequency and the 1°-per-step / 20-cycle schedule, and the model stack adds SA-DDES (detached-eddy hybrid) and an fl.Rotation model driven by the inlined roll-angle Expression.

In [ ]:
def make_static_params(project, alpha_deg):
    wall_surfaces = get_wall_surfaces(project.geometry)
    rotating_volume, farfield_volume = make_volumes(project.geometry, wall_surfaces)
    with fl.SI_unit_system:
        return fl.SimulationParams(
            meshing=make_common_meshing(
                project, rotating_volume, farfield_volume, wall_surfaces
            ),
            reference_geometry=make_reference_geometry(),
            operating_condition=make_operating_condition(
                mach=STATIC_MACH,
                reynolds_length_unit=STATIC_RE / MAC_IN,
                alpha_deg=alpha_deg,
            ),
            time_stepping=fl.Steady(
                CFL=fl.AdaptiveCFL(convergence_limiting_factor=0.7),
                max_steps=7500,
            ),
            models=[
                fl.Freestream(surfaces=project.geometry["Farfield"]),
                fl.Wall(surfaces=wall_surfaces),
                fl.Fluid(
                    navier_stokes_solver=fl.NavierStokesSolver(
                        absolute_tolerance=1e-10,
                        kappa_MUSCL=-1,
                        low_mach_preconditioner=True,
                    ),
                    turbulence_model_solver=fl.SpalartAllmaras(
                        absolute_tolerance=1e-8,
                        equation_evaluation_frequency=2,
                    ),
                ),
            ],
        )


def make_dynamic_params(project):
    wall_surfaces = get_wall_surfaces(project.geometry)
    rotating_volume, farfield_volume = make_volumes(project.geometry, wall_surfaces)
    rotating_volume.axis = (-1, 0, 0)
    rotating_volume.center = MOMENT_CENTER_IN * fl.u.inch

    omega = 2.0 * np.pi * DYNAMIC_FREQUENCY_HZ
    time_step = (np.pi * DYNAMIC_ANGLE_PER_STEP_DEG) / (omega * 180.0)
    total_steps = int((360.0 * DYNAMIC_TOTAL_CYCLES) / DYNAMIC_ANGLE_PER_STEP_DEG)

    with fl.SI_unit_system:
        return fl.SimulationParams(
            meshing=make_common_meshing(
                project, rotating_volume, farfield_volume, wall_surfaces
            ),
            reference_geometry=make_reference_geometry(),
            operating_condition=make_operating_condition(
                mach=DYNAMIC_MACH,
                reynolds_length_unit=DYNAMIC_RE / MAC_IN,
                alpha_deg=DYNAMIC_ALPHA,
            ),
            time_stepping=fl.Unsteady(
                CFL=fl.AdaptiveCFL(convergence_limiting_factor=0.7),
                step_size=time_step,
                steps=total_steps,
            ),
            models=[
                fl.Freestream(surfaces=project.geometry["Farfield"]),
                fl.Wall(surfaces=wall_surfaces),
                fl.Fluid(
                    navier_stokes_solver=fl.NavierStokesSolver(
                        absolute_tolerance=1e-10,
                        relative_tolerance=1e-2,
                        kappa_MUSCL=-1,
                        low_mach_preconditioner=True,
                    ),
                    turbulence_model_solver=fl.SpalartAllmaras(
                        absolute_tolerance=1e-8,
                        equation_evaluation_frequency=4,
                        relative_tolerance=1e-2,
                        hybrid_model=fl.DetachedEddySimulation(
                            shielding_function="DDES"
                        ),
                    ),
                ),
                fl.Rotation(
                    name="RollingFunction",
                    volumes=rotating_volume,
                    spec=fl.AngleExpression(Expression),
                ),
            ],
        )

Submit cases¶

Submit both branches to Flow360. The static sweep runs the 0° warm-up first and then each angle of attack in turn, forking every case from the previous one so the solver warm-starts. The dynamic branch submits the single rolling case. All submitted cases are collected in submitted_cases for the wait and post-processing steps.

In [ ]:
def run_static_cases(project):
    submitted_cases = []
    previous_case = None
    for alpha_deg in (STATIC_WARMUP_ALPHA, *STATIC_ALPHAS):
        case_name = f"scws_static_alpha_{alpha_deg:.1f}_{SOLVER_VERSION}"
        params = make_static_params(project, alpha_deg)
        case = project.run_case(
            name=case_name,
            params=params,
            solver_version=SOLVER_VERSION,
            use_beta_mesher=True,
            fork_from=previous_case,
        )
        print(f"Submitted static case: {case.name} ({case.id})")
        submitted_cases.append(case)
        previous_case = case
    return submitted_cases


def run_dynamic_cases(project):
    case_name = (
        f"scws_dynamic_alpha_{DYNAMIC_ALPHA:.1f}_phi_5.0_freq_0.0264_"
        f"cycles_{DYNAMIC_TOTAL_CYCLES}_{SOLVER_VERSION}"
    )
    case = project.run_case(
        name=case_name,
        params=make_dynamic_params(project),
        solver_version=SOLVER_VERSION,
        use_beta_mesher=True,
    )
    print(f"Submitted dynamic case: {case.name} ({case.id})")
    return [case]


if not RUN_STATIC and not RUN_DYNAMIC:
    raise RuntimeError("Both SCWS_RUN_STATIC and SCWS_RUN_DYNAMIC are disabled.")

submitted_cases = []
if RUN_STATIC:
    submitted_cases.extend(run_static_cases(project))
if RUN_DYNAMIC:
    submitted_cases.extend(run_dynamic_cases(project))

Wait for completion¶

Block until every submitted case has finished. This can take a long time, the static sweep is six steady solves and the dynamic case integrates 20 roll cycles at 1° per step (7200 physical steps).

In [ ]:
for case in submitted_cases:
    case.wait()

Postprocessing¶

The published results are four figures:

  1. Static rolling coefficients vs α: normal force, axial force, and pitching moment ($C_N$, $C_A$, $C_m$) against angle of attack.
  2. Static lift and drag vs α: $C_L$ and $C_D$ against angle of attack.
  3. Dynamic rolling time histories: the roll angle and the roll/yaw increments $\Delta C_l$, $\Delta C_n$ over time, with a $\pm\sigma$ phase-scatter band.
  4. Dynamic rolling coefficients vs roll phase: $\Delta C_l$ and $\Delta C_n$ plotted against the roll angle $\Delta\phi$ (hysteresis loops).

Each is compared against NASA Langley wind-tunnel data and the workshop participants' submissions. The cells below read results straight from the in-kernel case objects (no project is re-opened) and rebuild exactly these four figures. The unpublished intermediate CSVs written by the AE's script are omitted.

Setup¶

Styling colors, the static/dynamic case-name prefixes, the dynamic-analysis constants (phase binning, smoothing, cycles to skip while the wake develops), and the small numeric/plotting helpers that the figures share.

In [ ]:
FLOW360_COLOR = "#00A67E"   # Flow360 results
REFERENCE_COLOR = "black"   # experimental / reference data

STATIC_PREFIX = "scws_static_alpha_"
DYNAMIC_PREFIX = "scws_dynamic_alpha_"

SKIP_DYNAMIC_CYCLES = 3
DYNAMIC_SMOOTH_WINDOW = 7
DYNAMIC_PHASE_BINS = 180
DYNAMIC_ROLL_EXPRESSION = "0.08726646259971647*sin(1.23892180705839e-05*t)"

CASE_DIR = Path(".")


def parse_static_alpha(case_name):
    match = re.search(r"^scws_static_alpha_(-?\d+(?:\.\d+)?)_", case_name)
    if not match:
        raise ValueError(f"Could not parse static alpha from {case_name!r}")
    return float(match.group(1))


def extract_dynamic_amplitude(expression):
    match = re.search(r"^([^*]+)\*sin\(", expression.replace(" ", ""))
    if not match:
        raise ValueError("Could not extract dynamic roll amplitude.")
    return float(match.group(1))


def sum_coefficients(frame, suffix):
    columns = [column for column in frame.columns if str(column).endswith(f"_{suffix}")]
    if not columns:
        raise RuntimeError(f"No columns ending in _{suffix} were found.")
    return frame[columns].sum(axis=1).to_numpy(dtype=float)


def moving_average_circular(values, width):
    if width <= 1:
        return values
    if width % 2 == 0:
        width += 1
    pad = width // 2
    extended = np.r_[values[-pad:], values, values[:pad]]
    kernel = np.ones(width) / width
    return np.convolve(extended, kernel, mode="valid")


def phase_bin_statistics(theta, values, n_bins, lower_percentile, upper_percentile):
    edges = np.linspace(0.0, 2.0 * np.pi, n_bins + 1)
    centers = 0.5 * (edges[:-1] + edges[1:])
    bucket_index = np.digitize(theta, edges) - 1
    bucket_index = np.clip(bucket_index, 0, n_bins - 1)

    center = np.full(n_bins, np.nan)
    lower = np.full(n_bins, np.nan)
    upper = np.full(n_bins, np.nan)

    for bucket in range(n_bins):
        bucket_values = values[bucket_index == bucket]
        if len(bucket_values) == 0:
            continue
        center[bucket] = np.median(bucket_values)
        lower[bucket] = np.percentile(bucket_values, lower_percentile)
        upper[bucket] = np.percentile(bucket_values, upper_percentile)

    def fill_nan(values_1d):
        indices = np.arange(len(values_1d))
        valid = np.isfinite(values_1d)
        if valid.sum() < 2:
            return values_1d
        return np.interp(indices, indices[valid], values_1d[valid])

    center = moving_average_circular(fill_nan(center), DYNAMIC_SMOOTH_WINDOW)
    lower = moving_average_circular(fill_nan(lower), DYNAMIC_SMOOTH_WINDOW)
    upper = moving_average_circular(fill_nan(upper), DYNAMIC_SMOOTH_WINDOW)
    return centers, center, lower, upper


def periodic_interp(theta_query, theta_grid, value_grid):
    theta_wrapped = theta_query % (2.0 * np.pi)
    theta_extended = np.r_[theta_grid, theta_grid[0] + 2.0 * np.pi]
    values_extended = np.r_[value_grid, value_grid[0]]
    return np.interp(theta_wrapped, theta_extended, values_extended)


def nice_limits(*arrays, pad_frac=0.08):
    values = np.concatenate([np.asarray(array, dtype=float).ravel() for array in arrays])
    values = values[np.isfinite(values)]
    if len(values) == 0:
        return None
    ymin, ymax = float(values.min()), float(values.max())
    if np.isclose(ymin, ymax):
        return ymin - 1.0, ymax + 1.0
    pad = (ymax - ymin) * pad_frac
    return ymin - pad, ymax + pad


def fill_parametric_band(axis, x, y_lower, y_upper, **kwargs):
    axis.fill(
        np.r_[x, x[::-1]],
        np.r_[y_upper, y_lower[::-1]],
        linewidth=0.0,
        **kwargs,
    )


def configure_axes_grid(axes):
    for axis in axes:
        axis.minorticks_on()
        axis.grid(True, which="major", linewidth=0.8, alpha=0.7)
        axis.grid(True, which="minor", linestyle=":", linewidth=0.5, alpha=0.5)


def load_reference_curve(filename):
    return pd.read_csv(CASE_DIR / "ref_data" / filename, header=None, names=["x", "y"])


def load_participant_data():
    path = CASE_DIR / "ref_data" / "participants_data.csv"
    if not path.exists():
        path = CASE_DIR / "ref_data" / "participants.csv"
    return pd.read_csv(path)


def plot_participants(axis, participants, quantity, xy=False, linewidth=1.0):
    rows = participants[participants["quantity"].eq(quantity)]
    if rows.empty:
        return np.array([])

    if xy:
        column_pairs = [
            (column, f"{column[:-2]}_y")
            for column in participants.columns
            if column.endswith("_x") and f"{column[:-2]}_y" in participants.columns
        ]
    else:
        column_pairs = [
            ("x", column)
            for column in participants.columns
            if column not in ("x", "quantity") and not column.endswith(("_x", "_y"))
        ]

    y_values = []
    label = "Workshop Results"
    for x_column, y_column in column_pairs:
        data = rows[[x_column, y_column]].dropna()
        if data.empty:
            continue
        x = data[x_column].to_numpy(dtype=float)
        y = data[y_column].to_numpy(dtype=float)
        axis.plot(
            x,
            y,
            color="0.45",
            alpha=0.18 if xy else 0.55,
            linestyle="None" if xy else "-",
            linewidth=0 if xy else linewidth,
            marker="." if xy else None,
            markersize=2 if xy else None,
            label=label,
            zorder=2,
        )
        label = "_nolegend_"
        y_values.append(y)
    return np.concatenate(y_values) if y_values else np.array([])

Case records¶

Group the in-kernel submitted cases into the static and dynamic branches by name, and load the workshop-participant reference table used by every figure.

In [ ]:
records = [
    {"id": case.id, "name": case.name, "case": case} for case in submitted_cases
]
static_records = [r for r in records if r["name"].startswith(STATIC_PREFIX)]
dynamic_records = [r for r in records if r["name"].startswith(DYNAMIC_PREFIX)]

participants = load_participant_data()

Static coefficients¶

For each static case, read the surface-forces history (excluding the sting/TUBE mounting surfaces), average the last 10% of the converged history, and sum the per-surface contributions into the workshop coefficients $C_N$, $C_A$, $C_m$, $C_L$, and $C_D$.

In [ ]:
STATIC_COEFFICIENTS = {"CN": "CFz", "CA": "CFx", "Cm": "CMy", "CL": "CL", "CD": "CD"}

rows = []
for record in static_records:
    alpha_deg = parse_static_alpha(record["name"])
    surface_forces = record["case"].results.surface_forces.as_dataframe()
    surface_forces = surface_forces.loc[:, ~surface_forces.columns.str.contains("TUBE")]
    tail = surface_forces.tail(max(1, int(len(surface_forces) * 0.1)))

    row = {"case_name": record["name"], "alpha_deg": alpha_deg}
    for paper_name, suffix in STATIC_COEFFICIENTS.items():
        row[paper_name] = sum_coefficients(tail, suffix).mean()
    rows.append(row)

flow_summary = pd.DataFrame(rows).sort_values("alpha_deg").reset_index(drop=True)
flow_summary

Figure 1: Static rolling coefficients vs angle of attack¶

$C_N$, $C_A$, and $C_m$ versus α, comparing Flow360 (RANS) against the experiment and the workshop participants.

In [ ]:
panels = ["CN", "CA", "Cm"]
plot_rows = []
for panel in panels:
    for _, row in flow_summary.iterrows():
        plot_rows.append(
            {"series": "Flow360", "coefficient": panel,
             "alpha_deg": row["alpha_deg"], "value": row[panel]}
        )
    exp_curve = load_reference_curve(f"exp_{panel}.csv")
    for _, row in exp_curve.iterrows():
        plot_rows.append(
            {"series": "Experiment", "coefficient": panel,
             "alpha_deg": float(row["x"]), "value": float(row["y"])}
        )
static_cn_ca_cm = pd.DataFrame(plot_rows)

fig, axes = plt.subplots(
    1, 3, sharex=True, figsize=(10.5, 5.2),
    constrained_layout=True, gridspec_kw={"wspace": 0.25},
)
configure_axes_grid(axes)
ylabels = {"CN": r"$C_N$", "CA": r"$C_A$", "Cm": r"$C_m$"}

for axis, panel in zip(axes, panels):
    flow_values = static_cn_ca_cm[
        (static_cn_ca_cm["series"] == "Flow360")
        & (static_cn_ca_cm["coefficient"] == panel)
    ].sort_values("alpha_deg")
    exp_values = static_cn_ca_cm[
        (static_cn_ca_cm["series"] == "Experiment")
        & (static_cn_ca_cm["coefficient"] == panel)
    ].sort_values("alpha_deg")
    participant_y = plot_participants(axis, participants, panel)
    axis.plot(
        flow_values["alpha_deg"], flow_values["value"],
        color=FLOW360_COLOR, marker="o", linewidth=1.8, markersize=4.5,
        label="RANS Flow360",
    )
    axis.plot(
        exp_values["alpha_deg"], exp_values["value"],
        color=REFERENCE_COLOR, linestyle=":", marker="o", linewidth=1.8,
        markersize=4.5, label="Experiment",
    )
    axis.set_xlim(0, 20)
    axis.set_xticks([0, 5, 10, 15, 20])
    axis.set_ylim(*nice_limits(flow_values["value"], exp_values["value"], participant_y))
    axis.set_ylabel(ylabels[panel])
    axis.set_xlabel(r"$\alpha\ [deg]$")
    axis.legend()

os.makedirs("results", exist_ok=True)
fig.savefig("results/static_cn_ca_cm_vs_alpha.png", dpi=200, bbox_inches="tight")
plt.show()

Figure 2: Static lift and drag vs angle of attack¶

$C_L$ and $C_D$ versus α, again comparing Flow360 against the experiment and the workshop participants.

In [ ]:
polar_rows = []
for panel in ("CL", "CD"):
    for _, row in flow_summary.iterrows():
        polar_rows.append(
            {"series": "Flow360", "coefficient": panel,
             "alpha_deg": row["alpha_deg"], "value": row[panel]}
        )
    exp_curve = load_reference_curve(f"exp_{panel}.csv")
    for _, row in exp_curve.iterrows():
        polar_rows.append(
            {"series": "Experiment", "coefficient": panel,
             "alpha_deg": float(row["x"]), "value": float(row["y"])}
        )
static_cl_cd = pd.DataFrame(polar_rows)

fig, axes = plt.subplots(1, 2, figsize=(10.5, 5.2), constrained_layout=True)
configure_axes_grid(axes)

for axis, panel in zip(axes, ("CL", "CD")):
    flow_values = static_cl_cd[
        (static_cl_cd["series"] == "Flow360") & (static_cl_cd["coefficient"] == panel)
    ].sort_values("alpha_deg")
    exp_values = static_cl_cd[
        (static_cl_cd["series"] == "Experiment") & (static_cl_cd["coefficient"] == panel)
    ].sort_values("alpha_deg")
    participant_y = plot_participants(axis, participants, panel)
    axis.plot(
        flow_values["alpha_deg"], flow_values["value"],
        color=FLOW360_COLOR, marker="o", linewidth=1.8, markersize=4.5,
        label="RANS Flow360",
    )
    axis.plot(
        exp_values["alpha_deg"], exp_values["value"],
        color=REFERENCE_COLOR, linestyle=":", marker="x", linewidth=1.8,
        markersize=5.0, label="Experiment",
    )
    axis.set_xlim(0, 20)
    axis.set_xticks([0, 5, 10, 15, 20])
    axis.set_ylim(*nice_limits(flow_values["value"], exp_values["value"], participant_y))
    axis.set_xlabel(r"$\alpha\ [deg]$")
    axis.set_ylabel(panel)
    axis.legend()

fig.savefig("results/static_cl_cd_vs_alpha.png", dpi=200, bbox_inches="tight")
plt.show()

Dynamic roll analysis¶

Read the rolling case's surface-forces history and form the roll and yaw increments $\Delta C_l = -\sum C_{Mx}$ and $\Delta C_n = -\sum C_{Mz}$. Reconstruct the physical time from the roll frequency and time step, skip the first few cycles while the wake develops, and reduce the remaining samples into per-phase median curves with 16th/84th-percentile ($\pm\sigma$) bands.

In [ ]:
dynamic_case = dynamic_records[0]["case"]
frame = dynamic_case.results.surface_forces.as_dataframe()
frame.columns = frame.columns.str.strip()
frame = frame.loc[:, ~frame.columns.str.contains("TUBE")]

cl = -sum_coefficients(frame, "CMx")
cn = -sum_coefficients(frame, "CMz")

amplitude_rad = extract_dynamic_amplitude(DYNAMIC_ROLL_EXPRESSION)
omega = 2.0 * np.pi * DYNAMIC_FREQUENCY_HZ
time_step = (np.pi * 1.0) / (omega * 180.0)
period = 1.0 / DYNAMIC_FREQUENCY_HZ

time_series = np.arange(len(cl), dtype=float) * time_step
steps_per_cycle = int(round(period / time_step))
skip_steps = SKIP_DYNAMIC_CYCLES * steps_per_cycle
if skip_steps >= len(time_series):
    raise RuntimeError("Dynamic case does not contain enough samples after skipping.")

time_series = time_series[skip_steps:]
cl = cl[skip_steps:]
cn = cn[skip_steps:]
time_series = time_series - time_series[0]
theta = (omega * time_series) % (2.0 * np.pi)

theta_bins, cl_center, cl_lower, cl_upper = phase_bin_statistics(
    theta, cl, DYNAMIC_PHASE_BINS, 16, 84
)
_, cn_center, cn_lower, cn_upper = phase_bin_statistics(
    theta, cn, DYNAMIC_PHASE_BINS, 16, 84
)

Figure 3: Dynamic rolling time histories¶

Reconstruct two roll cycles from the per-phase curves and plot, top to bottom, the roll angle $\Delta\phi$ and the roll/yaw increments $\Delta C_l$ and $\Delta C_n$ over time, against the wind-tunnel measurements and the participant scatter band.

In [ ]:
recon_points = 720
recon_time = np.linspace(0.0, 2.0 * period, recon_points, endpoint=False)
recon_theta = (omega * recon_time) % (2.0 * np.pi)
recon_phi_deg = np.degrees(amplitude_rad * np.sin(omega * recon_time))

cl_recon = periodic_interp(recon_theta, theta_bins, cl_center)
cl_recon_lower = periodic_interp(recon_theta, theta_bins, cl_lower)
cl_recon_upper = periodic_interp(recon_theta, theta_bins, cl_upper)
cn_recon = periodic_interp(recon_theta, theta_bins, cn_center)
cn_recon_lower = periodic_interp(recon_theta, theta_bins, cn_lower)
cn_recon_upper = periodic_interp(recon_theta, theta_bins, cn_upper)

exp_cl = load_reference_curve("exp_dynamic_Cl.csv")
exp_cn = load_reference_curve("exp_dynamic_Cn.csv")
exp_cl["x"] = exp_cl["x"] - exp_cl["x"].iloc[0]
exp_cn["x"] = exp_cn["x"] - exp_cn["x"].iloc[0]

fig, axes = plt.subplots(
    3, 1, sharex=True, figsize=(11.5, 4.6),
    gridspec_kw={"hspace": 0.08}, constrained_layout=True,
)
configure_axes_grid(axes)

axes[0].plot(recon_time, recon_phi_deg, color=REFERENCE_COLOR, linestyle=":", linewidth=2.2)
axes[0].set_ylabel(r"$\Delta \phi\ [deg]$")
axes[0].set_ylim(-5, 5)
axes[0].set_yticks([-5, 0, 5])

axes[1].fill_between(
    recon_time, cl_recon_lower, cl_recon_upper,
    color=FLOW360_COLOR, alpha=0.18, label=r"Flow360 $\pm\sigma$ bounds",
)
participant_cl = plot_participants(axes[1], participants, "DeltaCl")
axes[1].plot(recon_time, cl_recon, color=FLOW360_COLOR, linewidth=2.0, label="Flow360")
axes[1].plot(
    exp_cl["x"], exp_cl["y"], color=REFERENCE_COLOR, linestyle=":",
    linewidth=1.5, label="Wind Tunnel",
)
axes[1].set_ylabel(r"$\Delta C_l$")
axes[1].set_ylim(*nice_limits(cl_recon_lower, cl_recon_upper, exp_cl["y"], participant_cl))
axes[1].set_yticks([-0.01, 0, 0.01])
axes[1].legend(loc="best", frameon=True)

axes[2].fill_between(
    recon_time, cn_recon_lower, cn_recon_upper,
    color=FLOW360_COLOR, alpha=0.18, label=r"Flow360 $\pm\sigma$ bounds",
)
participant_cn = plot_participants(axes[2], participants, "DeltaCn")
axes[2].plot(recon_time, cn_recon, color=FLOW360_COLOR, linewidth=2.0, label="Flow360")
axes[2].plot(
    exp_cn["x"], exp_cn["y"], color=REFERENCE_COLOR, linestyle=":",
    linewidth=1.5, label="Wind Tunnel",
)
axes[2].set_ylabel(r"$\Delta C_n$")
axes[2].set_xlabel(r"$t\ [s]$")
axes[2].set_ylim(*nice_limits(cn_recon_lower, cn_recon_upper, exp_cn["y"], participant_cn))
axes[2].set_yticks([-0.001, 0, 0.001])
axes[2].set_xlim(recon_time[0], recon_time[-1])
axes[2].legend(loc="best", frameon=True)

fig.savefig("results/dynamic_time_histories_one_sigma.png", dpi=200, bbox_inches="tight")
plt.show()

Figure 4: Dynamic rolling coefficients vs roll phase¶

Plot $\Delta C_l$ and $\Delta C_n$ against the roll angle $\Delta\phi$ over one oscillation to show the hysteresis loops, with the Flow360 $\pm\sigma$ band, the wind-tunnel points, and the tunnel $\pm\sigma$ markers.

In [ ]:
theta_plot = np.linspace(0.0, 2.0 * np.pi, 721)
phi_theta_deg = np.degrees(amplitude_rad * np.sin(theta_plot))
cl_theta = periodic_interp(theta_plot, theta_bins, cl_center)
cl_theta_lower = periodic_interp(theta_plot, theta_bins, cl_lower)
cl_theta_upper = periodic_interp(theta_plot, theta_bins, cl_upper)
cn_theta = periodic_interp(theta_plot, theta_bins, cn_center)
cn_theta_lower = periodic_interp(theta_plot, theta_bins, cn_lower)
cn_theta_upper = periodic_interp(theta_plot, theta_bins, cn_upper)

exp_cl_phi = load_reference_curve("exp_dynamic_Cl_Phi.csv")
exp_cn_phi = load_reference_curve("exp_dynamic_Cn_Phi.csv")
sigma1 = load_reference_curve("sigma1.csv")
sigma2 = load_reference_curve("sigma2.csv")

fig, axes = plt.subplots(
    1, 2, sharex=True, figsize=(10.5, 4.2),
    gridspec_kw={"wspace": 0.18}, constrained_layout=True,
)
configure_axes_grid(axes)

fill_parametric_band(
    axes[0], phi_theta_deg, cl_theta_lower, cl_theta_upper,
    color=FLOW360_COLOR, alpha=0.18, zorder=1, label=r"Flow360 $\pm\sigma$ bounds",
)
plot_participants(axes[0], participants, "DeltaClHysteresis", xy=True)
axes[0].plot(phi_theta_deg, cl_theta, color=FLOW360_COLOR, linewidth=2.4, zorder=3, label="Flow360")
axes[0].plot(
    exp_cl_phi["x"], exp_cl_phi["y"], color=REFERENCE_COLOR, marker="x",
    linestyle="None", markersize=6, zorder=4, label="Wind Tunnel",
)
axes[0].set_ylabel(r"$\Delta C_l$")
axes[0].set_xlabel(r"$\Delta \phi\ [deg]$")
axes[0].set_xlim(-5, 5)
axes[0].set_xticks([-5, 0, 5])
axes[0].legend(loc="best", frameon=True)

fill_parametric_band(
    axes[1], phi_theta_deg, cn_theta_lower, cn_theta_upper,
    color=FLOW360_COLOR, alpha=0.18, zorder=1, label=r"Flow360 $\pm\sigma$ bounds",
)
plot_participants(axes[1], participants, "DeltaCnHysteresis", xy=True)
axes[1].plot(phi_theta_deg, cn_theta, color=FLOW360_COLOR, linewidth=2.4, zorder=3, label="Flow360")
axes[1].plot(
    exp_cn_phi["x"], exp_cn_phi["y"], color=REFERENCE_COLOR, marker="x",
    linestyle="None", markersize=6, zorder=4, label="Wind Tunnel",
)
axes[1].plot(
    sigma1["x"], sigma1["y"], color=REFERENCE_COLOR, marker="o", linestyle="None",
    markersize=5, fillstyle="none", zorder=5, label=r"Wind Tunnel $\pm\sigma$ bounds",
)
axes[1].plot(
    sigma2["x"], sigma2["y"], color=REFERENCE_COLOR, marker="o", linestyle="None",
    markersize=5, fillstyle="none", zorder=5, label="_nolegend_",
)
axes[1].set_ylabel(r"$\Delta C_n$")
axes[1].set_xlabel(r"$\Delta \phi\ [deg]$")
axes[1].set_xlim(-5, 5)
axes[1].set_xticks([-5, 0, 5])
axes[1].set_yticks([-0.004, 0, 0.004])
axes[1].legend(loc="upper left", frameon=True)

fig.savefig("results/dynamic_cl_cn_vs_dphi_one_sigma.png", dpi=200, bbox_inches="tight")
plt.show()