2016-09-20 84 views
0

我想重现下面显示的两个图,它们是作为时间函数的温度和浓度分布图。我检查了我的方法和代码一百万次,但似乎无法找到它的错误,但我无法再现这些图。除了CA和T之外,所有的值都是常量。它可能是来自scipy的odeint的准确性问题吗?任何帮助将非常感激!用于ODE系统的scipy odeint的错误

这两个方程如下:

DCA/DT = Q *(CAI - CA)/ V - K * CA

的dT/dt的= W *(钛 - T)/(V p)+ d_HRķCA /(p C)+ UA *(TC - T)/(V p C)

的代码:

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() 

这里是图应该生产什么:

enter image description here

这里是我的代码的输出上面:

enter image description here enter image description here

+0

仔细检查此公式:'k = k0 * np.exp(8750 * 1/T)'。在一个完整的黑暗中,我改变了指数内的符号,所以'k = k0 * np.exp(-8750 * 1/T)',得到的结果看起来类似于你预期的情节显示。 –

+0

你是上帝......这是完全正确的!哇,真是一种解脱。如果您将此作为答案提交,我很乐意提交您的答案作为答案。谢谢! –

+0

好的,提交答案。 –

回答

1

显然有公式中的标志错误k = k0*np.exp(8750*1/T)。如果您将其更改为k = k0*np.exp(-8750*1/T),则会得到您所期望的情节。