Multi-zone diesel combustion and NOx¶

A packet-based, mixing-controlled DI-diesel model. Each injection packet seeds a zone; zones penetrate, entrain air, evaporate, mix and burn, and the heat-release shape is an outcome rather than an input.

The reason to want that is emissions. NO formation goes roughly as the fifth power of temperature, so a bulk-average cylinder temperature cannot predict it — you need the distribution. That is the whole argument for the model, and it is asserted directly in the test suite: a {2600, 1800} K zone pair must out-produce two zones at their mean of 2200 K, or the model has been quietly reduced to a mean.

Relationship to the single-zone predictive model

diesel_predictive_combustion predicts when combustion starts (Hardenberg–Hase) and how much burns premixed (Watson), then hands a prescribed Wiebe the shape. It is cheaper, well-tested and the right choice when you want heat release. This model is the right choice when you want the temperature distribution — that is, when you want NOx.

The zone state machine¶

Per packet, in this order. The order is physical and not interchangeable.

Seed. Liquid fuel mass m_f0 with axial velocity u₀.
Penetrate and entrain. X advances at u_spray; the entrainment factor GF follows the characteristic-length form below.
Evaporate, then mix. Liquid → vapour on τ_vap, vapour → mixed and burnable on τ_mix. Two stages, not one: a zone can be rich in vapour and still unable to burn.
Ignite. Once the spray clears the lift-off length and the ignition delay has elapsed.
Burn. dM_burnt = min(M_mix, M_air/AF_st) — mixing-controlled when fuel is ready, air-limited when it is not.
Update temperature. Mixing first, then the first law.

Two of those orderings matter enough to be pinned by tests:

Fuel cannot burn in the step it is injected. A fresh packet is entirely liquid. A model that lets it burn immediately is not mixing-controlled.
Mixing precedes the first law. Releasing heat into the pre-dilution mass overstates the peak, and T⁵ amplifies that fivefold in the one quantity the model exists to produce.

Correlations¶

Each is a named function with its source and validity range in the docstring. The exponents are verified by fitting them back out of two evaluations, not by inspection — see tests/validation/test_multizone_correlations_vs_literature.py.

Spray and entrainment¶

\[\tan(\alpha/2) = c\left[(\rho_g/\rho_l)^{0.19} - 0.0043\sqrt{\rho_l/\rho_g}\right]\]

\[x^{+} = \sqrt{\rho_l/\rho_g}\; d_{hole} / (a \tan(\alpha/2))\]

\[GF = \tfrac{1}{2}\left(\sqrt{1 + 16(x/x^{+})^2} - 1\right)\]

The two limits of GF are what make the form right rather than merely convenient: GF → 4(x/x⁺)² near the nozzle — quadratic, so a packet that has barely left entrains almost nothing — and GF → 2x/x⁺ in the far field, the classical jet result. The factor of 4 is easy to write as 2; getting it wrong misattributes an early-injection entrainment error to the spreading angle.

The geometry factor a = 0.66 is normally held fixed; calibrate through the spreading coefficient instead. Tuning both is double-counting, since a wider cone and a shorter x⁺ are the same knob.

Lift-off length¶

\[H = c\,T_a^{-3.74}\rho_a^{-0.85}d^{0.34}\,U\,(A/F)_{st}\]

SAE 2005-01-3843, with c = 7.04e8 in the correlation's own mixed units: H in mm, T in K, ρ in kg/m³, d in µm, U in m/s. Typical H is 6–8 mm — but only near 1000 K. It is steeply temperature-sensitive:

The shipped default is no longer the published constant

Since 2026-08-25 mz_c_liftoff defaults to 6.225e8, a calibrated value, which is 0.884 × the published 7.04e8 above. The table below is computed with the published constant, so read it as the shape of the correlation rather than as what the solver will produce; multiply by 0.884 for the shipped default (7.0 mm at 1000 K becomes 6.2 mm). The exponents are unchanged. See docs/changelog.md.

charge temperature	H
900 K	10.3 mm
1000 K	7.0 mm
1100 K	4.9 mm

Quoting an expected H without its temperature is how this correlation gets mis-calibrated.

Ignition delay¶

\[\tau_{id} = c\,p^{-1} X_{O_2}^{-1.2}\exp(E_a/RT), \qquad E_a/R = 3242\ \mathrm{K}\]

The oxygen exponent is the mechanism by which EGR retards combustion. The prefactor is anchored to give 1.0 ms at 60 bar / 900 K / air, which yields 0.52 ms at 80 bar / 1000 K and 1.53 ms at 30 % EGR.

The prefactor is the part that goes wrong silently

Every exponent above can be correct while the absolute scale is meaningless. An earlier draft of this model had a prefactor five orders of magnitude out — a 40 nanosecond ignition delay — and passed every sensitivity test, because those check ratios. Only an order-of-magnitude assertion caught it. The same happened to the lift-off constant.

Zone temperature¶

Two closures are implemented. The theory this follows specifies both and names neither as preferred; the reference implementation defaults to the 2013 one. This port defaults to the 2018 one.

Mode 2 — mixing then the first law (default). Entrained unburnt gas dilutes and cools the zone, then energy is added:

\[T \leftarrow \frac{M_{old}T_{old} + dM\,T_{ub}}{M_{new}}, \qquad dT = \frac{dQ - p\,dV}{c_v M_{new}}\]

Mode 1 — adiabatic-flame jump (parity only). An empirical ΔT_flame polynomial. Kept because parity against the reference needs it, and not default because its three coefficients have no cited source.

