Ablation Surface Multispecies

def _arrhenius_rate(T, alpha):
    """Arrhenius decomposition rate with numerical safeguards."""
    T_safe = np.clip(T, 200.0, 5000.0)
    remaining = np.clip(1.0 - alpha, 0.0, 1.0)
    exp_term = np.exp(-_EA_DECOMP / (_R_GAS * T_safe))
    return _A_DECOMP * exp_term * remaining ** _N_DECOMP

def _multispecies_rhs(t, y):
    dy = np.zeros(22)

    Y = np.clip(y[0:6], 0.0, 1.0)
    T = np.clip(y[6:14], 200.0, 5000.0)
    recession = y[14]
    rec_rate = y[15]
    char_thick = max(y[16], 0.0)
    alpha_sub = np.clip(y[17:22], 0.0, 1.0)

    T_surf = T[0]

    # --- Heterogeneous surface reactions ---
    r1 = _K1 * Y[_I_O2] * np.exp(-_E1 / (_R_GAS * T_surf))
    r2 = _K2 * Y[_I_CO2] * np.exp(-_E2 / (_R_GAS * T_surf))
    r3 = _K3 * Y[_I_H2O] * np.exp(-_E3 / (_R_GAS * T_surf))

    total_ablation_rate = r1 + r2 + r3
    dm_dt = total_ablation_rate  # mass loss rate per unit area

    # Molar masses: CO=28, CO₂=44, H₂=2, H₂O=18, N₂=28, O₂=32
    # Species production/consumption per unit area:
    #   R1: -O₂ (32g), +CO₂ (44g)
    #   R2: -CO₂ (44g), +2CO (56g)
    #   R3: -H₂O (18g), +CO (28g), +H₂ (2g)

    prod = np.zeros(6)
    prod[_I_CO2] = r1 * (44.0 / 32.0) - r2
    prod[_I_CO] = 2.0 * r2 * (28.0 / 44.0) + r3 * (28.0 / 18.0)
    prod[_I_H2] = r3 * (2.0 / 18.0)
    prod[_I_H2O] = -r3
    prod[_I_O2] = -r1
    prod[_I_N2] = 0.0

    # Species mass fraction evolution
    for i in range(6):
        dy[i] = (prod[i] - Y[i] * dm_dt) / _M_SURFACE

    # Freestream entrainment drives species back toward freestream composition
    _tau_mix = 0.5  # s, mixing timescale
    Y_free = np.array([0.0, 0.0, 0.0, 0.0, 0.77, 0.23])
    dy[0:6] += (Y_free - Y) / _tau_mix

    # --- Thermal (8-node 1D conduction) ---
    k_nodes = np.full(_N_THERMAL, 1.5)  # W/(m·K) char conductivity
    rho_cp = _RHO_CHAR2 * _CP
    dx2_inv = 1.0 / (_DX2 * _DX2)

    # Surface node: radiative + convective heating + reaction enthalpy
    q_rad = _EPSILON * _SIGMA_SB * (_T_RAD**4 - T[0]**4)
    q_conv = _H_CONV * (_T_RAD - T[0])
    q_react = -total_ablation_rate * 1.5e6  # net exothermic surface reactions (J/kg * rate)
    q_cond_0 = k_nodes[0] * (T[1] - T[0]) * dx2_inv
    dy[6] = (q_rad + q_conv + q_react) / (_DX2 * rho_cp) + q_cond_0 / rho_cp

    # Interior nodes
    for i in range(1, _N_THERMAL - 1):
        k_avg_l = 0.5 * (k_nodes[i - 1] + k_nodes[i])
        k_avg_r = 0.5 * (k_nodes[i] + k_nodes[i + 1])
        q_cond = (k_avg_l * (T[i - 1] - T[i]) + k_avg_r * (T[i + 1] - T[i])) * dx2_inv
        dy[6 + i] = q_cond / rho_cp

    # Back face: insulated
    k_avg = 0.5 * (k_nodes[_N_THERMAL - 2] + k_nodes[_N_THERMAL - 1])
    dy[6 + _N_THERMAL - 1] = k_avg * (T[_N_THERMAL - 2] - T[_N_THERMAL - 1]) * dx2_inv / rho_cp

    # --- Surface recession ---
    ds_dt = total_ablation_rate / _RHO_CHAR2
    dy[14] = ds_dt                        # total recession
    tau_rec = 1.0  # smoothing timescale
    dy[15] = (ds_dt - rec_rate) / tau_rec  # recession rate (smoothed)
    dy[16] = max(ds_dt * 0.3, 0.0)        # char thickness grows (fraction of recession)

    # --- Sublayer decomposition ---
    for i in range(5):
        T_sub = T[min(i + 1, _N_THERMAL - 1)]
        dy[17 + i] = _arrhenius_rate(T_sub, alpha_sub[i])

    return dy