PAW06 NASA 1507 Inlet: Supersonic Bleed Control Analysis¶
The NASA 1507 mixed-compression inlet is a benchmark for validating supersonic bleed control and shock-boundary-layer interaction in propulsion intakes. This case runs the PAW06 workshop Grid C configuration at Mach 3 with the Spalart-Allmaras turbulence model, AFT transition, and porous (Slater) bleed control across the four bleed regions.
Bleed mass flow is driven to the workshop targets in real time using Flow360 user-defined dynamics (a PI controller per bleed region), enabling true bleed convergence. This is a single steady solve; the post-processing compares Flow360 against workshop experiment and NASA Glenn Research Center (GRC) reference CFD for cowl static-pressure ratio and the throat-top boundary-layer rake total-pressure ratio.
The notebook runs top-to-bottom in a single kernel: load the mesh, build the simulation, submit the case, wait for completion, then post-process the in-kernel case object into the published figures.
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 matplotlib pandas pyvista
Run the cells top to bottom in a single kernel.
Imports¶
import os
import tarfile
import tempfile
from pathlib import Path
import flow360 as fl
import matplotlib.pyplot as plt
import pandas as pd
import pyvista as pv
from flow360.examples import download_benchmark_assets
Input data¶
The post-processing compares Flow360 against workshop and NASA reference data.
Fetch those files from the public benchmark bucket; this recreates a local
./ref_data/ directory that the plotting cells read from.
download_benchmark_assets("PAW06_1507", "ref_data")
REFERENCE_DIR = Path("ref_data")
Load project¶
The root asset for this case is the pre-generated workshop Grid C volume mesh. Download it from the public benchmark bucket and start a fresh project from it -- no private project id is needed. The mesh boundaries (walls, bleed regions, freestream, symmetry, outflow) are addressed by name throughout the setup.
root_asset_files = download_benchmark_assets("PAW06_1507", "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_asset_files if f.removesuffix(".zst").endswith((".cgns", ".ugrid"))
)
project = fl.Project.from_volume_mesh(volume_mesh_file, name="PAW06 NASA 1507 Inlet")
volume_mesh = project.volume_mesh
Physics setup¶
The operating point is Mach 3 with freestream density and temperature from the workshop, and a wedge sector spanning 27 of 360 degrees. The four bleed regions each receive a Slater porous-bleed wall, and their mass flow is driven to the workshop targets by a PI controller (user-defined dynamics), one per region. The constants below fix the operating condition, bleed geometry, and bleed targets.
WORKSHOP_MACH = 3.0
WEDGE_RATIO = 27 / 360
P_INF = 0.4084 * fl.u.psi
C_INF = 221.6278484 * fl.u.m / fl.u.s
INPUT_RHO = 0.080264886 * fl.u.kg / fl.u.m**3
INPUT_TEMP = 220 * fl.u.R
GRID_LENGTH = 0.0254 * fl.u.m
BLEED_PRESSURE_FACTORS = [1.5, 1.5, 1.5, 1.5]
BLEED_NAMES = ["I", "II", "III", "IV"]
TARGET_MASS_FLOW_RATES = [0.5071, 0.4351, 0.1924, 0.5777]
Each bleed region gets a UserDefinedDynamic PI controller that adjusts the
static-pressure ratio at the bleed patch to drive its mass-flow rate toward the
workshop target (non-dimensionalized by the sector density, speed of sound, and
grid length). The controller engages after 500 pseudo-steps.
def make_bleed_controllers(volume_mesh):
udds = []
start_step = 500
for bleed_name, target_mass_flow_rate in zip(BLEED_NAMES, TARGET_MASS_FLOW_RATES):
target = (
target_mass_flow_rate
* fl.u.lbm
/ fl.u.s
* WEDGE_RATIO
/ (INPUT_RHO * C_INF * GRID_LENGTH * GRID_LENGTH)
)
udds.append(
fl.UserDefinedDynamic(
name=f"massFlowRateController_{bleed_name}",
input_vars=["massFlowRate"],
constants={
"massFlowRateTarget": target,
"Kp": -target / 40,
"Ki": -target / 400,
"initialStaticPressureRatio": 0.58,
},
output_vars={
"staticPressureRatio": (
f"if (pseudoStep > {start_step}) state[0]; else initialStaticPressureRatio;"
)
},
state_vars_initial_value=["initialStaticPressureRatio", "0.0"],
update_law=[
(
f"if (pseudoStep > {start_step}) "
"state[0] + Kp * (massFlowRateTarget - massFlowRate) + Ki * state[1]; else state[0];"
),
(
f"if (pseudoStep > {start_step}) "
"state[1] + (massFlowRateTarget - massFlowRate); else state[1];"
),
],
input_boundary_patches=[volume_mesh[f"fluid/Bleed-{bleed_name}"]],
output_target=volume_mesh[f"fluid/Bleed-{bleed_name}"],
)
)
return udds
The boundary-condition and solver model list: viscous walls on the centerbody and cowl, a Slater porous bleed on each of the four bleed patches, freestream on the inflow/farfield, static-pressure outflow, slip walls on the symmetry planes, and the fluid solver itself (SA turbulence with AFT transition, a forced-transition trip box, and a MUSCL Navier-Stokes scheme).
def add_models(params, volume_mesh):
params.models = [
fl.Wall(
surfaces=[volume_mesh["fluid/CB-Walls"], volume_mesh["fluid/Cowl-Walls"]],
name="solid-walls",
),
fl.Wall(
entities=volume_mesh["fluid/Bleed-I"],
velocity=fl.SlaterPorousBleed(static_pressure=BLEED_PRESSURE_FACTORS[0] * P_INF, porosity=0.415),
name="Bleed-I",
),
fl.Wall(
entities=volume_mesh["fluid/Bleed-II"],
velocity=fl.SlaterPorousBleed(static_pressure=BLEED_PRESSURE_FACTORS[1] * P_INF, porosity=0.415),
name="Bleed-II",
),
fl.Wall(
entities=volume_mesh["fluid/Bleed-III"],
velocity=fl.SlaterPorousBleed(static_pressure=BLEED_PRESSURE_FACTORS[2] * P_INF, porosity=0.415),
name="Bleed-III",
),
fl.Wall(
entities=volume_mesh["fluid/Bleed-IV"],
velocity=fl.SlaterPorousBleed(static_pressure=BLEED_PRESSURE_FACTORS[3] * P_INF, porosity=0.415),
name="Bleed-IV",
),
fl.Freestream(
surfaces=[volume_mesh["fluid/Inflow"], volume_mesh["fluid/Farfield"]],
name="freestream",
turbulence_quantities=fl.TurbulenceQuantities(viscosity_ratio=1.0),
),
fl.Outflow(
surfaces=[volume_mesh["fluid/Exterior-Outflow"], volume_mesh["fluid/Nozzle-Outflow"]],
spec=fl.Pressure(0.1 * P_INF),
name="outflow",
),
fl.SlipWall(
surfaces=[volume_mesh["fluid/Symmetry1"], volume_mesh["fluid/Symmetry2"]],
name="symmetry",
),
fl.Fluid(
navier_stokes_solver=fl.NavierStokesSolver(
linear_solver=fl.LinearSolver(max_iterations=35),
limit_velocity=True,
limit_pressure_density=True,
update_jacobian_frequency=1,
kappa_MUSCL=-1,
),
turbulence_model_solver=fl.SpalartAllmaras(equation_evaluation_frequency=1),
transition_model_solver=fl.TransitionModelSolver(
linear_solver=fl.LinearSolver(max_iterations=25),
absolute_tolerance=1e-7,
update_jacobian_frequency=1,
equation_evaluation_frequency=1,
N_crit=7.0,
trip_regions=[
fl.Box.from_principal_axes(
name="TripRegion",
center=(66.4, 0.0, 0.0) * fl.u.inch,
size=(80, 80, 80) * fl.u.inch,
axes=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0)),
)
],
),
),
]
return params
Outputs setup¶
The case writes field/surface/slice outputs for visualization, plus the monitor
quantities the post-processing needs: surface integrals of the mass-flow metric
on each bleed patch and the nozzle outflow, and probe rakes carrying
total-pressure ratios through the boundary layers. Volume, surface and slice
outputs carry primitiveVars and Mach. The x-normal slice stations below are
given in inches along the inlet; two z-normal slices through the symmetry plane
are added for the wedge sector.
SLICE_LOCATIONS = [
("x_slice_throat", 36.72),
("x_slice_probes", 37.5),
("x_slice_duct", 61.0),
("x_slice_outflow", 82.0),
("x_41in_slice_vg", 41.0),
("x_41p25in_slice_vg", 41.25),
("x_41p58in_slice_vg", 41.58),
("x_42in_slice_vg", 42.0),
("x_43in_slice_vg", 43.0),
("x_44in_slice_vg", 44.0),
("x_45in_slice_vg", 45.0),
]
RAKES = [
("cb_rake_27_71", [27.71, 6.13, -0.1], [0.025, 0.054, 0.091, 0.144, 0.202, 0.301], 1),
("cb_rake_32_50", [32.50, 7.191, -0.1], [0.024, 0.053, 0.090, 0.143, 0.200, 0.294], 1),
("cowl_rake_32_50", [32.50, 9.8109, -0.1], [0.021, 0.049, 0.091, 0.140, 0.198, 0.293], -1),
(
"throat_rake_bottom_37_51",
[37.51, 7.6485, -0.1],
[0.054, 0.154, 0.304, 0.508, 0.729, 0.946, 1.146, 1.296, 1.396],
1,
),
(
"throat_rake_top_37_51",
[37.51, 9.1395, -0.1],
[0.054, 0.154, 0.304, 0.508, 0.729, 0.946, 1.146, 1.296, 1.396],
-1,
),
]
def make_locations(name, base_location, heights, direction):
return [
fl.Point(
name=f"{name}_{index}",
location=[base_location[0], base_location[1] + direction * height, base_location[2]] * fl.u.inch,
)
for index, height in enumerate(heights)
]
def add_main_outputs(params, volume_mesh):
params.outputs.append(
fl.VolumeOutput(output_format=["paraview"], output_fields=["primitiveVars", "Mach"])
)
params.outputs.append(
fl.SurfaceOutput(
surfaces=[volume_mesh["*"]],
output_fields=["primitiveVars", "Mach"],
output_format=["paraview"],
)
)
for name, location in SLICE_LOCATIONS:
params.outputs.append(
fl.SliceOutput(
entities=[fl.Slice(normal=(1, 0, 0), origin=(location, 0, 0) * fl.u.inch, name=name)],
output_format=["paraview"],
output_fields=["primitiveVars", "Mach"],
)
)
params.outputs.append(
fl.SliceOutput(
entities=[fl.Slice(normal=(0, 0, 1), origin=(0, 0, -1e-6) * fl.u.inch, name="Z=-1e6_slice")],
output_format=["paraview"],
output_fields=["primitiveVars", "Mach"],
)
)
params.outputs.append(
fl.SliceOutput(
entities=[fl.Slice(normal=(0, 0, 1), origin=(0, 0, 0) * fl.u.inch, name="Z=0_slice")],
output_format=["paraview"],
output_fields=["primitiveVars", "Mach"],
)
)
return params
The monitor quantities are defined as user variables: a total-pressure coefficient/ratio (referenced to the workshop Mach), its pitot form behind a normal shock, and a normal mass-flow metric. These feed the bleed/nozzle surface integrals and the boundary-layer rake and nozzle probes.
def add_user_variables_and_integrals(params, volume_mesh):
gamma = 1.4
rho_inf = params.operating_condition.thermal_state.density
speed_of_sound = params.operating_condition.thermal_state.speed_of_sound
mach_sq = fl.solution.Mach ** 2
mach_ref_sq = WORKSHOP_MACH ** 2
total_pressure_inf = 103421.0 * fl.u.Pa
total_pressure_coeff = fl.UserVariable(
name="total_pressure_coeff",
value=(
gamma
* fl.solution.pressure
/ (rho_inf * speed_of_sound * speed_of_sound)
* (1.0 + ((gamma - 1.0) / 2.0) * mach_sq) ** (gamma / (gamma - 1.0))
- (1.0 + ((gamma - 1.0) / 2.0) * mach_ref_sq) ** (gamma / (gamma - 1.0))
)
/ ((gamma / 2.0) * mach_ref_sq),
)
total_pressure_metric = fl.UserVariable(
name="total_pressure_metric",
value=total_pressure_coeff * 0.5 * rho_inf * (WORKSHOP_MACH * speed_of_sound) ** 2 + total_pressure_inf,
)
total_pressure_ratio = fl.UserVariable(
name="total_pressure_ratio",
value=total_pressure_metric / total_pressure_inf,
)
shock_mach_sq = fl.math.max(fl.solution.Mach, 1.03) ** 2
shock_wave_coeff = fl.UserVariable(
name="shock_wave_coeff",
value=(
((gamma + 1.0) * shock_mach_sq) / (((gamma - 1.0) * shock_mach_sq) + 2.0)
)
** (gamma / (gamma - 1.0))
* (((gamma + 1.0) / ((2.0 * gamma * shock_mach_sq) - (gamma - 1.0))) ** (1.0 / (gamma - 1.0))),
)
total_pressure_pitot_ratio = fl.UserVariable(
name="total_pressure_pitot_ratio",
value=total_pressure_ratio * shock_wave_coeff,
)
mass_flow_rate_normal_metric = fl.UserVariable(
name="mass_flow_rate_normal_metric",
value=fl.math.dot(fl.solution.velocity, fl.solution.node_unit_normal) * fl.solution.density,
)
probe_fields = [total_pressure_ratio, total_pressure_pitot_ratio]
integral_fields = [mass_flow_rate_normal_metric]
for bleed_name in BLEED_NAMES:
params.outputs.append(
fl.SurfaceIntegralOutput(
name=f"bleed{bleed_name}_mfr_metric",
entities=volume_mesh[f"fluid/Bleed-{bleed_name}"],
output_fields=integral_fields,
)
)
params.outputs.append(
fl.SurfaceIntegralOutput(
name="nozzle_outflow_mfr_metric",
entities=volume_mesh["fluid/Nozzle-Outflow"],
output_fields=integral_fields,
)
)
for name, base_location, heights, direction in RAKES:
params.outputs.append(
fl.ProbeOutput(
name=name,
entities=make_locations(name, base_location, heights, direction),
output_fields=probe_fields,
)
)
params.outputs.append(
fl.ProbeOutput(
name="nozzle_outflow",
entities=fl.Point(name="nozzle_outflow", location=[82, 5, -0.1] * fl.u.inch),
output_fields=probe_fields,
)
)
return params
Simulation Params¶
make_base_params sets the reference geometry, the Mach-3 aerospace operating
condition, steady time stepping (up to 6000 steps) on a tuned adaptive CFL ramp,
and the bleed controllers. make_params then layers on the boundary/solver
models, the monitor user-variables and integrals, and the field outputs to
produce the complete SimulationParams for the case.
def make_base_params(volume_mesh):
with fl.SI_unit_system:
return fl.SimulationParams(
reference_geometry=fl.ReferenceGeometry(
moment_center=(0, 0, 0) * fl.u.inch,
moment_length=1 * fl.u.inch,
area=1 * fl.u.inch * fl.u.inch,
),
operating_condition=fl.AerospaceCondition.from_mach(
mach=WORKSHOP_MACH,
thermal_state=fl.ThermalState(temperature=INPUT_TEMP, density=INPUT_RHO),
),
time_stepping=fl.Steady(
max_steps=6000,
CFL=fl.AdaptiveCFL(max=1e8, max_relative_change=10, convergence_limiting_factor=0.4),
),
models=[],
outputs=[],
user_defined_dynamics=make_bleed_controllers(volume_mesh),
)
def make_params(volume_mesh):
params = make_base_params(volume_mesh)
params = add_models(params, volume_mesh)
params = add_user_variables_and_integrals(params, volume_mesh)
params = add_main_outputs(params, volume_mesh)
return params
Submit case¶
Build the full parameter set from the loaded mesh and submit the single Grid C
case to Flow360. The solver version can be overridden via the
SOLVER_VERSION_OVERRIDE environment variable; it defaults to the benchmark's
pinned release.
SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.11")
case = project.run_case(
name=f"gridc_mach3_sa_aft_{SOLVER_VERSION}",
params=make_params(volume_mesh),
solver_version=SOLVER_VERSION,
)
Wait for completion¶
Block until the case has finished. This can take a long time -- it is a full steady solve of up to 6000 steps on a large mesh with bleed control.
case.wait()
Postprocessing¶
The published results are two figures:
- Cowl static pressure ratio (zoomed) -- p/p0 along the cowl wall over the throat region.
- Throat rake top total pressure ratio -- boundary-layer pitot profile at the
throat_rake_top_37_51station.
Each is compared against the workshop experiment and NASA Glenn Research Center
(GRC) reference CFD. The cells below read results straight from the in-kernel
case object (no project is re-opened) and rebuild exactly these two figures.
Setup¶
Plot styling colors, VTK helpers, and a persistent scratch directory. We download the case monitors and surface outputs once and extract them; both figures are built from this local copy.
os.makedirs("results", exist_ok=True)
# Plot styling (inlined from the shared benchmark style).
FLOW360_COLOR = "#00643c" # Flow360 results
REFERENCE_COLOR = "black" # experiment
NASA_GRC_CFD_COLOR = "#C00" # NASA GRC CFD reference
VTK_SUFFIXES = (".pvtu", ".pvtp", ".vtp", ".vtu", ".vtk")
def extract_tarball(path, destination):
with tarfile.open(path) as archive:
archive.extractall(destination)
download_root = Path(tempfile.mkdtemp(prefix="paw06_results_"))
monitors_tar = download_root / "monitors.tar.gz"
surfaces_tar = download_root / "surfaces.tar.gz"
case.results.monitors.download(str(monitors_tar), overwrite=True)
case.results.surfaces.download(str(surfaces_tar), overwrite=True)
extract_tarball(monitors_tar, download_root)
extract_tarball(surfaces_tar, download_root)
Helper functions¶
File lookup, monitor-CSV reading, probe-value extraction, and PyVista surface-mesh preparation used by the figures below. These preserve the exact numerical logic of the benchmark post-processing.
def find_downloaded_file(root, filename):
matches = sorted(root.rglob(filename))
if not matches:
raise FileNotFoundError(filename)
return matches[0]
def get_monitor_df(root, stem):
return pd.read_csv(find_downloaded_file(root, f"{stem}.csv"))
def extract_probe_index(column_name):
clean = column_name.strip()
if "_Point" in clean:
return int(clean.rsplit("_Point", 1)[1].split("_", 1)[0])
digits = "".join(ch for ch in clean.rsplit("_", 2)[-2] if ch.isdigit())
return int(digits) + 1 if digits else 999999
def get_rake_pitot_values(df, rake_name):
columns = []
for column in df.columns:
clean = column.strip()
if rake_name in clean and clean.endswith("_total_pressure_pitot_ratio"):
columns.append((extract_probe_index(clean), column))
values = []
for _, column in sorted(columns):
value = df[column].iloc[-1]
value_str = str(value).strip().lower()
if value_str in ("nan", ""):
values.append(float("nan"))
continue
values.append(round(float(value), 3))
return values
def load_reference_slice(sheet, row_range, column_pair):
start, stop = row_range
x_col, y_col = column_pair
return pd.DataFrame(
{
"total_pressure_ratio": pd.to_numeric(sheet.iloc[start - 2 : stop - 2, x_col], errors="coerce"),
"height_in": pd.to_numeric(sheet.iloc[start - 2 : stop - 2, y_col], errors="coerce"),
}
)
def prepare_mesh(mesh, required_arrays=()):
if isinstance(mesh, pv.MultiBlock):
blocks = [prepare_mesh(block, required_arrays) for block in mesh if block is not None]
blocks = [block for block in blocks if block is not None and getattr(block, "n_points", 0) > 0]
if not blocks:
return pv.PolyData()
if len(blocks) == 1:
return blocks[0]
return pv.merge(blocks, merge_points=False)
if mesh is None:
return None
if required_arrays and any(name in mesh.cell_data and name not in mesh.point_data for name in required_arrays):
return mesh.cell_data_to_point_data(pass_cell_data=True)
if not mesh.point_data and mesh.cell_data:
return mesh.cell_data_to_point_data(pass_cell_data=True)
return mesh
def load_meshes(paths, required_arrays=()):
meshes = [prepare_mesh(pv.read(str(path)), required_arrays) for path in paths]
meshes = [mesh for mesh in meshes if mesh is not None and getattr(mesh, "n_points", 0) > 0]
if not meshes:
return None
if len(meshes) == 1:
return meshes[0]
return pv.merge(meshes, merge_points=False)
def find_matching_vtk_files(root, stem):
root = Path(root)
matches = []
for suffix in VTK_SUFFIXES:
exact = sorted(root.rglob(f"{stem}{suffix}"))
if exact:
return exact
matches.extend(root.rglob(f"{stem}*{suffix}"))
return sorted(set(matches))
def find_array_name(mesh, candidates):
names = list(mesh.point_data.keys()) + list(mesh.cell_data.keys())
lower_map = {name.lower(): name for name in names}
for candidate in candidates:
if candidate.lower() in lower_map:
return lower_map[candidate.lower()]
for candidate in candidates:
lowered = candidate.lower()
for name in names:
if lowered in name.lower():
return name
raise ValueError(f"Missing array {candidates[0]}")
Cowl static pressure ratio (zoomed)¶
Slice the cowl-wall surface output along the mid-plane and normalize the static pressure (p/p0, with p0 = 1/1.4) and axial coordinate (x/Rc, Rc = 10). Points below p/p0 = 2.5 are dropped: they come from the outer side of the cowl and would otherwise clutter the plot. The Flow360 profile is compared against the workshop experiment and NASA GRC CFD, zoomed to the throat region (x/Rc in [3.0, 4.5]).
def build_pressure_profile(download_root, stem):
rc = 10.0
p0 = 1 / 1.4
mesh = load_meshes(find_matching_vtk_files(download_root, stem), required_arrays=("p",))
if mesh is None or mesh.n_points == 0:
raise ValueError(stem)
extracted = mesh.slice(origin=(0, 0, -1), normal=(0, 0, 1))
pressure_name = find_array_name(extracted, ["p"])
profile = pd.DataFrame(
{
"x_over_rc": extracted.points[:, 0] / rc,
"p_over_p0": extracted.point_data[pressure_name] / p0,
"source": "Flow360",
}
)
return profile.sort_values("x_over_rc")
def save_pressure_plot(data, stem, x_label, xlim):
plt.figure(figsize=(7, 5))
for source, color, style, marker in (
("Experiment", REFERENCE_COLOR, "-", "o"),
("NASA Glenn Research Center (GRC) CFD", NASA_GRC_CFD_COLOR, "--", None),
("Flow360", FLOW360_COLOR, "-", None),
):
subset = data[data["source"] == source]
if subset.empty:
continue
plt.plot(
subset["x_over_rc"],
subset["p_over_p0"],
color=color,
linestyle=style,
marker=marker,
label=source,
)
plt.xlabel(x_label)
plt.ylabel("Static Pressure Ratio p/p0")
plt.xlim(*xlim)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"results/{stem}.png", dpi=300)
plt.show()
# Points from the wrong (outer) side of the cowl sit well below the duct pressure;
# filtering them out keeps the comparison readable.
COWL_P_MIN = 2.5
zoom_limit = (3.0, 4.5)
experiment = pd.read_csv(REFERENCE_DIR / "cowl_pop0.csv").rename(
columns={"xoRc": "x_over_rc", "pop0": "p_over_p0"}
)
experiment["source"] = "Experiment"
grc = pd.read_csv(REFERENCE_DIR / "cowl_GRC_cfd_pop0.csv").rename(
columns={"xoRc": "x_over_rc", "pop0": "p_over_p0"}
)
grc["source"] = "NASA Glenn Research Center (GRC) CFD"
flow360 = build_pressure_profile(download_root, "surface_fluid_Cowl-Walls")
cowl_data = pd.concat([experiment, grc, flow360], ignore_index=True)
cowl_data = cowl_data[cowl_data["p_over_p0"] >= COWL_P_MIN]
save_pressure_plot(
cowl_data[
(cowl_data["x_over_rc"] >= zoom_limit[0]) & (cowl_data["x_over_rc"] <= zoom_limit[1])
],
"cowl_static_pressure_ratio_comparison_zoomed",
"Cowl x/Rc",
zoom_limit,
)
Throat rake top total pressure ratio¶
The rake plot compares the Flow360 pitot total-pressure ratio (taken as the last
recorded value at each probe height) against the workshop experiment and the NASA
GRC reference CFD solution, read from the workshop rake sheet. RAKE_SPEC selects
the row/column ranges in that sheet for the published throat_rake_top_37_51
station.
RAKE_SPEC = {
"name": "throat_rake_top_37_51",
"exp_rows": (85, 94),
"exp_cols": (4, 5),
"grc_cfd_rows": (42, 227),
"grc_cfd_cols": (49, 50),
}
def save_rake_plot(data, stem):
plt.figure(figsize=(7, 5))
for series, color, linestyle, marker in (
("Experiment", REFERENCE_COLOR, "None", "o"),
("NASA Glenn Research Center (GRC) CFD", NASA_GRC_CFD_COLOR, "--", None),
("Flow360 AFT transition + SA", FLOW360_COLOR, "-", "o"),
):
subset = data[data["series"] == series]
if subset.empty:
continue
plt.plot(
subset["total_pressure_ratio"],
subset["height_in"],
color=color,
linestyle=linestyle,
marker=marker,
label=series,
)
plt.xlabel("Total Pressure Ratio")
plt.ylabel("Height (in)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig(f"results/{stem}.png", dpi=300)
plt.show()
def write_rake_comparison(download_root, spec):
# Mirror the original pd.read_excel() behavior, which consumed the first row as the header.
sheet = pd.read_csv(REFERENCE_DIR / "bl_rakes_sheet.csv", header=0)
monitor_df = get_monitor_df(download_root, f"monitor_{spec['name']}_v2")
experiment = load_reference_slice(sheet, spec["exp_rows"], spec["exp_cols"]).reset_index(drop=True)
grc_cfd = load_reference_slice(sheet, spec["grc_cfd_rows"], spec["grc_cfd_cols"]).reset_index(drop=True)
flow360 = experiment[["height_in"]].copy()
flow360["total_pressure_ratio"] = pd.Series(
get_rake_pitot_values(monitor_df, spec["name"])
).reindex(flow360.index)
flow360["series"] = "Flow360 AFT transition + SA"
plot_df = pd.concat(
[
experiment.assign(series="Experiment"),
flow360,
grc_cfd.assign(series="NASA Glenn Research Center (GRC) CFD"),
],
ignore_index=True,
).dropna(subset=["height_in", "total_pressure_ratio"])
save_rake_plot(plot_df, f"{spec['name']}_total_pressure_ratio_comparison")
write_rake_comparison(download_root, RAKE_SPEC)