Finner Missile Pitch Damping Derivative: Transient Supersonic Analysis¶

This case reproduces Flow360's prediction of the dynamic pitch-damping derivative for the Army-Navy Basic Finner missile, a long-standing benchmark for supersonic dynamic stability that is backed by a wealth of experimental data. The flow is a forced small-amplitude pitching oscillation at Mach 1.96 and Reynolds number 187,000 (based on diameter), with a reduced frequency of 0.16 and a 0.32° amplitude, solved unsteady with SA-ZDES. The study sweeps 14 mean angles of attack from -10° to 85°, spanning attached supersonic flow through to separated, shock-dominated blunt-body behavior. At each angle a steady SA case is run first and used to fork the unsteady forced-pitch case(s); the pitching-moment history is then reduced to the mean pitching moment $C_m$ and the pitch-damping sum $C_{m_q}+C_{m_{\dot{\alpha}}}$, compared against AEDC experiment and published CFD (AIAA 2012-0691, FUN3D, KESTREL).

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 os
import math
import re
import warnings
from contextlib import nullcontext

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 case's forced-pitch motion is prescribed by an angle expression generated by a small local helper. It is inlined here so the notebook is self-contained. Given the Mach number, reduced frequency, pitch amplitude and the acoustic nondimensionalization lengths, it returns both the Flow360 angle-expression string (in nondimensional acoustic time t) and the physical oscillation frequency in Hz used to size the unsteady time step.

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
):
    """Return (angle-expression string, physical frequency [Hz]).

    k = l_omega * omega / V, with V = M * c, so omega = k*M*c/l_omega and
    f_real = omega / (2*pi). In nondimensional acoustic time t* = t_phys*c/l_time
    the phase becomes k*M*(l_time/l_omega)*t*.
    """
    amplitude_rad = amplitude_deg * math.pi / 180.0
    omega = (k * M * c) / l_omega              # [rad/s]
    f_real = omega / (2.0 * math.pi)           # [Hz]
    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 AEDC experimental data and published CFD (AIAA 2012-0691, FUN3D, KESTREL). Fetch those reference curves from the public benchmark bucket; this recreates a local ./ref_data/ directory that the plotting cells read from.

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

Load project¶

The root asset for this case is the missile geometry. Download it from the public benchmark bucket, start a fresh project from it, and group the surfaces by their faceName tag so the meshing refinements and boundary conditions can address them by name.

In [ ]:
root_asset_files = download_benchmark_assets("Finner_Missile_Pitching", "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="Finner Missile Pitch Damping Derivative"
)
project.geometry.group_faces_by_tag("faceName")

Case constants¶

Operating point, geometry references, and the sweep definition. The mean angle-of-attack sweep ALPHA_VALUES drives the whole study, and the tagging / time-stepping constants encode the high-angle treatment (first-order startup and a doubled temporal resolution above 70°/80°).

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

GAMMA = 1.4
GAS_CONSTANT = 287.05

DIAMETER_IN = 1.25
DIAMETER = DIAMETER_IN * fl.u.inch
MACH = 1.96
REYNOLDS_DIAMETER = 1.87e5
TOTAL_TEMPERATURE = 522 * fl.u.R
ALPHA_VALUES = [-10, -5, 0, 5, 10, 20, 30, 40, 45, 50, 60, 70, 80, 85]
CYLINDER_AXIS = (0, 1, 0)

CENTER_OF_MOTION_X_IN = 6.1 * DIAMETER_IN
LENGTH_IN = 10.0 * DIAMETER_IN
REFERENCE_AREA = np.pi * (DIAMETER / 2) ** 2
ROTATION_VOLUME_SPACING_MM = 25.0
PHYSICAL_STEPS_PER_CYCLE = 360
HIGH_ALPHA_ZDES_PHYSICAL_STEPS_PER_CYCLE = 2 * PHYSICAL_STEPS_PER_CYCLE
CYCLES = 3
HIGH_ALPHA_URANS_CYCLES = 1
HIGH_ALPHA_ZDES_CYCLES = 2
PITCH_AMPLITUDE_DEG = 0.32
REDUCED_FREQUENCY = 0.16
TAG_SPLIT_ALPHA = 70.0
ZDES_DISSIPATION_FACTORS = {80: 1.0, 85: 0.2}

