SDM Dynamic: Pitch Damping Derivative¶

This case benchmarks Flow360 on the Standard Dynamics Model (SDM) undergoing a prescribed forced-pitch oscillation, and extracts the dynamic pitch-damping derivative (C_{mq} + C_{m\dot{\alpha}}) across a transonic Mach sweep. The motion is a zero-mean, (\pm 5^{\circ}) sinusoidal pitch at a fixed reduced frequency (k = 0.0493), imposed with a rotating sub-volume so the model physically pitches inside the mesh.

The study sweeps 10 Mach numbers (0.40 to 1.10). Each case is a fully transient, unsteady RANS run using the k-omega SST turbulence model with a quadratic constitutive relation and DDES hybrid modeling, integrated over 5 oscillation cycles at 360 physical steps per cycle. For each case the pitching moment is harmonically analyzed over the settled cycles to recover the damping derivative, which is compared against wind-tunnel, free-flight, and published CFD/Euler reference data.

The notebook runs top-to-bottom in a single kernel: load the geometry, submit the Mach sweep, wait for completion, then post-process the in-kernel case objects into the published pitch-damping figure.

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

Run the cells top to bottom in a single kernel.

Imports¶

In [ ]:
import math
import os
import re
import warnings

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

import flow360 as fl

from flow360.examples import download_benchmark_assets

Helper modules¶

Inlined from the case's pitch_function_generation helper so the notebook is self-contained. pitching_expression builds the Flow360 angle-expression string for the prescribed sinusoidal pitch (using the solver's non-dimensional acoustic time) and returns the matching physical oscillation frequency in Hz, which sets the physical time-step size.

In [ ]:
def pitching_expression(
    amplitude_deg: float = 5.0,
    k: float = 0.0493,
    l_omega: float = 0.2299 / 2,   # [m] length used in k definition (e.g., MAC/2)
    M: float = 0.4,                # [-] Mach number (V/c)
    c: float = 340.074144758659,   # [m/s] speed of sound
    l_time: float = 1.0,           # [m] length used in time nondimensionalization
):
    """
    Returns:
      - exp:    expression string using nondimensional (acoustic) time t*:
                   amplitude_rad * sin( phase_factor * t )
               where t is t* and phase_factor = k*M*(l_time/l_omega)
      - f_real: physical oscillation frequency in Hz

    Definitions used:
      k = l_omega * omega / V
      V = M * c
      => omega = k * M * c / l_omega
      => f_real = omega / (2*pi)

    Nondimensional time convention (acoustic):
      t* = t_phys * c / l_time  ->  t_phys = (l_time/c) * t*
      => omega * t_phys = k*M*(l_time/l_omega) * t*
    """
    amplitude_rad = amplitude_deg * math.pi / 180.0

    # Physical angular frequency and frequency
    omega = (k * M * c) / l_omega              # [rad/s]
    f_real = omega / (2.0 * math.pi)           # [Hz]

    # Phase factor in terms of nondimensional acoustic time t*
    phase_factor = k * M * (l_time / l_omega)  # [-] multiplies t* in sin()

    exp = f"{amplitude_rad}*sin({phase_factor}*t)"
    return exp, f_real

Input data¶

The post-processing compares Flow360 against wind-tunnel, free-flight, and published CFD/Euler reference measurements. 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("SDM_Dynamic", "ref_data")

Load project¶

The root asset for this case is the SDM geometry. Download it from the public benchmark bucket and start a fresh project from it, no private project id is needed. The mesh is generated from this geometry by the meshing parameters set up below.

In [ ]:
root_asset_files = download_benchmark_assets("SDM_Dynamic", "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="SDM Dynamic Pitch Damping Derivative"
)

Configuration¶

Geometry reference dimensions, the surface faces that receive the finer surface refinement, the Mach sweep, and the motion parameters. LIMITER_MACHS marks the transonic/supersonic Mach numbers that additionally enable the pressure/density and velocity limiters for solver robustness. The motion is 5 cycles at 1 degree of pitch phase per physical step (360 steps per cycle).

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

