Catalysts, filters and surface chemistry¶
The aftertreatment closures: what a pipe's [pipe.catalyst], [pipe.filter],
[pipe.species] reaction, [pipe.thermal_wall] and [pipe.porous] blocks
turn into, what each closure adds to the residual (heat, species, momentum)
and to the slow solid states (catalyst solid, soot, cake, wall), which driver
reads which closure, and where each refusal happens. The contracts of every
function below are in their docstrings (the standard is
Function contracts and flow charts).
Units to keep in mind, because they are the classic traps here: a rate in
mol/(m^3 s) is per volume of the brick, dH is per mole of the indicator
reactant for the TWC, the mechanism and the SCR closures but per kg for
the soot regeneration (regen_dH, J per kg soot) and the CSF (csf_dH_co,
J per kg fuel); Y is a mass fraction and every dY/dt is per second on the
gas mass, so the mole rate is multiplied by the molar mass and divided by
rho.
From the case to the closures¶
What enters: the parsed PipeSpec. What leaves: a Pipe with
catalyst, filter, py_reaction_multi / py_reaction_surface /
py_reaction, py_heat and porous attached, plus the coolant and wall
lists and the DOC-sulfur controller.
flowchart TD
SPEC["PipeSpec: catalyst, filter, species.reaction, thermal_wall, porous"] --> AC{"_attach_closures"}
AC --> AP{"_attach_porous"}
AP --> PS["porous_source: Darcy-Forchheimer momentum sink, read by the residual kernels"]
AC --> CAT{"_attach_catalyst"}
CAT -- "cpsi > 0: monolith form" --> MR{"resolve_monolith_params / resolve_monolith_var_nu_params"}
MR -- "cpsi <= 0 or not finite, ofa outside (0, 1), unknown substrate" --> X1["ValueError"]
MR --> CS["CatalystState: kind, params, T_solid"]
CAT -- "explicit rho_s_cs, a_v, h" --> CS
AC --> FI{"_attach_filter -> _dpf_state"}
FI -- "frac_fast outside [0, 1]" --> X2["ValueError"]
FI --> FS["FilterState: trap, cake, regen, thermal, CSF"]
AC --> SP{"_attach_species"}
SP -- "PM species not one of the extras" --> X3["ValueError, in _cache_filter_species_columns"]
SP --> RX{"reaction kind"}
RX -- "hc_oxidation / doc_oxidation" --> R1["make_hc_oxidation, make_doc_oxidation, make_doc_oxidation_voltz"]
RX -- "doc_oxidation_S" --> R2["DOCSulfurOxidation + DOCSulfurController"]
RX -- "hnco_hydrolysis" --> R3["make_hnco_hydrolysis"]
RX -- "nh3_ads_des / nh3_scr_kinetics / nh3_slip_catalyst" --> R4["make_nh3_ads_des, make_nh3_scr_kinetics, make_nh3_slip_catalyst"]
RX -- "twc" --> R5{"_rx_twc"}
R5 -- "HC species is not C3H6" --> X4["ValueError"]
R5 --> R5b["make_twc_pipe_demo"]
RX -- "twc_brinkmeier" --> R6["make_brinkmeier_twc_pipe_demo"]
R6 -- "an array of the wrong shape" --> X5["ValueError"]
RX -- "mechanism" --> R7{"make_mechanism"}
R7 -- "duplicate species, no molar mass or one <= 0, bad site capacity, no reactions, undeclared species" --> X6["ValueError"]
R7 --> R7b["gas: py_reaction_multi; surface: py_reaction_surface + theta"]
AC --> TW{"_build_coolant_system"}
TW -- "wall without heat capacity or area" --> X7["ValueError, in _check_standalone_wall_spec / _check_coupled_wall_spec"]
TW --> WL["ThermalWall / MultiLayerWall, gas-side py_heat overlay"]The seams are PyReactionFn (T, Y, rho) -> (dY/dt, s_E) for a single
species (make_hc_oxidation), PyReactionMultiFn (T, Y, Y_extras, rho) ->
(dY/dt, dY_extras/dt, s_E) for the gas closures, and PyReactionSurfaceFn
(T, Y, Y_extras, rho, theta) -> (dY/dt, dY_extras/dt, s_E, dtheta/dt) for
the SCR storage, the Brinkmeier TWC and surface mechanisms. The species a
closure reads sit at fixed columns of Y_extras (SpeciesContract): the
build resolves them once and hands the closure plain indices, and a column
below zero means "species absent, read as zero".
What each closure adds to one step¶
Where the sources land: every closure adds to R[:, 2] (gas energy), a
cake or porous drag adds to R[:, 1] (momentum), species closures return
dY/dt for the species transport, and the slow solids return their own
rates for the driver to integrate.
flowchart TD
ST["pipe state: U, Y, Y_extras; slow states T_solid, soot, T_cake, T_wall"] --> CATR{"catalyst_rates_into"}
CATR -- "no catalyst" --> NC["None"]
CATR -- "ideal, per-cell gamma, real gas" --> CK["catalyst_apply_cells / _var_gamma / catalyst_apply_cells_eos"]
CK --> CO["R[:, 2] += h a_v (T_s - T_g); dT_solid/dt returned"]
ST --> CKI{"catalyst_K_into: quasi-steady only"}
CKI -- "real gas, or no catalyst" --> NC2["None: the explicit overlay is used instead"]
CKI --> CKO["K = h a_v per cell, T_solid (live array)"]
ST --> TRAP{"_overlay_filter_trap"}
TRAP -- "no PM column" --> NT["None"]
TRAP --> TR["filter_trap_rates: dY_PM/dt = -lambda abs(u) Y_PM; soot += the same mass; rho debited"]
TRAP -- "epsilon < 0" --> X1["ValueError"]
ST --> CAKE{"cake_into"}
CAKE -- "no filter, no soot array, resistance <= 0" --> NCK["untouched"]
CAKE --> CAKEO["R[:, 1] += -mu_ref u (soot + soot_fast) R_cake; energy untouched"]
ST --> REGEN{"regen_into"}
REGEN -- "no filter, no soot, direct path bare" --> NR["None"]
REGEN -- "T_cake size 0" --> RD["dpf_regen_apply_cells: R[:, 2] += dH rate"]
REGEN -- "T_cake size N" --> RT["dpf_regen_thermal_apply_cells: gas gets h_eff (T_cake - T_gas); dT_cake/dt returned"]
RD --> RO["d_soot_dt returned; fast rate left in FilterState._d_soot_dt_fast"]
RT --> RO
ST --> WALLF{"substrate_wall_into"}
WALLF -- "no wall array or h_wall <= 0" --> NW["None"]
WALLF --> WO["R[:, 2] += h_wall (T_wall - T_gas); dT_wall/dt returned"]
ST --> PYR{"py_reaction_multi / py_reaction / surface"}
PYR --> PYO["dY/dt, dY_extras/dt, s_E into the species transport and R[:, 2]"]
RO --> POST["after the step: regen_species_fold (carbon burned -> CO2, O2 / NO2 debited, rho up)"]
ST --> CSF{"csf_into"}
CSF -- "no filter, no CSF constants, both fuels untracked" --> NCS["None"]
CSF --> CSFO["heat into R[:, 2]; species rates returned"]
CSFO --> CSFS["apply_csf_transient: step scaled so no species is over-drawn"]The gas heat closures of a wall (make_wall_coupled_py_heat,
make_multilayer_wall_py_heat) are the py_heat overlay: they read the wall's
live temperature and return h_gas (T_wall - T_gas) A_wall / V_cell. The
walls themselves advance at the end of a macro step:
flowchart TD
A["end of a macro step: coolant triples, wall pairs, dt"] --> B{"advance_thermal_couplings"}
B --> C{"advance_coolant_system"}
C -- "N mismatch between gas pipe, wall and coolant pipe; wall mass or cp not > 0" --> X1["ValueError from wall_coolant_rates"]
C --> C1["sub-cycled at 1 / k_max: T_wall and T_cool advanced in place"]
B --> D{"advance_standalone_walls"}
D -- "ThermalWall" --> D1["forward Euler: gas heat, convection, radiation, axial conduction"]
D -- "MultiLayerWall, explicit" --> D2["_advance_multilayer_wall: T updated in place"]
D -- "MultiLayerWall, implicit" --> D3["_advance_multilayer_wall_implicit: Crank-Nicolson, T rebound to a new array"]Which driver reads which closure¶
| Closure | Explicit RK (_ExplicitNetwork) |
Implicit Newton | Quasi-steady |
|---|---|---|---|
| porous drag | porous_source inside every residual kernel |
porous_source, and the analytic porous_jacobian / porous_jacobian_sutherland_ideal |
the steady flow solve carries the same drag |
| catalyst heat | catalyst_rates_into in stage_rates |
_augment_residual_with_catalyst and _augment_jacobian_with_catalyst (ideal gas: catalyst_source_gas_energy, catalyst_jacobian); the solid by the Lie split |
catalyst_K_into folded semi-implicitly, else catalyst_rates_into |
| DPF trap | _overlay_filter_trap inside _species_rate |
the same _species_rate |
the quasi-steady species source |
| DPF cake drag | cake_into |
_augment_residual_with_dpf_cake and dpf_cake_jacobian; else cake_into |
in the pressure balance only when cake_qs_drag is on |
| DPF regeneration | regen_into, then regen_species_fold in post_rk |
regen_into in the Lie split; opt-in regen_gas_source_into folds the direct-path heat into the residual |
regen_into in _qs_advance_pipe, then _qs_advance_soot and _qs_advance_cake |
| substrate wall | substrate_wall_into |
_augment_residual_with_substrate_wall and substrate_wall_jacobian, the wall by the Lie split |
substrate_wall_into for the gas side, a closed-form step for the wall |
| CSF oxidation | apply_csf_transient in post_rk |
apply_csf_transient |
csf_into, then _qs_advance_csf_species |
gas reaction (py_reaction_multi) |
_overlay_species_reaction |
the explicit overlay, or the chemistry split | _qs_species_source |
surface reaction (py_reaction_surface) |
not in the residual | the chemistry split | _qs_surface_chem_step |
| walls and coolant | advance_thermal_couplings at the end of the step |
same | same |
The Strang driver runs the same transport stage as the explicit one with the
reaction overlay switched off, and advances every reaction closure in the
chemistry sub-step below (chemistry_step_on_pipes).
The chemistry sub-step¶
A pipe with a reaction closure has its chemistry advanced apart from the
transport, over the accumulated dt with the density and the kinetic energy
held fixed.
flowchart TD
IN["pipe with a reaction closure, dt, rtol, atol, chem_solver"] --> SEL{"chemistry_step_on_pipe: chem_solver"}
SEL -- "unknown name" --> X0["ValueError"]
SEL -- "bs_jit / auto with a plain Arrhenius pipe" --> BJ["_chem_step_bs_jit"]
SEL -- "bs_brinkmeier" --> BB{"_chem_step_bs_brinkmeier"}
BB -- "no _bs_jit_step on the seam" --> LS["falls through to LSODA"]
BB --> STEP{"the closure's _bs_jit_step"}
STEP -- "C3H6 oxidation, NH3 formation or the OSC heat split on" --> X1["NotImplementedError"]
STEP -- "no theta_Ce coverage" --> X2["RuntimeError"]
STEP --> KER["Bader-Deuflhard kernel: Y, extras, theta_Ce and E advanced in place"]
KER -- "a cell did not converge" --> X3["ChemistryConvergenceError; the caller restores Y, extras, E, theta"]
SEL -- "numbalsoda" --> NB{"_chem_step_numbalsoda"}
NB -- "not a compiled mechanism" --> LS
SEL -- "lsoda" --> LS
LS --> SURF{"py_reaction_surface?"}
SURF -- yes --> SL["chemistry_step_surface_lsoda: coverages clipped to [0, 1], over-draw refused"]
SURF -- no --> GL["_gas_lsoda_step_on_pipe: per cell, committed together after the loop"]Under the quasi-steady driver the Brinkmeier surface step runs through
_qs_surface_chem_step, which can call the fused @njit right-hand side
_qs_brinkmeier_rhs_fused instead of the Python seam wrapper; the two give
the same rates (test_the_fused_rhs_equals_the_seam_wrapper_in_canonical_order).
A coverage that the step over-draws outside [0, 1] is refused there, not
clipped.
Functions¶
src/aeolus1d/closures/catalyst.py::catalyst_ratessrc/aeolus1d/closures/catalyst.py::_catalyst_rates_at_Tsrc/aeolus1d/closures/catalyst.py::catalyst_jacobiansrc/aeolus1d/closures/catalyst.py::catalyst_source_gas_energysrc/aeolus1d/closures/catalyst.py::catalyst_apply_cellssrc/aeolus1d/closures/catalyst.py::catalyst_apply_cells_var_gammasrc/aeolus1d/closures/catalyst.py::catalyst_K_apply_cellssrc/aeolus1d/closures/catalyst.py::catalyst_rates_intosrc/aeolus1d/closures/catalyst.py::catalyst_K_intosrc/aeolus1d/closures/monolith.py::resolve_monolith_paramssrc/aeolus1d/closures/monolith.py::resolve_monolith_var_nu_paramssrc/aeolus1d/closures/porous.py::porous_sourcesrc/aeolus1d/closures/porous.py::porous_jacobiansrc/aeolus1d/closures/porous.py::porous_jacobian_sutherland_idealsrc/aeolus1d/closures/filter.py::filter_trap_ratessrc/aeolus1d/closures/filter.py::cake_intosrc/aeolus1d/closures/filter.py::cake_soot_totalsrc/aeolus1d/closures/filter.py::dpf_cake_apply_cellssrc/aeolus1d/closures/filter.py::dpf_cake_jacobiansrc/aeolus1d/closures/filter.py::regen_intosrc/aeolus1d/closures/filter.py::regen_gas_source_intosrc/aeolus1d/closures/filter.py::dpf_regen_apply_cellssrc/aeolus1d/closures/filter.py::dpf_regen_thermal_apply_cellssrc/aeolus1d/closures/filter.py::regen_species_foldsrc/aeolus1d/closures/filter.py::regen_species_fold_cellssrc/aeolus1d/closures/filter.py::substrate_wall_intosrc/aeolus1d/closures/filter.py::substrate_wall_jacobiansrc/aeolus1d/closures/filter.py::csf_intosrc/aeolus1d/closures/filter.py::apply_csf_transientsrc/aeolus1d/closures/hc_oxidation.py::make_hc_oxidationsrc/aeolus1d/closures/hc_oxidation.py::make_doc_oxidationsrc/aeolus1d/closures/hc_oxidation.py::make_doc_oxidation_voltzsrc/aeolus1d/closures/doc_sulfur.py::DOCSulfurOxidation.__call__src/aeolus1d/closures/doc_sulfur.py::DOCSulfurController.__call__src/aeolus1d/closures/hnco_hydrolysis.py::make_hnco_hydrolysissrc/aeolus1d/closures/urea_decomposition.py::make_urea_one_stepsrc/aeolus1d/closures/urea_decomposition.py::make_urea_thermolysissrc/aeolus1d/closures/scr_kinetics.py::make_nh3_ads_dessrc/aeolus1d/closures/scr_kinetics.py::make_nh3_scr_kineticssrc/aeolus1d/closures/scr_kinetics.py::make_nh3_slip_catalystsrc/aeolus1d/closures/scr_storage.py::make_storage_scr_demosrc/aeolus1d/closures/twc.py::make_twc_pipe_demosrc/aeolus1d/closures/twc_brinkmeier.py::make_brinkmeier_twc_pipe_demosrc/aeolus1d/closures/twc_brinkmeier.py::_qs_brinkmeier_rhs_fusedsrc/aeolus1d/closures/mechanism.py::make_mechanismsrc/aeolus1d/closures/wall_thermal.py::make_wall_coupled_py_heatsrc/aeolus1d/closures/wall_thermal.py::make_multilayer_wall_py_heatsrc/aeolus1d/closures/wall_thermal.py::advance_standalone_wallssrc/aeolus1d/closures/wall_thermal.py::_advance_multilayer_wallsrc/aeolus1d/closures/wall_thermal.py::_advance_multilayer_wall_implicitsrc/aeolus1d/closures/wall_coolant.py::advance_coolant_systemsrc/aeolus1d/io/build.py::_attach_closuressrc/aeolus1d/io/build.py::_attach_poroussrc/aeolus1d/io/build.py::_attach_catalystsrc/aeolus1d/io/build.py::_attach_filtersrc/aeolus1d/io/build.py::_attach_speciessrc/aeolus1d/io/build.py::_dpf_statesrc/aeolus1d/io/build.py::_cache_filter_species_columnssrc/aeolus1d/io/build.py::_rx_twcsrc/aeolus1d/io/build.py::_build_coolant_systemsrc/aeolus1d/io/build.py::_check_standalone_wall_specsrc/aeolus1d/io/build.py::_check_coupled_wall_specsrc/aeolus1d/solver/network_timeloop.py::_ExplicitNetwork.stage_ratessrc/aeolus1d/solver/network_timeloop.py::_ExplicitNetwork.post_rksrc/aeolus1d/solver/network_timeloop.py::_overlay_filter_trapsrc/aeolus1d/solver/network_timeloop.py::_overlay_species_reactionsrc/aeolus1d/solver/implicit_kernels.py::_augment_residual_with_catalystsrc/aeolus1d/solver/implicit_kernels.py::_augment_jacobian_with_catalystsrc/aeolus1d/solver/implicit_kernels.py::_augment_residual_with_dpf_cakesrc/aeolus1d/solver/implicit_kernels.py::_augment_residual_with_dpf_regensrc/aeolus1d/solver/implicit_kernels.py::_augment_residual_with_substrate_wallsrc/aeolus1d/solver/implicit_kernels.py::_augment_jacobian_with_dpf_cakesrc/aeolus1d/solver/implicit_kernels.py::_augment_jacobian_with_substrate_wallsrc/aeolus1d/solver/quasi_steady_network.py::_qs_advance_pipesrc/aeolus1d/solver/quasi_steady_network.py::_qs_catalyst_couplingsrc/aeolus1d/solver/quasi_steady_network.py::_qs_advance_sootsrc/aeolus1d/solver/quasi_steady_network.py::_qs_advance_cakesrc/aeolus1d/solver/quasi_steady_network.py::_qs_advance_csf_speciessrc/aeolus1d/solver/quasi_steady_network.py::_qs_species_sourcesrc/aeolus1d/solver/quasi_steady_network.py::_qs_surface_chem_stepsrc/aeolus1d/solver/solver_ops.py::advance_thermal_couplingssrc/aeolus1d/solver/chemistry_ops.py::chemistry_step_on_pipessrc/aeolus1d/solver/chemistry_ops.py::chemistry_step_on_pipesrc/aeolus1d/solver/chemistry_ops.py::_chem_step_bs_brinkmeiersrc/aeolus1d/solver/chemistry_ops.py::_chem_step_bs_jitsrc/aeolus1d/solver/chemistry_ops.py::_chem_step_numbalsodasrc/aeolus1d/solver/chemistry_ops.py::surface_mechanism_step_on_pipesrc/aeolus1d/solver/chemistry_ops.py::chemistry_step_surface_lsodasrc/aeolus1d/solver/chemistry_ops.py::_gas_lsoda_step_on_pipe