Meshing setup¶

build_meshing assembles the beta-mesher parameters: a fine anisotropic boundary layer (40 layers, 1.5e-6 m first layer, halved on the nose and leading edges), the automated farfield plus a RotationVolume around the pitching cylinder, near-field and interface box refinements, and mesh slice outputs on the x/y/z planes.

In [ ]:
def build_meshing(farfield, rotating_volume, cylinder, box, slice_x, slice_y, slice_z, geometry):
    with fl.SI_unit_system:
        return fl.MeshingParams(
            outputs=[
                fl.MeshSliceOutput(slices=[slice_x, slice_y, slice_z]),
            ],
            defaults=fl.MeshingDefaults(
                surface_max_edge_length=0.025 * fl.u.inch,
                boundary_layer_first_layer_thickness=1.5e-6 * fl.u.m,
                curvature_resolution_angle=5 * fl.u.deg,
                boundary_layer_growth_rate=1.125,
                number_of_boundary_layers=40,
            ),
            volume_zones=[farfield, rotating_volume],
            refinements=[
                fl.UniformRefinement(
                    name="Interface",
                    entities=cylinder,
                    spacing=ROTATION_VOLUME_SPACING_MM * fl.u.mm,
                ),
                fl.UniformRefinement(
                    name="Near-Field",
                    entities=box,
                    spacing=6.25 * fl.u.mm,
                ),
                fl.SurfaceRefinement(
                    name="LE",
                    faces=[geometry["LE*"]],
                    curvature_resolution_angle=2 * fl.u.deg,
                    max_edge_length=0.001 * 15 * fl.u.inch,
                ),
                fl.SurfaceRefinement(
                    name="Fin",
                    faces=geometry["fin"],
                    curvature_resolution_angle=2 * fl.u.deg,
                    max_edge_length=0.01 * 2 * fl.u.inch,
                ),
                fl.BoundaryLayer(
                    faces=[
                        geometry["*Nose"],
                        geometry["LE*"],
                    ],
                    first_layer_thickness=0.45 * 1.5e-6 * fl.u.m,
                ),
            ],
        )

Physics setup¶

The flow physics: a freestream / viscous-wall boundary pair and a compressible Navier-Stokes solver with the Spalart-Allmaras turbulence model (rotation correction on). The small thermodynamic helpers derive the static temperature and speed of sound from the total-temperature operating point; unsteady_time_stepping sizes the physical time step from the pitching frequency; and enable_zdes switches the turbulence model into hybrid ZDES mode (with a tuned Roe dissipation factor at the highest angles). warning_bypass suppresses the length-scale-mismatch advisory that is expected for this geometry.

In [ ]:
def warning_bypass(warning_id):
    bypass = getattr(fl, "warning_bypass", None)
    return bypass(warning_id) if bypass else nullcontext()


def static_temperature():
    return TOTAL_TEMPERATURE / (1 + ((GAMMA - 1) / 2) * MACH**2)


def speed_of_sound_mps():
    return float(np.sqrt(GAMMA * GAS_CONSTANT * static_temperature().to("K")))


def unsteady_time_stepping(cycles, frequency_hz, convergence_limiter, physical_steps_per_cycle=PHYSICAL_STEPS_PER_CYCLE):
    return fl.Unsteady(
        CFL=fl.AdaptiveCFL(convergence_limiting_factor=convergence_limiter),
        step_size=(1.0 / (physical_steps_per_cycle * frequency_hz)) * fl.u.s,
        steps=physical_steps_per_cycle * cycles,
        max_pseudo_steps=30,
    )


def enable_zdes(params, alpha):
    fluid = params.models[2]
    fluid.turbulence_model_solver.hybrid_model = fl.DetachedEddySimulation(
        shielding_function="ZDES",
        grid_size_for_LES="maxEdgeLength",
    )
    if alpha in ZDES_DISSIPATION_FACTORS:
        fluid.navier_stokes_solver.riemann_solver = fl.RoeFlux(
            numerical_dissipation_factor=ZDES_DISSIPATION_FACTORS[alpha]
        )