REF_FACE_IDS = [
    "body00001_face00023",
    "body00001_face00024",
    "body00001_face00028",
    "body00001_face00029",
    "body00001_face00030",
    "body00001_face00063",
    "body00001_face00064",
    "body00001_face00069",
    "body00001_face00070",
    "body00001_face00071",
    "body00001_face00072",
    "body00001_face00073",
    "body00001_face00074",
    "body00001_face00075",
    "body00001_face00076",
    "body00001_face00077",
    "body00001_face00078",
    "body00001_face00091",
    "body00001_face00098",
    "body00001_face00099",
    "body00001_face00100",
    "body00001_face00101",
    "body00001_face00102",
    "body00001_face00103",
    "body00001_face00104",
    "body00001_face00105",
    "body00001_face00106",
    "body00001_face00107",
    "body00001_face00108",
    "body00001_face00109",
    "body00001_face00115",
    "body00001_face00117",
    "body00001_face00118",
    "body00001_face00119",
    "body00001_face00120",
    "body00001_face00121",
    "body00001_face00122",
    "body00001_face00123",
    "body00001_face00124",
    "body00001_face00125",
    "body00001_face00148",
    "body00001_face00149",
    "body00001_face00153",
    "body00001_face00154",
    "body00001_face00155",
    "body00001_face00188",
    "body00001_face00189",
    "body00001_face00194",
    "body00001_face00195",
    "body00001_face00196",
    "body00001_face00197",
    "body00001_face00198",
    "body00001_face00199",
    "body00001_face00200",
    "body00001_face00201",
    "body00001_face00202",
    "body00001_face00203",
    "body00001_face00216",
    "body00001_face00223",
    "body00001_face00224",
    "body00001_face00225",
    "body00001_face00226",
    "body00001_face00227",
    "body00001_face00228",
    "body00001_face00229",
    "body00001_face00231",
    "body00001_face00232",
    "body00001_face00233",
    "body00001_face00234",
    "body00001_face00240",
    "body00001_face00242",
    "body00001_face00243",
    "body00001_face00244",
    "body00001_face00245",
    "body00001_face00246",
    "body00001_face00247",
    "body00001_face00248",
    "body00001_face00249",
    "body00001_face00250",
    "body00001_face00002",
    "body00001_face00127",
]
MACH_VALUES = [0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.94, 0.96, 1.0, 1.1]
LIMITER_MACHS = {0.8, 0.96, 1.0, 1.1}

SPAN_M = 0.6096
LENGTH_M = 0.9429
MEAN_AERO_CHORD_M = 0.2299
REFERENCE_AREA_M2 = 0.1238
CENTER_OF_MOTION_X_M = LENGTH_M * 0.55
ROTATION_VOLUME_SPACING_MM = 25.0
ANGLE_PER_STEP_DEG = 1.0
CYCLES = 5

Meshing setup¶

Build the meshing parameters. The automated farfield encloses a rotating sub-volume (the pitching cylinder) around the model. Surface refinement is applied to the listed body faces, and a uniform refinement matches the rotating-interface spacing so the sliding interface is well resolved.

In [ ]:
def build_meshing(farfield, rotating_volume, cylinder, geometry):
    with fl.SI_unit_system:
        return fl.MeshingParams(
            defaults=fl.MeshingDefaults(
                surface_max_edge_length=12.5 * fl.u.mm,
                boundary_layer_first_layer_thickness=2e-6 * fl.u.m,
                curvature_resolution_angle=5 * fl.u.deg,
            ),
            volume_zones=[farfield, rotating_volume],
            refinements=[
                fl.SurfaceRefinement(
                    faces=[geometry[face] for face in REF_FACE_IDS],
                    max_edge_length=6 * fl.u.mm,
                    curvature_resolution_angle=1 * fl.u.deg,
                ),
                fl.UniformRefinement(
                    entities=cylinder,
                    spacing=ROTATION_VOLUME_SPACING_MM * fl.u.mm,
                ),
            ],
        )

Physics setup¶

Build the boundary conditions and flow solvers. The farfield is a freestream, the model surfaces are viscous walls, and the fluid uses the k-omega SST turbulence model with the quadratic constitutive relation and a DDES hybrid model. The prescribed pitch is imposed as a Rotation on the cylinder sub-volume using the angle expression built by pitching_expression.

