DARPA SUBOFF: Submarine Resistance Validation¶

DARPA's DTRC 5470 is a standard marine CFD benchmark backed by high-quality towing-tank data. This case uses it to validate Flow360's liquid operating condition on an external-flow problem where the total resistance can be compared directly against experimental measurements.

The study sweeps six speeds from 5.93 to 17.79 knots. At each speed it runs a steady Spalart-Allmaras RANS solve and then a Spalart-Allmaras DDES (detached-eddy) solve forked from the converged RANS state, giving 12 cases in total. The published result compares the predicted total resistance against the towing-tank data across the speed range.

The notebook runs top-to-bottom in a single kernel: load the mesh, submit the sweep, wait for completion, then post-process the in-kernel case objects into the resistance-versus-speed 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 couple of plotting/data packages:

pip install pandas matplotlib

Run the cells top to bottom in a single kernel.

Imports¶

In [ ]:
import os

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

from flow360.examples import download_benchmark_assets

Input data¶

The post-processing compares Flow360 against the towing-tank measurements. Fetch the reference data from the public benchmark bucket; this recreates a local ./ref_data/ directory that the plotting cell reads from.

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

Load project¶

The root asset for this case is a pre-generated volume mesh of the fully appended DARPA SUBOFF model. Download it from the public benchmark bucket and start a fresh project from it, no private project id is needed.

In [ ]:
root = download_benchmark_assets("DARPA_Suboff", "root_assets")
# The snapshot bundles the mesh with solver metadata and logs; select the
# volume mesh file (.cgns or .ugrid, possibly .zst-compressed) to load.
volume_mesh_file = next(
    f for f in root if f.removesuffix(".zst").endswith((".cgns", ".ugrid"))
)
project = fl.Project.from_volume_mesh(volume_mesh_file, name="DARPA SUBOFF")

Simulation setup¶

Every case shares the same geometry and boundary conditions; only the speed, time-stepping scheme, and turbulence/fluid configuration change. create_params assembles the full SimulationParams for one run: a liquid operating condition for water at the requested velocity, the hull reference geometry, a no-slip wall on the hull, a freestream farfield, and three streamwise slices (at 0.978, 1.04, and 1.2 hull-lengths) writing Cp and the primitive variables.

The two fluid builders define the turbulence approach: make_rans_fluid is a steady Spalart-Allmaras RANS with the low-Mach preconditioner, and make_ddes_fluid is Spalart-Allmaras run as a Detached-Eddy Simulation (DDES shielding) with looser per-step tolerances for the unsteady solve.

In [ ]:
def create_params(velocity, project, time_stepping, fluid):
    with fl.SI_unit_system:
        length = 4356
        slice_ratios = [0.978, 1.04, 1.2]
        scaled_slices = [length * ratio for ratio in slice_ratios]
        volume_mesh = project.volume_mesh

        params = fl.SimulationParams(
            operating_condition=fl.LiquidOperatingCondition(
                velocity_magnitude=velocity,
                material=fl.Water(
                    name="Water",
                    density=1000 * fl.u.kg / fl.u.m**3,
                    dynamic_viscosity=0.001002 * fl.u.kg / fl.u.m / fl.u.s,
                ),
                reference_velocity_magnitude=velocity,
            ),
            reference_geometry=fl.ReferenceGeometry(
                moment_center=(2009.47, 0, 2.17) * fl.u.mm,
                moment_length=(4.261, 1, 1) * fl.u.m,
                area=6.34949082 * fl.u.m**2,
            ),
            time_stepping=time_stepping,
            models=[
                fl.Wall(
                    surfaces=volume_mesh["blk-1/fluid:wall"],
                ),
                fl.Freestream(
                    surfaces=volume_mesh["blk-1/fluid:farfield"],
                    turbulence_quantities=fl.TurbulenceQuantities(
                        viscosity_ratio=10,
                    ),
                ),
                fluid,
            ],
            outputs=[
                fl.SliceOutput(
                    slices=[
                        fl.Slice(
                            name=f"{scaled_slices[0]}",
                            normal=(1, 0, 0),
                            origin=(scaled_slices[0], 0, 0) * fl.u.mm,
                        ),
                        fl.Slice(
                            name=f"{scaled_slices[1]}",
                            normal=(1, 0, 0),
                            origin=(scaled_slices[1], 0, 0) * fl.u.mm,
                        ),
                        fl.Slice(
                            name=f"{scaled_slices[2]}",
                            normal=(1, 0, 0),
                            origin=(scaled_slices[2], 0, 0) * fl.u.mm,
                        ),
                    ],
                    output_format="both",
                    output_fields=[
                        "primitiveVars",
                        "Cp",
                    ],
                ),
            ],
        )
        return params