def build_models(farfield, geometry):
    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,
                    limit_pressure_density=True,
                    limit_velocity=True,
                ),
                turbulence_model_solver=fl.SpalartAllmaras(
                    absolute_tolerance=1e-8,
                    equation_evaluation_frequency=1,
                    relative_tolerance=1e-2,
                    rotation_correction=True,
                ),
            ),
        ]

Simulation Params¶

create_params assembles the full SimulationParams for one mean angle of attack: the mesh slices, farfield and rotation volume, the reference geometry (moment center at the missile's center of motion), the Mach/Reynolds operating condition at the requested alpha, and the meshing/physics blocks defined above.

In [ ]:
def create_params(project, cylinder, box, alpha):
    geometry = project.geometry
    with fl.SI_unit_system:
        slice_y = fl.Slice(origin=(0, 0, 0), normal=(0, 1, 0), name="y_slice")
        slice_z = fl.Slice(origin=(0, 0, 0), normal=(0, 0, 1), name="z_slice")
        slice_x = fl.Slice(
            origin=(11.5, 0, 0) * fl.u.inch,
            normal=(1, 0, 0),
            name="x_slice",
        )
        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["*"],
        )
        with warning_bypass("potential_length_scale_mismatch"):
            return fl.SimulationParams(
                meshing=build_meshing(
                    farfield, rotating_volume, cylinder, box, slice_x, slice_y, slice_z, geometry
                ),
                reference_geometry=fl.ReferenceGeometry(
                    area=REFERENCE_AREA,
                    moment_center=(CENTER_OF_MOTION_X_IN, 0.0, 0.0) * fl.u.inch,
                    moment_length=(LENGTH_IN, DIAMETER_IN, LENGTH_IN) * fl.u.inch,
                ),
                operating_condition=fl.AerospaceCondition.from_mach_reynolds(
                    mach=MACH,
                    reynolds_mesh_unit=REYNOLDS_DIAMETER / DIAMETER_IN,
                    project_length_unit=1 * fl.u.inch,
                    alpha=alpha * fl.u.deg,
                    temperature=static_temperature(),
                ),
                models=build_models(farfield, geometry),
            )

Submit cases¶

This is the full study. For every mean angle of attack the notebook builds one RotationVolume cylinder and near-field box, generates the pitching angle expression, then submits:

  1. a steady SA startup case (first-order and more conservative above 70°),
  2. the unsteady forced-pitch case forked from that steady solution, switched to second-order SA-ZDES; above 80° an intermediate URANS case is inserted and the ZDES case uses a doubled temporal resolution over two cycles.

All submitted cases are collected in submitted_cases. The launches are non-blocking; completion is awaited in the next section.

In [ ]:
cylinder = fl.Cylinder(
    name="Cylinder",
    center=(CENTER_OF_MOTION_X_IN, 0, 0) * fl.u.inch,
    axis=CYLINDER_AXIS,
    height=LENGTH_IN * fl.u.inch,
    outer_radius=(2.5 * LENGTH_IN) * fl.u.inch,
)
box = fl.Box(
    name="Volume Ref",
    center=(6.625, 0, 5) * fl.u.inch,
    size=(14.5, 6.25, 20) * fl.u.inch,
)

pitch_expression_text, frequency_hz = pitching_expression(
    M=MACH,
    l_omega=float(DIAMETER.to("m") / 2),
    k=REDUCED_FREQUENCY,
    amplitude_deg=PITCH_AMPLITUDE_DEG,
    l_time=0.0254,
    c=speed_of_sound_mps(),
)


def slug_number(value):
    return f"{value:g}".replace(".", "p").replace("-", "m")


def build_case_name(alpha, label):
    return f"alpha{slug_number(alpha)}_{label}_{SOLVER_VERSION}"


