DTU 10 MW Reference Wind Turbine

Flow360 wind-turbine rotor-aerodynamics validation against DTU CFD reference data. The DTU 10 MW Reference Wind Turbine is a three-bladed upwind rotor (178.3 m diameter, 119 m hub height), the open 10 MW successor to the NREL 5 MW reference, widely used to validate wind-turbine aerodynamics codes.

This case runs a below-rated velocity sweep (U = 8, 9, 10, 11 m/s) using the Spalart-Allmaras DDES turbulence model on a sliding-mesh rotating zone, plus a mesh-convergence study at U = 11 m/s over four mesh levels. Each operating point is a two-stage unsteady run: stage 1 settles the flow (10°/step, 16 revs), stage 2 forks from it and gathers statistics (3°/step, 4 revs, low-dissipation numerics). Results are compared against DTU CFD for shaft power, rotor thrust, and spanwise sectional loads (Fx, Fz, Ct, Cp).

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

Everything needed is either flow360 or a standard scientific-Python library.

import os
import re
import flow360 as fl
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

from flow360.examples import download_benchmark_assets

Input data

Fetch the DTU CFD reference data (turb/ power curve and spanwise-load files) that the post-processing compares against.

download_benchmark_assets("DTU_WindTurbine", "turb")

Load project

Fetch the published turbine geometry and start a fresh Flow360 project from it, then group the surface faces by their faceName tag so the hub, blade and tip walls can be referenced by name.

