H2/O2 combustion 9-species (Mueller subset)

def _h2_mean_mw(Y):
    """Mean molecular weight from mass fractions (Y for 8 reactive species)."""
    Y_N2 = max(1.0 - np.sum(Y), 0.0)
    inv_W = np.sum(Y / _H2_MW) + Y_N2 / _MW["N2"]
    return 1.0 / max(inv_W, 1e-30)

def h2_combustion_9sp_rhs(t, y):
    T = max(y[0], 200.0)
    Y = np.maximum(y[1:], 0.0)

    W_mix = _h2_mean_mw(Y)
    rho = _P_ATM * W_mix / (_R_U * T) * 1e-3   # kg/m^3 (MW in g/mol -> /1000)

    # Molar concentrations [mol/cm^3] (cgs for rate constants)
    rho_cgs = rho * 1e-3   # g/cm^3
    C = np.empty(8)
    for i in range(8):
        C[i] = max(rho_cgs * Y[i] / _H2_MW[i], 0.0)
    C_N2 = max(rho_cgs * max(1.0 - np.sum(Y), 0.0) / _MW["N2"], 0.0)
    C_M = np.sum(C) + C_N2   # third-body concentration

    inv_T = 1.0 / T

    # Arrhenius rates [mol/cm^3/s]
    k1 = 3.52e16 * T**(-0.7) * np.exp(-17070.0 * inv_T)
    k2 = 5.06e4 * T**2.67 * np.exp(-6290.0 * inv_T)
    k3 = 2.65e19 * T**(-1.3)                           # third-body
    k4 = 1.17e9 * T**1.3 * np.exp(-3626.0 * inv_T)
    k5 = 4.20e14 * np.exp(-12000.0 * inv_T)
    k6 = 2.95e14 * np.exp(-24370.0 * inv_T)            # unimolecular + M

    # Reaction rates  (indices: H2=0 O2=1 H2O=2 H=3 O=4 OH=5 HO2=6 H2O2=7)
    r1 = k1 * C[3] * C[1]          # H + O2 -> OH + O
    r2 = k2 * C[4] * C[0]          # O + H2 -> OH + H
    r3 = k3 * C[3] * C[1] * C_M    # H + O2 + M -> HO2 + M
    r4 = k4 * C[5] * C[0]          # OH + H2 -> H2O + H
    r5 = k5 * C[6] * C[6]          # HO2 + HO2 -> H2O2 + O2
    r6 = k6 * C[7] * C_M           # H2O2 + M -> 2 OH + M

    # Species production rates omega_dot [mol/cm^3/s]
    omega = np.zeros(8)
    omega[0] += -r2 - r4                  # H2
    omega[1] += -r1 - r3 + r5            # O2
    omega[2] += r4                        # H2O
    omega[3] += -r1 - r3 + r2 + r4       # H
    omega[4] += r1 - r2                   # O
    omega[5] += r1 + r2 - r4 + 2.0 * r6  # OH
    omega[6] += r3 - 2.0 * r5            # HO2
    omega[7] += r5 - r6                   # H2O2

    # Convert to SI: omega [mol/cm^3/s] -> [mol/m^3/s] (* 1e6)
    omega_si = omega * 1e6

    dy = np.empty(9)
    # dY_k/dt = W_k * omega_dot_k / rho  (W_k in kg/mol = g/mol * 1e-3)
    for i in range(8):
        dy[i + 1] = (_H2_MW[i] * 1e-3) * omega_si[i] / max(rho, 1e-30)

    # dT/dt = -1/(rho*cp) * sum(hf_k * omega_dot_k)  [hf in J/mol, omega in mol/m^3/s]
    q_dot = np.dot(_H2_HF, omega_si)
    dy[0] = -q_dot / max(rho * _CP, 1e-30)

    return dy