2017-08-08 83 views
1

我正在使用Spark燃燒引擎模型,並且因爲我使用python來模擬燃燒的一些原因。我試圖使用ODE的解算器,但收益完全超出了現實。我發現圓柱體積的集成是錯誤的。我已經嘗試過使用「odeint」和「ode」解算器,但結果是一樣的。 該代碼顯示了Volume與theta的導數並整合以查找音量。我把分析方程進行比較。Python中的ODE求解器錯誤

OBS:我有一個類似的問題,使用Matlab,但是當我嘗試在三角函數中使用度。當我改變弧度時,問題解決了。 代碼如下:

from scipy.integrate import odeint 
from scipy.integrate import ode 
from scipy import integrate 
import math 
import sympy 
from sympy import sqrt, sin, cos, tan, atan 
from pylab import * 
from RatesComp import * 
V_real=np.zeros((100)) 

def Volume(V,theta): 
    V_sol = V[0] 
    dVdtheta = Vtdc*(r-1)/2 *(sin(theta) + eps/2*sin(2*theta)/sqrt(1-(eps**2)*sin(theta)**2)) 
    return [dVdtheta] 

#Geometry 
eps = 0.25;  # half stroke to rod ratio, s/2l 
r = 10;   # compression ratio 
Vtdc = 6.9813e-05 # volume at TDC 

# Initial Conditions 
theta0 = - pi 
V_init = 0.0006283 
theta = linspace(-pi,pi,100) 
solve = odeint(Volume, V_init, theta) 

# Analytical Result 
Size = len(theta) 

for i in range(0, Size,1): 
    V_real[i] = Vtdc*(1+(r-1)/2*(1-cos(theta[i])+ 1/eps*(1-(1-(eps**2)*sin(theta[i])**2)**0.5))) 

figure(1) 
plot(theta, solve[:,0],label="Comput") 
plot(theta, V_real[0:Size],label="Real") 
ylabel('Volume [m^3]') 
xlabel('CA [Rad]') 
legend() 
grid(True) 
show() 

我展示的是圓柱體的體積的圖。結果實和計算

https://i.stack.imgur.com/MIoEV.png

能有人爲什麼這個問題發生的信息幫助嗎?

+0

如果增加矢量theta的大小,你會得到相同的行爲嗎? – jmoon

回答

2

顯然你使用python2。 r=10的聲明給出rinteger這種類型導致在'真正'解決方案中的(r-1)/2中出現不需要的整數除法。在派生函數中,浮點值爲Vtdc作爲產品中的第一個因子,之後整個產品評估處於浮動狀態。

因此更改爲r=10.0或使用(r-1.0)/20.5*(r-1)


而且你還要設置V_init = r*Vtdc爲是V_real(-pi)值。

+0

謝謝您的幫助。代碼現在運行。 –