Battery Pack 6-Cell Echem

def _butler_volmer(i0, eta, T):
    f = _F / (_R_GAS * T)
    arg_a = np.clip(_ALPHA_BV * f * eta, -_EXP_CLAMP, _EXP_CLAMP)
    arg_c = np.clip(-_ALPHA_BV * f * eta, -_EXP_CLAMP, _EXP_CLAMP)
    return i0 * (np.exp(arg_a) - np.exp(arg_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 _radial_diffusion_3node(c, D_s, r_p):
    dr = r_p / 2.0
    dc = np.zeros(3)
    dc[0] = 6.0 * D_s * (c[1] - c[0]) / (dr * dr)
    r_m = dr
    flux_out = D_s * ((r_m + 0.5 * dr) ** 2) * (c[2] - c[1]) / dr
    flux_in = D_s * ((r_m - 0.5 * dr) ** 2) * (c[1] - c[0]) / dr
    dc[1] = (flux_out - flux_in) / (r_m * r_m * dr)
    r_s = r_p
    flux_in_s = D_s * ((r_s - 0.5 * dr) ** 2) * (c[2] - c[1]) / dr
    dc[2] = -flux_in_s / (r_s * r_s * dr)
    return dc

def _single_cell_echem(state):
    """Full echem cell dynamics. Returns dy[16] and 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)
    eta_a = state[5]
    eta_c = state[6]
    c_a = np.clip(state[7:10], 0.0, _CS_MAX_A)
    c_c = np.clip(state[10:13], 0.0, _CS_MAX_C)
    Q_total = state[13]
    P_gas = state[14]
    R_int = state[15]

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

    i_bv_a = _butler_volmer(_I0_ANODE, eta_a, T)
    i_bv_c = _butler_volmer(_I0_CATHODE, eta_c, T)

    tau_relax = 1.0
    theta_a = c_a[2] / _CS_MAX_A
    theta_c = c_c[2] / _CS_MAX_C
    U_a = 0.6 - 0.5 * theta_a
    U_c = 4.2 - 0.8 * theta_c
    V_cell = U_c - U_a
    d_eta_a = (-eta_a + (V_cell * 0.5 - U_a)) / tau_relax
    d_eta_c = (-eta_c + (U_c - V_cell * 0.5)) / tau_relax

    dc_a = _radial_diffusion_3node(c_a, _DS_ANODE, _RP_ANODE)
    dc_c = _radial_diffusion_3node(c_c, _DS_CATHODE, _RP_CATHODE)
    dc_a[2] += i_bv_a / _F
    dc_c[2] -= i_bv_c / _F

    q_echem = abs(i_bv_a * eta_a) + abs(i_bv_c * eta_c)
    q_gen = q_dot_decomp * _M_CELL + q_echem * _A_SURF
    q_cool = _H_CONV * _A_SURF * (T - _T_AMB)
    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(16)
    dy[0] = d_sei
    dy[1] = d_ae
    dy[2] = d_ca
    dy[3] = d_el
    dy[4] = dT
    dy[5] = d_eta_a
    dy[6] = d_eta_c
    dy[7:10] = dc_a
    dy[10:13] = dc_c
    dy[13] = dQ
    dy[14] = dP
    dy[15] = dR
    return dy, q_gen

def _bp6e_rhs(t, y):
    cells = [y[i * _CELL_DIM_ECHEM:(i + 1) * _CELL_DIM_ECHEM]
             for i in range(_BP6E_N_CELLS)]
    T = np.array([np.clip(cells[i][4], 250.0, 2000.0)
                  for i in range(_BP6E_N_CELLS)])

    flux_offset = _BP6E_CELL_BLOCK
    contact_offset = flux_offset + _BP6E_N_FLUX
    q_tracked = y[flux_offset:contact_offset]
    Tc = y[contact_offset:contact_offset + _BP6E_N_FLUX]

    dy = np.zeros(_BP6E_DIM)

    cell_dy = []
    for i in range(_BP6E_N_CELLS):
        cdy, _ = _single_cell_echem(cells[i])
        cell_dy.append(cdy)

    for i in range(_BP6E_N_FLUX):
        q_cond = _K_CONTACT * _A_CONTACT * (T[i] - T[i + 1]) / _D_GAP
        q_rad = _SIGMA * _EMISSIVITY * _A_CONTACT * (T[i]**4 - T[i + 1]**4)
        q_total = q_cond + q_rad

        cell_dy[i][4] -= q_total / (_M_CELL * _CP)
        cell_dy[i + 1][4] += q_total / (_M_CELL * _CP)

        dy[flux_offset + i] = (q_total - q_tracked[i]) / 0.1
        dy[contact_offset + i] = (0.5 * (T[i] + T[i + 1]) - Tc[i]) / 1.0

    for i in range(_BP6E_N_CELLS):
        dy[i * _CELL_DIM_ECHEM:(i + 1) * _CELL_DIM_ECHEM] = cell_dy[i]

    return dy