def run_tag(alpha):
    return "first_run_better" if alpha > TAG_SPLIT_ALPHA else "first_run"


def run_case(project, **kwargs):
    with warning_bypass("potential_length_scale_mismatch"):
        return project.run_case(**kwargs)


submitted_cases = []
for alpha in ALPHA_VALUES:
    params = create_params(project=project, cylinder=cylinder, box=box, alpha=alpha)
    params.time_stepping = fl.Steady(
        max_steps=3000,
        CFL=fl.AdaptiveCFL(convergence_limiting_factor=0.5),
    )
    convergence_limiter = 0.75
    if alpha > TAG_SPLIT_ALPHA:
        params.models[2].turbulence_model_solver.order_of_accuracy = 1
        params.models[2].navier_stokes_solver.order_of_accuracy = 1
        convergence_limiter = 0.15
    tag = run_tag(alpha)

    steady = run_case(
        project,
        name=build_case_name(
            alpha,
            "steady_sa_first_order" if alpha > TAG_SPLIT_ALPHA else "steady_sa",
        ),
        params=params,
        use_beta_mesher=True,
        solver_version=SOLVER_VERSION,
        tags=[CASE_TAG, SOLVER_VERSION, "steady", tag],
    )

    params.models[2].turbulence_model_solver.order_of_accuracy = 2
    params.models[2].navier_stokes_solver.order_of_accuracy = 2

    params.models.append(
        fl.Rotation(
            name="Pitching Function",
            volumes=cylinder,
            spec=fl.AngleExpression(pitch_expression_text),
        )
    )
    if alpha >= 80:
        params.time_stepping = unsteady_time_stepping(
            HIGH_ALPHA_URANS_CYCLES, frequency_hz, convergence_limiter
        )
        urans = run_case(
            project,
            name=build_case_name(alpha, "pitch_urans_sa"),
            params=params,
            use_beta_mesher=True,
            solver_version=SOLVER_VERSION,
            tags=[CASE_TAG, SOLVER_VERSION, tag, "pitch_urans_sa"],
            fork_from=steady,
        )
        enable_zdes(params, alpha)
        params.time_stepping = unsteady_time_stepping(
            HIGH_ALPHA_ZDES_CYCLES,
            frequency_hz,
            convergence_limiter,
            HIGH_ALPHA_ZDES_PHYSICAL_STEPS_PER_CYCLE,
        )
        case = run_case(
            project,
            name=build_case_name(alpha, "pitch_zdes_sa"),
            params=params,
            use_beta_mesher=True,
            solver_version=SOLVER_VERSION,
            tags=[CASE_TAG, SOLVER_VERSION, tag, "pitch_zdes_sa"],
            fork_from=urans,
        )
    else:
        params.time_stepping = unsteady_time_stepping(
            CYCLES, frequency_hz, convergence_limiter
        )
        enable_zdes(params, alpha)
        case = run_case(
            project,
            name=build_case_name(alpha, "pitch_zdes_sa"),
            params=params,
            use_beta_mesher=True,
            solver_version=SOLVER_VERSION,
            tags=[CASE_TAG, SOLVER_VERSION, tag, "pitch_zdes_sa"],
            fork_from=steady,
        )
    submitted_cases.append(case)

Wait for completion¶

Block until every submitted case has finished. This can take a long time (hours), since each mean angle runs a steady startup plus one or more unsteady forced-pitch solves.

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

Postprocessing¶

The published results are two curves versus mean angle of attack:

  1. Pitch damping derivative: $C_{m_q}+C_{m_{\dot{\alpha}}}$ vs $\alpha$.
  2. Mean pitching moment: $C_m$ vs $\alpha$.

Both are extracted from the same reduction of each case's pitching-moment history and compared against the AEDC experiment and published CFD references. The cells below read the results straight from the in-kernel submitted_cases objects (no project is re-opened).

Setup¶

Styling colors (inlined from the shared benchmark plot style), reference-data locations and markers, and the harmonic-reduction constants: the pitch amplitude and reduced frequency, the per-angle averaging window, and the moment-plot y-limits.

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

