※Canteraバージョン：3.0.0
Canteraでエンジン筒内計算を行ってみる。ゼロ次元の反応器でピストンを動かし、完全予混合のn-Heptaneの自着火の解析を行う。今回は簡単のために、吸排気や熱逃げはない閉じた系とする。
計算モデル
筒内は、ゼロ次元の理想気体反応器IdealGasReactorを使う。ピストンとして移動できるWallを定義し、反応器の体積を変化させることで、エンジン筒内をモデル化する。
計算プログラム
# Input Parameters
rpm = 600.0 # engine speed [rpm]
bore = 82.55 # bore diameter [mm]
stroke = 114.3 # stroke [mm]
cratio = 10.0 # compression ratio [-]
conrod = 200.0 # connecting rod [mm]エンジン緒言を設定する。エンジン回転数rpm、ボア径bore、ストロークstroke、圧縮比cratio、コンロッド長conrodを入力する。
# initial temperature, pressure, and equivalence ratio
T_ini = 350.0 # [K]
p_ini = 1.0e5 # [Pa]
phi = 0.33
# outer temperature, pressure, and composition
T_out = 300.0 # [K]
p_out = 1.0e5 # [Pa]
c_out = 'O2:1.0, N2:3.76'
# Reaction mechanism name
reaction_mechanism = 'reduced_247.yaml'筒内の初期温度T_ini、初期圧力p_init、当量比phiを入力する。外部の温度T_out、圧力p_out、組成c_outを入力する（今回は閉じた系なので関係ない）。反応メカニズムreaction_mechanismを入力する（以前リダクションしたn-Heptaneのメカニズムを使う）。
# Simulation time
ca_start = -144.0 # start CA [deg]
ca_end = 180.0 # end CA [deg]
ca_step = 0.01 # step CA [deg]
ca_out = 0.2 # output CA [deg]計算開始ca_start、終了ca_end、時間刻みca_step、結果出力刻みca_outをクランク角[deg]で入力する。ここでは、上死点が0[deg]である。
# load reaction mechanism
gas = ct.Solution(reaction_mechanism)
# define initial state
gas.TP = T_ini, p_ini
gas.set_equivalence_ratio(phi, 'NC7H16', 'O2:1.0, N2:3.76')
r = ct.IdealGasReactor(gas)
sim = ct.ReactorNet([r])
gas.TPX = T_out, p_out, c_out
outer = ct.Reservoir(gas)gasを定義し初期条件を設定する。IdealGasReactorで理想気体の反応器を定義する。outerは、外部の条件で後でピストンの設定に使う。
# convert time to crank angle [rad]
rps = rpm / 60.0
def crank_angle(t):
return 2.0 * np.pi * rps * t + ca_start * np.pi / 180.0時間を入れると、クランク角[rad]を返す関数crank_angle(t)を定義する。
# set up IC engine parameters
stroke *= 0.001
bore *= 0.001
conrod *= 0.001
area = 0.25 * np.pi * bore * bore
vol_h = stroke * area # volume cylinder
vol_c = vol_h / (cratio - 1.0) # volume combustion dome
ca = crank_angle(0.0) # initial CA
r_ca = stroke * 0.5 # crank radius
vol_ini= (r_ca + conrod - (r_ca * np.cos(ca) + np.sqrt(conrod**2 - r_ca**2 * np.sin(ca)**2))) * area + vol_c
r.volume = vol_ini # initial volumeエンジン緒言を設定する。areaはピストン頂面積、vol_hはシリンダ容積、vol_cは燃焼室容積、r_caはクランク半径、vol_iniは初期の体積を計算している。r.volumeで初期体積を設定する。
# set up piston
piston = ct.Wall(outer, r)
piston.area = area # piston area
def piston_speed(t):
ca = crank_angle(t)
return -2.0 * np.pi * rps * (r_ca * np.sin(ca) + r_ca**2 * np.sin(2.0 * ca) / 2.0 / np.sqrt(conrod**2 - r_ca**2 * np.sin(ca)**2))
piston.velocity = piston_speed # piston speedピストンの設定をする。pitonを壁Wallで定義する。Wallは移動することができる。piston.areaでピストン頂面積を設定する。piston_speedは時間を入れるとピストン速度を返す関数で、piston.velocityでピストン速度を設定している。
$$ V_{piston} = -\omega \left(r \sin \theta + \frac{r^2 \sin 2 \theta }{2 \sqrt{l^2 - r^2 \sin^2 \theta}} \right)$$
# set up time
t_sim = (ca_end - ca_start) / rps / 360.0 # simulation time
t_step = ca_step / rps / 360.0 # simulation time step
t_out = ca_out / rps / 360.0 # simulation output time
ttt = 0.0時間を設定する。t_simは全計算時間、t_stepは時間刻み、t_outは結果出力時間刻みを[秒]で計算している。
# set up output data arrays
states = ct.SolutionArray(r.thermo)
t = []
heat_release_rate = []
# output file
outfile = open('ic_engine.csv', 'w', newline="")
csvfile = csv.writer(outfile)
csvfile.writerow(['ca[deg]','P[bar]','T[K]','HR[J/deg]'])結果出力用の配列を用意して、結果出力用のCSVファイルを設定する。
# do simulation
for t_i in np.arange(0, t_sim, t_step):
sim.advance(t_i)
# write output data
if t_i >= ttt:
ca = crank_angle(t_i) * 180.0 / np.pi
t.append(ca)
states.append(r.thermo.state)
hr = -r.volume * ct.gas_constant * r.T * np.sum(gas.standard_enthalpies_RT * r.thermo.net_production_rates, 0)
hr = hr * t_step / ca_step
heat_release_rate.append(hr)
csvfile.writerow([ca, r.thermo.P / 1.0e5, r.T, hr])
ttt += t_out時間ループにより計算を実行する。sim.advanceで各時刻まで計算を進める。if文は、設定した結果出力刻み毎に結果を出力している。tにクランク角[deg]、statesに熱力学データ、heat_release_rateに熱発生率[J/deg]を出力する。CSVファイルにも結果を出力する。
熱発生率は、生成物の熱発生から計算している。
最後に圧力、温度、熱発生率をグラフに出力している。上死点付近で着火し圧力、温度が上昇しているのがわかる。
