Battery Multicell Propagation

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 _single_cell_dynamics(state, T_amb_eff):
    """4-stage dynamics for one cell with effective ambient temperature.

    ``T_amb_eff`` replaces the global _T_AMB so coupling heat can be
    injected through the cooling term.

    Returns (dy[8], q_gen) where q_gen is the raw heat generation rate (W).
    """
    alpha_sei = np.clip(state[0], 0.0, 1.0)
    alpha_ae  = np.clip(state[1], 0.0, 1.0)
    alpha_ca  = np.clip(state[2], 0.0, 1.0)
    alpha_el  = np.clip(state[3], 0.0, 1.0)
    T         = np.clip(state[4], 250.0, 2000.0)

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

    q_gen = q_dot * _M_CELL
    q_cool = _H_CONV * _A_SURF * (T - T_amb_eff)
    dT = (q_gen - q_cool) / (_M_CELL * _CP)

    dQ = q_gen

    n_gas_max = 0.01
    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

    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, q_gen

def _multicell_rhs(t, y):
    # Unpack per-cell states
    cells = [y[i * _CELL_DIM:(i + 1) * _CELL_DIM] for i in range(_N_CELLS)]
    T = np.array([np.clip(cells[i][4], 250.0, 2000.0) for i in range(_N_CELLS)])

    # Inter-cell heat fluxes (state y[24:27])
    q_12 = y[24]   # heat flux cell 1→2
    q_23 = y[25]   # heat flux cell 2→3
    q_13 = y[26]   # heat flux cell 1→3 (radiation)

    # Contact temperatures (state y[27:30])
    Tc_12 = y[27]
    Tc_23 = y[28]
    Tc_13 = y[29]

    dy = np.zeros(30)

    # Compute single-cell dynamics with ambient cooling
    cell_dy = []
    for i in range(_N_CELLS):
        cdy, _ = _single_cell_dynamics(cells[i], _T_AMB)
        cell_dy.append(cdy)

    # Inter-cell conduction: q = k * A * ΔT / d
    q_cond_12 = _K_CONTACT * _A_CONTACT * (T[0] - T[1]) / _D_GAP
    q_cond_23 = _K_CONTACT * _A_CONTACT * (T[1] - T[2]) / _D_GAP

    # Inter-cell radiation: q = σ * ε * A * (T_i⁴ - T_j⁴)
    q_rad_12 = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[0]**4 - T[1]**4)
    q_rad_23 = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[1]**4 - T[2]**4)
    q_rad_13 = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[0]**4 - T[2]**4)

    # Total heat exchange
    q_total_12 = q_cond_12 + q_rad_12
    q_total_23 = q_cond_23 + q_rad_23

    # Inject coupling heat into cell temperature derivatives
    # Cell 1 loses heat to cells 2 and 3
    cell_dy[0][4] -= (q_total_12 + q_rad_13) / (_M_CELL * _CP)
    # Cell 2 gains from cell 1, loses to cell 3
    cell_dy[1][4] += (q_total_12 - q_total_23) / (_M_CELL * _CP)
    # Cell 3 gains from cell 2 and from cell 1 (radiation)
    cell_dy[2][4] += (q_total_23 + q_rad_13) / (_M_CELL * _CP)

    # Pack per-cell derivatives
    for i in range(_N_CELLS):
        dy[i * _CELL_DIM:(i + 1) * _CELL_DIM] = cell_dy[i]

    # Heat flux state derivatives (track actual fluxes for diagnostics)
    tau_flux = 0.1  # relaxation time for flux tracking, s
    dy[24] = (q_total_12 - q_12) / tau_flux
    dy[25] = (q_total_23 - q_23) / tau_flux
    dy[26] = (q_rad_13 - q_13) / tau_flux

    # Contact temperature dynamics (thermal mass of contact interface)
    # Thin interface: relaxes toward average of adjacent cell temps
    tau_contact = 1.0  # s
    dy[27] = (0.5 * (T[0] + T[1]) - Tc_12) / tau_contact
    dy[28] = (0.5 * (T[1] + T[2]) - Tc_23) / tau_contact
    dy[29] = (0.5 * (T[0] + T[2]) - Tc_13) / tau_contact

    return dy