Four further modes (3, 4, 11, 22) existed upstream and are not ported — undocumented, no equations, undated labels. Each is recorded with its reason in UNPORTED_TZONE_MODES, so nobody re-derives them from the source.

NOx¶

Per zone, the reduced extended Zeldovich rate (Heywood eq. 11.6):

\[\frac{d[NO]}{dt} = \frac{2R_1(1-\beta^2)}{1 + \beta R_1/(R_2+R_3)}, \qquad \beta = [NO]/[NO]_{eq}\]

The equilibrium radicals [O]ₑ, [N]ₑ, [OH]ₑ come from the twelve-species equilibrium solve, not from partial-equilibrium correlations. That is what the {H, O, N, NO₂} extension to chemistry.equilibrium was for.

The (1 − β²) factor makes this relax rather than only grow: a zone that cooled after forming NO at a higher temperature has a negative rate and decays back toward equilibrium.

Rate constants are cm³/(mol·s)

Heywood's constants are quoted in cm³/(mol·s). Feeding them SI concentrations overstates every rate by 10⁶ and pins NO at whatever guard ceiling is in place — which reads as clipped, not as wrong units. A ceiling-valued result is indistinguishable from a saturated-but-correct one by inspection, so the test asserts NO is well below the ceiling.

NO carries almost no information about K_N

Correcting K_N from a wrong pair to a NASA-7-derived one changed [N] by 367× at 1800 K and the NO rate by under 0.02 %. [N] enters only through β·R₁/(R₂+R₃), and R₂, R₃ are both linear in it while R₁ is not — so the term goes as 1/[N] and is already negligible. K_N must come from thermodynamics; NO data cannot constrain it.

Validation against measurement¶

Tier 1 checks each correlation against its published source. Tier 3 checks the model against a test bed: 247 operating points of measured engine-out NOx, speed, torque, lambda and crank-resolved cylinder pressure.

With generic constants and no fitting, this port's NO driver rank-correlates with measured engine-out NOx at r = +0.50 over 229 usable points. That is a claim about structure — the operating points that made the engine produce more NO make this model produce more too — not about magnitude, which would require a calibration this package deliberately does not ship.

Validate the oracle before comparing anything to it

Tier 2, parity against the reference implementation's own results, is skipped — and that is a finding rather than an omission. Its noxc channel rank-correlates with lambda at exactly +1.0000 over 223 cases, so it is a deterministic function of lambda rather than an independent NOx result. A model "agreeing" with it would only have been agreeing with lambda, and the agreement would have been reported as validation.

The give-away was checking against peak combustion temperature: −0.98. A Zeldovich-driven quantity cannot fall that cleanly with temperature.

The same check run against the test-bed data passes: −0.93 with lambda, +0.84 with torque, +0.74 with peak cylinder pressure — strong, physical, and nowhere near ±1.000.

The measurement data is proprietary and stays out of the repository; only the aggregate statistics above are recorded. See examples/validation/DATA_MANIFEST.md.

Traceability¶

Quantity	Source	Implements	Validity
Spreading angle	Siebers, SAE 980809 / 1999-01-0528	`spreading_angle_tangent`	collapses at low `ρ_g`
Characteristic length `x⁺`	Siebers, SAE 1999-01-0528	`characteristic_length`	mm scale in-cylinder
Entrainment `GF`	Siebers mixing-limited jet	`entrainment_factor`	4(x/x⁺)² → 2x/x⁺
Lift-off length	Siebers & Higgins, SAE 2005-01-3843	`liftoff_length`	6–8 mm near 1000 K
Ignition delay	SAE 2005-01-3837 basis	`ignition_delay_seconds`	0.5–2 ms at TDC
Zone T (2018)	theory report, `TzoneMode = 2`	`zone_temperature_energy_balance`	no fitted coefficients
Zone T (2013)	theory report, `TzoneMode = 1`	`zone_temperature_adiabatic_flame`	uncited coefficients
Extended Zeldovich	Heywood eq. 11.6	`multizone_nox`	needs `[O]ₑ, [N]ₑ, [OH]ₑ`
Series restrictions	injector report	`equivalent_area`	`(ΣA⁻²)^(−½)`
Common rail	injector report	`rail_pressure_derivative`	mean droop only

What is deliberately absent¶

Stated so the boundary is visible rather than inferred:

Swirl. Carried as a parameter, defaults to zero, as in the reference.
Radiative loss from the zone. The theory carries it as optional and gives no coefficient.
Wall impingement beyond the k_mix mixing slowdown.
Water-hammer and pulse-to-pulse rail dynamics. The injector model is lumped; the source is explicit that these need a distributed model and can materially change a main injection after a pilot. A fitted discharge coefficient is absorbing some of what they would have done.

Configuration¶

[cylinder.combustion]
kind = "multizone_diesel"
mz_c_xchar = 0.66          # held fixed; calibrate mz_c_tan_alpha instead
mz_c_tau_vap = 0.15e-3
mz_c_tau_mix = 0.8e-3
mz_tzone_mode = 2          # the only mode; 1 is refused (it had no code path)
inj_a_holes_m2 = 7.0e-8
inj_rail_pressure_pa = 1.8e8

A worked case is examples/89_dd13_multizone.toml.

Defaults are generic, never engine-fitted

Every shipped default is a literature value or a stated generic. Values calibrated to a particular engine are calibration, not physics, and are kept out of the package — test_proprietary_data_stays_local enforces that mechanically, because a fitted float is indistinguishable from a chosen one by inspection.