Theory¶
aeolus1d solves the one-dimensional compressible-Euler equations on networks of pipes coupled through zero-dimensional volumes, with optional source-term closures (friction, wall heat transfer, porous-media drag, two-phase catalyst, multi-species reactive chemistry, surface kinetics) and an extensive boundary-condition catalogue (characteristic, valve, transient stagnation, non-reflecting LODI). This document collects the equations and references behind every model the solver implements.
The numbered citations at the bottom ([1] … [40]) are the
literature aeolus1d actually draws from — every named correlation,
every numerical scheme, every closure points back to a reference
the implementation comments cite.
Looking for a specific component?
The Component index maps every component to its model, the theory section here, a runnable example, and its tutorial — all in one table.
1. Conservation laws¶
1.1 Quasi-1D Euler¶
Let \(A(x)\) be the pipe cross-sectional area and \(\mathbf{U} = (\rho,\ \rho u,\ \rho E)^\mathsf{T}\) the conservative state. The quasi-one-dimensional Euler system in conservative-variable form is
with the Euler flux
the pressure-on-wall geometric source
and (optional) per-cell closure sources for wall friction, heat transfer, porous drag, catalyst heat exchange, and reactive-chemistry energy release (see §3). The total energy is \(E = e + \tfrac12 u^2\) and the specific internal energy \(e\) is recovered from the EOS (§1.2).
1.2 Equation of state¶
aeolus1d supports four EOSs through a kind-dispatched protocol. Each
exposes the same three primitives — pressure(T, ρ),
internal_energy(T, ρ), and sound_speed(T, ρ) — so the HLLC
flux and the per-cell pressure recovery work identically across
the choices.
Ideal gas (default):
Redlich-Kwong (RK) [22, §3.4]:
with departure-function
Peng-Robinson (PR) [23]:
In both cubic EOSs the speed of sound includes the residual correction
which is recovered by analytical differentiation of the pressure
expression — see aeolus1d.closures.eos. Here \(c_v\) is the
EOS-corrected real specific heat
\(c_{v,\text{real}} = (\partial e/\partial T)_\rho
= c_{v,\text{ideal}} + \partial e_\text{departure}/\partial T\), not
the ideal-gas reference \(c_{v,\text{ideal}}\): at dense /
near-critical states the departure term is comparable to
\(c_{v,\text{ideal}}\), so using the ideal value biases the sound
speed (audit finding #12, 2026-07). The per-pipe [pipe.eos] block
may also override \(R_\text{gas}\) independently of the bulk gas.
Tabulated (CoolProp) is loaded as a 2-D NumPy interpolation table from a saved NPZ; the protocol stays the same so the HLLC kernel is unchanged. Used for refrigerants and CO\(_2\) where the cubic EOSs are too coarse near the critical point.
For real-gas pipes the inverse \((\rho, e) \to (p, T)\) step is a 1-D Newton iteration at fixed \(\rho\) solving \(e_\text{eos}(T, \rho) = e_\text{target}\). This is the EOS-aware \((m, E_\text{internal}) \to (p, T)\) recovery used both inside the network volume update (Phase 21b) and inside the HLLC pressure recovery on each face.
2. Numerical methods¶
2.1 Finite-volume discretisation¶
The pipe is discretised into \(N\) cells of width \(\Delta x = L/N\) with \(N_g\) ghost cells on each side (\(N_g = 1\) for first-order, 2 for MUSCL, 3 for MP5). Conservative variables are stored cell-centred; fluxes live on the \(N+1\) interior faces.
The discretised quasi-1D residual per cell \(i\) is
where \(A_L, A_R\) are the bracketing face areas, \(A_c\) the cell area, the pressure-on-wall term injects only into the momentum component (\(\mathbf{e}_2\)), and the closure sources are per-unit-volume (already in \(d\mathbf{U}/dt\) form). A uniform-area pipe collapses exactly to the standard \(-(\mathbf{F}_R - \mathbf{F}_L)/\Delta x + \mathbf{S}\) residual.
2.2 Spatial reconstruction¶
Three reconstruction kinds are dispatched per pipe.
First-order (Godunov) — piecewise-constant. Used by default for any pipe with \(N_g = 1\). The interior flux is \(\mathbf{F}_{i+1/2} = \mathrm{HLLC}(\mathbf{U}_i, \mathbf{U}_{i+1})\).
MUSCL with limiter [27, Ch. 14] — second-order TVD reconstruction in primitive variables \(\mathbf{W} = (\rho, u, p)\). The limited slope is
and the face states
Two limiters are available:
- minmod — \(\mathrm{minmod}(a, b) = \mathrm{sgn}(a)\, \min(|a|, |b|)\) for sign-consistent \(a, b\); zero on opposite signs. Most diffusive but rock-solid on shocks.
- van Leer — \(2 a b / (a + b)\) (sign-consistent), zero on opposite signs. Sharper near contacts.
Achieves observed order \(\approx 2.0\) on smooth flow and drops \(\mathcal{L}_1\) error on Sod from \(\sim 1.5\,\%\) (1st order, 400 cells) to \(\sim 0.3\,\%\) (MUSCL + van Leer, 400 cells).
MP5 (Suresh-Huynh, 1997) [25] — fifth-order monotonicity- preserving reconstruction. The unlimited 5th-order upwind interpolant is
The monotonicity-preserving limiter computes
and accepts \(w_\text{orig}\) when \((w_\text{orig} - w_i)(w_\text{orig} - w_\text{MP}) \le \varepsilon\); otherwise it falls back to a TVD median value computed from the 4-argument minmod of second differences. Requires \(N_g = 3\) and gives \(\sim 50\times\) lower \(\mathcal{L}_1\) error than MUSCL at the same resolution on smooth flow while still preserving monotonicity at shocks.
2.3 Riemann solver¶
HLLC [27, Ch. 10] — three-wave approximate Riemann solver that resolves the contact discontinuity (which the simpler HLL smears). aeolus1d uses Davis wave-speed estimates:
with star-state densities and energies
and the upwind flux selected by which wave the interface lies
on. Per-cell \(\gamma\) is supported (Phase-28e) by evaluating the
sound speed and pressure recovery on each side with its own
\(\gamma\); the wave-speed estimates remain \(\gamma\)-agnostic and
the solver collapses bit-exactly to standard HLLC for uniform
\(\gamma\). Per-cell \(\gamma\) is opt-in via
[pipe.species] per_cell_gamma = true; the per-cell mapping then
defaults to the energy-consistent gamma_T_from_energy form
(\(\gamma\) recovered from the local internal energy each step),
which suppresses the Karni-style oscillation the naive
temperature-only mapping produced at composition/temperature
contacts.
HLLC-LM [40] — Thornber et al. low-Mach correction. At low Mach the standard HLLC velocity reconstruction is overly dissipative; the fix is to rescale the reconstructed face velocities toward their average:
Pressure and density pass unchanged. At \(M \ge 0.1\) the flux is
bit-exactly standard HLLC. At \(M \to 0\) the velocity jump goes to
zero, suppressing the spurious acoustic dissipation. Opt-in via
[case] hllc_low_mach = true.
2026-09-13 — the flag was inert on single-pipe runs until this date
The selector above was read and announced by the startup banner on every
driver, but only the explicit network driver's residual acted on it. A
single-pipe case setting hllc_low_mach = true ran standard HLLC while
being told it was running the corrected flux. Both explicit drivers now
honour it; the implicit driver still does not, and its banner omits the
row. See docs/changelog.md.
2.4 Time integration¶
Four explicit Runge-Kutta variants are wired through the same
stage interface (three fully-implicit drivers — implicit Euler,
BDF2, and Radau IIA(3) — are also available for CFL > 1 runs; see
aeolus1d.solver.implicit_network). The CFL-limited timestep is
with \(\mathrm{CFL} \approx 0.5\) for SSP-RK2/3 + first-order or MUSCL on shocked flow; \(\mathrm{CFL} \approx 0.3\) for MP5.
SSP-RK2 (Heun), 2nd-order TVD [28]:
SSP-RK3 (Shu-Osher), 3rd-order TVD [28]:
RK4 (classical), 4-stage 4th-order, not TVD — for smooth-flow benchmarks where the temporal-error floor is the bottleneck.
RK5 (Cash-Karp), 6-stage 5th-order, also non-TVD. Pairs naturally with MP5 reconstruction.
3. Source-term closures¶
3.1 Wall friction¶
All friction kinds expose the same Darcy-Weisbach momentum source plus mechanical-to-internal energy conversion:
Two friction-factor closures are shipped:
darcy_const — user-supplied constant \(f\).
colebrook — Colebrook-White implicit equation, evaluated
through Haaland's explicit approximation [4]:
with a linear blend in the transition window \(2300 < \mathrm{Re} < 4000\). Hydraulic diameter \(D_h\), roughness ratio \(\varepsilon/D\), and gas kinematic viscosity \(\nu = \mu/\rho\) feed \(\mathrm{Re} = |u|\,D_h / \nu\).
3.2 Wall heat transfer¶
Constant wall-temperature exchange through a \(4/D_h\) per-volume exchange area for a circular pipe:
3.2.1 Dynamic wall + coolant pipe (Phase Coolant)¶
For engine-warmup or jacket-coupled studies the wall temperature is not known a-priori — it's an ODE state driven by the gas heat flux on one side and a coolant loop on the other. Three state layers per cell:
Energy is conserved per cell: every joule the gas loses across
the wall is gained by the wall + coolant. The coolant pipe is
1-D incompressible (pump-prescribed \(u_c, \rho_c, c_c\); no
pressure / momentum state — engine-cooling regime). Implementation:
aeolus1d.components.coolant_pipe.CoolantPipe /
ThermalWall, rate functions in
aeolus1d.closures.wall_coolant. The coupling is registered via
top-level [[coolant_pipe]] + [[coolant_coupling]] blocks
in TOML; advance_network(coolant_couplings=...) advances them
forward-Euler-style once per outer gas step (coolant time constants
≫ acoustic dt, so first order is well within the gas-side
discretisation error).
The 0-D lumped equivalent (aeolus1d.coolant.CoolantState from
Phase E3) is also available for bulk-block warmup studies where
spatial resolution along the jacket isn't needed; pick the 1-D
CoolantPipe when temperature gradients along the cylinder head
matter.
3.3 Porous media¶
Darcy-Forchheimer two-term drag [12]:
The first term is linear (Darcy regime, \(\mathrm{Re}_K \to 0\)); the second is the inertial Forchheimer correction that dominates at higher flow speeds.
3.4 Two-phase catalyst¶
This two-phase (gas + solid) channel model is the standard 1+1D flow-through monolith formulation; Millo et al. [41, §11.2] give the canonical derivation (gas-phase species/energy, solid-phase energy with axial conduction, and washcoat species balances), which aeolus1d's catalyst closures follow.
A solid-phase temperature \(T_s(x, t)\) is integrated alongside the gas state, coupled by Newton-cooling heat exchange:
with \(\rho_s c_s\) the volumetric heat capacity of the substrate and \(a_v\) the gas-solid surface area per unit reactor volume. Three calibration kinds are exposed:
Constant-h (Phase 15): user supplies \((\rho_s c_s,\ a_v,\ h)\) directly.
Monolith resolver (Phase 45a) [13, §4]: derives \((a_v, h)\)
from honeycomb geometry. With cpsi cells per square inch and
open-frontal-area ofa,
with a constant Nusselt number \(\mathrm{Nu} = 2.98\) (laminar fully-developed square channel).
Variable-Nu (Phase 45b): replaces the constant Nu with the Hawthorn correlation
so the heat-transfer coefficient tracks local Reynolds number.
Useful for transient cold-start where the thermal entry length
matters. A sieder_tate correlation
\(\mathrm{Nu} = 0.027\,\mathrm{Re}^{0.8}\,\mathrm{Pr}^{1/3}\) is also
available. Unlike Hawthorn (whose additive \(2.98\) baseline is a built-in
floor), the Sieder-Tate form vanishes as \(\mathrm{Re}\to 0\), severing
the gas↔solid heat path at idle / soak / cool-down. The opt-in
catalyst_nu_laminar_floor case flag floors it at the channel's
fully-developed-laminar asymptote — the Graetz-Nusselt limit, \(2.98\) for
a square channel (Millo et al. [41, Table 11.1]: \(2.977\) constant-\(T\) /
\(3.089\) constant-flux; the constant-\(T\) limit applies after light-off,
the constant-flux limit before). Default off (validation-sensitive): see
Run settings.
Mechanism-driven catalyst (Phase 46d): the catalyst kind
mechanism reads a CHEMKIN-Surface or Cantera YAML file and
runs a full surface-mechanism step on every cell. Surface
species coverages \(\theta_k\) evolve via mass-action rates;
the gas-phase species feel both adsorption fluxes and
surface-reaction products. The mechanism-load path is built on
the existing CHEMKIN parser (Phase 27 [33]).
3.5 Species transport¶
A primary advected species \(Y(x, t)\) (mass fraction) is transported with the gas mass flux:
where \(r\) is the per-mass reaction rate from the chosen reaction kind.
For multi-species cases (Phase 22), an additional column array \(\mathbf{Y}_\text{extras}(x, t)\) of shape \((N, n_\text{extra})\) advects passively alongside \(Y\). Species are upwind-only (HLLC mass-flux direction selects which side's \(Y\) is transported) so the species-transport kernel is bit-exact mass-conservative in tandem with the gas-mass update.
Scheme-consistent species mass flux. The per-face mass flux that upwinds the species must be the same flux the density \(\rho\) was advanced with — the spatial scheme's reconstructed-state HLLC flux, not a first-order scalar-\(\gamma\) HLLC of the raw cell centres. On a MUSCL / MP5 / variable-\(\gamma\) / real-gas-EOS pipe the two differ, and a species advected on the wrong flux drifts even when the composition is uniform (it sees a different \(\dot m\) than the \(\rho\) it normalises against). The scheme flux is recovered without re-running the Riemann solve by telescoping the mass residual the scheme already produced: on a constant-area pipe
so one cumulative sum, anchored at the boundary-face flux \(\dot m_0\),
reconstructs \(\dot m\) at every interior face. The anchor must itself be
the scheme's reconstructed face-0 flux: a first-order anchor telescopes
the interior faces off the wrong constant and biases every face by
\(\dot m_0^\text{true} - \dot m_0^\text{anchor}\) — a DC shift that fakes a
uniform advective drift of any non-uniform composition at a
prescribed-inflow boundary (it is identically zero for uniform
composition, for walls where \(\dot m = 0\), and for first_order pipes,
where the raw-cell HLLC already is the scheme flux). The boundary-face
flux is rebuilt with the scheme's own reconstruction for every MUSCL and
MP5 variant — scalar-\(\gamma\), variable-\(\gamma\), and real-gas EOS — so
no spatial scheme is left on the residual-inconsistent first-order anchor;
first_order pipes keep the raw-cell anchor, where it already is the
scheme flux. The matching boundary kernels
(boundary_face_flux_{muscl,mp5}{,_var_gamma}_jit,
boundary_face_flux_eos_{muscl,mp5}_jit) each mirror their residual
kernel's face reconstruction, and the same helpers anchor the
flux-conservative junction port-rate (junction._pipe_port_flux).
(aeolus1d.solver.network_timeloop._species_rate, _left_boundary_mass_flux.)
Inflow-boundary composition. At a boundary with a prescribed
species inflow (a mass / velocity / stagnation inlet, or a junction
coupling) the inlet ghost is the feed, not a zero-gradient copy of
the edge cell. This matters for a reaction product the feed omits
(e.g. CO₂/H₂O over an oxidation brick): a zero-gradient ghost mirrors
cell 0, so the advective in-flux \(\dot m\,Y_0\) cancels the out-flux
\(\dot m\,Y_0\) and the product can never be flushed at the inflow
face — chemistry then injects it every step with no advective removal,
\(\sum Y\) drifts above 1, and mass is no longer conserved. A
complete feed (the specified mass fractions sum to \(\approx 1\)) is
treated as the full incoming composition, so omitted species are
absent (ghost \(0\)) and products are flushed; a partial trickle
feed (sum \(\ll 1\), a dilute-tracer pipe whose carrier gas is not in
the feed) keeps the zero-gradient ghost for the unspecified remainder.
The split is the shared _inflow_ghost_extras helper, applied
identically by the explicit, implicit-Newton, and quasi-steady species
transports. See ADR-0009.
Five reaction kinds are wired through the species seam.
arrhenius_scr_demo (Phase 15b) — single-species
first-order Arrhenius for NH3 → products:
arrhenius_scr_multi_demo (Phase 22b) — two-species standard
SCR with 1:1 stoichiometry [14, 15]:
nh3_ads_des (Phase 49a) — Langmuir adsorption / desorption
on a single storage site [16]. The catalyst surface coverage
\(\theta_\text{NH3} \in [0, 1]\) is integrated alongside the gas
state by the chemistry sub-step (§7):
with \(W_\text{cat}\) the total NH3-site density (mol-sites / m\(^3\)).
nh3_scr_kinetics (Phase 49a) — adds NH3 oxidation and
standard SCR on top of ads/des. Both new steps consume adsorbed
ammonia (Eley-Rideal):
Stoichiometry: \(4\,\text{NH}_3(s) + 3\,\text{O}_2 \to 2\,\text{N}_2 + 6\,\text{H}_2\text{O}\), \(4\,\text{NH}_3(s) + 4\,\text{NO} + \text{O}_2 \to 4\,\text{N}_2 + 6\,\text{H}_2\text{O}\).
Pollutant emissions (Phase 25m): two engine-relevant post-flame closures, both fed by the multi-species advection infrastructure.
- Extended Zeldovich NOx [11, §11.3.1]:
$$ \frac{\mathrm{d}[\text{NO}]}{\mathrm{d}t} = 2\,k_{1f}(T)\,[\text{N}_2]\,[\text{O}] \ [\text{mol}/(\text{m}^3 \cdot \text{s})] $$
$$ k_{1f}(T) = A_{1f}\,\exp(-E_{a,1f}/RT), \quad [\text{O}] = K_O(T)\,\sqrt{[\text{O}_2]} $$
$$ K_O(T) = K_{O,A}\,\exp(-E_{a,K_O}/RT). $$
The O-atom concentration is taken from the partial equilibrium
\(\tfrac12\text{O}_2 \rightleftharpoons \text{O}\), whose
equilibrium constant is
\(K_O = [\text{O}]/[\text{O}_2]^{1/2}\), so \(K_O\) multiplies
\(\sqrt{[\text{O}_2]}\) directly (the Heywood correlation form; a
\(\sqrt{K_O\,[\text{O}_2]}\) grouping would instead treat \(K_O\)
as the full \(\text{O}_2\rightleftharpoons 2\text{O}\) constant and
halve the effective activation temperature — corrected in the
2026-07 chemistry audit). Since 2026-09-29 \(K_O\) is a van't Hoff
fit to the repo's NASA-7 constant (\(K_{O,A} = 4.507\times10^3\)
(mol/m³)^½, \(E_{a,K_O} = 246.52\) kJ/mol, within 1 % over
1600-3000 K) and \(k_{1f}\) is Heywood's Table 11.1
\(7.6\times10^{13}\exp(-38000/T)\) cm³/(mol·s). The previous pair
(258.5 kJ/mol with \(K_{O,A} = 3.6\times10^3\), and the Bowman/Baulch
\(k_{1f}\), ~2.4× Heywood's) put [O] at 0.36-0.50× thermodynamic and
compensated in the thermal rate only. Quasi-steady on N atoms,
reverse reactions neglected — valid for lean post-flame gases. A
full extended-Zeldovich variant (make_zeldovich_nox_extended,
N₂O path + relaxation toward equilibrium) and a per-cycle
post-processing NO integrator (aeolus1d.analysis.nox_post, whose
stiff solve_ivp step escalates BDF → LSODA → Radau and falls
back to an equilibrium-limited envelope) build on the same rate.
The N₂O pair (\(k_4\), \(K_{\text{N}_2\text{O}}\)) was re-derived on
2026-08-26 — see the closure docstring; the integrator caps
max_step at the caller's own sample spacing, because a rate that
is exactly zero outside the active window otherwise lets BDF step
over a ~1 ms burn and converge to an all-zero history.
- Hiroyasu-Kadota soot formation/oxidation [17, 11, §11.3.2]:
$$ \frac{\mathrm{d}\rho_\text{soot}}{\mathrm{d}t} = R_f - R_o $$
$$ R_f = A_f\,\rho_\text{fv}\,p^{0.5}\,\exp(-E_{a,f}/RT) $$
$$ R_o = A_o\,\rho_\text{soot}\,X_{\text{O}2}\,\frac{p}{\sqrt{T}} \,\exp(-E{a,o}/RT). $$
3.6 Junction-port loss¶
Optional Bassett-style port loss for Corberán junctions (§5.2):
with \(K_\text{loss}\) calibrated from datasheet measurements (typical range 0–2 for sharp-edged ports).
4. Boundary conditions¶
All boundary conditions are ghost-cell fillers with signature
bc(pipe, end, t) -> None. The interior cells then sample
U_pad[Ng-1] / U_pad[N+Ng] as if they were ordinary neighbours,
and HLLC at the boundary face produces the correct numerical
flux for that BC kind.
4.1 Wall and transmissive¶
Wall — exact reflection. \(\rho\) and \(E\) mirrored across the boundary; momentum flipped: \((\rho u)_\text{ghost} = -(\rho u)_\text{int}\).
Transmissive — zero-gradient extrapolation. Crude non-reflecting; lets waves leave the domain at the cost of some reflection on strong shocks.
4.2 Stagnation inlet¶
Subsonic-inlet driven by stagnation \((p_0, T_0)\). The outgoing characteristic \(R^- = u - 2 a/(\gamma - 1)\) carries information from the interior; the incoming reservoir provides the isentropic stagnation closure. Solving the quadratic in boundary sound speed \(a_b\):
gives the boundary state, after which
Outflow through the inlet. When the interior pushes gas out (backflow), the end is an open end discharging into the reservoir: static \(p_b = p_0\), the interior entropy and the outgoing characteristic — the §4.3 construction with \(p_\text{back} = p_0\) — or plain extrapolation if that outflow is supersonic. The switch is taken when that outflow state points out of the pipe. Before 2026-09-26 the inflow relation above was applied to outflow as well, which made the exit a lossless diffuser: its ghost pressure fell as the backflow sped up, and the backflow ran away. When the quadratic has no positive root the boundary is the reservoir state at rest.
The valve_inlet BC (Phase 25 [29]) wraps this with a
crank-angle-windowed lift schedule:
\(f(t) = \mathrm{lift}(t) / \mathrm{lift}_\text{max}\) modulates
the flow from full stagnation inlet down to a wall when
\(f \to 0\).
4.3 Pressure outlet¶
Subsonic-outlet imposing static \(p_\text{back}\). The entropy and the \(R^+\) characteristic carry from the interior:
A supersonic-outlet branch falls back to extrapolation.
4.4 Non-reflecting (LODI)¶
Local one-dimensional inviscid characteristic BC [21]. The incoming Riemann invariant is pinned to a user-supplied freestream \((\rho_\infty, u_\infty, p_\infty)\) and the outgoing invariant is read from the interior:
The same construction is used as a port boundary on junctions
when port.lodi = true (Phase 23): the pipe-side outgoing
invariant pins the boundary state and the plenum's \(p_V\)
provides the pressure reference.
4.5 Transient stagnation / pressure outlet¶
The same characteristic-based BCs as §4.2/4.3 but with
\(p_0(t),\ T_0(t),\) or \(p_\text{back}(t)\) interpolated from a
user-supplied piecewise-linear schedule (p0_t, p0_v). The
solver clock \(t\) advances with the timeloop and feeds into the
BC's lookup.
4.6 Valve orifice (Winterbone-Pearson)¶
Phase 25b orifice-flow valve [30, Ch. 6]. The effective area modulates with the lift schedule:
The boundary mass flow then follows compressible-orifice flow:
- Choked (\(\mathrm{PR} < \mathrm{PR}_\text{crit}\)):
$$ \dot m = A_\text{eff}\,p_0\,\sqrt{\frac{\gamma}{R\,T_0}} \cdot \left(\frac{2}{\gamma + 1}\right)^{!(\gamma+1)/(2(\gamma-1))}. $$
- Subsonic (\(\mathrm{PR} \ge \mathrm{PR}_\text{crit}\)):
$$ \dot m = A_\text{eff}\,p_0\, \sqrt{\frac{2\gamma}{(\gamma-1)\,R\,T_0}} \cdot \sqrt{\mathrm{PR}^{2/\gamma} - \mathrm{PR}^{(\gamma+1)/\gamma}}. $$
The boundary velocity recovers from \(\dot m\) and the geometric \(A_\text{pipe}\); temperature follows the isentropic relation \(T_b = T_0\,(p_b/p_0)^{(\gamma-1)/\gamma}\).
4.7 Mass-flow and velocity inlets¶
Two inlet variants for cases where the upstream actuator imposes a flow quantity rather than a pressure reservoir.
Mass-flow inlet — given \((\dot m,\ T_0)\), reuse the same \(R^-\) characteristic + isentropic-stagnation closure as §4.2 to obtain \((u_b,\ a_b,\ T_b)\), then close \(\rho_b\) from the mass-flow constraint at the boundary face:
The sign of \(u_b\) is fixed by the characteristic so the BC works on either pipe end. As \(|u_b| \to 0\) the closure degenerates; the implementation clamps \(|u_b|\) to a small floor to keep the ghost finite — review the run if it triggers, the upstream actuator is asking for more mass than the interior can accept.
Velocity inlet — given \((u,\ T_0,\ p_0)\), pure Dirichlet ghost via total-enthalpy + isentropic-stagnation:
No characteristic from the interior — this can introduce mild wave reflections compared to §4.2 / mass-flow inlet. Acceptable trade for piston / blower / fan boundaries where \(u\) is the controlled quantity.
For real-gas pipes both BCs route the conservative ghost energy
through eos_internal_energy(T_b, ρ_b, …) rather than the ideal-
gas closed form, so the ghost stays consistent with the EOS-aware
interior (mirrors the Phase-21 junction ghost-fill pattern).
4.8 Transient inlets with species schedules¶
The transient variants (§4.5) of the stagnation, mass-flow, and
velocity inlets all accept optional species schedules so the
inlet's gas composition can vary in time. For a primary species
(e.g. \(Y_{\mathrm{NH_3}}\)) the user supplies
\((Y_t,\ Y_v)\) — piecewise-linear in seconds — and the BC
interpolates the value at each timestep. For each extra column on
Pipe.Y_extras (e.g. \(Y_{\mathrm{NO}}\)) a
BCExtraSpeciesSchedule(name, t, v) entry attaches its own
schedule.
The BC writes the interpolated values to
Pipe.species_inflow_left (or _right) — a small dict keyed by
species name — and the species advection consults it before
applying the upwind selector at the inlet face. Reverse flow is
handled cleanly: the upwinder picks the interior \(Y\) over the
ghost-imposed \(Y\) when the boundary mass flux runs outward, so
the schedule has no effect during back-flow, matching the
characteristic-BC convention.
5. Network coupling¶
Pipes are joined through 0-D plena. Each junction kind exposes
the same port_flux(pipe_state) interface so the time loop can
mix and match.
5.1 Plenum (mass + energy ODE)¶
A Volume carries a conservative \((m, E_\text{internal})\)
state. The time derivatives are
where the summed contributions are signed mass flows and total specific enthalpies on each port (positive = into the volume). The per-port flux is computed by filling the pipe-side ghost and calling HLLC at the port face.
Port coupling: the plenum is a reservoir mouth¶
Until 2026-08-02 the ghost was the plenum rest state \((\rho_V, 0, p_V)\), and this section claimed that was "exact for mean-flow mass and energy". It is not, and the error is large.
Writing \(u = 0\) puts a velocity jump on the port face. A Godunov face between a state at rest and a moving state does not resolve a contraction — it resolves a piston, whose pressure response is the acoustic impedance \(\rho a \Delta u\), not the dynamic head \(\tfrac12 \rho u^2\). The two differ by \(2/M\), so at the \(M \approx 0.13\) of an engine manifold the port invented roughly 15× the correct loss. Feeding the port the closed-form steady reservoir-discharge state and asking whether it is a fixed point measures it directly: the old ghost returned a mass flux 50 % low and a momentum residual of 3.6–12.8 dynamic heads, at every subsonic pressure ratio.
The port is now treated as what it physically is, with the flow direction deciding whether the dynamic head is recovered or destroyed. Feed (volume → pipe) — the volume is a stagnation reservoir and gas accelerates from rest isentropically, so the port static pressure sits below \(p_V\):
with \(a_0^2 = \gamma p_V/\rho_V\), clamped at \(a^2 = a_0^2\,\tfrac{2}{\gamma+1}\) — the port chokes at \(M = 1\). Dump (pipe → volume) — a sudden expansion, whose dynamic head is dissipated inside the volume rather than at the face, so \(p = p_V\) and the gas keeps its own entropy, \(\rho = \rho_i (p_V/p_i)^{1/\gamma}\).
Both branches carry \(u = u_\text{int}\): no velocity jump, hence no spurious \(\rho a \Delta u\). Both make the closed-form steady state an exact fixed point (residual \(\le 8\times10^{-16}\) relative). At \(u_\text{int} = 0\) the feed branch returns \((\rho_V, 0, p_V)\) bit-identically, so a network at rest is unchanged.
This completes ADR-0008, which found the same zero-velocity ghost on the cylinder ports ("~60 % of the correct orifice flow"), replaced it there with a compressible-orifice model, and left it in place on plenum ports. The two couplings now assume the same physics. See ADR-0012.
For real-gas plena (Phase 21b), \(p_V, T_V\) are recovered from \((m, E_\text{int})\) via the EOS Newton inverse; the ghost fill uses those. The isentropic relations above are the ideal-gas forms evaluated at the pipe's \(\gamma\) — the same approximation class as the \(p/(\gamma-1)\) they replaced — while the ghost energy still routes through the pipe's EOS.
A neck-pipe correction (Phase 14) inserts a short internal FV pipe at the port; its length defaults to the Rayleigh end correction \(L_\text{neck} = 0.85\,\sqrt{A_\text{port} / \pi}\) for a flanged opening [1, §6.6].
5.2 Corberán CP+h junction¶
Algebraic constant-pressure + enthalpy-balance junction with no chamber state [3]. For an \(N\)-port junction, the closure system is
\(p\) at every port equal — the CP closure — solved by Riemann-invariant matching across each port. Faster than plenum on short merges where chamber dynamics are artificial. Per-port Bassett \(K_\text{loss}\) (§3.6) is applied as a pressure-drop adjustment on each port's invariant.
5.3 Extended Helmholtz Resonator (EHR)¶
A plenum with an additional per-port lumped neck momentum ODE (Phase 14d) [2]:
where \(L_\text{eff}\) is the effective neck length (geometric + Rayleigh end correction). The integrated \(u_n\) then drives the port mass flow \(\dot m = \rho_\text{up}\,A_\text{port}\,u_n\). This is the only junction kind that gives the correct Helmholtz mode on a 0-D chamber.
5.4 Lumped Helmholtz resonator (textbook 0-D primitive)¶
For design-stage sizing only (HelmholtzResonator class, used
in benches and examples), aeolus1d exposes the textbook 2-DOF
resonator with the analytical natural frequency [1, §6.6]:
This is not a solver component — the time-domain network solver uses the EHR or quasi-2D Helmholtz path instead.
6. Engine components¶
6.1 Reciprocating cylinder (slider-crank)¶
Phase 25g cylinder volume [11, §2.3]. Piston displacement and volume are
with \(R = \text{stroke}/2\), \(L = \text{con\_rod\_length}\), \(D = \text{bore}\). The cylinder is treated as a 0-D volume with its own \((m, E_\text{int})\) state plus the moving boundary; the \(-p\,\mathrm{d}V/\mathrm{d}t\) work term is carried into the energy ODE.
6.2 Wiebe combustion¶
Phase 25h heat-release [11, §9.2.2]:
The Phase 25j two-zone variant tracks burnt and unburnt zone masses + temperatures separately; per-zone \(\gamma\) tracking (Phase 25k) lets the variable-\(\gamma\) stack feel the composition shift across the flame.
Wall heat loss uses the Woschni correlation [38]:
with mean piston speed \(\bar S_p\). The Hohenberg correlation [39] is offered as an alternative. Both feed \(\dot Q_\text{wall} = h\,A_\text{wall}(\theta)\,(T_\text{wall} - T_\text{gas})\).
Two-zone dissociation cap (opt-in). A frozen-composition
burnt zone overpredicts the peak temperature because it ignores
the endothermic dissociation of the combustion products. With
WiebeCombustion.burnt_dissociation = true the burnt-zone energy
balance solves \(c_{v,b}\,T_b + \Delta e_\text{dissoc}(T_b) =
e_{b,\text{target}}\), where \(\Delta e_\text{dissoc}\) is a
NASA-7-derived sensible-energy augmentation (CO/CO₂ dissociation
≈ 5 % at 2400 K, ≈ 21 % at 2800 K) that absorbs energy at high
\(T\) and caps the burnt-zone temperature. Default off
(bit-identical when unset).
Multi-stage Wiebe. MultiStageWiebe (aeolus1d.closures.multi_stage_wiebe)
superposes two or three WiebeStage sub-profiles — each with its own
\((a, m, \theta_\text{start}, \Delta\theta)\) and a mass_fraction
weight (summing to 1.0) — for HCCI / RCCI low-temperature +
main-burn (+ optional late soot-oxidation) heat-release shapes. The
composite \(\mathrm{d}x_b/\mathrm{d}\theta\) is the mass-weighted sum
of the per-stage derivatives.
Knock (Livengood-Wu). A knock-prediction integral accumulates the inverse ignition delay over the compression window \([\text{IVC},\ \theta_\text{start}]\) each cycle:
with Douaud-Eyzat defaults \(A_k = 17.68\,\text{ms}\), \(n_k = 1.7\),
\(B_k = 3800\,\text{K}\). Knock is predicted when \(I_\text{knock}\)
reaches 1.0 ahead of combustion start (the end-gas would auto-ignite
before the spark / injection event). The integral resets on IVC and
only accumulates over the pre-ignition window (integrating through the
burned gas would fire spuriously every cycle). By default the event
feeds a closed-loop spark-retard controller
(knock_retard_deg_per_event = 2.0 CAD/event, capped at
knock_retard_max_deg = 10.0, recovering knock_advance_deg_per_cycle
= 0.5 CAD/cycle); set both retard fields to 0 for a diagnostic-only
integral with no phasing feedback.
Chen-Flynn FMEP. When a [cylinder.fmep] block attaches
FMEPCoeffs, the engine-metrics post-processor evaluates the classical
Chen-Flynn friction correlation
so brake work follows from indicated work. With no coefficients attached (the legacy default) friction is unmodelled (\(\text{FMEP} = 0\), \(\text{BMEP} = \text{IMEP}\)).
Compressible-orifice valve coupling (ADR-0008). With the
default-on cylinder_orifice_valve_flow case flag, each open
[[cylinder.port]] couples to its runner pipe through a compressible
isentropic-orifice mass/energy flow (aeolus1d.components.valve_flow)
rather than the legacy zero-velocity-ghost + \(A_\text{eff}\)-credit
coupling that leaks mass. The pipe face is set to a wall and the
orifice \(\dot m / \dot E\) is applied symmetrically (added to the
cylinder, removed from the runner boundary cell by the timeloop's
valve-source step), which recovered engine volumetric efficiency from
~6 % to ~74 % on the EPW diesel harness. Inflow-ghost species
conservation follows ADR-0009 (§3.5).
Charge composition. Beyond \((m, E_\text{int})\) the cylinder carries a well-mixed mass inventory \(m_k\) of every species column its ports touch, driven by the same port mass flux that drives \(\mathrm{d}m/\mathrm{d}t\) with an upwind composition,
so gas leaving a valve carries the charge composition into the runner instead of the runner's own (the zero-gradient ghost default, which mirrors the pipe back at itself). Combustion converts the burnt fuel mass into complete-combustion products of a \(\mathrm{CH}_y\) fuel,
with \(y\) the fuel's H/C atom ratio (fuel_H_C_ratio), giving per-kilogram
mass yields \((1 + y/4)M_{\mathrm{O_2}}/M_{\mathrm{CH}_y}\) consumed and
\(M_{\mathrm{CO_2}}/M_{\mathrm{CH}_y}\), \((y/2)M_{\mathrm{H_2O}}/M_{\mathrm{CH}_y}\)
produced. The stoichiometric air-fuel ratio is then an output,
\(\text{AFR}_\text{st} = Y_{\mathrm{O_2}}/0.232\); at the diesel default
\(y = 1.85\) it returns 14.54, matching the 14.5 used by the predictive
diesel closure. Single-zone and complete-combustion only — see
scope-and-limitations for what that excludes.
The inventory is inert on a case whose pipes carry no non-fuel column.
6.3 Compressors and turbines¶
Phase 25c/d quasi-steady map-driven turbomachinery [31, 32].
The map-physics layer reads corrected mass flow and isentropic efficiency from a \((\mathrm{PR}, N_\text{corr})\) lookup:
Three map kinds: synthetic_centrifugal (analytical PR ↔ ṁ
factory), synthetic_turbine (Stodola-style swallowing
capacity), synthetic_vgt_turbine /
synthetic_igv_centrifugal (3-D variants with a vane_open
axis added on top of the (PR, N_corr) plane), from_csv /
from_npz (load a tabulated 2-D or 3-D map from disk — see
Turbo map file format for the
expected columns / arrays / units), and learned_npz
(LearnedTurboMap MLP surrogate).
When two turbo components share a Shaft, Phase 25d couples
them through the inertia ODE:
with a floor \(\omega_\text{min}\) (default 50 rad/s) to avoid
the singularity at zero speed. N independent shafts with N×N
turbomachines compose without code change — see
examples/60_twin_sequential_turbo.toml for a HP+LP two-stage
turbocharger built as two shafts in series.
Variable-geometry vane. When the map carries a vane axis
(map.has_vane_axis = True), the trilinear lookup
map.lookup(PR_inv, N_corr, vane=v) returns the vane-corrected
\((\dot m_\text{corr}, \eta_\text{is})\) — the map is the single
source of truth for both swallowing capacity and efficiency vs
vane position, replacing the legacy linear ṁ *= vane scaling.
2-D maps without a vane axis fall back to the linear scaling on
Turbine (Compressor never had it).
Wastegate. Turbines carry an additional bypass_open
fraction \(\beta \in [0, 1]\) that scales extracted shaft power
and blends the post-expansion stagnation temperature
\(T_{0,\text{out}}^\text{full}\) toward the inlet stagnation
temperature \(T_{0,\text{in}}\):
so \(\beta = 0\) is wastegate-closed (full expansion) and
\(\beta = 1\) is full bypass (no shaft work). A
WastegateController (velocity-form PI in aeolus1d.control)
drives \(\beta\) to a target manifold pressure with sign-
inverted action vs the VGT boost controller.
Surge / choke regime detection. Compressor.last_regime
classifies each call as "normal" / "idle" /
"backflow" / "surge" / "choke" based on the requested
\(\mathrm{PR}\) relative to the map's PR-grid edges (and on
\(\omega \le 0\) / \(\mathrm{PR} < 1\)). last_surge_margin
and last_choke_margin are normalised fractional distances
from the top and bottom edges respectively (positive when
inside the envelope, \(\le 0\) at or past the edge). The
wall-close fall-through behaviour is unchanged; the new fields
are pure diagnostics that controllers and post-processors can
read.
6.4 Variable-Valve-Timing (Phase E5)¶
Two actuation hooks layered on top of the continuous-lift schedule
exposed by CylinderPort:
- Cam phaser — a crank-angle offset \(\Delta\theta\) applied to
the schedule lookup. Positive = advance (valve opens earlier in
the cycle):
$$
\theta_\text{eff} = (\theta_\text{crank} + \Delta\theta)
\bmod 720^\circ.
$$
Both the binary open / close window and the piecewise-linear
lift_y(lift_cad)interpolation use \(\theta_\text{eff}\). - Lift scale — a dimensionless multiplier \(s_l \in [0, 1+]\) on the final effective-area fraction: $$ \frac{A_\text{eff}(\theta)}{A_\text{port}} = C_d \cdot \frac{\mathrm{lift}(\theta_\text{eff})}{\mathrm{lift}_\text{max}} \cdot s_l. $$ \(s_l = 1\) is the nominal cam profile; \(s_l < 1\) models a variable-lift mechanism (Honda VTEC low-cam, GM AFM half-lift); \(s_l = 0\) is effective valve closure (cylinder deactivation — Mercedes ZAS, GM AFM, Honda VCM).
Both default to identity (\(\Delta\theta = 0^\circ,\ s_l = 1\)) so
every pre-E5 cylinder case is bit-equal. Both are TOML-configurable
on each [[cylinder.port]] block AND mutable at runtime — a
controller can dial them per step (closed-loop knock-margin
cam-phaser, load-following cylinder deactivation, etc.).
The cylinder's slider-crank angle itself is not part of the control state — it advances kinematically as \(\theta = \theta_0 + \omega_\text{crank} \cdot t\). VVT modifies the mapping from \(\theta\) to port area, not the kinematics themselves.
7. Strang-split chemistry¶
Stiff reactive-source terms are integrated by Strang splitting [26]:
where \(\mathcal{T}\) is hyperbolic transport (the standard SSP-RK
timeloop with include_reaction=False) and \(\mathcal{C}\) is a
per-cell stiff ODE solve at frozen \((\rho, u)\). Strang gives
second-order accuracy on the composite step and is unconditionally
stable for the stiff part regardless of the transport CFL.
Two ODE solvers are wired through the same chemistry sub-step:
lsoda — scipy.solve_ivp(method="LSODA") per cell. Calls
the user's Python py_reaction(_multi) callable. Robust; works
with any reaction shape including the surface-coverage variant
(Phase 49) where the per-cell state vector is
\((Y,\ \mathbf{Y}_\text{extras},\ \boldsymbol{\theta},\ E)\).
bs_jit — Phase-19d Bader-Deuflhard semi-implicit
Bulirsch-Stoer (Numerical Recipes §17.5 [6, 24]). The
semi-implicit midpoint rule is, with frozen Jacobian \(\mathbf{J}\):
The 2×2 analytical Jacobian for the bundled Arrhenius reaction is closed-form (Cramer-rule solve, no LU). Step refinement uses the semi-implicit BS sequence \(\{2, 6, 10, 14, 22, 34, 50\}\) and a Gustafsson step controller [9]:
clamped to \([0.2 H,\ 5 H]\). After Numba warm-up the kernel runs \(\sim\)50–100\(\times\) faster than the LSODA path on stiff Arrhenius problems.
8. Acoustic mode analysis¶
Two complementary tools answer "what natural frequencies does this geometry support?" without running the time-domain solver.
8.1 Helmholtz lumped-parameter¶
The textbook 2-DOF resonator (§5.4) is exposed as a class with the analytical formula \(f_H = (a / 2\pi) \cdot \sqrt{A_n / (V\,L_\text{eff})}\) — useful at design time.
8.2 FEM continuous-Galerkin eigensolver¶
aeolus1d.modes.solve_modes runs a high-order continuous-Galerkin
solve on the variable-area duct eigenproblem [37]:
with Lagrange basis polynomials of arbitrary order \(p\).
Spectral convergence in \(p\), algebraic in n_elements. Returns
natural frequencies + mode shapes; defaults p=4, n_elements=8
are 1e-3 % accurate on the first 5 modes of a uniform pipe.
The Phase-42 multi-pipe variant assembles a connected duct network with junction degrees of freedom, so a chain of non-uniform pipes (e.g. an SCR + intake + plenum) is solved as one block-sparse generalised eigenvalue problem rather than piecewise.
8.3 State-space linearization (Phase L)¶
aeolus1d.analysis.linearize builds a first-order Taylor expansion
around an operating point \((\mathbf{x}^\star, \mathbf{u}^\star)\):
with
The resulting StateSpace plugs straight into
scipy.signal.StateSpace / python-control for Bode plots, pole
analysis, LQR / H∞ controller synthesis, and Kalman observer
design.
Phase L1 — single pipe. State vector is one pipe's interior
conservative variables \((\rho, \rho u, \rho E)\) per cell,
flattened cell-major. \(A\) via central FD on
pipe_residual_first_order. Cost: \(2n\) residual evaluations.
Phase L2 — multi-pipe + junction-coupled. State vector is the
concatenation of every pipe's interior cells; junctions enter
through apply_port_bcs() between FD perturbations so the
junction-driven ghost coupling is captured. Volume / port-u_n
state initially frozen at the operating point.
Phase L2b — live volume dynamics. Each junction's
\((m, E_\text{internal})\) joins the packed state vector; rates
come from Junction.accumulate_rates(). The plenum's mass /
energy dynamics now show up as eigenvalues alongside the acoustic
modes. Opt back to L2 frozen-volume mode with
frozen_volumes=True when the slow plenum dynamics are not the
target of the study, or when the network contains a component
whose state beyond \((m, E)\) (cylinder slider-crank, turbo shaft
\(\omega\)) would need an L2c / L3 treatment.
Phase L2c — port-level kinematic state. Turbo shaft
\(\omega\) (Shaft.ports[0].u_n) and extended-Helmholtz-resonator
slug velocity (EHRJunction.ports[i].u_n) join the packed state.
Detection is duck-typed on the network timeloop's
port_u_rates: list[float] + ports[i].u_n: float contract —
every component matching it works without registration. Cylinder
crank angle is deliberately not packed (kinematic clock, not
an integrable ODE; cycle dynamics belong in a harmonic / cycle-
averaged framework, not in state-space).
Phase L3 — analytic Jacobian. \(A\) is pulled directly from
the implicit Newton solver's closed-form block-tridiagonal
assembly (_assemble_jacobian_first_order) — the same one
Phases B16 → B22 already ship for every closure kind. One
assembly call replaces \(2n\) residual sweeps; machine-precision-
exact rather than FD-truncated; ~5-50× wall-clock speedup on
50+ cell pipes. Pure pipe-and-far-BC networks (no junctions).
Phase L3b — hybrid for junction networks. Drops L3's
"junctions empty" restriction: pipe-on-pipe diagonal blocks are
overlaid analytically onto the FD-derived \(A\); cross-blocks
(pipe ↔ volume, pipe ↔ u_n, pipe-pipe-through-junction) remain
FD. Cost equals FD; accuracy on the dominant block improves to
analytic precision. Pure-analytic cross-blocks (the per-junction
ghost-gradient stitch \(\partial U_\text{ghost} /
\partial(m, E_\text{internal}) \cdot \partial F / \partial U_R\))
are out of scope here — Phase L3c when needed.
Non-normality caveat. The Euler operator is strongly non- normal. Two \(A\) matrices agreeing element-wise to \(10^{-7}\) can produce eigenvalues that disagree by ~10 % — see Trefethen's pseudospectra. The unit tests assert element-wise agreement (strict) and spectral abscissa \(\max \mathrm{Re}(\lambda)\) agreement (the physically meaningful stability metric) rather than pointwise eigenvalue agreement.
9. Inverse problems / optimization¶
Phase 44 ships aeolus1d.optimize for case-driven inverse
problems: fit one or more case parameters to match measured
time-series traces.
A ParamSpec.path is a dotted/bracketed walk into the Case
tree (pipes[0].init.p, pipes["tube"].species.A_ads, …). A
ProbeObjective carries a reference \((t, v)\) array and the
\((\text{pipe}, x_\text{frac}, \text{var})\) triple of a virtual
sensor; the objective runs the case with tracing on, interpolates
each probe series to the reference timestamps, and returns the
weighted-sum metric \(\text{metric}(\text{probe},\ \text{reference})\)
— \(\mathcal{L}_1 / \mathcal{L}_2 / \text{rmse} / \text{max}\).
Optimizers wrap scipy.optimize.minimize with method strings
"nelder-mead", "powell", "l-bfgs-b", "slsqp", …; bounded
methods require bounds on every parameter. The [optimize] TOML
block exposes the same shape end-to-end; aeolus1d optimize
case.toml runs the fit and writes the tuned parameters to a
sibling *.tuned.toml.
The Phase 49 SCR rate-fitting examples
(examples/optimize/scr/) drive this seam end-to-end: fit
adsorption / desorption (and full-SCR oxidation + standard-SCR)
rate parameters from a measured \(Y_\text{NH3}(t)\) (and
optionally \(Y_\text{NO}(t)\)) trace, recovering the truth to
0.001 % on synthesised data.
References¶
[1] L. E. Kinsler, A. R. Frey, A. B. Coppens, J. V. Sanders, Fundamentals of Acoustics, 4th ed., Wiley 2000. Helmholtz resonator + end correction.
[2] M. L. Munjal, Acoustics of Ducts and Mufflers, 2nd ed., Wiley 2014. Extended Helmholtz Resonator (EHR) port-impedance formulation.
[3] J. M. Corberán, M. Fernández-Pacheco, F. Payri, Modeling of unsteady flow in branches of intake and exhaust systems of reciprocating internal combustion engines. Proc. IMechE Part D 217 (2003). CP+h algebraic junction.
[4] S. E. Haaland, Simple and explicit formulas for the friction factor in turbulent pipe flow. ASME J. Fluids Eng. 105 (1983) 89-90. Explicit Colebrook approximation.
[5] Numerical Recipes 3rd ed., §17.5 (Bader-Deuflhard semi-implicit Bulirsch-Stoer).
[6] G. Bader, P. Deuflhard, A semi-implicit mid-point rule for stiff systems of ordinary differential equations. Numer. Math. 41 (1983) 373-398. Semi-implicit BS midpoint step.
[9] K. Gustafsson, Control-theoretic techniques for stepsize selection in implicit Runge-Kutta methods. ACM TOMS 20 (1994) 496-517. Step controller for the BS extrapolation.
[11] J. B. Heywood, Internal Combustion Engine Fundamentals, 2nd ed., McGraw-Hill 2018. Slider-crank kinematics (§2.3), Wiebe heat-release (§9.2.2), Zeldovich NOx (§11.3.1), Hiroyasu soot (§11.3.2), Woschni heat-loss (§12.4).
[12] D. A. Nield, A. Bejan, Convection in Porous Media, 5th ed., Springer 2017. Darcy-Forchheimer two-term drag law.
[13] R. M. Heck, R. J. Farrauto, S. T. Gulati, Catalytic Air Pollution Control: Commercial Technology, 3rd ed., Wiley 2009. Substrate technology + monolith geometry (§4).
[14] M. Ciardelli, I. Nova, E. Tronconi, D. Chatterjee, B. Bandl-Konrad, NH3 SCR of NOx for diesel exhausts aftertreatment: role of NO2 in catalytic mechanism, unsteady kinetics and monolith size effects. Chem. Eng. Sci. 62 (2007) 5001-5006. Standard-SCR rate law parametrisation.
[15] I. Nova, D. Bounechada, R. Maestri, E. Tronconi, A. Heibel, T. Collins, T. Boger, Influence of the substrate properties on the performances of NH3-SCR catalysts. Catal. Today 114 (2006) 3-12. Engineering parameter ranges for Cu-zeolite SCR.
[16] L. Olsson, H. Sjövall, R. J. Blint, A kinetic model for ammonia selective catalytic reduction over Cu-ZSM-5. Appl. Catal. B 87 (2009) 200-210. Langmuir NH3 ads/des on Cu-zeolite.
[17] H. Hiroyasu, T. Kadota, M. Arai, Development and use of a spray combustion modeling to predict diesel engine efficiency and pollutant emissions. JSME Bulletin 26 (1983) 569-583. Soot formation and oxidation rate laws.
[21] T. J. Poinsot, S. K. Lele, Boundary conditions for direct simulations of compressible viscous flows. J. Comp. Phys. 101 (1992) 104-129. LODI characteristic boundary conditions.
[22] O. Redlich, J. N. S. Kwong, On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions. Chem. Rev. 44 (1949) 233-244. Redlich-Kwong EOS.
[23] D.-Y. Peng, D. B. Robinson, A new two-constant equation of state. Ind. Eng. Chem. Fund. 15 (1976) 59-64. Peng-Robinson EOS.
[24] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer 1996, §IV.9. Stiff ODE methods, BS extrapolation, Bader-Deuflhard derivation.
[25] A. Suresh, H. T. Huynh, Accurate monotonicity-preserving schemes with Runge-Kutta time stepping. J. Comp. Phys. 136 (1997) 83-99. MP5 5th-order monotonicity-preserving reconstruction.
[26] G. Strang, On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5 (1968) 506-517. Operator splitting.
[27] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer 2009. HLLC (Ch. 10), MUSCL (Ch. 14), Roe (Ch. 11).
[28] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong stability- preserving high-order time discretization methods. SIAM Review 43 (2001) 89-112. SSP-RK2/RK3 + TVD analysis.
[29] J.-Y. Chalet, J.-P. Chesse, J. Hetet, Detailed analysis of the steady-flow performance of intake and exhaust valves. Proc. IMechE Part D 219 (2005). Effective-area valve BC.
[30] D. E. Winterbone, R. J. Pearson, Theory of Engine Manifold Design, Professional Engineering Publishing 2000. Compressible- orifice valve BC.
[31] J. R. Serrano, F. Arnau, A. Gil, P. Olmeda, F. Payri, Methodology to characterize a turbocharger turbine map under choked flow. SAE Technical Paper 2013-01-2475 (2013). Turbine-map measurement and synthesis.
[32] H. P. Hodson, T. P. Hynes, A simple semi-actuator-disc representation of axial-compressor performance. ASME GT2010-23133 (2010). Quasi-steady centrifugal-compressor PR(N, ṁ) factory.
[33] R. J. Kee, F. M. Rupley, J. A. Miller, CHEMKIN-II, Sandia Report SAND89-8009B (1989). CHEMKIN mechanism format.
[37] B. Cockburn, G. E. Karniadakis, C.-W. Shu (eds.), Discontinuous Galerkin Methods, Springer 2000; J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover 2001. High-order continuous-Galerkin eigensolver background.
[38] G. Woschni, A universally applicable equation for the instantaneous heat transfer coefficient in the internal combustion engine. SAE Technical Paper 670931 (1967). Cylinder wall heat-loss correlation.
[39] G. Hohenberg, Advanced approaches for heat transfer calculations. SAE Technical Paper 790825 (1979). Cylinder wall heat-loss correlation (alternative).
[40] B. Thornber, A. Mosedale, D. Drikakis, D. Youngs, R. J. R. Williams, An improved reconstruction method for compressible flows with low Mach number features. J. Comp. Phys. 227 (2008) 4873-4894. HLLC-LM low-Mach correction.
[41] F. Millo, S. Gundlapally, W. Wang, S. Wahiduzzaman, 1D Simulation Models for Aftertreatment Components, Ch. 11 in A. Onorati, G. Montenegro (eds.), 1D and Multi-D Modeling Techniques for IC Engine Simulation, SAE International (2020), pp. 387-416. ISBN 978-0-7680-9352-0. Canonical reference for the 1+1D flow-through monolith model (gas/solid/washcoat conservation equations, Graetz-Nusselt Nu/Sh asymptotics [Table 11.1], parallel- pore effective diffusivity, the QSS limit-point/concentration-jump robustness problem, and the inner/outer adaptive-mesh decomposition).