Ablation Charring Pyrolysis

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 _permeability(alpha):
    """Darcy permeability increasing with char fraction."""
    return 1.0e-12 * (1.0 + 10.0 * np.clip(alpha, 0.0, 1.0))

def _thermal_conductivity(alpha):
    """Effective conductivity: 0.5 (virgin) → 2.0 (char) W/(m·K)."""
    return 0.5 + 1.5 * np.clip(alpha, 0.0, 1.0)

def _charring_pyrolysis_rhs(t, y):
    dy = np.zeros(18)

    alpha = np.clip(y[0:5], 0.0, 1.0)
    T = np.clip(y[5:10], 200.0, 5000.0)
    gas_flux = y[10:13]
    P = np.maximum(y[13:18], 1000.0)

    # --- Decomposition (Arrhenius) ---
    rate = _arrhenius_rate(T, alpha)
    dy[0:5] = rate

    # --- Gas generation from pyrolysis ---
    gas_gen = np.zeros(_NL1)
    for i in range(_NL1):
        gas_gen[i] = _RHO_VIRGIN * rate[i] * _V_LAYER

    # --- Gas pressure (ideal-gas-like response to generation + flow) ---
    K_perm = _permeability(alpha)
    dx_inv = 1.0 / _DX1

    for i in range(_NL1):
        dP_dx_in = 0.0
        dP_dx_out = 0.0
        if i > 0:
            dP_dx_in = (P[i - 1] - P[i]) * dx_inv
        if i < _NL1 - 1:
            dP_dx_out = (P[i] - P[i + 1]) * dx_inv

        flux_in = (K_perm[i] / _MU_GAS) * dP_dx_in if i > 0 else 0.0
        flux_out = (K_perm[i] / _MU_GAS) * dP_dx_out if i < _NL1 - 1 else 0.0

        # Layer 0 (surface) vents to atmosphere
        if i == 0:
            flux_out_vent = (K_perm[0] / _MU_GAS) * (P[0] - _P_ATM) * dx_inv
            net_flux = flux_in - flux_out_vent
        else:
            net_flux = flux_in - flux_out

        # P evolves from generation and net flux divergence
        # dp/dt ~ (gas_gen * R_specific * T / V - net_flux_divergence * P) / (rho * V)
        dp = gas_gen[i] * _R_GAS * T[i] / (0.029 * _V_LAYER) + net_flux * dx_inv * 1e3
        dy[13 + i] = dp

    # --- Gas mass flux at 3 interior interfaces (between layers 0-1, 1-2, 2-3) ---
    for j in range(3):
        i = j  # interface between layer j and j+1
        K_avg = 0.5 * (K_perm[i] + K_perm[i + 1])
        dP = P[i] - P[i + 1]
        flux = -(K_avg / _MU_GAS) * dP * dx_inv
        tau_flux = 0.01  # relaxation timescale for flux
        dy[10 + j] = (flux - gas_flux[j]) / tau_flux

    # --- Energy equation ---
    k_eff = _thermal_conductivity(alpha)
    rho_cp = _RHO_VIRGIN * _CP  # simplified (should interpolate virgin→char)

    for i in range(_NL1):
        # Conduction (finite differences)
        if i == 0:
            # Surface: radiative + convective heating from environment
            q_rad = _EPSILON * _SIGMA_SB * (_T_RAD**4 - T[0]**4)
            q_conv = _H_CONV * (_T_RAD - T[0])
            q_in = q_rad + q_conv
            q_cond = k_eff[0] * (T[1] - T[0]) / (_DX1 * _DX1)
            dT = (q_in / _DX1 + q_cond) / rho_cp
        elif i == _NL1 - 1:
            # Insulated back face
            q_cond = k_eff[i] * (T[i - 1] - T[i]) / (_DX1 * _DX1)
            dT = q_cond / rho_cp
        else:
            # Interior: central difference
            q_cond = k_eff[i] * (T[i - 1] - 2.0 * T[i] + T[i + 1]) / (_DX1 * _DX1)
            dT = q_cond / rho_cp

        # Pyrolysis heat source (endothermic)
        dT += _Q_PYROLYSIS * rate[i] / _CP

        # Gas convective cooling within pores
        if i < 3:
            gas_cooling = _H_GAS * gas_flux[min(i, 2)] * (T[i] - _T_GAS_IN)
            dT -= gas_cooling / rho_cp

        dy[5 + i] = dT

    return dy