Battery Thermal Runaway 4stage

def _arrhenius_rate(A, E, T):
    """Safe Arrhenius rate: A * exp(-E/(R*T)) with clamped exponent."""
    arg = np.clip(E / (_R_GAS * T), 0.0, _EXP_CLAMP)
    return A * np.exp(-arg)

def _four_stage_decomposition(alpha_sei, alpha_ae, alpha_ca, alpha_el, T):
    """Compute the four decomposition rates and total volumetric heat.

    Returns (d_sei, d_ae, d_ca, d_el, q_dot) where q_dot is W/kg (mass-specific).
    """
    k_sei = _arrhenius_rate(_A_SEI, _E_SEI, T)
    k_ae = _arrhenius_rate(_A_AE, _E_AE, T)
    k_ca = _arrhenius_rate(_A_CA, _E_CA, T)
    k_el = _arrhenius_rate(_A_EL, _E_EL, T)

    # SEI: consumed (α decreases)
    d_sei = -k_sei * alpha_sei

    # Anode-electrolyte: consumed fraction increases
    d_ae = k_ae * alpha_ae * (1.0 - alpha_ae)
    # At α_ae=0 this would give zero rate, but the seed is the
    # SEI decomposition products — use a small baseline nucleation
    # once SEI has started decomposing.
    if alpha_ae < 1e-12 and alpha_sei < 0.15 - 1e-6:
        d_ae = k_ae * 1e-6

    # Cathode: decomposed fraction increases
    d_ca = k_ca * (1.0 - alpha_ca)

    # Electrolyte: decomposed fraction increases
    d_el = k_el * (1.0 - alpha_el)

    q_dot = (_Q_SEI * _W_SEI * abs(d_sei)
             + _Q_AE * _W_AE * d_ae
             + _Q_CA * _W_CA * d_ca
             + _Q_EL * _W_EL * d_el)

    return d_sei, d_ae, d_ca, d_el, q_dot

def _btr_4stage_rhs(t, y):
    alpha_sei = np.clip(y[0], 0.0, 1.0)
    alpha_ae  = np.clip(y[1], 0.0, 1.0)
    alpha_ca  = np.clip(y[2], 0.0, 1.0)
    alpha_el  = np.clip(y[3], 0.0, 1.0)
    T         = np.clip(y[4], 250.0, 2000.0)
    Q_total   = y[5]
    P_gas     = y[6]
    R_int     = y[7]

    d_sei, d_ae, d_ca, d_el, q_dot = _four_stage_decomposition(
        alpha_sei, alpha_ae, alpha_ca, alpha_el, T,
    )

    # Temperature
    q_gen = q_dot * _M_CELL
    q_cool = _H_CONV * _A_SURF * (T - _T_AMB)
    dT = (q_gen - q_cool) / (_M_CELL * _CP)

    # Cumulative heat
    dQ = q_gen

    # Gas pressure: gas moles proportional to electrolyte decomposition
    # n_gas ≈ α_el * n_gas_max;  dn/dt = n_gas_max * d_el
    n_gas_max = 0.01  # mol of gas at full electrolyte decomposition
    n_gas = alpha_el * n_gas_max
    dn_gas = n_gas_max * d_el
    dP = (dn_gas * _R_GAS * T + n_gas * _R_GAS * dT) / _V_HEAD

    # Internal resistance growth
    dR = _R0 * (0.5 * abs(d_sei) + 2.0 * d_ae)

    dy = np.empty(8)
    dy[0] = d_sei
    dy[1] = d_ae
    dy[2] = d_ca
    dy[3] = d_el
    dy[4] = dT
    dy[5] = dQ
    dy[6] = dP
    dy[7] = dR
    return dy