In [ ]:
def build_models(cylinder, farfield, geometry, pitch_expression_text):
    with fl.SI_unit_system:
        return [
            fl.Freestream(surfaces=farfield.farfield),
            fl.Wall(surfaces=geometry["*"]),
            fl.Fluid(
                navier_stokes_solver=fl.NavierStokesSolver(
                    absolute_tolerance=1e-10,
                    relative_tolerance=1e-2,
                ),
                turbulence_model_solver=fl.KOmegaSST(
                    absolute_tolerance=1e-8,
                    equation_evaluation_frequency=1,
                    relative_tolerance=1e-2,
                    quadratic_constitutive_relation=True,
                    hybrid_model=fl.DetachedEddySimulation(
                        shielding_function="DDES"
                    ),
                ),
            ),
            fl.Rotation(
                name="Pitching Function",
                volumes=cylinder,
                spec=fl.AngleExpression(pitch_expression_text),
            ),
        ]

Simulation Params¶

Assemble the full SimulationParams for one Mach number: the meshing, the reference geometry (area, moment center at the pitch axis, and moment lengths), the Mach operating condition, the unsteady time stepping (adaptive CFL, prescribed step size and step count), and the physics models with the pitch motion.

In [ ]:
def create_params(project, mach, cylinder, step_size_s, steps, pitch_expression_text):
    geometry = project.geometry
    with fl.SI_unit_system:
        farfield = fl.AutomatedFarfield()
        rotating_volume = fl.RotationVolume(
            spacing_axial=ROTATION_VOLUME_SPACING_MM * fl.u.mm,
            spacing_circumferential=ROTATION_VOLUME_SPACING_MM * fl.u.mm,
            spacing_radial=ROTATION_VOLUME_SPACING_MM * fl.u.mm,
            entities=cylinder,
            enclosed_entities=geometry["*"],
        )
        return fl.SimulationParams(
            meshing=build_meshing(farfield, rotating_volume, cylinder, geometry),
            reference_geometry=fl.ReferenceGeometry(
                area=REFERENCE_AREA_M2 * fl.u.m**2,
                moment_center=(CENTER_OF_MOTION_X_M, 0.0, 0.0) * fl.u.m,
                moment_length=(SPAN_M, MEAN_AERO_CHORD_M, SPAN_M) * fl.u.m,
            ),
            operating_condition=fl.AerospaceCondition.from_mach(mach=mach),
            time_stepping=fl.Unsteady(
                CFL=fl.AdaptiveCFL(convergence_limiting_factor=0.25),
                step_size=step_size_s * fl.u.s,
                steps=steps,
                max_pseudo_steps=30,
            ),
            models=build_models(cylinder, farfield, geometry, pitch_expression_text),
        )

Submit cases¶

Sweep every Mach number and submit its transient forced-pitch case to Flow360. For each Mach the pitch angle expression and its physical frequency are derived, the physical time-step size is set so one physical step advances the pitch phase by 1 degree, and the run covers 5 full cycles. The transonic/supersonic Mach numbers in LIMITER_MACHS additionally enable the Navier-Stokes limiters. Each submitted case object is kept in submitted_cases so the post-processing can read its results directly from this kernel.

In [ ]:
cylinder = fl.Cylinder(
    name="Cylinder",
    center=(CENTER_OF_MOTION_X_M, 0.0, 0.0) * fl.u.m,
    axis=(0, 1, 0),
    height=2.0 * SPAN_M * fl.u.m,
    outer_radius=2.0 * LENGTH_M * fl.u.m,
)

submitted_cases = []
steps = int(360 * CYCLES / ANGLE_PER_STEP_DEG)
for mach in MACH_VALUES:
    pitch_expression_text, frequency_hz = pitching_expression(
        M=mach, l_omega=MEAN_AERO_CHORD_M / 2.0
    )
    step_size_s = ANGLE_PER_STEP_DEG / (360.0 * frequency_hz)
    params = create_params(
        project=project,
        mach=mach,
        cylinder=cylinder,
        step_size_s=step_size_s,
        steps=steps,
        pitch_expression_text=pitch_expression_text,
    )
    if mach in LIMITER_MACHS:
        params.models[2].navier_stokes_solver.limit_pressure_density = True
        params.models[2].navier_stokes_solver.limit_velocity = True
    case = project.run_case(
        name=f"m{mach:.2f}_pitch_ddes_sstqcr_{SOLVER_VERSION}",
        params=params,
        solver_version=SOLVER_VERSION,
        use_beta_mesher=True,
    )
    submitted_cases.append(case)