FLOW360_COLOR = "#00643c"   # Flow360 results
REFERENCE_COLOR = "black"   # experimental / reference CFD

RESULTS_DIR = "results"
REF_DATA_DIR = "ref_data"

POINTS_PER_CYCLE = 360
DEFAULT_CYCLES_TO_AVERAGE = 1.0
USE_ALPHA_AVERAGING_WINDOW = True
AMPLITUDE_DEG = 0.32
REDUCED_FREQUENCY = 0.16
PHASE_OFFSET_DEG = 0.0
DT_CV_WARN = 0.02
CM_YLIM = (-10, 8)

ALPHA_RE = re.compile(r"alpha\s*=?\s*(?P<alpha>[-+]?\d*\.?\d+)", re.IGNORECASE)
SLUG_ALPHA_RE = re.compile(r"alpha(?P<alpha>m?\d+(?:p\d+)?)", re.IGNORECASE)

CMQ_REFERENCES = (
    ("AEDC-TR-75-43", "exp_data_Cmq.csv", "x", REFERENCE_COLOR),
    ("AIAA 2012-0691", "AIAA2012_Cmq.csv", "s", "#4d4d4d"),
    ("FUN3D", "FUN3D_Cmq.csv", "^", "#ff7f0e"),
    ("KESTREL", "KESTREL_Cmq.csv", "D", "#9467bd"),
)
CM_REFERENCES = (
    ("AEDC-TR-75-43", "exp_data_Cm.csv", "x", REFERENCE_COLOR),
    ("AIAA 2012-0691", "AIAA2012_Cm.csv", "s", "#4d4d4d"),
    ("FUN3D", "FUN3D_Cm.csv", "^", "#ff7f0e"),
    ("KESTREL", "KESTREL_Cm.csv", "D", "#9467bd"),
)

Reduction and plotting routines¶

For each case the pitching-moment coefficient CMy is extracted per physical step from the total-forces history, split into complete pitching cycles, and reduced by a periodic cubic-Lagrange harmonic fit: the mean gives $C_m$, and the in-/out-of-phase components give $C_{m_\alpha}$ and the damping sum $C_{m_q}+C_{m_{\dot{\alpha}}}$. The averaging window narrows at the highest angles (2 cycles up to 75°, 1 at 80°, half a cycle at 85°, with a doubled 720-point cycle at 80°/85°). load_reference reads a two-column comparison curve and make_plot renders one series vs angle of attack against the references.

In [ ]:
def parse_case_alpha(case, case_name):
    alpha_with_units = getattr(getattr(case, "params", None), "operating_condition", None)
    if alpha_with_units is not None:
        alpha_with_units = getattr(alpha_with_units, "alpha", None)
    if alpha_with_units is not None:
        match = ALPHA_RE.search(str(alpha_with_units))
        if match:
            return float(match.group("alpha"))
    for regex in (ALPHA_RE, SLUG_ALPHA_RE):
        match = regex.search(str(case_name))
        if match:
            return float(match.group("alpha").replace("m", "-", 1).replace("p", "."))
    raise RuntimeError(f"Could not parse alpha from case '{case_name}'.")


def cycles_to_average_for_alpha(alpha):
    if not USE_ALPHA_AVERAGING_WINDOW:
        return DEFAULT_CYCLES_TO_AVERAGE
    if alpha <= 75.0:
        return 2.0
    if math.isclose(alpha, 80.0):
        return 1.0
    if math.isclose(alpha, 85.0):
        return 0.5
    return DEFAULT_CYCLES_TO_AVERAGE


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_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_column(df, {"physical_step", "physicalstep", "step", "step_index", "physical_step_index"})
    pseudo_col = find_column(df, {"pseudo_step", "pseudostep", "pseudo", "substep"})
    time_col = find_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):
    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}' (std(dt)/mean(dt)={cv:.2%}).")


def build_cycle_time_axis(cycle_table, n_points):
    time_values = pd.to_numeric(cycle_table["time"], errors="coerce").to_numpy(dtype=float)
    if len(time_values) == n_points 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
    return np.arange(n_points, dtype=float) / float(n_points), 1.0