root_asset_files = download_benchmark_assets("DTU_WindTurbine", "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="DTU 10 MW Reference Wind Turbine")

geo = project.geometry
geo.group_faces_by_tag("faceName")

Configuration

Solver version, the study definition (time-stepping stages, mesh levels, turbulence mode), physical constants, and the four operating conditions (wind speed / rotor RPM).

SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.8")

FORK_FROM = None  # set to a case id to fork stage 1 from an existing run

# (dt_deg, num_revs) per stage, the last stage uses the low-dissipation solver
STAGES = [
    (10, 16),
    (3, 4),
]

# Mesh refinement levels (spacing multiplier, label). 1.5x between levels ~= 3.4x nodes.
MESH_LEVELS = [
    (1.5**2,     "Coarse"),     # 2.25 - coarsest
    (1.5**1,     "Medium"),     # 1.50 - reference for the velocity sweep
    (1.5**0,     "Fine"),       # 1.00
    (1.5**(-1),  "Very Fine"),  # 0.67 - finest
]
SWEEP_MESH_REFINE = 1.5**1  # Medium mesh used for the velocity sweep

TURBULENCE_MODES = ["DDES"]  # "DDES" = S-A + DetachedEddySimulation hybrid, low-dissipation NS on last stage

# Physical / geometric constants
R_TIP       = 89.166      # m  - rotor tip radius
RHO         = 1.225       # kg/m^3
TEMPERATURE = 288.15      # K
GAMMA       = 1.4
R_GAS       = 287.0       # J/(kg K)
DEG_PER_REV = 360
REF_AREA    = 24977.47    # m^2

# Operating conditions (wind speed and rotor angular velocity)
_all_conditions = [
    {"u_inf": 8.0,  "omega_rpm": 6.426},
    {"u_inf": 9.0,  "omega_rpm": 7.229},
    {"u_inf": 10.0, "omega_rpm": 8.032},
    {"u_inf": 11.0, "omega_rpm": 8.836},
]

Geometry and zones

Resolve the wall faces (hub, blade, tip), define the rotating-zone cylinder that encloses the rotor, and set up the automated farfield.

hub_faces   = geo["hub"]
blade_faces = geo["blade"]
tip_faces   = geo["tip"]
wall_faces  = [hub_faces, blade_faces, tip_faces]

rotation_cylinder = fl.Cylinder(
    name="rotation_cylinder",
    axis=(0, 0, 1),
    center=(0, 0, 0) * fl.u.m,
    height=40 * fl.u.m,
    inner_radius=0 * fl.u.m,
    outer_radius=110 * fl.u.m,
)
farfield = fl.AutomatedFarfield()

speed_of_sound = np.sqrt(GAMMA * R_GAS * TEMPERATURE)  # m/s
wall           = fl.Wall(name="Wall", entities=wall_faces)

Meshing setup

Build the meshing parameters for a given refinement level. All spacing targets scale linearly with mesh_refine: boundary-layer and surface spacing, a uniform refinement region above the rotor, extra surface refinement on the tip, and the rotating volume zone enclosing the blade walls.

def make_meshing(mesh_refine):
    return fl.MeshingParams(
        defaults=fl.MeshingDefaults(
            boundary_layer_first_layer_thickness=1e-6 * fl.u.m,
            surface_max_edge_length=mesh_refine * 1000 * fl.u.mm,
        ),
        refinements=[
            fl.UniformRefinement(
                name="Uniform refinement",
                spacing=mesh_refine * 2000 * fl.u.mm,
                entities=[
                    fl.Cylinder(
                        name="uniform_refinement",
                        axis=(0, 0, 1),
                        center=(0, 0, 55) * fl.u.m,
                        height=200 * fl.u.m,
                        inner_radius=0 * fl.u.m,
                        outer_radius=200 * fl.u.m,
                    )
                ],
            ),
            fl.SurfaceRefinement(
                name="Surface refinement",
                max_edge_length=mesh_refine * 10 * fl.u.mm,
                entities=tip_faces,
            ),
            fl.UniformRefinement(
                name="Uniform refinement",
                entities=[rotation_cylinder],
                spacing=mesh_refine * 1000 * fl.u.mm,
            ),
        ],
        volume_zones=[
            farfield,
            fl.RotationVolume(
                name="rotating_zone",
                spacing_axial=mesh_refine * 1000 * fl.u.mm,
                spacing_radial=mesh_refine * 1000 * fl.u.mm,
                spacing_circumferential=mesh_refine * 1000 * fl.u.mm,
                entities=[rotation_cylinder],
                enclosed_entities=wall_faces,
            ),
        ],
    )

Physics setup

Build the fluid model. In DDES mode the flow uses Spalart-Allmaras with a DetachedEddySimulation hybrid; the final time-stepping stage additionally switches the Navier-Stokes solver to low-dissipation numerics. (URANS mode would use plain S-A with no hybrid and no low-dissipation.)

def make_fluid(urans, low_dissipation=False):
    if urans:
        return fl.Fluid(
            navier_stokes_solver=fl.NavierStokesSolver(),
            turbulence_model_solver=fl.SpalartAllmaras(),
        )
    ns = (
        fl.NavierStokesSolver(numerical_dissipation_factor=0.2, kappa_MUSCL=-0.33)
        if low_dissipation
        else fl.NavierStokesSolver()
    )
    return fl.Fluid(
        navier_stokes_solver=ns,
        turbulence_model_solver=fl.SpalartAllmaras(hybrid_model=fl.DetachedEddySimulation()),
    )

Outputs setup

Request the solver outputs: a time-averaged spanwise force distribution (for the sectional-load plots) and time-averaged rotor force/moment coefficients CFz, CMz (for integrated thrust and power), averaged over a moving window of two revolutions.

def make_outputs(dt_deg_eff):
    return [
        fl.TimeAverageForceDistributionOutput(
            name="y_force_distribution",
            distribution_direction=[0, 1, 0],
            distribution_type="incremental",
        ),
        fl.ForceOutput(
            name="rotor_forces",
            models=[wall],
            output_fields=["CFz", "CMz"],
            moving_statistic=fl.MovingStatistic(
                method="mean",
                moving_window_size=int(DEG_PER_REV * 2 / dt_deg_eff),
            ),
        ),
    ]

Case builder

Assemble and submit all time-stepping stages for one (wind speed, mesh level, turbulence mode) triple. Stage 1 is a fresh run; each later stage forks from the previous one, swaps in the appropriate fluid model (low-dissipation on the last stage), and re-times the run. The operating condition is set at 90° incidence (rotor plane normal to the freestream) with a rotating-zone angular velocity.

def run_operating_condition(op, mesh_refine, urans, fork_id=None):
    """Submit all STAGES for one (wind speed, mesh_refine, urans) triple. Returns list of Case objects."""
    cases = []
    with fl.SI_unit_system:
        u_inf     = op["u_inf"]
        omega_rpm = op["omega_rpm"]
        omega     = omega_rpm * fl.u.rpm

        mach           = u_inf / speed_of_sound
        reference_mach = ((omega_rpm * fl.u.rpm).to("rad/s") * R_TIP * fl.u.m).value / speed_of_sound
        meshing        = make_meshing(mesh_refine)

        def make_params(dt_deg, num_revs, low_dissipation=False):
            dt_deg_eff = dt_deg
            dt    = ((dt_deg_eff * fl.u.deg) / omega).to("s")
            steps = int(num_revs * DEG_PER_REV / dt_deg_eff)
            return fl.SimulationParams(
                operating_condition=fl.AerospaceCondition.from_mach(
                    mach=mach,
                    alpha=90 * fl.u.deg,
                    beta=0 * fl.u.deg,
                    thermal_state=fl.ThermalState(
                        temperature=TEMPERATURE * fl.u.K,
                        density=RHO * fl.u.kg / fl.u.m**3,
                    ),
                    reference_mach=reference_mach,
                ),
                models=[
                    make_fluid(urans, low_dissipation),
                    wall,
                    fl.Freestream(name="Freestream", surfaces=farfield.farfield),
                    fl.Rotation(name="Rotation", spec=fl.AngularVelocity(omega_rpm * fl.u.rpm), entities=[rotation_cylinder]),
                ],
                time_stepping=fl.Unsteady(order_of_accuracy=2, steps=steps, step_size=dt, CFL=fl.AdaptiveCFL()),
                reference_geometry=fl.ReferenceGeometry(
                    area=REF_AREA * fl.u.m**2,
                    moment_center=(0, 0, 0) * fl.u.m,
                    moment_length=(R_TIP, R_TIP, R_TIP) * fl.u.m,
                ),
                meshing=meshing,
                outputs=make_outputs(dt_deg_eff),
            )

        turb_tag = "URANS" if urans else "DDES"
        print(f"\n== {turb_tag}  u_inf={u_inf} m/s, omega={omega_rpm} rpm, MR={mesh_refine:.2f} ==")
        for i, (dt_deg, num_revs) in enumerate(STAGES):
            low_diss   = (i == len(STAGES) - 1) and not urans
            dt_deg_eff = dt_deg
            stage_name = f"{turb_tag}_u{u_inf:.0f}_omega{omega_rpm:.3f}_MR{mesh_refine:.2f}_S{i+1}_{dt_deg_eff:.2f}deg_{SOLVER_VERSION}"

            if fork_id:
                parent     = fl.Case(fork_id)
                dt         = ((dt_deg_eff * fl.u.deg) / omega).to("s")
                steps      = int(num_revs * DEG_PER_REV / dt_deg_eff)
                new_models = parent.params.models.copy()
                new_models[0] = make_fluid(urans, low_diss)
                stage_params = parent.params.copy(update=dict(
                    time_stepping=fl.Unsteady(order_of_accuracy=2, steps=steps, step_size=dt, CFL=fl.AdaptiveCFL()),
                    models=new_models,
                ))
                case = project.run_case(stage_params, name=stage_name, solver_version=SOLVER_VERSION,
                                        use_beta_mesher=True, fork_from=parent)
            else:
                stage_params = make_params(dt_deg, num_revs, low_diss)
                case = project.run_case(stage_params, name=stage_name, solver_version=SOLVER_VERSION,
                                        use_beta_mesher=True)

            print(f"  Submitted stage {i+1}: {case.name}  ({dt_deg_eff:.2f} deg/step, {num_revs} revs{', low-diss' if low_diss else ''})")
            cases.append(case)
            fork_id = case.id

    return cases

Submit cases

Submit the full study for each turbulence mode: first the mesh-convergence run at U = 11 m/s over all four mesh levels, then the velocity sweep at U = 8, 9, 10 m/s on the Medium mesh. (U = 11 m/s on the Medium mesh is already covered by the convergence run, so it is not duplicated.) All submitted cases are collected for later post-processing.

submitted_cases = []

for turb_tag in TURBULENCE_MODES:
    urans = (turb_tag == "URANS")
    print(f"\n{'#'*60}\n#  Turbulence mode: {turb_tag}\n{'#'*60}")

    # Mesh convergence: U = 11 m/s at all four mesh levels
    op11 = _all_conditions[-1]  # {"u_inf": 11.0, "omega_rpm": 8.836}
    for mesh_refine, mesh_label in MESH_LEVELS:
        print(f"\n{'='*60}")
        print(f"  {turb_tag} mesh convergence - {mesh_label}  (MR={mesh_refine:.2f})")
        print(f"{'='*60}")
        submitted_cases.extend(run_operating_condition(op11, mesh_refine, urans, FORK_FROM))

    # Velocity sweep: U = 8, 9, 10 m/s on the Medium mesh
    sweep_ops = [op for op in _all_conditions if op["u_inf"] != 11.0]
    print(f"\n{'='*60}")
    print(f"  {turb_tag} velocity sweep - Medium mesh  (MR={SWEEP_MESH_REFINE:.2f})")
    print(f"{'='*60}")
    for op in sweep_ops:
        submitted_cases.extend(run_operating_condition(op, SWEEP_MESH_REFINE, urans, FORK_FROM))

print(f"\nSubmitted {len(submitted_cases)} cases.")

Wait for completion

Block until every submitted case has finished. This can take a long time (unsteady rotor runs over several mesh levels and wind speeds).

for case in submitted_cases:
    case.wait()

Postprocessing

Reuse the case objects submitted above (one kernel, no project reload) to reproduce the three published figures: the power/thrust velocity sweep, the spanwise sectional loads at U = 11 m/s across mesh levels, and the spanwise loading by wind speed. All are compared against the DTU CFD reference data fetched into turb/.

Setup

Constants, plot colors, DTU reference data, and selection of the final-stage cases from the in-kernel sweep, the convergence cases (U = 11 m/s, all mesh levels) and the Medium-mesh velocity-sweep cases (all wind speeds).

# Flow360 = brand green, reference/experiment = black
FLOW360_COLOR   = "#00643c"
REFERENCE_COLOR = "black"

L_CHORD_REF             = 3.4
U_INF_CONVERGENCE_M_S   = 11
OMEGA_CONVERGENCE_RAD_S = (8.836 * fl.u.rpm).to("rad/s") / fl.u.rad

# DTU CFD reference: power curve and the U = 11 m/s spanwise loads
power_curve = pd.read_csv("turb/power_curve.dat", sep=r"\s+", comment="#",
                          names=["winf", "power", "thrust", "Cp", "Ct", "pitch", "omega_ref"]
                          ).set_index("winf")
df_turb = pd.read_csv(f"turb/wsp_{U_INF_CONVERGENCE_M_S}_spanwise_loads.dat", sep=r"\s+", comment="#",
                      names=["radius", "fx", "fz", "localcp", "localct"])

# Final time-stepping stage number (run.py names stages S1, S2, ...)
_stage_re    = re.compile(r"_S(\d+)_")
_final_stage = max(int(m.group(1)) for c in submitted_cases
                   for m in [_stage_re.search(c.name)] if m)
_final_tag   = f"_S{_final_stage}_"
_sweep_mr_tag = f"_MR{SWEEP_MESH_REFINE:.2f}_"

_u_inf_re = re.compile(r"(?:^|_)u(\d+)_")
def _parse_u_inf(name):
    m = _u_inf_re.search(name)
    return int(m.group(1)) if m else None

# Convergence cases: U = 11 m/s at each mesh level (final stage), coarse -> fine
_MR_LABELS = [("_MR2.25_", "Coarse"), ("_MR1.50_", "Medium"),
              ("_MR1.00_", "Fine"), ("_MR0.67_", "Very Fine")]
CONVERGENCE_CASES = []
for mr_tag, label in _MR_LABELS:
    for c in submitted_cases:
        if c.name.startswith("DDES_") and "_u11_" in c.name and mr_tag in c.name and _final_tag in c.name:
            CONVERGENCE_CASES.append((c, label))
            break

# Velocity-sweep cases: Medium mesh, final stage, all wind speeds
VELOCITY_CASES = [
    (c, _parse_u_inf(c.name))
    for c in submitted_cases
    if c.name.startswith("DDES_") and _final_tag in c.name and _sweep_mr_tag in c.name
]
VELOCITY_CASES = [(c, u) for c, u in VELOCITY_CASES if u is not None]

Integrated loads at U = 11 m/s (medium mesh)

Integrate the Medium-mesh U = 11 m/s rotor force/moment coefficients into total thrust (kN) and shaft power (MW, using the reference RPM), and compare against the DTU CFD power curve. Written to results/DDES/integrated_loads.csv and shown below.

# Integrated loads at U = 11 m/s (Medium mesh) vs DTU CFD reference.
# Reuse the in-kernel Medium convergence case (U = 11 m/s, medium mesh); the
# force/moment integration matches the sectional-load extraction below.
os.makedirs("results/DDES", exist_ok=True)

_medium_case = next(c for c, label in CONVERGENCE_CASES if label == "Medium")
_medium_case.wait()

u_ref = _medium_case.params.operating_condition.reference_velocity_magnitude
rho   = _medium_case.params.operating_condition.thermal_state.density
q     = 0.5 * rho * u_ref**2
A     = _medium_case.params.reference_geometry.area
L     = _medium_case.params.reference_geometry.moment_length

CF = _medium_case.results.total_forces.get_averages(1/12)
Fz = float((CF["CFz"] * q * A).to("kN").to_value())
Mz = (CF["CMz"] * q * A * L[2]).to("N*m")
P  = float((Mz * OMEGA_CONVERGENCE_RAD_S).to("MW").to_value())

# DTU CFD reference at U = 11 m/s (power curve).
ref_row       = power_curve.loc[U_INF_CONVERGENCE_M_S]
ref_thrust_kN = ref_row["thrust"] / 1e3
ref_power_MW  = ref_row["power"] / 1e6

def _err_pct(v_f360, v_ref):
    return (v_f360 - v_ref) / abs(v_ref) * 100 if v_ref != 0 else 0.0

integrated_loads = pd.DataFrame([
    {"quantity": "Thrust (U = 11 m/s, medium mesh)", "units": "kN",
     "flow360": Fz, "reference": ref_thrust_kN, "error_pct": _err_pct(Fz, ref_thrust_kN)},
    {"quantity": "Shaft Power (U = 11 m/s, medium mesh)", "units": "MW",
     "flow360": P,  "reference": ref_power_MW,  "error_pct": _err_pct(P,  ref_power_MW)},
])

integrated_loads.to_csv("results/DDES/integrated_loads.csv", index=False)
print(integrated_loads.to_string(index=False))
integrated_loads

Velocity-sweep data

For each Medium-mesh case, integrate the rotor force/moment coefficients into thrust (kN) and shaft power (MW, using the reference RPM), and extract the spanwise sectional loads (Fx, Fz, Ct, Cp) used by the plots below.

sweep = []
for case, u_inf_v in VELOCITY_CASES:
    print(f"Velocity sweep - U = {u_inf_v} m/s ...")
    case.wait()
    u_ref = case.params.operating_condition.reference_velocity_magnitude
    rho   = case.params.operating_condition.thermal_state.density
    q     = 0.5 * rho * u_ref**2
    A     = case.params.reference_geometry.area
    L     = case.params.reference_geometry.moment_length

    CF = case.results.total_forces.get_averages(1/12)
    Fz = float((CF["CFz"] * q * A).to("kN").to_value())
    Mz = float((CF["CMz"] * q * A * L[2]).to("N*m").to_value())
    om = power_curve.loc[u_inf_v, "omega_ref"] if u_inf_v in power_curve.index else float("nan")
    P  = Mz * om  # W

    # Spanwise sectional loads
    df_sl  = case.results.y_slicing_force_distribution.as_dataframe()
    y_sl   = df_sl["Y"].values
    cfx_sl = np.zeros(len(y_sl))
    cfz_sl = np.zeros(len(y_sl))
    for comp in ["hub", "blade", "tip"]:
        for key, arr in [(f"rotation_cylinder/{comp}_CFx_per_span", cfx_sl),
                         (f"rotation_cylinder/{comp}_CFz_per_span", cfz_sl)]:
            if key in df_sl.columns:
                arr += df_sl[key].values
    idx_sl = np.argsort(y_sl)
    y_sl   = y_sl[idx_sl]
    cfx_sl = -cfx_sl[idx_sl]
    cfz_sl = cfz_sl[idx_sl]
    fx_sl  = (cfx_sl * q * A).to("N")
    fz_sl  = (cfz_sl * q * A).to("N")
    xi_sl  = y_sl / R_TIP
    ct_sl  = fz_sl / (0.5 * rho * (om * y_sl)**2 * L_CHORD_REF) * xi_sl
    cp_sl  = (fx_sl * y_sl) / (0.5 * rho * (om * y_sl)**2 * L_CHORD_REF) * xi_sl * (om / u_inf_v)
    mask   = y_sl > 4

    sweep.append(dict(u_inf=u_inf_v, Fz=Fz, Power=P / 1e6,
                      y=y_sl[mask], fx=fx_sl[mask], fz=fz_sl[mask],
                      ct=ct_sl[mask], cp=cp_sl[mask]))

sweep = sorted(sweep, key=lambda r: r["u_inf"])
df_sweep = pd.DataFrame([{k: r[k] for k in ("u_inf", "Fz", "Power")} for r in sweep])

Velocity sweep: Power and thrust

Flow360 shaft power and rotor thrust versus wind speed (Medium mesh) against the DTU CFD power curve.

simulated_u = df_sweep["u_inf"].values
ref_at_sim  = power_curve.reindex(simulated_u).dropna()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

ax1.plot(ref_at_sim.index, ref_at_sim["power"] / 1e6, color=REFERENCE_COLOR, linestyle="--", marker="s", label="DTU CFD")
ax1.plot(df_sweep["u_inf"], df_sweep["Power"], color=FLOW360_COLOR, marker="o", label="Flow360 Medium")
ax1.set(xlabel="U_inf (m/s)", ylabel="Power (MW)", title="Power curve")
ax1.legend(); ax1.grid(True, alpha=0.3)

ax2.plot(ref_at_sim.index, ref_at_sim["thrust"] / 1e3, color=REFERENCE_COLOR, linestyle="--", marker="s", label="DTU CFD")
ax2.plot(df_sweep["u_inf"], df_sweep["Fz"], color=FLOW360_COLOR, marker="o", label="Flow360 Medium")
ax2.set(xlabel="U_inf (m/s)", ylabel="Thrust (kN)", title="Thrust curve")
ax2.legend(); ax2.grid(True, alpha=0.3)

plt.tight_layout()
os.makedirs("results/DDES", exist_ok=True)
fig.savefig("results/DDES/velocity_sweep.png", dpi=150, bbox_inches="tight")
plt.show()

Sectional-load data at U = 11 m/s

Extract the spanwise sectional loads (Fx, Fz, Ct, Cp) at U = 11 m/s for each mesh level, using the convergence-study RPM.

conv_results = []
for case, label in CONVERGENCE_CASES:
    print(f"Sectional loads ({label}) ...")
    case.wait()
    u_ref = case.params.operating_condition.reference_velocity_magnitude
    rho   = case.params.operating_condition.thermal_state.density
    q     = 0.5 * rho * u_ref**2
    A     = case.params.reference_geometry.area

    df  = case.results.y_slicing_force_distribution.as_dataframe()
    y   = df["Y"].values
    cfx = np.zeros(len(y))
    cfz = np.zeros(len(y))
    for comp in ["hub", "blade", "tip"]:
        for key, arr in [(f"rotation_cylinder/{comp}_CFx_per_span", cfx),
                         (f"rotation_cylinder/{comp}_CFz_per_span", cfz)]:
            if key in df.columns:
                arr += df[key].values
    idx = np.argsort(y)
    y   = y[idx]; cfx = -cfx[idx]; cfz = cfz[idx]
    fx  = (cfx * q * A).to("N")
    fz  = (cfz * q * A).to("N")
    xi  = y / R_TIP
    ct  = fz / (0.5 * rho * (OMEGA_CONVERGENCE_RAD_S * y)**2 * L_CHORD_REF) * xi
    cp  = (fx * y) / (0.5 * rho * (OMEGA_CONVERGENCE_RAD_S * y)**2 * L_CHORD_REF) * xi * (OMEGA_CONVERGENCE_RAD_S / U_INF_CONVERGENCE_M_S)
    mask = y > 4
    conv_results.append(dict(label=label, y=y[mask], fx=fx[mask], fz=fz[mask], ct=ct[mask], cp=cp[mask]))

Sectional loads at U = 11 m/s

Spanwise Fx, Fz, local Ct and local Cp along the blade at U = 11 m/s for each mesh level, compared against the DTU CFD spanwise loads.

_MESH_COLORS = ["tab:blue", "tab:orange", "tab:green", "tab:red",
                "tab:purple", "tab:brown", "tab:pink", "tab:cyan"]

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))

