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.