def make_rans_fluid():
    return fl.Fluid(
        turbulence_model_solver=fl.SpalartAllmaras(),
        navier_stokes_solver=fl.NavierStokesSolver(
            low_mach_preconditioner=True,
        ),
    )


def make_ddes_fluid():
    return fl.Fluid(
        turbulence_model_solver=fl.SpalartAllmaras(
            relative_tolerance=1e-2,
            hybrid_model=fl.DetachedEddySimulation(shielding_function="DDES"),
        ),
        navier_stokes_solver=fl.NavierStokesSolver(
            relative_tolerance=1e-2,
            low_mach_preconditioner=True,
        ),
    )

Submit cases¶

Sweep the six speeds. At each speed we submit a steady RANS case (up to 5000 steps, adaptive CFL) and then an unsteady DDES case (2500 steps of 5e-5 s) forked from the converged RANS case so it starts from a settled field. Cases are named u<speed>kt_<rans|ddes>_<solver_version>. We keep each submitted case alongside its speed and model label so the post-processing can use the in-kernel objects directly.

In [ ]:
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")
VELOCITIES_KNOTS = [5.93, 10.0, 11.85, 13.92, 16.0, 17.79]
KNOT_TO_MPS = 0.514444

submitted_cases = []
for velocity_knots in VELOCITIES_KNOTS:
    velocity = velocity_knots * KNOT_TO_MPS

    steady = fl.Steady(
        max_steps=5000,
        CFL=fl.AdaptiveCFL(
            convergence_limiting_factor=0.7,
        ),
    )
    rans_params = create_params(
        velocity=velocity,
        project=project,
        time_stepping=steady,
        fluid=make_rans_fluid(),
    )
    rans_name = f"u{velocity_knots:g}kt_rans_{SOLVER_VERSION}"
    rans_case = project.run_case(
        params=rans_params,
        name=rans_name,
        solver_version=SOLVER_VERSION,
    )
    submitted_cases.append({"velocity_knots": velocity_knots, "model": "rans", "case": rans_case})

    unsteady = fl.Unsteady(
        CFL=fl.AdaptiveCFL(
            convergence_limiting_factor=0.75,
        ),
        step_size=0.00005 * fl.u.s,
        steps=2500,
    )
    ddes_params = create_params(
        velocity=velocity,
        project=project,
        time_stepping=unsteady,
        fluid=make_ddes_fluid(),
    )
    ddes_name = f"u{velocity_knots:g}kt_ddes_{SOLVER_VERSION}"
    ddes_case = project.run_case(
        params=ddes_params,
        name=ddes_name,
        fork_from=rans_case,
        solver_version=SOLVER_VERSION,
    )
    submitted_cases.append({"velocity_knots": velocity_knots, "model": "ddes", "case": ddes_case})

Wait for completion¶

Block until every submitted case has finished. This can take a long time (minutes to hours): each speed runs a steady RANS solve followed by a time-accurate DDES solve.

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

Postprocessing¶

The published result is a single figure: total resistance versus speed, Flow360 RANS and DDES against the towing-tank experiment.

For each case we take the settled drag coefficient (the mean of the last 10% of the force history) and convert it to a total resistance, R = C_D · ½ρV² · A_ref, using the water density and the hull reference area. We then read the experimental resistance from ref_data/exp_data.csv, line up the RANS and DDES predictions against it per speed, and plot all three. The cells read straight from the in-kernel case objects, no project is re-opened.