for res, color in zip(conv_results, _MESH_COLORS):
    lbl = f"Flow360 ({res['label']})"
    ax1.plot(res["y"], res["fx"], color=color, label=lbl)
    ax2.plot(res["y"], res["fz"], color=color, label=lbl)
    ax3.plot(res["y"], res["ct"], color=color, label=lbl)
    ax4.plot(res["y"], res["cp"], color=color, label=lbl)

ax1.scatter(df_turb["radius"], df_turb["fx"],      color=REFERENCE_COLOR, s=20, zorder=5, label="DTU CFD")
ax2.scatter(df_turb["radius"], df_turb["fz"],      color=REFERENCE_COLOR, s=20, zorder=5, label="DTU CFD")
ax3.scatter(df_turb["radius"], df_turb["localct"], color=REFERENCE_COLOR, s=20, zorder=5, label="DTU CFD")
ax4.scatter(df_turb["radius"], df_turb["localcp"], color=REFERENCE_COLOR, s=20, zorder=5, label="DTU CFD")

for ax, ylabel, title in [
    (ax1, "Fx (N/m)", "Fx sectional"),
    (ax2, "Fz (N/m)", "Fz sectional (thrust)"),
    (ax3, "Ct (-)",   "Local Ct"),
    (ax4, "Cp (-)",   "Local Cp"),
]:
    ax.set(xlabel="Blade radius (m)", ylabel=ylabel, title=title)
    ax.legend(fontsize=8); ax.grid(True, alpha=0.3)

