如何停止积分时间相关的 ODE，但时间跨度的上限取决于另一个变量

Question

我正在求解一些与时间相关的方程式，这些方程式处理充电和放电时离子浓度的变化。 （c2、c3、c4、c5 离子的浓度） 例如：dc2dt = ((I/96485.0) - (k2c2S) - (2k5c5S) - (k4c4*S))/V -- c2 浓度的变化（V 是体积）

充放电终点用

标示

• 电压截止（充电终点：2V，放电终点：0.2V）或
• 按最高离子浓度（c2、c3、c4、c5 浓度不应高于 2.0 或低于 0）

电压也随着 c2 c3 c4 c5 的浓度变化为： E = E0 + ((8.314Tnp.log((c2c5)/(c3c4)))/96485.0) -- 能斯特方程

我正在使用 Scipy.integrate 的 odeint 函数，但我不知道充电和放电的终点“时间”应该是多少。我希望时间由电压限制驱动。

这是我尝试过的方法，我设置的时间跨度使浓度在限制范围内。但是我想要一些方法让时间向量依赖于上面提到的电压函数。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

k2 = 6.90*(10**-7)  # dms-1
k3 = 2.54*(10**-7)  # dms-1
k4 = 5.37*(10**-7)  # dms-1
k5 = 4.64*(10**-7)  # dms-1

c2 = 0.01
c5 = 0.01
c4 = 1.99
c3 = 1.99

E_upper = 1.7   # V
E_lower = 1.1   # V
E_0 = 1.4       # V
A = 3000.0      # cm2 electrode area
S = 30.0        # dm2 membrane area
d = 0.00127     # dm membrane length
V = 1.0         # Litre
T = 298.15      # K
R = 2.0         # Cell Resistance ohm*cm2
I = 180.0       # Current A



def charging_odes(c, t):
    c2 = c[0]
    c3 = c[1]
    c4 = c[2]
    c5 = c[3]

    dc2dt = ((I/96485.0) - (k2*c2*S) - (2*k5*c5*S) - (k4*c4*S))/V
    dc3dt = (-(I/96485.0) - (k3*c3*S) + (3*k5*c5*S) + (2*k4*c4*S))/V
    dc4dt = (-(I/96485.0) - (k4*c4*S) + (3*k2*c2*S) + (2*k3*c3*S))/V
    dc5dt = ((I/96485.0) - (k5*c5*S) - (2*k2*c2*S) - (k3*c3*S))/V
    return [dc2dt, dc3dt, dc4dt, dc5dt]


c0 = [0.01, 1.99, 1.99, 0.01]
cc = [0.01, 1.99, 1.99, 0.01]


def discharging_odes(c, t):

    c2 = c[0]
    c3 = c[1]
    c4 = c[2]
    c5 = c[3]

    dc2dt = (-(I/96485.0) - (k2*c2*S) - (2*k5*c5*S) - (k4*c4*S))/V
    dc3dt = ((I/96485.0) - (k3*c3*S) + (3*k5*c5*S) + (2*k4*c4*S))/V
    dc4dt = ((I/96485.0) - (k4*c4*S) + (3*k2*c2*S) + (2*k3*c3*S))/V
    dc5dt = (-(I/96485.0) - (k5*c5*S) - (2*k2*c2*S) - (k3*c3*S))/V

    return [dc2dt, dc3dt, dc4dt, dc5dt]


# print(charging_odes(c0, 0))
# declare a time vector

def single_charge_plot():
    t = np.linspace(0, 1020, 50)

    conc = odeint(charging_odes, c0, t)

    c2 = conc[:, 0]
    c3 = conc[:, 1]
    c4 = conc[:, 2]
    c5 = conc[:, 3]

    # plot the results

    plt.semilogy(t, c2)
    plt.semilogy(t, c3)
    plt.semilogy(t, c4)
    plt.semilogy(t, c5)
    plt.show()


def single_discharge_plot():
    t = np.linspace(0, 1020, 50)

    conc = odeint(discharging_odes, c0, t)

    c2 = conc[:, 0]
    c3 = conc[:, 1]
    c4 = conc[:, 2]
    c5 = conc[:, 3]

    print(c2)
    # plot the results

    plt.semilogy(t, c2)
    plt.semilogy(t, c3)
    plt.semilogy(t, c4)
    plt.semilogy(t, c5)
    plt.show()


def single_charge(cc):
    t = np.linspace(0, 1080, 50)

    conc = odeint(charging_odes, cc, t)

    c2 = conc[:, 0]
    c3 = conc[:, 1]
    c4 = conc[:, 2]
    c5 = conc[:, 3]

    conc_for_discharge = [c2[-1], c3[-1], c4[-1], c5[-1]]

    return conc_for_discharge


def single_discharge(cd):
    t = np.linspace(0, 950, 50)

    conc = odeint(discharging_odes, cd, t)

    c2 = conc[:, 0]
    c3 = conc[:, 1]
    c4 = conc[:, 2]
    c5 = conc[:, 3]

    conc_for_charge = [c2[-1], c3[-1], c4[-1], c5[-1]]

    return conc_for_charge


def single_cycle(cc):
    cd = single_charge(cc)
    print(cd)
    cc = single_discharge(cd)
    print(cc)


single_cycle(cc)

我已经查看了 odeint 文档，但无法弄明白。任何建议或线索都会有所帮助。