Wait for completion¶

Block until every submitted case has finished. This can take a long time (minutes to hours), since each case is a fully transient solve over 5 oscillation cycles.

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

Postprocessing¶

The published result is a single figure: the pitch-damping derivative (C_{mq} + C_{m\dot{\alpha}}) versus Mach, with Flow360 compared against wind-tunnel, free-flight, and published CFD/Euler references.

For each case the pitching-moment history (C_{My}) is reduced per physical step, and the settled oscillation cycles (skipping the first two of five) are analyzed harmonically: the in-phase and quadrature components of the moment relative to the prescribed motion give (C_{m\dot{\alpha}}) and (C_{mq} + C_{m\dot{\alpha}}). The cells below read the results straight from the in-kernel case objects (no project is re-opened).

Setup¶

Styling color, the harmonic-analysis window (5 cycles total, first 2 skipped), the motion parameters used in the reduction (5 deg amplitude, reduced frequency 0.0493), and the reference-data file names.

In [ ]:
os.makedirs("results", exist_ok=True)

FLOW360_COLOR = "#00643c"   # Flow360 results
WIND_TUNNEL_COLOR = "black"
FREE_FLIGHT_COLOR = "#ff7f0e"

REF_DATA_DIR = "ref_data"

# --- harmonic-analysis window and motion parameters ---
POINTS_PER_CYCLE = 360
TOTAL_CYCLES = 5
SKIP_CYCLES = 2
MIN_CYCLES_REQUIRED = 3
AMPLITUDE_DEG = 5.0
REDUCED_FREQUENCY = 0.0493
PHASE_OFFSET_DEG = 0.0
DT_CV_WARN = 0.02

# --- reference data ---
WIND_TUNNEL_FILE = "WindTunnel_AEDC_3.csv"
RANGE_FILE = "RangeData_AEDC_3.csv"
ADDITIONAL_FILES = ("PMB-C.csv", "da_ronch_fine.csv", "URANS_Tatar.csv")

MACH_RE = re.compile(r"\bm\s*=?\s*(?P<mach>\d+(?:\.\d+)?)", re.IGNORECASE)
RELEASE_SUFFIX_RE = re.compile(r"_release[-\d.]+$", re.IGNORECASE)

Reduction routines¶

Helpers that turn each case's total-forces history into the harmonic pitch-damping coefficients, and that load the reference curves. extract_cm_by_step pulls the per-physical-step (C_{My}); compute_harmonic_cycles integrates the moment against (\sin\theta)/(\cos\theta) over each settled cycle (periodic cubic Lagrange quadrature) to recover the derivatives and their cycle-to-cycle spread.

In [ ]:
def strip_release_suffix(name):
    return RELEASE_SUFFIX_RE.sub("", name)


def parse_mach(case_name):
    match = MACH_RE.search(case_name)
    if not match:
        raise RuntimeError(f"Could not parse Mach number from case name: {case_name}")
    return float(match.group("mach"))


def additional_dataset_label(filename):
    labels = {
        "PMB-C.csv": "Da Ronch Euler 0.7M",
        "da_ronch_fine.csv": "Da Ronch Euler 2.4M",
        "URANS_Tatar.csv": "Tatar URANS",
    }
    return labels.get(filename, os.path.splitext(filename)[0])


def marker_and_color(index):
    markers = ("D", "^", "P")
    colors = ("#d62728", "#9467bd", "#8c564b")
    return markers[index % len(markers)], colors[index % len(colors)]


def extract_total_coeff(df, coeff):
    columns = [str(column).strip() for column in df.columns]
    coeff_lower = coeff.lower()
    for column in columns:
        if column == coeff or column.lower() == coeff_lower:
            return df[column].to_numpy(dtype=float)
    matches = []
    for column in columns:
        column_lower = column.lower()
        if column_lower.endswith(f"_{coeff_lower}") or column_lower.endswith(f".{coeff_lower}"):
            matches.append(column)
            continue
        tokens = [token for token in re.split(r"[^a-z0-9]+", column_lower) if token]
        if tokens and tokens[-1] == coeff_lower:
            matches.append(column)
    if not matches:
        raise RuntimeError(
            f"Could not find '{coeff}' in total_forces columns. Available: {columns}"
        )
    if len(matches) > 1:
        total_matches = [column for column in matches if "total" in column.lower()]
        if len(total_matches) == 1:
            return df[total_matches[0]].to_numpy(dtype=float)
        raise RuntimeError(f"Ambiguous columns for '{coeff}' in total_forces: {matches}")
    return df[matches[0]].to_numpy(dtype=float)