def estimate_window_duration(window, n_points):
    time_values = pd.to_numeric(window["time"], errors="coerce").to_numpy(dtype=float)
    if len(time_values) > 1 and np.all(np.isfinite(time_values)):
        dt = np.diff(time_values)
        if np.all(dt > 0.0):
            return float((time_values[-1] - time_values[0]) + np.median(dt))
    return float(len(window)) / float(n_points)


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)
        x_nodes = []
        y_nodes = []
        for node in (i - 1, i, i + 1, i + 2):
            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 fit_harmonic_window(window, n_points):
    amplitude_rad = math.radians(AMPLITUDE_DEG)
    phase_offset = math.radians(PHASE_OFFSET_DEG)
    theta = (2.0 * np.pi * window["phase_idx"].to_numpy(dtype=float) / float(n_points)) + phase_offset
    design = np.column_stack([np.ones(len(window), dtype=float), np.sin(theta), np.cos(theta)])
    cm_mean, sin_coeff, cos_coeff = np.linalg.lstsq(design, window["cm"].to_numpy(dtype=float), rcond=None)[0]
    return {
        "cm_mean_global": float(cm_mean),
        "cm_alpha_bar_global": float(sin_coeff / amplitude_rad),
        "cmq_plus_cmadot_global": float(cos_coeff / (REDUCED_FREQUENCY * amplitude_rad)),
    }


def compute_harmonic_cycles(table, n_points, cycles_to_average):
    if cycles_to_average <= 0.0:
        raise RuntimeError("cycles_to_average must be positive.")
    amplitude_rad = math.radians(AMPLITUDE_DEG)
    step_rel = table["step"].to_numpy(dtype=int) - int(table["step"].min())
    cycle_id = step_rel // n_points
    phase_idx = step_rel % n_points
    work = pd.DataFrame(
        {
            "cycle_id": cycle_id,
            "phase_idx": phase_idx,
            "cm": table["cm"].to_numpy(dtype=float),
            "time": table["time"].to_numpy(dtype=float),
        }
    )
    complete_cycles = []
    for cycle, group in work.groupby("cycle_id", sort=True):
        group = group.sort_values("phase_idx")
        if group["phase_idx"].nunique() == n_points:
            complete_cycles.append((int(cycle), group))
    full_cycle_count = int(round(cycles_to_average))
    if not math.isclose(cycles_to_average, full_cycle_count):
        window_points = int(round(cycles_to_average * n_points))
        if window_points < 4:
            raise RuntimeError("Harmonic window must contain at least 4 points.")
        if len(complete_cycles) * n_points < window_points:
            raise RuntimeError(
                f"Not enough full cycles for damping estimate (need >= {cycles_to_average:g}, got {len(complete_cycles)})."
            )
        window = pd.concat([group for _, group in complete_cycles], ignore_index=True).tail(window_points)
        coeffs = fit_harmonic_window(window, n_points)
        return [
            {
                "cycle_id": int(window["cycle_id"].iloc[-1]),
                "period": estimate_window_duration(window, n_points),
                "cm_mean_cycle": coeffs["cm_mean_global"],
                "cm_alpha_bar_cycle": coeffs["cm_alpha_bar_global"],
                "cmq_plus_cmadot_cycle": coeffs["cmq_plus_cmadot_global"],
                "points_per_cycle": n_points,
                "points_averaged": window_points,
                "cycles_averaged": cycles_to_average,
            }
        ], coeffs
    if len(complete_cycles) < full_cycle_count:
        raise RuntimeError(
            f"Not enough full cycles for damping estimate (need >= {full_cycle_count}, got {len(complete_cycles)})."
        )
    complete_cycles = complete_cycles[-full_cycle_count:]
    cm_retained = np.concatenate([group["cm"].to_numpy(dtype=float) for _, group in complete_cycles])
    cm_global_mean = float(np.mean(cm_retained))
    phase_offset = math.radians(PHASE_OFFSET_DEG)
    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, n_points)
        theta = (2.0 * np.pi * phase / float(n_points)) + 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),
                "points_per_cycle": n_points,
            }
        )
    return cycle_values, {
        "cm_mean_global": cm_global_mean,
        "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_name, case):
    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)
    alpha = parse_case_alpha(case, case_name)
    points_per_cycle = 720 if alpha in (80.0, 85.0) else POINTS_PER_CYCLE
    cycles_to_average = cycles_to_average_for_alpha(alpha)
    _, global_coeffs = compute_harmonic_cycles(table, points_per_cycle, cycles_to_average)
    return {
        "case_id": case.id,
        "case_name": case_name,
        "alpha": alpha,
        "cm_mean_global": float(global_coeffs["cm_mean_global"]),
        "cm_alpha_bar_global": float(global_coeffs["cm_alpha_bar_global"]),
        "cmq_plus_cmadot_global": float(global_coeffs["cmq_plus_cmadot_global"]),
    }


