I am trying to reproduce the two graphs shown below, which are temperature and concentration profiles as a function of time. I have checked my method and code a million times but cannot seem to find an error in it but I cannot reproduce those graphs. All values are constants except CA and T. Could it be an issue with the accuracy of odeint from scipy? Any help would be much appreciated!
The two equations are as follows:
dCA/dt = q*(CAi - CA)/V - k*CA
dT/dt = w*(Ti - T)/(Vp) + d_HRkCA/(pC) + UA*(Tc - T)/(VpC)
The code is:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def ODESolve(y, t, q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc):
CA, T = y
k = k0*np.exp(8750*1/T)
dydt = [q*(CAi - CA)/V - k*CA, w*(Ti - T)/(V*p) + \
dH_R*k*CA/(p*C) + UA*(Tc - T)/(V*p*C)]
return dydt
q = 100
CAi = 1.0
V = 100
p = 1000
C = .239
dH_R = 5*(10**4)
k0 = 7.2*(10**10)
UA = 5*10**4
CA0 = .5
T0 = 350
Ti = T0
w = p*q
y0 = [CA0, T0]
t = np.linspace(0, 20, 100)
Tc = 305
sol1 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))
Tc = 300
sol2 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))
Tc = 290
sol3 = odeint(ODESolve, y0, t, args = (q, CAi, V, k0, w, Ti, p, dH_R, C, UA, Tc))
plt.figure(1)
plt.plot(t, sol1[:,0], label = 'Tc = 305')
plt.plot(t, sol2[:,0], label = 'Tc = 300')
plt.plot(t, sol3[:,0], label = 'Tc = 290')
plt.ylim(ymax = 1, ymin = 0)
plt.title ('CA(t)')
plt.legend(loc = 'best')
plt.figure(2)
plt.plot(t, sol1[:,1], label = 'Tc = 305')
plt.plot(t, sol2[:,1], label = 'Tc = 300')
plt.plot(t, sol3[:,1], label = 'Tc = 290')
plt.ylim(ymax = 450, ymin = 300)
plt.legend(loc = 'best')
plt.title ('T(t)')
plt.show()
Here is what the graphs are supposed to produce:
And here is the output of my code above:
Apparently there is a sign error in the formula k = k0*np.exp(8750*1/T)
. If you change that to k = k0*np.exp(-8750*1/T)
, you'll get plots like those you expect.