211 オットーサイクル¶
Example: Spark-Ignition. Four-Stroke Engine
読み込み¶
In [1]:
import math
import sys
sys.path.append(r'~\Lib\site-packages') # Path to library directory
import CoolProp
import CoolProp.CoolProp as CP
from CoolProp.Plots import PropertyPlot
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
パラメータ定義¶
In [2]:
WF = 'air'
K = 273.15
bore = 0.1 # m
stroke = 0.09 #m
CR = 8.3 # Compression ratio
N_cyl = 4 # number of cylinder
N = 60 # engine speed 1/s
AF = 16 # air-fuel ratio
HC = 44*10**6 # heat of combusion J/kg
T_amb = 32 + K # outdoor air temperature
P_atm = 100*10**3 # atomospheric air temperature Pa
Vol_dis_cyl = math.pi*bore**2 *stroke/4 # Displacement of each cylinder
Vol_cl = Vol_dis_cyl/(CR-1) # clearance volume
Vol_BDC = Vol_dis_cyl+Vol_cl #bottom dead center volume
Vol_TDC = Vol_cl # top dead center volume
データフレーム用意¶
In [3]:
index = list(range(1,4+1))
states = pd.DataFrame(columns=['T', 'P', 's', 'v', 'Vol', 'u', 'm'], index=index)#K, Pa, J/kg/K, m3/kg, m3/mol?, J/kg, kg/mol
Point 1¶
In [4]:
states.ix[1, 'T'] = T_amb
states.ix[1, 'P'] = P_atm
states.ix[1, 's'] = CP.PropsSI('S', 'T', states.ix[1, 'T'], 'P', states.ix[1, 'P'], WF)
states.ix[1, 'u'] = CP.PropsSI('U', 'T', states.ix[1, 'T'], 'P', states.ix[1, 'P'], WF)
states.ix[1, 'v'] = 1/CP.PropsSI('D', 'T', states.ix[1, 'T'], 'P', states.ix[1, 'P'], WF)
states.ix[1, 'Vol'] = Vol_BDC
states.ix[1, 'm'] = states.ix[1, 'Vol'] / states.ix[1, 'v']
Point 2¶
In [5]:
states.ix[2, 'm'] = states.ix[1, 'm']
states.ix[2, 'Vol'] = Vol_TDC
states.ix[2, 'v'] = states.ix[2, 'Vol'] / states.ix[2, 'm']
states.ix[2, 's'] = states.ix[1, 's']
states.ix[2, 'u'] = CP.PropsSI('U', 'D', 1/states.ix[2, 'v'], 'S', states.ix[2, 's'], WF)
states.ix[2, 'T'] = CP.PropsSI('T', 'D', 1/states.ix[2, 'v'], 'S', states.ix[2, 's'], WF)
states.ix[2, 'P'] = CP.PropsSI('P', 'D', 1/states.ix[2, 'v'], 'S', states.ix[2, 's'], WF)
Work etc.¶
In [6]:
W_comp = states.ix[2, 'm'] * states.ix[2, 'u'] - states.ix[1, 'm'] * states.ix[1, 'u'] #Compression work
m_in = (Vol_BDC - Vol_TDC) / (1/CP.PropsSI('D', 'T', T_amb, 'P', P_atm, WF)) #mass of incoming air
m_f = m_in / (AF+1) # mass of fuel
print('m_in:{:.02E}'.format(m_in))
print('m_f:{:.02E}'.format(m_f))
m_in:8.07E-04
m_f:4.75E-05
State 3¶
In [7]:
states.ix[3, 'm'] = states.ix[2, 'm']
states.ix[3, 'Vol'] = Vol_TDC
states.ix[3, 'v'] = states.ix[3, 'Vol'] / states.ix[3, 'm']
states.ix[3, 'u'] = (m_f * HC - states.ix[2, 'm'] * states.ix[2, 'u']) / states.ix[3, 'm'] # from Energy balance
states.ix[3, 'T'] = CP.PropsSI('T', 'U', states.ix[3, 'u'], 'D', 1/states.ix[3, 'v'], WF)
states.ix[3, 'P'] = CP.PropsSI('P', 'U', states.ix[3, 'u'], 'D', 1/states.ix[3, 'v'], WF)
states.ix[3, 's'] = CP.PropsSI('S', 'U', states.ix[3, 'u'], 'D', 1/states.ix[3, 'v'], WF)
State 4¶
In [8]:
def find_P_from_v_s(WF, v, s, Pstart, Pstop, eps=1e-6):
v1 = 1/CP.PropsSI('D','P',Pstart,'S', s, WF)
v2 = 1/CP.PropsSI('D','P',Pstop,'S', s, WF)
while abs(v2-v1)/v > eps:
Pm = (Pstart+Pstop)/2.0
vm = 1/CP.PropsSI('D','P',Pm,'S',s,WF)
if (vm-v)*(v2-v) < 0: Pstart,v1 = Pm,vm
else: Pstop,v2 = Pm,vm
return Pm
states.ix[4, 'm'] = states.ix[3, 'm']
states.ix[4, 'Vol'] = Vol_BDC
states.ix[4, 'v'] = states.ix[4, 'Vol'] / states.ix[4, 'm']
states.ix[4, 's'] = states.ix[3, 's']
#states.ix[4, 'P'] = CP.PropsSI('P', 'D', 1/states.ix[4, 'v'], 'S', 1/states.ix[4, 's'], WF)
states.ix[4, 'P'] = find_P_from_v_s(WF, states.ix[4, 'v'], states.ix[4, 's'], states.ix[1, 'P'], states.ix[3, 'P'])
#states.ix[4, 'T'] = CP.PropsSI('T', 'D', 1/states.ix[4, 'v'], 'S', 1/states.ix[4, 's'], WF)
states.ix[4, 'T'] = CP.PropsSI('T', 'D', 1/states.ix[4, 'v'], 'P', states.ix[4, 'P'], WF)
#states.ix[4, 'u'] = CP.PropsSI('U', 'D', 1/states.ix[4, 'v'], 'S', 1/states.ix[4, 's'], WF)
states.ix[4, 'u'] = states.ix[3, 'u'] - states.ix[2, 'u']
結果テーブルの表示¶
In [9]:
W_exp = states.ix[4, 'm'] * states.ix[4, 'u'] - states.ix[3, 'm'] * states.ix[3, 'u'] # expansion work
W_net =W_exp - W_comp # net work per power stroke
eta = W_net / (m_f * HC)
print('W_net {:.02f} = W_exp {:.02f} - W_comp {:.02f}'.format(W_net, W_exp, W_comp, ))
print('eta: {:.02f}'.format(eta))
display(states)
W_net -848.25 = W_exp -581.94 - W_comp 266.31
eta: -0.41
| T | P | s | v | Vol | u | m | |
|---|---|---|---|---|---|---|---|
| 1 | 305.15 | 100000 | 3907.63 | 0.875698 | 0.000803688 | 343914 | 0.000917768 |
| 2 | 695.502 | 1.90557e+06 | 3907.63 | 0.105506 | 9.68299e-05 | 634085 | 0.000917768 |
| 3 | 1831.19 | 5.02999e+06 | 4754.13 | 0.105506 | 9.68299e-05 | 1.64231e+06 | 0.000917768 |
| 4 | 929.61 | 305023 | 4754.13 | 0.875698 | 0.000803688 | 1.00823e+06 | 0.000917768 |
PVグラフ¶
In [10]:
nb_points = 1000
P12 = np.linspace(states.ix[1, 'P'], states.ix[2, 'P'], nb_points)
P23 = np.linspace(states.ix[2, 'P'], states.ix[3, 'P'], nb_points)
P34 = np.linspace(states.ix[3, 'P'], states.ix[4, 'P'], nb_points)
P41 = np.linspace(states.ix[4, 'P'], states.ix[1, 'P'], nb_points)
v12 = 1/CP.PropsSI('D', 'P', P12, 'S', states.ix[1, 's'], WF)
v23 = np.linspace(states.ix[2, 'v'], states.ix[3, 'v'], nb_points)
v34 = 1/CP.PropsSI('D', 'P', P34, 'S', states.ix[3, 's'], WF)
v41 = np.linspace(states.ix[4, 'v'], states.ix[1, 'v'], nb_points)
L_P = [P12, P23, P34, P41]
L_v = [v12, v23, v34, v41]
for i in range(len(L_v)):
plt.plot(L_v[i]*1e3,L_P[i]/1e5,label='{}$\\to${}'.format(i+1,(i+1)%5+1))
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)
plt.xlabel('Volume $v$ m$^3/$kg')
plt.ylabel('Pressure bar (=x$10^2$ kPa)')
plt.grid()
plt.show()
Tsグラフ¶
In [11]:
nb_points = 1000
T12 = np.linspace(states.ix[1, 'T'], states.ix[2, 'T'], nb_points)
T23 = np.linspace(states.ix[2, 'T'], states.ix[3, 'T'], nb_points)
T34 = np.linspace(states.ix[3, 'T'], states.ix[4, 'T'], nb_points)
T41 = np.linspace(states.ix[4, 'T'], states.ix[1, 'T'], nb_points)
s12 = np.linspace(states.ix[1, 's'], states.ix[2, 's'], nb_points)
s23 = CP.PropsSI('S', 'T', T23, 'D', 1/states.ix[2, 'v'], WF)
s34 = np.linspace(states.ix[3, 's'], states.ix[4, 's'], nb_points)
s41 = CP.PropsSI('S', 'T', T41, 'D', 1/states.ix[1, 'v'], WF)
L_T = [T12, T23, T34, T41]
L_s = [s12, s23, s34, s41]
for i in range(len(L_v)):
plt.plot(L_s[i]*1e3,L_T[i]/1e5,label='{}$\\to${}'.format(i+1,(i+1)%5+1))
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)
plt.xlabel('Entropy $s$ J/kg/K')
plt.ylabel('Temperature $T$ K')
plt.grid()
plt.show()
