Battery Pack 24-Cell Module

def _bp24_cell_idx(r, c):
    return r * _BP24_NC + c

def _four_stage_decomposition(alpha_sei, alpha_ae, alpha_ca, alpha_el, T):
    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)

    d_sei = -k_sei * alpha_sei
    d_ae = k_ae * alpha_ae * (1.0 - alpha_ae)
    if alpha_ae < 1e-12 and alpha_sei < 0.15 - 1e-6:
        d_ae = k_ae * 1e-6
    d_ca = k_ca * (1.0 - alpha_ca)
    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_4stage(state, T_amb_eff):
    """4-stage dynamics for one cell. Returns (dy[8], q_gen)."""
    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 _bp24_rhs(t, y):
    cells = [y[i * _CELL_DIM_4STAGE:(i + 1) * _CELL_DIM_4STAGE]
             for i in range(_BP24_N_CELLS)]
    T = np.array([np.clip(cells[i][4], 250.0, 2000.0)
                  for i in range(_BP24_N_CELLS)])

    flux_offset = _BP24_CELL_BLOCK
    loss_offset = flux_offset + _BP24_N_EDGES
    q_edge = y[flux_offset:loss_offset]
    q_loss = y[loss_offset:]

    dy = np.zeros(_BP24_DIM)

    cell_dy = []
    for i in range(_BP24_N_CELLS):
        cdy, _ = _single_cell_4stage(cells[i], _T_AMB)
        cell_dy.append(cdy)

    # Horizontal edges
    edge_idx = 0
    for r in range(_BP24_NR):
        for c in range(_BP24_NC - 1):
            ci = _bp24_cell_idx(r, c)
            cj = _bp24_cell_idx(r, c + 1)
            q_cond = _K_CONTACT * _A_CONTACT * (T[ci] - T[cj]) / _D_GAP
            q_rad = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[ci]**4 - T[cj]**4)
            q_total = q_cond + q_rad

            cell_dy[ci][4] -= q_total / (_M_CELL * _CP)
            cell_dy[cj][4] += q_total / (_M_CELL * _CP)

            dy[flux_offset + edge_idx] = (q_total - q_edge[edge_idx]) / 0.1
            edge_idx += 1

    # Vertical edges
    for r in range(_BP24_NR - 1):
        for c in range(_BP24_NC):
            ci = _bp24_cell_idx(r, c)
            cj = _bp24_cell_idx(r + 1, c)
            q_cond = _K_CONTACT * _A_CONTACT * (T[ci] - T[cj]) / _D_GAP
            q_rad = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[ci]**4 - T[cj]**4)
            q_total = q_cond + q_rad

            cell_dy[ci][4] -= q_total / (_M_CELL * _CP)
            cell_dy[cj][4] += q_total / (_M_CELL * _CP)

            dy[flux_offset + edge_idx] = (q_total - q_edge[edge_idx]) / 0.1
            edge_idx += 1

    # Pack cell derivatives and loss trackers
    for i in range(_BP24_N_CELLS):
        dy[i * _CELL_DIM_4STAGE:(i + 1) * _CELL_DIM_4STAGE] = cell_dy[i]
        q_loss_actual = _H_CONV * _A_SURF * (T[i] - _T_AMB)
        dy[loss_offset + i] = (q_loss_actual - q_loss[i]) / 1.0

    return dy