Settled resistance per case¶

Constants and two small helpers: mean_last_10_percent averages the tail of the convergence history to get a settled coefficient, and resistance_from_cd converts a drag coefficient to a total resistance in newtons. The DDES history is ordered by physical step and the RANS history by pseudo step before averaging.

In [ ]:
FLOW360_COLOR = "#00643c"   # Flow360 results
REFERENCE_COLOR = "black"   # experimental reference

KNOT_TO_MPS = 0.514444
WATER_DENSITY = 1000.0
REFERENCE_AREA = 6.34949082


def mean_last_10_percent(values):
    n_values = len(values)
    if n_values == 0:
        return None
    start = max(int(n_values * 0.9), 0)
    return float(values.iloc[start:].mean())


def resistance_from_cd(cd_value, velocity_knots):
    velocity_mps = velocity_knots * KNOT_TO_MPS
    dynamic_pressure = 0.5 * WATER_DENSITY * velocity_mps**2
    return float(cd_value) * dynamic_pressure * REFERENCE_AREA


case_rows = []
for record in submitted_cases:
    case = record["case"]
    model = record["model"]
    velocity_knots = record["velocity_knots"]

    force_data = case.results.total_forces.as_dataframe()
    if model == "ddes" and "physical_step" in force_data.columns:
        force_data = force_data.sort_values("physical_step")
    elif model == "rans" and "pseudo_step" in force_data.columns:
        force_data = force_data.sort_values("pseudo_step")

    cd_mean = mean_last_10_percent(force_data["CD"])
    case_rows.append(
        {
            "model": model,
            "velocity_knots": velocity_knots,
            "cd_mean": cd_mean,
            "resistance_n": resistance_from_cd(cd_mean, velocity_knots),
        }
    )

all_cases_df = pd.DataFrame(case_rows).sort_values(["velocity_knots", "model"]).reset_index(drop=True)

Total resistance vs speed¶

Merge the RANS and DDES resistances onto the experimental curve by speed, then plot Flow360 DDES and RANS against the towing-tank measurements.

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

reference_df = pd.read_csv("ref_data/exp_data.csv").rename(
    columns={"V_knots": "velocity_knots", "Resistance": "reference_resistance_n"}
)
reference_df["velocity_knots"] = pd.to_numeric(reference_df["velocity_knots"])
reference_df["reference_resistance_n"] = pd.to_numeric(reference_df["reference_resistance_n"])

rans_df = (
    all_cases_df[all_cases_df["model"] == "rans"][["velocity_knots", "resistance_n"]]
    .rename(columns={"resistance_n": "rans_resistance_n"})
)
ddes_df = (
    all_cases_df[all_cases_df["model"] == "ddes"][["velocity_knots", "resistance_n"]]
    .rename(columns={"resistance_n": "ddes_resistance_n"})
)

comparison_df = reference_df.merge(rans_df, on="velocity_knots", how="left").merge(
    ddes_df, on="velocity_knots", how="left"
)

plt.figure(figsize=(7, 4.5))
plt.plot(
    comparison_df["velocity_knots"],
    comparison_df["ddes_resistance_n"],
    label="Flow360 DDES",
    color=FLOW360_COLOR,
    linewidth=2,
    marker="o",
    markersize=4,
)
plt.plot(
    comparison_df["velocity_knots"],
    comparison_df["rans_resistance_n"],
    label="Flow360 RANS",
    color=FLOW360_COLOR,
    linewidth=2,
    linestyle="--",
    marker="s",
    markersize=4,
    alpha=0.8,
)
plt.scatter(
    comparison_df["velocity_knots"],
    comparison_df["reference_resistance_n"],
    marker="x",
    color=REFERENCE_COLOR,
    label="Experiment",
)
plt.xlabel("Velocity (knots)")
plt.ylabel("Resistance (N)")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("results/resistance_vs_v.png", dpi=300)
plt.show()