I am solving some time dependent equations that deal with change in concentration of ions while charging and discharging. (concentrations of c2 ,c3,c4, c5 ions)

ex: dc2dt = ((I/96485.0) - (k2*c2*S) - (2*k5*c5*S) - (k4*c4*S))/V -- change in c2 concentration (V is volume)

Charging and discharging end points are marked by

  • voltage cutoffs(Charging end point : 2V, discharging end point :0.2V) or
  • by highest concentration of ions (c2,c3,c4,c5 concentration shouldn't go above 2.0 and below 0)

Voltage also varies with the concentration of c2 c3 c4 c5 as: E = E0 + ((8.314*T*np.log((c2*c5)/(c3*c4)))/96485.0) -- Nernst Equation

I am using the odeint function from Scipy.integrate but I don't know what the end point 'time' of charging and discharging should be. And I want the time to be driven by the voltage limits.

This is what I have tried, I have put the time spans such that the concentrations are within the limits. But I want some method to make the time vector be dependent on the voltage function mentioned above.


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)

I have looked through odeint documentation but could'nt figure it out. Any suggestion or leads will be helpful.

0

There are 0 best solutions below