Engine-Cycle Simulation¶
This tutorial goes over a complete engine cycle, using an example of a dual-expander methane/oxygen cycle. It assumes that you are already comfortable with the thrust-chamber, cooling-circuit, and regenerative heat-transfer objects introduced in Minimal Simulation.
The complete example consists of two scripts:
examples/advanced/
├── sizer_sim.py
└── post_process.py
RUN_MODE = "regen_only" is useful while developing the thrust chamber: it runs the regenerative circuits with explicit inlet boundary conditions and does not construct the complete cycle. This tutorial concentrates on RUN_MODE = "full_cycle", where the cooling circuits are embedded in a coupled pump–regenerator–turbine–injector network.
Run the simulation and report generator from the repository root:
python examples/advanced/sizer_sim.py
python examples/advanced/post_process.py
The simulation makes a result object, and saves it in results.pkl. Post-processing reads that file and writes methane_engine_report.html, as well as the plots shown in this tutorial. All the plots you see in this tutorial can be found embedded in the report html. The html reports collects the data generated by pyskyfire in an easy to navigate and flexible format. The last part of this tutorial goes through how to set up and save a report.
Define the design point¶
make_params() contains the editable numerical and physical inputs. It holds four kinds of data:
chamber design point, propellants, contour parameters, wall stack, and cooling geometry;
tank pressures, pump efficiencies and speeds, turbine efficiencies, and duct pressure ratios;
small recirculation and turbine-bypass fractions;
independent inlet conditions for the standalone regeneration mode.
def make_params():
"""Return the editable input set for this example."""
params = dict(
# Core engine parameters
p_c=100e5,
p_e=0.8e5,
MR=2.8,
AR_c=2.0,
F=100e3,
# Fuel/oxidizer parameters
cea_fu=psf.common.Fluid(type="fuel", propellants=["CH4"], fractions=[1.0]),
cea_ox=psf.common.Fluid(type="oxidizer", propellants=["O2"], fractions=[1.0]),
coolprop_fu=psf.common.Fluid(type="fuel", propellants=["methane"], fractions=[1.0],),
coolprop_ox=psf.common.Fluid(type="oxidizer", propellants=["oxygen"], fractions=[1.0],),
T_gas_fu_in=300, # Adjusted to fall within NASA CEA tables
T_gas_ox_in=300,
# Propellant tanks
p_tank_ox=5e5,
p_tank_fu=5e5,
# Chamber/nozzle parameters
theta_conv=35,
R_1f=0.5,
R_2f=1.0,
R_3f=0.5,
length_fraction=1.0, # 80 % nozzle
L_star=1.1,
# Cooling channels
copper_roughness_height=0.0030e-3,
fuel_channel_count=140,
nozzle_channel_count=120,
channel_blockage_ratio=0.1,
fuel_hot_gas_surface_area_multiplier=1.0,
ox_outbound_hot_gas_surface_area_multiplier=1.2,
ox_return_hot_gas_surface_area_multiplier=1.2,
fuel_throat_channel_height=2.0e-3,
fuel_channel_height_slope=6.0e-3,
ox_channel_height=6e-3,
ox_single_pass_channel_height=2.5e-3,
# All cooling circuits use the same uncoated copper wall.
copper_thickness=0.6e-3,
# Materials
copper=psf.common.solids.GRCop42,
# Pump parameters
eta_pump_fu=0.72,
n_fu=50000, # rpm
eta_pump_ox_stage1=0.76,
eta_pump_ox_stage2=0.72,
n_ox=25000, # rpm
# Turbine efficiencies
eta_turbine_fu=0.75,
eta_turbine_ox=0.72,
# Duct pressure ratios: fuel side
eta_pump_regen_fu=0.95,
eta_regen_turbine_fu=0.98,
eta_turbine_injector_fu=0.94,
eta_fu_injector=0.88,
# Duct pressure ratios: oxidizer side
eta_pump_regen_ox=0.98,
eta_regen_turbine_ox=0.98,
eta_turbine_injector_ox=0.88,
eta_ox_injector=0.88,
# Mass-flow leakage fractions
zeta_fu_recirc=0.005,
zeta_ox_recirc=0.005,
ox_regen_flow_fraction=0.4,
zeta_fu_turbine_bypass=0.005,
zeta_ox_turbine_bypass=0.005,
# Standalone regeneration-analysis boundary conditions, used in regen_only mode.
T_regen_fu_in=111.0,
p_regen_fu_in=150e5,
T_regen_ox_in=93.0,
p_regen_ox_in=150e5,
)
# Derived tank properties used by the full-cycle initial guesses.
params["T_tank_ox"] = CP.PropsSI("T", "P", params["p_tank_ox"], "Q", 0, params["coolprop_ox"].propellants[0],)
params["T_tank_fu"] = CP.PropsSI("T", "P", params["p_tank_fu"], "Q", 0, params["coolprop_fu"].propellants[0],)
return params
The combustion propellants and coolant propellants are both defined because the hot-gas model and coolant-property model use different backends. The parameter dictionary serves as the single source of design assumptions.
Build the thrust chamber¶
The next thing to do is to build a thrust chamber. This can be a bit of an iterative process, getting acceptable regenerative cooling while keeping pressure drops to a minimum. Use regen_only mode to get the thrust chamber right. Since this is an expander cycle engine, heat pickup in the chamber is also important. Notice how a high chamber aspect ratio increases the surface area of the chamber, and thereby the power extracted for the cycle.
Since this is a dual expander cycle, several cooling circuits are defined around the chamber for both the fuel and the oxidizer to flow through. The block below computes the combustion-gas transport model, generates the contour, defines walls and cooling circuits, and combines them into one ThrustChamber.
def setup_thrust_chamber(params):
"""Build the chamber geometry, walls, channels, and gas-side transport."""
aerothermodynamics = psf.skycea.Aerothermodynamics.from_F_pe_Lstar(
fu=params["cea_fu"],
ox=params["cea_ox"],
T_fu_in=params["T_gas_fu_in"],
T_ox_in=params["T_gas_ox_in"],
MR=params["MR"],
p_c=params["p_c"],
F=params["F"],
p_e=params["p_e"],
L_star=params["L_star"],
p_amb=1e5,
minimum_cea_temperature=305.0,
)
xs, rs = psf.regen.contour.get_contour(
V_c=aerothermodynamics.V_c,
r_t=aerothermodynamics.r_t,
area_ratio=aerothermodynamics.eps,
AR_c=params["AR_c"],
theta_conv=params["theta_conv"],
nozzle="rao",
R_1f=params["R_1f"],
R_2f=params["R_2f"],
R_3f=params["R_3f"],
length_fraction=params["length_fraction"],
)
contour = psf.regen.Contour(xs, rs, name="Methane Engine")
copper_wall = psf.regen.Wall(
material=params["copper"],
thickness=params["copper_thickness"],
)
def fuel_channel_height(x):
return (
params["fuel_throat_channel_height"]
+ params["fuel_channel_height_slope"] * abs(x)
)
ox_outbound_span = (0.11, 1.0)
ox_return_span = (1.0, 0.21)
def span_fraction_to_x(span_fraction):
if span_fraction >= 0.0:
return span_fraction * float(contour.xs[-1])
return span_fraction * -float(contour.xs[0])
ox_outbound_bounds = tuple(map(span_fraction_to_x, ox_outbound_span))
ox_return_bounds = tuple(map(span_fraction_to_x, ox_return_span))
ox_overlap_start = max(min(ox_outbound_bounds), min(ox_return_bounds))
def ox_channel_height(x):
if x < ox_overlap_start:
return params["ox_single_pass_channel_height"]
return params["ox_channel_height"]
cross_section = psf.regen.CrossSectionSquared(
blockage_ratio=params["channel_blockage_ratio"],
)
fuel_transport = psf.skycea.CoolantTransport(params["coolprop_fu"])
ox_transport = psf.skycea.CoolantTransport(params["coolprop_ox"])
fuel_channel_placement = psf.regen.SurfacePlacement(n_channel_positions=params["fuel_channel_count"],)
nozzle_channel_placement = psf.regen.SurfacePlacement(n_channel_positions=params["nozzle_channel_count"],)
fuel_chamber_pass = psf.regen.CoolingCircuit(
name="Fuel Chamber Pass",
contour=contour,
coolant_transport=fuel_transport,
cross_section=cross_section,
span=[0.1, -1.0],
placement=fuel_channel_placement,
walls=[copper_wall],
roughness=params["copper_roughness_height"],
channel_height=fuel_channel_height,
hot_gas_surface_area_multiplier=params[
"fuel_hot_gas_surface_area_multiplier"
],
)
ox_nozzle_pass_outbound = psf.regen.CoolingCircuit(
name="Oxidizer Nozzle Pass Outbound",
contour=contour,
coolant_transport=ox_transport,
cross_section=cross_section,
span=ox_outbound_span,
placement=nozzle_channel_placement,
walls=[copper_wall],
roughness=params["copper_roughness_height"],
channel_height=ox_channel_height,
hot_gas_surface_area_multiplier=params[
"ox_outbound_hot_gas_surface_area_multiplier"
],
)
ox_nozzle_pass_return = psf.regen.CoolingCircuit(
name="Oxidizer Nozzle Pass Return",
contour=contour,
coolant_transport=ox_transport,
cross_section=cross_section,
span=ox_return_span,
placement=nozzle_channel_placement,
walls=[copper_wall],
roughness=params["copper_roughness_height"],
channel_height=ox_channel_height,
hot_gas_surface_area_multiplier=params[
"ox_return_hot_gas_surface_area_multiplier"
],
)
return psf.regen.ThrustChamber(
contour=contour,
combustion_transport=aerothermodynamics,
cooling_circuits=[fuel_chamber_pass, ox_nozzle_pass_outbound, ox_nozzle_pass_return],
h_gas_corr=1.0,
h_cold_corr=1.0,
n_nodes=100,
)
In principle, the thrust chamber aerothermodynamic properties are dependent on the condition of the propellants at the inlet of the engine. It could therefore be a part of the simulation loop. However, the computational cost to this is great, and the added precision it gives is marginal. It is therefore more practical to guess the inlet conditions to the thrust chamber, and perhaps update them once the inlet conditions to the thrust chamber has been established by running the simulation.
Looking at the thrust chamber contour, one can notice the high aspect ratio in the combustion chamber section, contributing to increased heat pickup:
This example contains one fuel chamber pass and two serial oxidizer nozzle passes. The fuel enters just downstream of the throat and travels toward the injector. The oxidizer enters just downstream of the throat, travels to the nozzle exit, and then returns to its starting axial position before entering the turbine. All three circuits use uncoated copper walls. Inspect the 3D engine below:
Initial station guesses¶
A full engine cycle has fluid states at every component interface. Pyskyfire represents each state as a Station(p, T, mdot), containing pressure in Pa, temperature in K, and mass flow in kg/s.
The network needs an initial station dictionary before it can converge:
def setup_initial_stations(params, thrust_chamber):
"""Create the initial station guesses used by the full-cycle solver."""
stations = {}
mdot_fu_est = thrust_chamber.combustion_transport.mdot_fu
stations["fu_engine_in"] = psf.common.Station(params["p_tank_fu"], params["T_tank_fu"], mdot_fu_est*0.6,)
stations["fu_pump_in"] = psf.common.Station(stations["fu_engine_in"].p, stations["fu_engine_in"].T, mdot_fu_est*0.6,)
stations["fu_pump_out"] = psf.common.Station(params["p_c"] * 1.4, params["T_tank_fu"] + 10, mdot_fu_est*0.6,)
stations["fu_shaft_recirc"] = psf.common.Station(params["p_c"] * 1.4, params["T_tank_fu"] + 10, 0.1,)
stations["fu_regen_duct_in"] = psf.common.Station(params["p_c"] * 1.4, params["T_tank_fu"] + 10, mdot_fu_est*0.6,)
stations["fu_regen_in"] = psf.common.Station(params["p_c"] * 1.5, params["T_tank_fu"] + 10, mdot_fu_est*0.6,)
stations["fu_regen_out"] = psf.common.Station(params["p_c"] * 1.4, params["T_tank_fu"] + 250, mdot_fu_est*0.6,)
stations["fu_turbine_inlet_split"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_fu"] + 250, mdot_fu_est,)
stations["fu_bypass_valve"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_fu"] + 250, mdot_fu_est,)
stations["fu_turbine_in"] = psf.common.Station(params["p_c"] * 1.2, params["T_tank_fu"] + 250, mdot_fu_est,)
stations["fu_turbine_out"] = psf.common.Station(params["p_c"] * 1.1, params["T_tank_fu"] + 200, mdot_fu_est,)
stations["fu_turbine_outlet_merge"] = psf.common.Station(params["p_c"] * 1.1, params["T_tank_fu"] + 200, mdot_fu_est,)
stations["fu_injector_plenum"] = psf.common.Station(params["p_c"] * 1.1, params["T_tank_fu"] + 200, mdot_fu_est,)
stations["fu_chamber_in"] = psf.common.Station(params["p_c"], params["T_tank_fu"] + 200, mdot_fu_est,)
ox_regen_flow_fraction = float(params["ox_regen_flow_fraction"])
if not 0.0 < ox_regen_flow_fraction < 1.0:
raise ValueError("ox_regen_flow_fraction must be strictly between 0 and 1")
mdot_ox_est = thrust_chamber.combustion_transport.mdot_ox
mdot_ox_main_est = mdot_ox_est * 0.8
mdot_ox_regen_est = mdot_ox_main_est * ox_regen_flow_fraction
mdot_ox_direct_est = mdot_ox_main_est - mdot_ox_regen_est
mdot_ox_bypass_est = mdot_ox_regen_est * params["zeta_ox_turbine_bypass"]
mdot_ox_turbine_est = mdot_ox_regen_est - mdot_ox_bypass_est
stations["ox_engine_in"] = psf.common.Station(params["p_tank_ox"], params["T_tank_ox"], mdot_ox_est*0.8,)
stations["ox_pump_in"] = psf.common.Station(stations["ox_engine_in"].p, stations["ox_engine_in"].T, mdot_ox_est*0.8,)
stations["ox_stage1_pump_out"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 10, mdot_ox_est*0.8,)
stations["ox_shaft_recirc"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 10, 0.1,)
stations["ox_main_flow"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 10, mdot_ox_main_est,)
stations["ox_direct_branch"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 10, mdot_ox_direct_est,)
stations["ox_stage2_pump_in"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 10, mdot_ox_regen_est,)
stations["ox_regen_duct_in"] = psf.common.Station(params["p_c"] * 1.8, params["T_tank_ox"] + 15, mdot_ox_regen_est,)
stations["ox_regen_in"] = psf.common.Station(params["p_c"] * 1.7, params["T_tank_ox"] + 15, mdot=mdot_ox_regen_est,)
stations["ox_regen_interstage"] = psf.common.Station(params["p_c"] * 1.7, params["T_tank_ox"] + 130, mdot_ox_regen_est,)
stations["ox_regen_out"] = psf.common.Station(params["p_c"] * 1.7, params["T_tank_ox"] + 250, mdot_ox_regen_est,)
stations["ox_turbine_inlet_split"] = psf.common.Station(params["p_c"] * 1.6, params["T_tank_ox"] + 250, mdot_ox_regen_est,)
stations["ox_bypass_valve"] = psf.common.Station(params["p_c"] * 1.6, params["T_tank_ox"] + 250, mdot_ox_bypass_est,)
stations["ox_turbine_in"] = psf.common.Station(params["p_c"] * 1.6, params["T_tank_ox"] + 250, mdot_ox_turbine_est,)
stations["ox_turbine_out"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 200, mdot_ox_turbine_est,)
stations["ox_heated_branch"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 200, mdot_ox_regen_est,)
stations["ox_main_flow_merge"] = psf.common.Station(params["p_c"] * 1.3, params["T_tank_ox"] + 100, mdot_ox_main_est,)
stations["ox_injector_plenum"] = psf.common.Station(params["p_c"] * 1.1, params["T_tank_ox"] + 100, mdot_ox_main_est,)
stations["ox_chamber_in"] = psf.common.Station(params["p_c"], params["T_tank_ox"] + 100, mdot_ox_main_est,)
return stations
The station names are not merely labels. They are the links between blocks: a block reads its named inlet stations and writes its named outlet stations. For example, the fuel pump writes fu_pump_out; the recirculation splitter then reads fu_pump_out and writes fu_regen_duct_in and fu_shaft_recirc.
These values are initial guesses, not prescribed results. They should nevertheless be physically credible. A useful starting point follows the expected trend through the engine:
pressure rises across a pump and decreases across ducts, cooling circuits, turbines, and injectors;
temperature rises through regenerative cooling and decreases across a turbine;
mass flow follows the intended split, merge, recirculation, and bypass fractions.
A fixed-point solver is not a global optimiser that can reliably recover from arbitrary guesses. If an early sweep sends a fluid to an invalid thermodynamic state, makes a turbine outlet temperature non-physical, or generates an enormous pressure change, property calls and the regenerative solver may fail before the network has a chance to settle. Use rough but reasonable values; they do not need to be accurate final-cycle predictions.
Initial scalar signals¶
Stations represent flowing propellant. Signals represent scalar quantities shared between blocks but not carried by a single fluid stream: pump power, turbine power requirement, pressure drops, targets, or future control variables.
def setup_initial_signals(params):
"""Create the scalar initial guesses for the full-cycle solver."""
return {
"p_c": params["p_c"],
"P_fuel_turbine_required": 2.8e5,
"P_ox_turbine_required": 2.0e5,
}
This example begins with estimates for the fuel and oxidizer turbine power requirements. The transmission blocks later replace them with the pump loads calculated by the network. Pressure-drop signals are added separately after the blocks are created because each loss-producing block declares its own dp_key.
Create the engine network¶
The network setup first gathers the station guesses, scalar signals, and component blocks.
stations = setup_initial_stations(params, thrust_chamber)
signals = setup_initial_signals(params)
blocks = []
regen_nodes = make_regen_node_spacings(thrust_chamber)
fuel_medium = params["coolprop_fu"].propellants[0]
ox_medium = params["coolprop_ox"].propellants[0]
Each block advertises the station and signal keys it consumes and emits. The EngineNetwork executes the blocks in the order supplied by the list, merges each block’s outputs into its station and signal dictionaries, and measures the maximum relative change in the updated quantities.
At present, Pyskyfire preserves the supplied block order. It does not yet perform a dependency-aware topological sort. Arrange the blocks in a valid flow order: merge or source, pump, splitter, ducts and cooling passages, turbine, downstream ducts, injector, then transmission coupling.
Fluid and signal blocks¶
The example uses the following block types.
Block |
Role in the network |
|---|---|
|
Combines same-fluid branches by mass and enthalpy. |
|
Raises pressure to meet its required load and emits the pump shaft-power signal. |
|
Divides a stream into fixed fractions; the branch pressure and temperature are unchanged by the ideal split itself. |
|
Applies a fixed pressure ratio and an adiabatic, constant-enthalpy pressure loss. |
|
Runs a regenerative-cooling circuit using its inlet station, writes its outlet station, and emits its pressure drop. |
|
Expands the inlet stream enough to deliver the shaft-power signal requested by the transmission. |
|
Sums shaft power demands and writes the power signal consumed by the turbine. |
Fuel-side blocks and the pump load¶
The fuel-side block sequence reads much like an engine flow schematic:
# Fuel-side blocks
blocks.append(
psf.common.MassFlowMergerBlock(
name="Fuel Inlet Merge",
st_in=["fu_engine_in", "fu_shaft_recirc"],
st_out="fu_pump_in",
medium=fuel_medium,
)
)
blocks.append(
psf.common.PumpBlock(
name="Fuel Pump",
st_in="fu_pump_in",
st_out="fu_pump_out",
overcome=[
"Duct Pump-Regen Fuel",
"Regen Fuel Chamber Pass",
"Duct Regen-Turbine Fuel",
"Fuel Turbine",
"Duct Turbine-Injector Fuel",
"Fu Injector",
],
load_fraction=1.0,
p_base=params["p_c"],
input_p=params["p_tank_fu"],
eta=params["eta_pump_fu"],
n=params["n_fu"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.MassFlowSplitterBlock(
name="Fuel Recirc. Split",
st_in="fu_pump_out",
st_out=["fu_regen_duct_in", "fu_shaft_recirc"],
fractions=[
1 - params["zeta_fu_recirc"],
params["zeta_fu_recirc"],
],
medium=fuel_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Pump-Regen Fuel",
st_in="fu_regen_duct_in",
st_out="fu_regen_in",
pressure_ratio=params["eta_pump_regen_fu"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.RegenBlock(
name="Regen Fuel Chamber Pass",
st_in="fu_regen_in",
st_out="fu_regen_out",
circuit_index=0,
thrust_chamber=thrust_chamber,
medium=fuel_medium,
nodes=regen_nodes[0],
post_process_nodes=regen_nodes[0],
heat_curvature_correction=False,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Regen-Turbine Fuel",
st_in="fu_regen_out",
st_out="fu_turbine_inlet_split",
pressure_ratio=params["eta_regen_turbine_fu"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.MassFlowSplitterBlock(
name="Fuel Split Turbine Bypass",
st_in="fu_turbine_inlet_split",
st_out=["fu_turbine_in", "fu_bypass_valve"],
fractions=[
1 - params["zeta_fu_turbine_bypass"],
params["zeta_fu_turbine_bypass"],
],
medium=fuel_medium,
)
)
blocks.append(
psf.common.TurbineBlock(
name="Fuel Turbine",
st_in="fu_turbine_in",
st_out="fu_turbine_out",
P_req_key="P_fuel_turbine_required",
eta=params["eta_turbine_fu"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.MassFlowMergerBlock(
name="Merge Turbine Bypass",
st_in=["fu_turbine_out", "fu_bypass_valve"],
st_out="fu_turbine_outlet_merge",
medium=fuel_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Turbine-Injector Fuel",
st_in="fu_turbine_outlet_merge",
st_out="fu_injector_plenum",
pressure_ratio=params["eta_turbine_injector_fu"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Fu Injector",
st_in="fu_injector_plenum",
st_out="fu_chamber_in",
pressure_ratio=params["eta_fu_injector"],
medium=fuel_medium,
)
)
blocks.append(
psf.common.TransmissionBlock(
name="Fuel Shaft",
sink_keys=["P_Fuel Pump"],
source_keys=["P_fuel_turbine_required"],
)
)
PumpBlock.overcome deserves special attention. It lists the downstream blocks whose pressure losses the pump must overcome:
overcome=[
"Duct Pump-Regen Fuel",
"Regen Fuel Chamber Pass",
"Duct Regen-Turbine Fuel",
"Fuel Turbine",
"Duct Turbine-Injector Fuel",
"Fu Injector",
]
A pressure-loss block named Regen Fuel Chamber Pass writes a signal named dp_Regen Fuel Chamber Pass. The pump collects dp_<block name> for every listed item, adds the baseline chamber-side pressure requirement, and determines the target pump outlet pressure. Its required shaft power is then emitted as P_Fuel Pump.
This is somewhat clunky. The pressure path is not inferred automatically from the network graph, so modifying the cycle topology requires manually updating overcome. It remains the current solution because it keeps the block models local and the fixed-point solver simple: the pump only needs scalar pressure-drop signals, rather than a general graph traversal and algebraic loop formulation. Treat each overcome list as part of the cycle definition and review it whenever a component is added, removed, bypassed, or moved.
Oxidizer-side blocks¶
The oxidizer side uses two pump stages. Stage 1 pumps the complete oxidizer flow only high enough to supply the common injector path. After the recirculation split, ox_regen_flow_fraction selects the fraction sent through stage 2; the example uses 0.5. The remaining oxidizer bypasses the regenerative branch and proceeds directly toward the chamber.
Stage 2 pumps only the selected fraction to the higher pressure needed by the two serial nozzle cooling passes and oxidizer turbine. Its overcome list consequently contains only components in that high-pressure branch. After turbine expansion, the heated branch recombines with the direct branch before the shared downstream duct and injector. The oxidizer transmission sums the shaft powers of both pump stages when setting the turbine power requirement:
# Oxidizer-side blocks
blocks.append(
psf.common.MassFlowMergerBlock(
name="Merge Ox Recirc",
st_in=["ox_engine_in", "ox_shaft_recirc"],
st_out="ox_pump_in",
medium=ox_medium,
)
)
blocks.append(
psf.common.PumpBlock(
name="Ox Pump Stage 1",
st_in="ox_pump_in",
st_out="ox_stage1_pump_out",
overcome=[
"Duct Turbine-Injector Ox",
"Ox Injector",
],
load_fraction=1.0,
p_base=params["p_c"],
input_p=params["p_tank_ox"],
n=params["n_ox"],
eta=params["eta_pump_ox_stage1"],
medium=ox_medium,
)
)
blocks.append(
psf.common.MassFlowSplitterBlock(
name="Split Ox Recirc",
st_in="ox_stage1_pump_out",
st_out=["ox_main_flow", "ox_shaft_recirc"],
fractions=[
1 - params["zeta_ox_recirc"],
params["zeta_ox_recirc"],
],
medium=ox_medium,
)
)
blocks.append(
psf.common.MassFlowSplitterBlock(
name="Split Ox Main Flow",
st_in="ox_main_flow",
st_out=["ox_stage2_pump_in", "ox_direct_branch"],
fractions=[
params["ox_regen_flow_fraction"],
1 - params["ox_regen_flow_fraction"],
],
medium=ox_medium,
)
)
blocks.append(
psf.common.PumpBlock(
name="Ox Pump Stage 2",
st_in="ox_stage2_pump_in",
st_out="ox_regen_duct_in",
overcome=[
"Duct Pump-Regen Ox",
"Regen Ox Nozzle Outbound",
"Regen Ox Nozzle Return",
"Duct Regen-Turbine Ox",
"Ox Turbine",
],
load_fraction=1.0,
# Stage 2 supplies only the additional pressure required by the
# regenerative-cooling and turbine branch.
p_base=0.0,
input_p=0.0,
n=params["n_ox"],
eta=params["eta_pump_ox_stage2"],
medium=ox_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Pump-Regen Ox",
st_in="ox_regen_duct_in",
st_out="ox_regen_in",
pressure_ratio=params["eta_pump_regen_ox"],
medium=ox_medium,
)
)
blocks.append(
psf.common.RegenBlock(
name="Regen Ox Nozzle Outbound",
st_in="ox_regen_in",
st_out="ox_regen_interstage",
circuit_index=1,
thrust_chamber=thrust_chamber,
medium=ox_medium,
nodes=regen_nodes[1],
post_process_nodes=regen_nodes[1],
heat_curvature_correction=False,
)
)
blocks.append(
psf.common.RegenBlock(
name="Regen Ox Nozzle Return",
st_in="ox_regen_interstage",
st_out="ox_regen_out",
circuit_index=2,
thrust_chamber=thrust_chamber,
medium=ox_medium,
nodes=regen_nodes[2],
post_process_nodes=regen_nodes[2],
heat_curvature_correction=False,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Regen-Turbine Ox",
st_in="ox_regen_out",
st_out="ox_turbine_inlet_split",
pressure_ratio=params["eta_regen_turbine_ox"],
medium=ox_medium,
)
)
blocks.append(
psf.common.MassFlowSplitterBlock(
name="Ox Split Turbine Bypass",
st_in="ox_turbine_inlet_split",
st_out=["ox_turbine_in", "ox_bypass_valve"],
fractions=[
1 - params["zeta_ox_turbine_bypass"],
params["zeta_ox_turbine_bypass"],
],
medium=ox_medium,
)
)
blocks.append(
psf.common.TurbineBlock(
name="Ox Turbine",
st_in="ox_turbine_in",
st_out="ox_turbine_out",
P_req_key="P_ox_turbine_required",
eta=params["eta_turbine_ox"],
medium=ox_medium,
)
)
blocks.append(
psf.common.MassFlowMergerBlock(
name="Merge Ox Bypass",
st_in=["ox_turbine_out", "ox_bypass_valve"],
st_out="ox_heated_branch",
medium=ox_medium,
)
)
blocks.append(
psf.common.MassFlowMergerBlock(
name="Merge Ox Main Flow",
st_in=["ox_direct_branch", "ox_heated_branch"],
st_out="ox_main_flow_merge",
medium=ox_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Duct Turbine-Injector Ox",
st_in="ox_main_flow_merge",
st_out="ox_injector_plenum",
pressure_ratio=params["eta_turbine_injector_ox"],
medium=ox_medium,
)
)
blocks.append(
psf.common.SimpleDuctBlock(
name="Ox Injector",
st_in="ox_injector_plenum",
st_out="ox_chamber_in",
pressure_ratio=params["eta_ox_injector"],
medium=ox_medium,
)
)
blocks.append(
psf.common.TransmissionBlock(
name="Ox Shaft",
sink_keys=["P_Ox Pump Stage 1", "P_Ox Pump Stage 2"],
source_keys=["P_ox_turbine_required"],
)
)
This split avoids raising the entire oxidizer flow to the regenerative-circuit pressure. Change ox_regen_flow_fraction to explore the tradeoff between coolant/turbine mass flow and second-stage pump demand.
Keeping track of which block is connected where can be difficult in the script view. Therefore, pyskyfire implements visualisation to view the engine network you have created. A visualisation of the above network is shown below:
Seed pressure-drop signals and solve¶
Before the first iteration, the script creates initial pressure-drop signals from the station guesses.
# Initialise pressure-drop signals from the station guesses. # TODO internalise this process
for block in blocks:
if not hasattr(block, "dp_key"):
continue
st_in = block.station_inputs[0] if block.station_inputs else None
st_out = block.station_outputs[0] if block.station_outputs else None
if st_in in stations and st_out in stations:
dp0 = max(stations[st_in].p - stations[st_out].p, 0.0)
else:
dp0 = 0.0
signals.setdefault(block.dp_key, dp0)
This bootstrap step exists because the pumps read downstream pressure-drop signals on their first sweep, while those signals are only calculated by the loss-producing blocks during the sweep. The initial values do not need to be exact, but they should be non-negative and consistent with the station guesses.
The network is then constructed and solved:
net = psf.common.EngineNetwork(stations, signals, blocks)
net.run_fixed_point(tol=1e-3, max_iter=200)
return {
"net": net,
"stations": net.stations,
"signals": net.signals,
"residuals": net.residuals,
"block_results": net.block_results,
}
During one fixed-point sweep, every block receives the current station and signal dictionaries, computes its outputs, and overwrites the corresponding dictionary entries. Pyskyfire records the largest relative update among pressure, temperature, mass flow, and scalar signals. It repeats complete sweeps until that residual falls below tol or the iteration limit is reached.
A converged network performs one additional post-process sweep. Every block has a post_process(stations, signals) hook. Most blocks return an empty dictionary because their ordinary station result is already sufficient. RegenBlock uses this hook to rerun its cooling calculation on a detailed axial grid using the final converged inlet condition. The resulting temperature, pressure, heat-flux, and wall-temperature profiles are collected in net.block_results, keyed by block name.
This split is important: a network iteration only needs enough information to update the coupled cycle. Detailed axial arrays are more expensive and are therefore generated after convergence, once rather than once per fixed-point sweep.
Normalize cooling results and save the simulation¶
The standalone regeneration path returns named cooling data directly. The full-cycle path obtains equivalent cooling data from the regenerative blocks’ post-process outputs:
def cooling_data_from_full_cycle(block_results):
"""Normalize full-cycle block results to the common cooling-data contract."""
return {
"fuel_chamber": block_results["Regen Fuel Chamber Pass"],
"oxidizer_nozzle_outbound": block_results["Regen Ox Nozzle Outbound"],
"oxidizer_nozzle_return": block_results["Regen Ox Nozzle Return"],
}
The main function selects the mode, runs the requested calculation, and saves a portable Results object.
if RUN_MODE not in {"regen_only", "full_cycle"}:
raise ValueError(
f"RUN_MODE must be 'regen_only' or 'full_cycle', not {RUN_MODE!r}."
)
params = make_params()
thrust_chamber = setup_thrust_chamber(params)
start_time = time.time()
if RUN_MODE == "regen_only":
cooling_data = run_regen_only(params, thrust_chamber)
cycle_results = None
else:
cycle_results = engine_sizer(params, thrust_chamber)
cooling_data = cooling_data_from_full_cycle(
cycle_results["block_results"],
)
duration = time.time() - start_time
print(f"{RUN_MODE} simulation completed in {duration:.2f} s")
results = psf.common.Results()
results.add(name="mode", obj=RUN_MODE)
results.add(name="params", obj=params)
results.add(name="thrust_chamber", obj=thrust_chamber)
results.add(name="cooling_data", obj=cooling_data)
if cycle_results is not None:
for name, value in cycle_results.items():
results.add(name=name, obj=value)
if output_dir is None:
output_dir = Path(__file__).resolve().parent
output_dir.mkdir(parents=True, exist_ok=True)
output_path = output_dir / RESULTS_FILENAME
results.save(output_path)
print(f"Results saved to {output_path}")
Both modes save the common result contract:
mode
params
thrust_chamber
cooling_data
A full-cycle result also stores:
net
stations
signals
residuals
block_results
Saving results separates the expensive calculation from visualisation. You can run the cycle once, inspect the residual history, then adjust graph selection, captions, plot order, or report layout repeatedly without rerunning the thermal solver and fixed-point cycle.
Generate the report¶
post_process.py loads results.pkl and first generates report tabs that are meaningful for either run mode: input parameters, thrust-chamber geometry, cooling data, combustion properties, and through-wall thermal gradients.
def add_common_report_content(output_dir, report, params, thrust_chamber, cooling_data):
"""Add report tabs available for both regen-only and full-cycle results."""
# Parameters
tab_params = report.add_tab("Parameters")
tab_params.add_table(
params,
caption="Input Parameters",
key_title="Parameter",
value_title="Value",
precision=3,
)
tab_params.add_table(
thrust_chamber.combustion_transport.optimum,
caption="Optimal Values",
key_title="Parameter",
value_title="Value",
precision=3,
)
# Engine overview
tab_overview = report.add_tab("Engine Overview")
engine_viewer = psf.viz.make_engine_3d(
thrust_chamber,
stride=3,
show=False,
)
engine_viewer.save_html(output_dir / "engine-3d.html")
tab_overview.add_iframe(engine_viewer.data_url, caption="Engine 3D")
engine_viewer.close()
contour_plot = psf.viz.PlotContour(thrust_chamber.contour)
tab_overview.add_figure(contour_plot)
contour_plot.save_html(output_dir / "contour.html")
# Cooling data
ordered_cooling_data = list(cooling_data.values())
tab_cooling = report.add_tab("Cooling Data")
tab_cooling.add_figure(
psf.viz.PlotWallTemperature(
*ordered_cooling_data,
plot_hot=True,
plot_coolant_wall=True,
plot_interfaces=True,
)
)
tab_cooling.add_figure(psf.viz.PlotCoolantTemperature(*ordered_cooling_data))
tab_cooling.add_figure(
psf.viz.PlotCoolantPressure(cooling_data["fuel_chamber"]),
caption="Fuel coolant pressure.",
)
tab_cooling.add_figure(
psf.viz.PlotCoolantPressure(
cooling_data["oxidizer_nozzle_outbound"],
cooling_data["oxidizer_nozzle_return"],
),
caption="Oxidizer coolant pressure.",
)
tab_cooling.add_figure(psf.viz.PlotHeatFlux(*ordered_cooling_data))
tab_cooling.add_figure(psf.viz.PlotVelocity(*ordered_cooling_data))
# Thrust-chamber properties
tab_chamber = report.add_tab("Thrust Chamber Properties")
tab_chamber.add_figure(psf.viz.PlotCoolantArea(thrust_chamber))
tab_chamber.add_figure(psf.viz.PlotHydraulicDiameter(thrust_chamber))
tab_chamber.add_figure(
psf.viz.PlotRadiusOfCurvature(thrust_chamber),
caption="Radius-of-curvature computation is still experimental.",
)
tab_chamber.add_figure(psf.viz.PlotdAdxThermalHotGas(thrust_chamber))
tab_chamber.add_figure(psf.viz.PlotdAdxThermalCoolant(thrust_chamber))
tab_chamber.add_figure(psf.viz.PlotdAdxCoolantArea(thrust_chamber))
# Combustion transport
tab_combustion = report.add_tab("Combustion")
transport = thrust_chamber.combustion_transport
for prop in ("M", "gamma", "T", "p", "h", "cp", "k", "mu", "Pr", "rho", "a"):
tab_combustion.add_figure(
psf.viz.PlotTransportProperty(
transport, prop=prop, results=ordered_cooling_data
)
)
# Through-wall temperatures at three axial locations
tab_gradient = report.add_tab("Thermal Gradient")
fuel_chamber_data = cooling_data["fuel_chamber"]
profile_x = sorted(float(x) for x in fuel_chamber_data.x)
for x in (profile_x[0], 0.0, profile_x[-1]):
tab_gradient.add_figure(
psf.viz.PlotTemperatureProfile(
fuel_chamber_data,
thrust_chamber,
0,
x,
)
)
For a full-cycle result, the post-processing script adds a dedicated cycle tab containing the convergence history, station pressure/temperature/mass-flow plots, and fuel/oxidizer pressure-temperature paths.
def add_full_cycle_report_content(output_dir, report, results):
"""Add network and station plots that require a full-cycle result."""
stations = results["stations"]
residuals = results["residuals"]
net = results["net"]
tab_cycle = report.add_tab("Engine Cycle")
tab_cycle.add_figure(
psf.viz.PlotResidualHistory(residuals),
caption="Maximum relative residual per fixed-point iteration.",
)
fuel_stations = [
"fu_engine_in",
"fu_pump_in",
"fu_pump_out",
"fu_regen_in",
"fu_regen_out",
"fu_turbine_in",
"fu_turbine_out",
"fu_injector_plenum",
"fu_chamber_in",
]
oxidizer_stations = [
"ox_engine_in",
"ox_pump_in",
"ox_stage1_pump_out",
"ox_direct_branch",
"ox_stage2_pump_in",
"ox_regen_duct_in",
"ox_regen_in",
"ox_regen_out",
"ox_turbine_in",
"ox_turbine_out",
"ox_main_flow_merge",
"ox_injector_plenum",
"ox_chamber_in",
]
for property_name, caption in (
("p", "Fuel-side pressure"),
("T", "Fuel-side temperature"),
):
tab_cycle.add_figure(
psf.viz.PlotStationProperty(
station_dicts=stations,
station_list=fuel_stations,
property_name=property_name,
),
caption=caption,
)
for property_name, caption in (
("p", "Oxidizer-side pressure"),
("T", "Oxidizer-side temperature"),
):
tab_cycle.add_figure(
psf.viz.PlotStationProperty(
station_dicts=stations,
station_list=oxidizer_stations,
property_name=property_name,
),
caption=caption,
)
sankey = psf.viz.PlotMassFlowSankey(engine_network=net, title="Methane Engine Mass Flow",)
tab_cycle.add_figure(sankey)
tab_cycle.add_figure(
psf.viz.PlotPTDiagram(
station_dicts=[stations],
station_list=fuel_stations,
fluid_name="methane",
title="Fuel-side P-T path",
scale="linear",
)
)
tab_cycle.add_figure(
psf.viz.PlotPTDiagram(
station_dicts=[stations],
station_list=oxidizer_stations,
fluid_name="oxygen",
title="Oxidizer-side P-T path",
scale="linear",
)
)
NETWORK_LAYOUT_PATH = (
Path(__file__).resolve()
.with_name("methane-engine-cycle.layout.json")
)
tab_network = report.add_tab("Engine Network")
viewer = psf.viz.make_network_viz(
results["net"],
title="Methane engine cycle",
layout=NETWORK_LAYOUT_PATH,
)
viewer.save_html(path=output_dir / "network.html")
tab_network.add_iframe(
viewer.data_url,
caption="Editable engine-cycle schematic",
height="900px",
)
Finally, the script validates the saved result contract, chooses the appropriate report content from mode, and writes the HTML report.
def main(
output_dir: Path | None = None,
input_path: Path | None = None,
):
script_dir = Path(__file__).resolve().parent
if output_dir is None:
output_dir = script_dir
output_dir.mkdir(parents=True, exist_ok=True)
output_path = output_dir / REPORT_FILENAME
if input_path is None:
input_path = script_dir / RESULTS_FILENAME
results = psf.common.Results.load(input_path)
required_keys = {"mode", "params", "thrust_chamber", "cooling_data"}
missing_keys = required_keys.difference(results)
if missing_keys:
missing = ", ".join(sorted(missing_keys))
raise ValueError(
f"{input_path.name} does not use the expected result format. "
f"Missing: {missing}."
)
mode = results["mode"]
if mode not in {"regen_only", "full_cycle"}:
raise ValueError(f"Unknown result mode: {mode!r}")
print(f"Generating {mode} report from {input_path}")
report = psf.viz.Report("Methane Engine")
add_common_report_content(
output_dir=output_dir,
report=report,
params=results["params"],
thrust_chamber=results["thrust_chamber"],
cooling_data=results["cooling_data"],
)
if mode == "full_cycle":
required_cycle_keys = {"stations", "residuals"}
missing_cycle_keys = required_cycle_keys.difference(results)
if missing_cycle_keys:
missing = ", ".join(sorted(missing_cycle_keys))
raise ValueError(
f"Full-cycle result is missing required data: {missing}."
)
add_full_cycle_report_content(output_dir, report, results)
report.save_html(output_path)
print(f"Report saved to {output_path}")
A regeneration-only result produces a thrust-chamber and cooling report without an engine-cycle tab. A full-cycle result includes both the chamber analysis and the coupled-cycle diagnostics. You can view the report here: Engine Cycle Report