Skip to content

Cylinders and combustion

A reciprocating cylinder is a 0-D variable-volume junction: the slider-crank prescribes its volume, and the network integrates only its mass and internal energy. Each RK stage the cylinder sets the ends of its runner pipes (apply_port_bcs), then reports how fast its mass and energy change (accumulate_rates); once per step, after the last stage, it advances its inventories (apply_fuel_step). The contracts of every function below are in their docstrings (the standard is Function contracts and flow charts); the cylinder's own are in src/aeolus1d/components/cylinder.py.

Bases and signs

These are the places this feature has shipped, or nearly shipped, a wrong number. Each is stated once in a docstring and pinned by a test named there.

Quantity Convention
crank angle degrees CAD in the user's frame, reduced into [0, 720); TDC is at cad_at_tdc (default 0), so a firing TDC at 360 needs cad_at_tdc = 360 when the case counts that way. Radians appear only inside _volume_from_theta and _kinematics
valve window half-open, [open, close), mod 720; a window of 720 deg or more is always open; phase_offset_deg shifts window and lift lookup together (positive = earlier)
flow sign mdot and edot are the amounts into the cylinder: positive for intake (runner to cylinder), negative for blowback. On the legacy path the pipe end sets the sign (+1 at a RIGHT end, -1 at a LEFT end)
enthalpy carried the upstream side's stagnation enthalpy: the runner's cp T0 on intake, the cylinder's own specific enthalpy h on blowback
heat release the Wiebe burn is per degree; Q_dot [W] = m_fuel LHV (dx_b/dtheta_deg) omega_deg, omega_deg = omega 180/pi. A cycle releases m_fuel LHV (1 - exp(-a)), 99.3 % of nominal at a = 5
injection mdot [kg/s] = m_fuel_per_cycle shape(theta) / integral(shape) omega_deg; the integral is in shape x degrees, cached per injector
wall heat h A_wall (T_wall - T_gas) [W]: negative for a hot gas. S_p = stroke omega / pi is the mean piston speed
piston work -p dV/dt [W] on the internal energy: expansion removes it
species and fuel inventories are masses; the fraction a runner sees is m_species / volume.m

One stage: what enters, what leaves

What enters: clock.t (the stage time), the cylinder's volume.m and volume.E_internal, and the runner pipes' boundary cells. What leaves: the ghost rows or the valve stash the runner sees, and (dm, dE) for the network to integrate.

flowchart TD
    T["clock.t, or 0 with no clock: the crank angle then never moves"] --> K["_kinematics: theta, V, dV/dt (user CAD, mod 720)"]
    K --> S{"first call?"}
    S -- yes --> SEED["_seed_species_inventory, then the inert-products warning"]
    SEED -- "species_init names a column no port carries" --> X1["ValueError; the seeded flag is put back, so a retry refuses again"]
    SEED --> V
    S -- no --> V["volume.V := V(theta)  (prescribed, not integrated)"]
    V --> C{"orifice_valve_flow?"}
    C -- no --> L["legacy: per port _port_area_fraction (window, lift, VVT, Cd)"]
    L --> LO{"area fraction > 0?"}
    LO -- yes --> LG["_fill_plenum_ghost: the ghost is the cylinder state; fuel and species ghosts written"]
    LO -- no --> LW["wall_bc; fuel and species keys cleared"]
    C -- yes --> O["_apply_port_bcs_orifice (second chart)"]
    LG --> R
    LW --> R
    O --> R["accumulate_rates"]
    R --> P{"coupling"}
    P -- legacy --> PL["_pipe_port_flux x port area x area fraction, signed by the pipe end"]
    P -- orifice --> PO["_port_rates_orifice: sum of the stashed mdot, edot"]
    PL --> F
    PO --> F["_accumulate_rates_finish (third chart)"]
    F --> OUT["(dm, dE): kg/s and W, positive into the cylinder"]
    OUT --> NET["the explicit network driver: volume.m, volume.E_internal advance (Euler / SSP); orifice flows are debited from the runner cell by _valve_sources_into"]
    NET --> LAST{"last stage of the step?"}
    LAST -- no --> T
    LAST -- yes --> FS["apply_fuel_step (fourth chart)"]

