Two-temperature ZND detonation on Cantera 3: evaluating rates at Park temperature and keeping the density solve stable

Goal: a steady ZND (detonation structure) model where O2, N2, H2 and CO keep their own vibrational temperature behind the shock (Landau-Teller relaxation), with Park-style vibration-chemistry coupling (k evaluated at T^s * Tv^(1-s) for reactions that break a vibrating molecule).

Tricks that worked: - Rates at a different temperature without changing concentrations: keep a second ct.Solution and set g2.TDY = Ta, rho, Y. Cantera concentrations depend only on density and Y, so g2.net_rates_of_progress gives k(Ta) times the true concentrations. Splice those rates into the selected reactions, then form species rates with nu @ q, where nu = product_stoich_coeffs - reactant_stoich_coeffs (dense in 3.2). Use one extra Solution per vibrating species group. Skip reaction types containing three-body, falloff or chemically-activated, so recombinations are not coupled. - Nonequilibrium enthalpy: h = sum Y_k h_k(T) - sum_s Y_s/W_s [e_v(T) - e_v(Tv_s)], with harmonic-oscillator e_v. Get h_k(T) quickly with gh.TP = T, p; gh.standard_enthalpies_RT*R*T/W, which needs no composition and only one Cantera call per T. Tv inverts in closed form: Tv = theta/ln(1 + R*theta/e_v). - Steady 1-D conservation (rho u = m, p + rho u^2 = P0, h + u^2/2 = H0) reduces to one equation in density at each state. Use brentq in a narrow bracket around the previous density (continuation on the subsonic branch). For the frozen shock jump, scan density from high to low for the first sign change, so the trivial root at rho1 is skipped. - Stiff BDF trial steps can release more heat than the flow can take (past the sonic point), so the density root vanishes and the run dies with "root lost". Fix: inside the RHS, return the density that minimises |residual| (the choked state) and let error control reject the step. In post-processing, re-solve strictly and stop at the first failure. Also clamp vibrational energies to be non-negative before inverting for Tv, or NaNs appear. Stopping the integration at about 60% of heat release (an event on T) avoids the sonic singularity, and the steepest-temperature-rise point used as the induction length comes before it. - Millikan-White per collision pair, as 1/tau_s = sum_r X_r/tau_sr plus Park's high-temperature limit, is a few lines of numpy. It is poorly calibrated for N2 and O2 against hydrocarbons and H2O, so scale tau by 0.3 and 3 to show sensitivity.

Cost: about 5 to 30 s per mixture with GRI-Mech 3.0, and an 8-process Pool ran 30 mixtures times 4 models in under a minute.