def find_step_column(df, names):
    preferred = set(names)
    for column in df.columns:
        normalized = re.sub(r"[^a-z0-9]+", "_", str(column).strip().lower()).strip("_")
        if normalized in preferred:
            return str(column)
    return None


def extract_cm_by_step(df):
    cm = extract_total_coeff(df, "CMy")
    step_col = find_step_column(
        df, {"physical_step", "physicalstep", "step", "step_index", "physical_step_index"}
    )
    pseudo_col = find_step_column(df, {"pseudo_step", "pseudostep", "pseudo", "substep"})
    time_col = find_step_column(df, {"physical_time", "time", "simulation_time"})

    if step_col is None:
        return pd.DataFrame({"step": np.arange(len(df), dtype=int), "cm": cm, "time": np.nan})

    step_values = pd.to_numeric(df[step_col], errors="coerce")
    pseudo_values = (
        pd.to_numeric(df[pseudo_col], errors="coerce")
        if pseudo_col is not None
        else pd.Series(np.arange(len(df), dtype=float), index=df.index)
    )
    time_values = (
        pd.to_numeric(df[time_col], errors="coerce")
        if time_col is not None
        else pd.Series(np.nan, index=df.index)
    )
    table = pd.DataFrame(
        {
            "step": step_values,
            "pseudo": pseudo_values,
            "cm": pd.to_numeric(cm, errors="coerce"),
            "time": time_values,
        }
    ).dropna(subset=["step", "pseudo", "cm"])
    if table.empty:
        raise RuntimeError("No valid numeric rows found for step/pseudo/cm.")
    table["step"] = np.rint(table["step"]).astype(int)
    table["pseudo"] = np.rint(table["pseudo"]).astype(int)
    return table.sort_values(["step", "pseudo"]).groupby("step", as_index=False).last()


def warn_if_nonuniform_dt(table, case_name):
    if "time" not in table.columns:
        return
    time_values = pd.to_numeric(table["time"], errors="coerce").to_numpy(dtype=float)
    time_values = time_values[np.isfinite(time_values)]
    if time_values.size < 4:
        return
    dt = np.diff(time_values)
    dt = dt[dt > 0.0]
    if dt.size < 3:
        return
    cv = float(np.std(dt) / np.mean(dt))
    if cv > DT_CV_WARN:
        warnings.warn(
            f"Non-uniform timestep detected for '{case_name}' "
            f"(std(dt)/mean(dt)={cv:.2%})."
        )


def build_cycle_time_axis(cycle_table):
    if "time" in cycle_table.columns:
        time_values = pd.to_numeric(cycle_table["time"], errors="coerce").to_numpy(dtype=float)
        if len(time_values) == POINTS_PER_CYCLE and np.all(np.isfinite(time_values)):
            dt = np.diff(time_values)
            if np.all(dt > 0.0):
                dt_med = float(np.median(dt))
                period = float((time_values[-1] - time_values[0]) + dt_med)
                if period > 0.0:
                    return time_values, period
    synthetic = np.arange(POINTS_PER_CYCLE, dtype=float) / float(POINTS_PER_CYCLE)
    return synthetic, 1.0


def integrate_periodic_cubic_lagrange(x, y, period):
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if x.size != y.size or x.size < 4:
        raise RuntimeError("Periodic cubic Lagrange integration requires matching 1D arrays.")
    if not np.all(np.diff(x) > 0.0):
        raise RuntimeError("x must be strictly increasing for cubic Lagrange integration.")

    integral = 0.0
    for i in range(x.size):
        a = float(x[i])
        b = float(x[i + 1]) if i + 1 < x.size else float(x[0] + period)
        nominal = [i - 1, i, i + 1, i + 2]
        x_nodes = []
        y_nodes = []
        for node in nominal:
            index = node % x.size
            shift = node // x.size
            x_nodes.append(float(x[index] + shift * period))
            y_nodes.append(float(y[index]))
        x_nodes = np.array(x_nodes, dtype=float)
        y_nodes = np.array(y_nodes, dtype=float)
        interval = 0.0
        for m in range(4):
            others = [x_nodes[q] for q in range(4) if q != m]
            basis = np.poly(others)
            denom = 1.0
            for q in range(4):
                if q != m:
                    denom *= x_nodes[m] - x_nodes[q]
            weights = np.polyint(basis / denom)
            interval += float(np.polyval(weights, b) - np.polyval(weights, a)) * y_nodes[m]
        integral += interval
    return float(integral)