The implicit network driver applies no valve sources, so an orifice cylinder with an implicit rk would gain mass and energy the walled runner never loses; the build refuses that combination (_resolve_orifice_valve_flow).

Since 2026-10-11 the cylinder refuses a dead state of its own: apply_port_bcs calls _refuse_dead_state before any ghost is written, and m <= 0, a non-finite m or E, or p <= 0 raises a RuntimeError naming the cylinder, the time and the crank angle (until then it was caught, if at all, by the runner pipes' own checks). Its construction refuses what cannot be a cylinder (the last chart).

The orifice valve coupling

With cylinder_orifice_valve_flow (the default in a built case) the runner always sees a wall at the port, and the valve flow is a compressible isentropic orifice between the runner's boundary cell and the cylinder. The same mdot and edot are added to the cylinder here and removed from the runner cell by the network, so mass and energy cross the valve exactly.

flowchart TD
    IN["theta, volume p T h gamma_eff R_eff Y_fuel"] --> PT["per port: _port_area_fraction, then wall_bc on the runner"]
    PT --> OPEN{"area fraction > 0?"}
    OPEN -- no --> CL["no stash row; the runner fuel key cleared"]
    OPEN -- yes --> RHO{"runner rho > 0?"}
    RHO -- no --> VAC["port marked closed so the stash and last_port_open stay in step"]
    RHO -- yes --> DEC["_runner_bc_index, _runner_gas: decode p, T with the runner's own gamma and R"]
    DEC --> VF["valve_flow: direction from p (or total head), orifice_mass_flux: subsonic or choked"]
    VF --> SG{"sign of mdot"}
    SG -- "positive: intake" --> UP1["upstream is the runner: edot = mdot cp T0; Y_up from the runner cell"]
    SG -- "negative: blowback" --> UP2["upstream is the cylinder: edot = mdot h; Y_up = cylinder Y_fuel"]
    UP1 --> ROW
    UP2 --> ROW["stash row: pipe, cell, V_cell from the cell-centre area, mdot, edot, Y_up, fuel column, species upwind"]
    ROW --> BK["_port_rates_orifice: dm, dE and the fuel, burnt-gas, species and pocket bookkeeping by sign"]
    BK --> REC["_record_port_transport: latest values, and the stage means apply_fuel_step will integrate"]

The stage tail: work, heat release, injection, wall heat

The four terms _accumulate_rates_finish adds to the port flow. The combustion models are exclusive, in this order; the cylinder does not warn when two are set.

flowchart TD
    IN["dm, dE from the ports; last_theta_deg, last_V, last_dVdt"] --> W["dE += -p dV/dt"]
    W --> CM{"combustion model"}
    CM -- "diesel_combustion" --> QD["Q_dot = the stored diesel value"]
    CM -- "si_combustion" --> QS["Q_dot = the stored SI value"]
    CM -- "Wiebe combustion" --> QW["_wiebe_heat_release: per degree, times omega_deg; injectors' fuel replaces m_fuel"]
    QW --> ZD["_update_two_zone_diagnostics: _burn_progress, ignition snapshot, and with two_zone the zone temperatures"]
    ZD --> ZU["_unburnt_zone_temperature, _burnt_zone_temperature (frozen, dissociation or caloric split), _bounded_T_burnt"]
    CM -- none --> Q0["no heat release"]
    QD --> QF["dE += Q_dot; last_fuel_burn_rate = Q_dot / LHV"]
    QS --> QF
    ZU --> QF
    QW --> QF
    Q0 --> INJ
    QF --> INJ{"injection schedule?"}
    INJ -- yes --> IJ["_fuel_injection_rates: dm += mdot, dE += mdot (cv T_fuel - h_vap)"]
    INJ -- no --> HT
    IJ --> HT{"heat_transfer model?"}
    HT -- yes --> WH["heat_transfer_coefficient (Woschni, Hohenberg, Chang), then dE += h A_wall (T_wall - T_gas); 3-zone split if a zone temperature is set"]
    HT -- no --> Z["wall heat = 0"]
    WH --> OUT["(dm, dE)"]
    Z --> OUT

Where a refusal happens here: WoschniHeatTransfer refuses a non-positive bore; the Wiebe path asserts a model is set; a zero duration_deg divides by zero when the angle sits exactly on the start (the case load does not check it).

Once per step: apply_fuel_step

Called by the explicit and implicit drivers, never by the quasi-steady one (which refuses cylinders). It reads the last stage's angle and the stage means the accumulators hold, so a valve that reverses between stages does not debit the cylinder for a transfer the runner was not credited with.

flowchart TD
    IN["dt, last stage state, the stage-mean accumulators"] --> PM{"predictive model"}
    PM -- diesel --> AD["_advance_diesel_combustion: Q_dot, burn rate, progress"]
    PM -- SI --> AS["_advance_si_combustion: spark latch, cycle re-arm, Q_dot, progress"]
    AD --> KN
    AS --> KN
    PM -- "none / Wiebe" --> KN["_update_knock_integral: dt / tau_ind in its window; reset at intake valve close; retard on the first knock"]
    KN --> SP["_advance_species_inventory: m_species += dt x stage-mean rate, floored at 0"]
    SP --> BI["_advance_burnt_inventory: progress increment x unburnt mass, plus port transport; volume.x_burnt"]
    BI --> CA["_sync_caloric_composition: table built once; T kept across the first switch"]
    CA --> DR["drain every stage accumulator"]
    DR --> Z{"every fuel rate zero?"}
    Z -- yes --> END["return"]
    Z -- no --> PIN["pockets: gas returns at the pocket's mean fuel fraction, capped by what entered"]
    PIN --> CAP["sinks (burn + blowback) scaled together to the fuel available"]
    CAP --> MF["volume.m_fuel += injected + flowed_in - burnt - blown_out, clamped at 0; m_fuel_burnt += burnt"]
    MF --> DP{"burnt > 0 and tracked columns?"}
    DP -- yes --> DC["_deposit_combustion_products: O2 debited, CO2 and H2O credited, scaled together if O2 is short"]
    DP -- no --> END
    DC --> END

Construction, resume and where the refusals are

flowchart TD
    A["Cylinder(...)"] --> G{"bore, stroke, clearance > 0 and rod > stroke / 2?"}
    G -- no --> X1["ValueError naming the value"]
    G -- yes --> VOL["Volume at V(theta_init), 1 bar 300 K"]
    VOL --> INJ{"injectors, or m_fuel_per_cycle > 0"}
    INJ -- "fuel but no schedule with a positive integral" --> X2["ValueError: no injection schedule"]
    INJ -- "table unusable" --> X3["NonMonotoneLiftTableError or ValueError (CyclicCadTable.build, unwrap_cyclic_cad)"]
    INJ --> OK["injector dicts completed in place; _norm cached even at 0 kg"]
    OK --> SI["set_initial_state(p, T)"]
    SI -- "p or T not > 0" --> X4["ValueError"]
    SI --> RUN["network stages (first chart)"]
    R["resume: restore_burnt_state(x_burnt, m, progress)"] -- "non-finite" --> X5["ValueError"]
    R --> RUN
    LT["CylinderPort.lift_table(), first use"] -- "length mismatch or bad abscissa" --> X6["ValueError / NonMonotoneLiftTableError, thousands of steps in unless validate.py caught it at load"]

Two known gaps were held for the maintainer on 2026-09-30 (the fix moves a bundled example's diagnostic) and fixed on 2026-10-11, their strict expected failures un-marked: with burnt_dissociation the burnt-zone temperature floored at 1200 K for an all-burnt charge and read 6000 K above the dissociation grid (_burnt_zone_temperature, now _dissociated_burnt_T), and a multi-stage Wiebe window that crosses 720 deg had heat release but no burn progress after the wrap (_burn_progress, now the single-stage reduction per stage).

Functions

  • src/aeolus1d/components/cylinder.py::Cylinder.__post_init__
  • src/aeolus1d/components/cylinder.py::Cylinder.set_initial_state
  • src/aeolus1d/components/cylinder.py::Cylinder._volume_from_theta
  • src/aeolus1d/components/cylinder.py::Cylinder._kinematics
  • src/aeolus1d/components/cylinder.py::Cylinder._is_open
  • src/aeolus1d/components/cylinder.py::Cylinder._port_area_fraction
  • src/aeolus1d/components/cylinder.py::Cylinder._seed_species_inventory
  • src/aeolus1d/components/cylinder.py::Cylinder._warn_if_combustion_products_is_inert
  • src/aeolus1d/components/cylinder.py::Cylinder._write_species_ghosts
  • src/aeolus1d/components/cylinder.py::Cylinder.apply_port_bcs
  • src/aeolus1d/components/cylinder.py::Cylinder._refuse_dead_state
  • src/aeolus1d/components/cylinder.py::Cylinder._runner_bc_index
  • src/aeolus1d/components/cylinder.py::Cylinder._runner_gas
  • src/aeolus1d/components/cylinder.py::Cylinder._apply_port_bcs_orifice
  • src/aeolus1d/components/cylinder.py::Cylinder._orifice_species_upwind
  • src/aeolus1d/components/cylinder.py::Cylinder._port_rates_orifice
  • src/aeolus1d/components/cylinder.py::Cylinder._record_port_transport
  • src/aeolus1d/components/cylinder.py::Cylinder.accumulate_rates
  • src/aeolus1d/components/cylinder.py::Cylinder._accumulate_rates_finish
  • src/aeolus1d/components/cylinder.py::Cylinder._update_two_zone_diagnostics
  • src/aeolus1d/components/cylinder.py::Cylinder._burn_progress
  • src/aeolus1d/components/cylinder.py::Cylinder._unburnt_zone_temperature
  • src/aeolus1d/components/cylinder.py::Cylinder._burnt_zone_temperature
  • src/aeolus1d/components/cylinder.py::_dissociated_burnt_T
  • src/aeolus1d/components/cylinder.py::Cylinder._burnt_T_caloric
  • src/aeolus1d/components/cylinder.py::Cylinder._bounded_T_burnt
  • src/aeolus1d/components/cylinder.py::Cylinder._update_knock_integral
  • src/aeolus1d/components/cylinder.py::Cylinder._advance_si_combustion
  • src/aeolus1d/components/cylinder.py::Cylinder._advance_diesel_combustion
  • src/aeolus1d/components/cylinder.py::Cylinder._sync_caloric_composition
  • src/aeolus1d/components/cylinder.py::Cylinder.apply_fuel_step
  • src/aeolus1d/components/cylinder.py::Cylinder.restore_burnt_state
  • src/aeolus1d/components/cylinder.py::Cylinder._advance_burnt_inventory
  • src/aeolus1d/components/cylinder.py::Cylinder._advance_species_inventory
  • src/aeolus1d/components/cylinder.py::Cylinder._deposit_combustion_products
  • src/aeolus1d/components/cylinder.py::Cylinder._fuel_injection_rates
  • src/aeolus1d/components/cylinder.py::Cylinder._wiebe_heat_release
  • src/aeolus1d/components/cylinder.py::CylinderPort.lift_table
  • src/aeolus1d/components/cylinder.py::CyclicCadTable.build
  • src/aeolus1d/components/cylinder.py::unwrap_cyclic_cad
  • src/aeolus1d/components/cylinder.py::WoschniHeatTransfer.heat_transfer_coefficient
  • src/aeolus1d/components/cylinder.py::HohenbergHeatTransfer.heat_transfer_coefficient
  • src/aeolus1d/components/cylinder.py::ChangHeatTransfer.heat_transfer_coefficient
  • src/aeolus1d/components/valve_flow.py::valve_flow
  • src/aeolus1d/components/valve_flow.py::orifice_mass_flux
  • src/aeolus1d/components/junction.py::_fill_plenum_ghost
  • src/aeolus1d/components/junction.py::_pipe_port_flux
  • src/aeolus1d/bc/base.py::wall_bc
  • src/aeolus1d/io/build.py::_resolve_orifice_valve_flow
  • src/aeolus1d/bc/base.py::fill_network_ghosts
  • src/aeolus1d/solver/network_timeloop.py::_accumulate_at_weight
  • src/aeolus1d/solver/network_timeloop.py::_valve_sources_into
  • src/aeolus1d/solver/network_timeloop.py::_ExplicitNetwork.post_rk
  • src/aeolus1d/closures/si_turbulent_combustion.py::SITurbulentCombustion.advance
  • src/aeolus1d/closures/diesel_predictive_combustion.py::DieselPredictiveCombustion.advance