def load_reference(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 pd.DataFrame(columns=["alpha", "value"])
    df = pd.read_csv(path, header=None)
    if df.shape[1] < 2:
        warnings.warn(f"Reference data file must have at least two columns: {path}")
        return pd.DataFrame(columns=["alpha", "value"])
    return pd.DataFrame(
        {
            "alpha": pd.to_numeric(df.iloc[:, 0], errors="coerce"),
            "value": pd.to_numeric(df.iloc[:, 1], errors="coerce"),
        }
    ).dropna(subset=["alpha", "value"])


def build_reference_sets(reference_defs):
    return [
        (label, load_reference(filename), marker, color)
        for label, filename, marker, color in reference_defs
    ]


def make_plot(results, references, value_key, ylabel, ylim, out_path):
    ordered = sorted(results, key=lambda row: row["alpha"])
    x = [row["alpha"] for row in ordered]
    y = [row[value_key] 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.plot(x, y, "o-", linewidth=2.0, markersize=6, color=FLOW360_COLOR, label="Flow360 SA-ZDES")
    for label, data, marker, color in references:
        if data.empty:
            continue
        ax.plot(
            data["alpha"].to_numpy(dtype=float),
            data["value"].to_numpy(dtype=float),
            linestyle="None",
            marker=marker,
            markersize=7,
            color=color,
            label=label,
            zorder=5,
        )
    ax.set_xlabel("Alpha (deg)")
    ax.set_ylabel(ylabel)
    if ylim is not None:
        ax.set_ylim(*ylim)
    ax.legend(loc="best", frameon=True)
    fig.savefig(out_path, dpi=220, bbox_inches="tight")
    plt.show()

Reduce the swept-case results¶

Run the harmonic reduction on every completed case (using the in-kernel objects directly) to obtain the mean pitching moment and pitch-damping sum at each mean angle of attack.

In [ ]:
results = []
for case in submitted_cases:
    print(f"Processing case: {case.name} ({case.id})")
    results.append(process_case(case.name, case))

if not results:
    raise RuntimeError("No cases were successfully processed.")

Pitch damping derivative versus angle of attack¶

The published pitch-damping curve, $C_{m_q}+C_{m_{\dot{\alpha}}}$ vs $\alpha$, against the AEDC experiment and the AIAA 2012, FUN3D and KESTREL references (results/pitch_damping_summary.png).

In [ ]:
cmq_references = build_reference_sets(CMQ_REFERENCES)
make_plot(
    results,
    cmq_references,
    "cmq_plus_cmadot_global",
    r"$C_{m_q}+C_{m_{\dot{\alpha}}}$",
    None,
    os.path.join(RESULTS_DIR, "pitch_damping_summary.png"),
)

Mean pitching moment versus angle of attack¶

The published mean pitching-moment curve, $C_m$ vs $\alpha$, against the same reference set (results/pitch_moment_summary.png).

In [ ]:
cm_references = build_reference_sets(CM_REFERENCES)
make_plot(
    results,
    cm_references,
    "cm_mean_global",
    r"$C_m$",
    CM_YLIM,
    os.path.join(RESULTS_DIR, "pitch_moment_summary.png"),
)