def compute_harmonic_cycles(table):
    amplitude_rad = math.radians(AMPLITUDE_DEG)
    step_rel = table["step"].to_numpy(dtype=int) - int(table["step"].min())
    cycle_id = step_rel // POINTS_PER_CYCLE
    phase_idx = step_rel % POINTS_PER_CYCLE
    keep = (cycle_id >= SKIP_CYCLES) & (cycle_id < TOTAL_CYCLES)
    if not np.any(keep):
        raise RuntimeError("No samples left after the configured cycle window.")

    work = pd.DataFrame(
        {
            "cycle_id": cycle_id[keep],
            "phase_idx": phase_idx[keep],
            "cm": table["cm"].to_numpy(dtype=float)[keep],
            "time": table["time"].to_numpy(dtype=float)[keep]
            if "time" in table.columns
            else np.nan,
        }
    )

    phase_offset = math.radians(PHASE_OFFSET_DEG)
    complete_cycles = []
    for cycle, group in work.groupby("cycle_id", sort=True):
        group = group.sort_values("phase_idx")
        if group["phase_idx"].nunique() == POINTS_PER_CYCLE:
            complete_cycles.append((int(cycle), group))
    if len(complete_cycles) < MIN_CYCLES_REQUIRED:
        raise RuntimeError(
            f"Not enough full cycles for damping estimate "
            f"(need >= {MIN_CYCLES_REQUIRED}, got {len(complete_cycles)})."
        )

    cm_retained = np.concatenate(
        [group["cm"].to_numpy(dtype=float) for _, group in complete_cycles]
    )
    cm_global_mean = float(np.mean(cm_retained))
    cycle_values = []
    integral_sin_total = 0.0
    integral_cos_total = 0.0
    total_period = 0.0

    for cycle, group in complete_cycles:
        cm_cycle = group["cm"].to_numpy(dtype=float)
        phase = group["phase_idx"].to_numpy(dtype=float)
        time_axis, period = build_cycle_time_axis(group)
        theta = (2.0 * np.pi * phase / float(POINTS_PER_CYCLE)) + phase_offset
        delta_cm = cm_cycle - cm_global_mean
        integral_sin = integrate_periodic_cubic_lagrange(
            time_axis, delta_cm * np.sin(theta), period
        )
        integral_cos = integrate_periodic_cubic_lagrange(
            time_axis, delta_cm * np.cos(theta), period
        )
        cm_alpha_bar = (2.0 / (amplitude_rad * period)) * integral_sin
        cmq_plus_cmadot = (2.0 / (REDUCED_FREQUENCY * amplitude_rad * period)) * integral_cos
        integral_sin_total += integral_sin
        integral_cos_total += integral_cos
        total_period += period

        cycle_values.append(
            {
                "cycle_id": cycle,
                "period": float(period),
                "cm_mean_cycle": float(np.mean(cm_cycle)),
                "cm_alpha_bar_cycle": float(cm_alpha_bar),
                "cmq_plus_cmadot_cycle": float(cmq_plus_cmadot),
            }
        )
    if total_period <= 0.0:
        raise RuntimeError("Total retained period is non-positive.")

    return cycle_values, {
        "cm_alpha_bar_global": float(
            (2.0 / (amplitude_rad * total_period)) * integral_sin_total
        ),
        "cmq_plus_cmadot_global": float(
            (2.0 / (REDUCED_FREQUENCY * amplitude_rad * total_period))
            * integral_cos_total
        ),
    }


