Cantera:エンジン筒内燃焼解析

Canteraでエンジン筒内計算を行ってみる。ゼロ次元の反応器でピストンを動かし、完全予混合のn-Heptaneの自着火の解析を行う。今回は簡単のために、吸排気や熱逃げはない閉じた系とする。

計算モデル

筒内は、ゼロ次元の理想気体反応器IdealGasReactorを使う。ピストンとして移動できるWallを定義し、反応器の体積を変化させることで、エンジン筒内をモデル化する。

計算プログラム

ソースコード(GitHub)

# 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.cti'

筒内の初期温度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.set_velocity(piston_speed)  # piston speed

ピストンの設定をする。pitonを壁Wallで定義する。Wallは移動することができる。piston.areaでピストン頂面積を設定する。piston_speedは時間を入れるとピストン速度を返す関数で、piston.set_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ファイルにも結果を出力する。

熱発生率は、生成物の熱発生から計算している。

最後に圧力、温度、熱発生率をグラフに出力している。上死点付近で着火し圧力、温度が上昇しているのがわかる。

スポンサーリンク
科学技術計算のご相談は「キャットテックラボ」へ

科学技術計算やCAEに関するご相談、計算用プログラムの開発などお困りのことは「株式会社キャットテックラボ」へお問い合わせください。

お問い合わせはこちら

フォローする

スポンサーリンク