plt.tight_layout()
os.makedirs("results/DDES", exist_ok=True)
fig.savefig("results/DDES/sectional_loads.png", dpi=150, bbox_inches="tight")
plt.show()

Spanwise loading by wind speed

Spanwise Fx, Fz, local Ct and local Cp along the blade for each wind speed (Medium mesh), with DTU CFD reference points where available.

_TAB_COLORS = ["tab:blue", "tab:orange", "tab:green", "tab:red",
               "tab:purple", "tab:brown", "tab:pink", "tab:cyan"]
colors = [_TAB_COLORS[i % len(_TAB_COLORS)] for i in range(len(sweep))]

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))

legend_handles = []
for res, color in zip(sweep, colors):
    u = res["u_inf"]
    for ax, key in [(ax1, "fx"), (ax2, "fz"), (ax3, "ct"), (ax4, "cp")]:
        ax.plot(res["y"], res[key], color=color, linewidth=1.8)

    ref_file = f"turb/wsp_{int(u)}_spanwise_loads.dat"
    try:
        df_r = pd.read_csv(ref_file, sep=r"\s+", comment="#",
                           names=["radius", "fx", "fz", "localcp", "localct"])
        for ax, col in [(ax1, "fx"), (ax2, "fz"), (ax3, "localct"), (ax4, "localcp")]:
            ax.scatter(df_r["radius"], df_r[col], color=REFERENCE_COLOR, s=25, zorder=5,
                       marker="^", edgecolors="none")
    except FileNotFoundError:
        pass

    legend_handles.append(Line2D([0], [0], color=color, linewidth=1.8, label=f"Flow360 Medium (U={u} m/s)"))

dtu_handle = Line2D([0], [0], color=REFERENCE_COLOR, marker="^", markersize=6, linestyle="none", label="DTU CFD")

for ax, ylabel, title in [
    (ax1, "Fx (N/m)", "Fx sectional"),
    (ax2, "Fz (N/m)", "Fz sectional (thrust)"),
    (ax3, "Ct (-)",   "Local Ct"),
    (ax4, "Cp (-)",   "Local Cp"),
]:
    ax.set(xlabel="Blade radius (m)", ylabel=ylabel, title=title)
    ax.legend(handles=legend_handles + [dtu_handle], fontsize=7)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
os.makedirs("results/DDES", exist_ok=True)
fig.savefig("results/DDES/velocity_sweep_sectional.png", dpi=150, bbox_inches="tight")
plt.show()