def process_case(case):
    case_name = case.name
    case.wait()
    df = case.results.total_forces.as_dataframe()
    df.columns = df.columns.str.strip()
    table = extract_cm_by_step(df)
    warn_if_nonuniform_dt(table, case_name)
    cycle_values, global_coeffs = compute_harmonic_cycles(table)

    cycle_mq = np.array(
        [row["cmq_plus_cmadot_cycle"] for row in cycle_values], dtype=float
    )
    return {
        "case_id": case.id,
        "case_name": case_name,
        "mach": parse_mach(strip_release_suffix(case_name)),
        "cmq_plus_cmadot_global": float(global_coeffs["cmq_plus_cmadot_global"]),
        "cmq_plus_cmadot_cycle_std": float(np.std(cycle_mq, ddof=1))
        if len(cycle_mq) > 1
        else 0.0,
    }


def load_reference_points(filename):
    path = os.path.join(REF_DATA_DIR, filename)
    if not os.path.isfile(path):
        warnings.warn(f"Reference data file not found: {path}")
        return None
    df = pd.read_csv(path, header=None)
    if df.shape[1] < 2:
        warnings.warn(f"Reference data file must have at least 2 columns: {path}")
        return None
    out = pd.DataFrame(
        {
            "mach": pd.to_numeric(df.iloc[:, 0], errors="coerce"),
            "value": pd.to_numeric(df.iloc[:, 1], errors="coerce"),
        }
    ).dropna(subset=["mach", "value"])
    return out if not out.empty else None

Pitch damping derivative versus Mach¶

Reduce every case to its global (C_{mq} + C_{m\dot{\alpha}}) (with the cycle-to-cycle standard deviation as an error bar), load the reference datasets, and build the published figure comparing Flow360 against the wind-tunnel, free-flight, and CFD/Euler references.

In [ ]:
results = [process_case(case) for case in submitted_cases]

reference_sets = {
    "Wind Tunnel": load_reference_points(WIND_TUNNEL_FILE),
    "Free Flight": load_reference_points(RANGE_FILE),
}
for filename in ADDITIONAL_FILES:
    reference_sets[filename] = load_reference_points(filename)

ordered = sorted(results, key=lambda row: row["mach"])
x = [row["mach"] for row in ordered]
y = [row["cmq_plus_cmadot_global"] for row in ordered]
yerr = [row["cmq_plus_cmadot_cycle_std"] for row in ordered]

fig, ax = plt.subplots(figsize=(8.8, 4.8), constrained_layout=True)
ax.minorticks_on()
ax.grid(True, which="major", linewidth=0.8, alpha=0.7)
ax.grid(True, which="minor", linestyle=":", linewidth=0.5, alpha=0.5)

ax.errorbar(
    x, y, yerr=yerr,
    fmt="o-", capsize=4, linewidth=2.0, markersize=6,
    color=FLOW360_COLOR, label="Flow360 k-omega SST DDES",
)

if reference_sets["Wind Tunnel"] is not None:
    ref = reference_sets["Wind Tunnel"]
    ax.plot(
        ref["mach"].to_numpy(dtype=float), ref["value"].to_numpy(dtype=float),
        linestyle="None", marker="x", markersize=7,
        color=WIND_TUNNEL_COLOR, label="Wind Tunnel AEDC", zorder=5,
    )
if reference_sets["Free Flight"] is not None:
    ref = reference_sets["Free Flight"]
    ax.plot(
        ref["mach"].to_numpy(dtype=float), ref["value"].to_numpy(dtype=float),
        linestyle="None", marker="s", markersize=6,
        markerfacecolor="none", markeredgewidth=1.3,
        color=FREE_FLIGHT_COLOR, label="Free Flight", zorder=5,
    )
for index, filename in enumerate(ADDITIONAL_FILES):
    ref = reference_sets.get(filename)
    if ref is None:
        continue
    marker, color = marker_and_color(index)
    ax.plot(
        ref["mach"].to_numpy(dtype=float), ref["value"].to_numpy(dtype=float),
        linestyle="None", marker=marker, markersize=6,
        color=color, label=additional_dataset_label(filename), zorder=6,
    )

ax.set_xlabel("Mach")
ax.set_ylabel(r"$C_{mq} + C_{m\dot{\alpha}}$")
ax.set_xlim(0.35, 1.2)
ax.set_ylim(-20.0, 0.0)
ax.legend(loc="best", frameon=True)
fig.savefig("results/pitch_damping_summary.png", dpi=220, bbox_inches="tight")
plt.show()