Ablation Recession Coupled

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 _recession_coupled_rhs(t, y):
    dy = np.zeros(15)

    T = np.clip(y[0:6], 200.0, 5000.0)
    alpha = np.clip(y[6:9], 0.0, 1.0)
    gas_flux = y[9:12]
    s = max(y[12], 0.0)          # recession (non-negative)
    T_surf = np.clip(y[13], 200.0, 5000.0)
    char_thick = max(y[14], 0.0)

    L_eff = max(_L_TPS - s, 1.0e-4)  # remaining TPS thickness, bounded away from zero
    L_inv = 1.0 / L_eff
    L_inv2 = L_inv * L_inv

    # --- Surface recession rate (B' formulation) ---
    B_prime = 0.5 * np.exp(-_EA_ABLATION / (_R_GAS * T_surf))
    m_dot_abl = B_prime * _RHO_E * _U_E * _C_H
    ds_dt = m_dot_abl / _RHO_CHAR

    # --- Thermal conductivity and diffusivity ---
    # Char fraction at each node: nodes 0,1 assumed fully charred near surface;
    # nodes 2,3,4 use the interior decomposition state; node 5 is virgin.
    alpha_full = np.zeros(_NT3)
    alpha_full[0] = 1.0
    alpha_full[1] = 1.0
    alpha_full[2:5] = alpha
    alpha_full[5] = 0.0

    k_eff = _thermal_conductivity(alpha_full)
    rho_cp = _RHO_VIRGIN * _CP
    kappa = k_eff / rho_cp  # thermal diffusivity per node

    # --- Energy equation in moving frame ---
    # dT/dt = κ * d²T/dξ² / (L-s)² + (ds/dt)*ξ/(L-s)*dT/dξ + Q_pyro/(ρ*cp)*dα/dt
    dxi = _DXI

    for i in range(_NT3):
        xi = i * dxi  # transformed coordinate

        # d²T/dξ² via finite differences
        if i == 0:
            # Surface node: radiative + convective + ablation enthalpy
            q_rad = _EPSILON * _SIGMA_SB * (_T_RAD**4 - T[0]**4)
            q_conv = _H_CONV * (_T_RAD - T[0])
            q_abl = -m_dot_abl * 3.0e6  # ablation enthalpy (J/kg)
            # One-sided second derivative
            d2T = (T[1] - 2.0 * T[0] + T[0]) / (dxi * dxi)  # ghost = T[0] (Neumann-like)
            dT_dxi = (T[1] - T[0]) / dxi
            dT = kappa[0] * d2T * L_inv2 + ds_dt * xi * L_inv * dT_dxi
            dT += (q_rad + q_conv + q_abl) / (_DXI * L_eff * rho_cp)
        elif i == _NT3 - 1:
            # Back face: insulated (dT/dξ = 0 at ξ=1)
            d2T = (T[i - 1] - T[i]) / (dxi * dxi)  # ghost T[N] = T[N-1]
            dT_dxi = 0.0
            dT = kappa[i] * d2T * L_inv2
        else:
            d2T = (T[i - 1] - 2.0 * T[i] + T[i + 1]) / (dxi * dxi)
            dT_dxi = (T[i + 1] - T[i - 1]) / (2.0 * dxi)
            dT = kappa[i] * d2T * L_inv2 + ds_dt * xi * L_inv * dT_dxi

        # Pyrolysis source at interior nodes
        if 2 <= i <= 4:
            j = i - 2
            pyro_rate = _arrhenius_rate(T[i], alpha[j])
            dT += _Q_PYROLYSIS * pyro_rate / _CP

        dy[i] = dT

    # --- Interior decomposition (nodes 2, 3, 4 → indices 0, 1, 2 in alpha) ---
    for j in range(3):
        T_node = T[j + 2]
        dy[6 + j] = _arrhenius_rate(T_node, alpha[j])

    # --- Gas mass flux at 3 interior points ---
    K_perm_avg = _permeability(alpha)
    P_grad_scale = 5000.0  # Pa/m characteristic pressure gradient from pyrolysis
    for j in range(3):
        target_flux = -(K_perm_avg[j] / _MU_GAS) * P_grad_scale * dy[6 + j]
        tau_f = 0.05
        dy[9 + j] = (target_flux - gas_flux[j]) / tau_f

    # --- Recession and char thickness ---
    dy[12] = ds_dt

    # Surface temperature tracks grid node 0 with slight lag (separate ODE for stiffness)
    tau_surf = 0.1
    dy[13] = (T[0] - T_surf) / tau_surf

    # Char thickness grows as decomposition front advances
    mean_alpha = np.mean(alpha)
    dy[14] = ds_dt * 0.5 + _arrhenius_rate(T[2], mean_alpha) * _DXI * L_eff * 0.1

    return dy