Zach Thomas Zach Thomas - 2 months ago 13
Python Question

Error in scipy odeint use for system of ODE's

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:

enter image description here

And here is the output of my code above:

enter image description here
enter image description here

Answer

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.

Comments