1.燃焼の基礎(その2)

以下記事の続き。

1.7 Canteraによる断熱燃焼温度の計算

　断熱火炎温度を化学反応解析ソフトcanteraを使用して計算する。

各混合比における完全燃焼を想定した理論断熱火炎温度と、完全燃焼は起こらず熱解離を伴う化学平衡状態で反応を想定した平衡断熱燃焼温度を計算し、その結果を比較する。
対象とする化学種は、CH4として”gri30.yaml”というcanteraの反応機構モデルを使用して計算を行う。
反応機構モデルとは、燃焼反応の素過程とそのアレニウスパラメータを集めたものであり、"gri30.yamal"はメタンを中心とした燃焼反応メカニズム(NOx 生成と再燃焼を含む)として53 化学種，325反応素過程により構築され，77種の実験データ (ターゲット) に対して 31 のパラメータが最適化されたものであり，通常の実験で得られるマクロな物理量 (温度や圧力など) を再現するモデルとしての信頼度は相当に高い。（gri30.yamalはGRI-Mech3.0という反応機構モデルをcanteraに読み込み可能な形式で記述したものである。）

canteraを用いて、理論断熱火炎温度と平衡断熱火炎温度を算出するプログラムを以下に示す。canteraではpythonによるプログラムが可能なため、pythonにて記述している。

flame_temp.py

"""
Adiabatic flame temperature and equilibrium composition for a fuel/air mixture
as a function of equivalence ratio, including formation of solid carbon.

Requires: cantera >= 2.5.0, matplotlib >= 2.0
Keywords: equilibrium, combustion, multiphase
"""

import cantera as ct
import numpy as np
import sys
import csv

##############################################################################
# Edit these parameters to change the initial temperature, the pressure, and
# the phases in the mixture.

T = 300.0
P = 101325.0

# phases
gas = ct.Solution('gri30.yaml')
carbon = ct.Solution('graphite.yaml') #dummy
species = {S.name: S for S in ct.Species.list_from_file('gri30.yaml')}
complete_species = [species[S] for S in ('CH4','O2','N2','CO2','H2O')]
gas2 = ct.Solution(thermo='ideal-gas', species=complete_species)

# the phases that will be included in the calculation, and their initial moles
mix_phases = [(gas, 1.0), (carbon, 0.0)]
mix_phases2 = [(gas2, 1.0), (carbon, 0.0)]

# gaseous fuel species
fuel_species = 'CH4'

# equivalence ratio range
npoints = 50
phi = np.linspace(0.3, 3.5, npoints)

##############################################################################

mix = ct.Mixture(mix_phases)
mix2 = ct.Mixture(mix_phases2)

# create some arrays to hold the data
tad = np.zeros(npoints)
tad2 = np.zeros(npoints)
xeq = np.zeros((mix.n_species, npoints))
xeq2 = np.zeros((mix2.n_species, npoints))

for i in range(npoints):
    # set the gas state
    gas.set_equivalence_ratio(phi[i], fuel_species, 'O2:1.0, N2:3.76')
    gas2.set_equivalence_ratio(phi[i], fuel_species, 'O2:1.0, N2:3.76')

    # create a mixture of 1 mole of gas, and 0 moles of solid carbon.
    mix = ct.Mixture(mix_phases)
    mix.T = T
    mix.P = P
    mix2 = ct.Mixture(mix_phases2)
    mix2.T = T
    mix2.P = P

    # equilibrate the mixture adiabatically at constant P
    mix.equilibrate('HP', solver='gibbs', max_steps=1000)
    mix2.equilibrate('HP', solver='gibbs', max_steps=1000)


    tad[i] = mix.T
    tad2[i] = mix2.T
    print('At phi = {0:12.4g}, Tad = {1:12.4g}, Tad2 = {1:12.4g}'.format(phi[i], tad[i], tad2[i]))
    xeq[:, i] = mix.species_moles
    xeq2[:, i] = mix2.species_moles

# write output CSV file for importing into Excel
csv_file = 'adiabatic.csv'
with open(csv_file, 'w', newline='') as outfile:
    writer = csv.writer(outfile)
    writer.writerow(['phi', 'T (K)'] + mix.species_names)
    for i in range(npoints):
        writer.writerow([phi[i], tad[i]] + list(xeq[:, i]))
print('Output written to {0}'.format(csv_file))

csv_file2 = 'theoretical_adiabatic.csv'
with open(csv_file2, 'w', newline='') as outfile2:
    writer = csv.writer(outfile2)
    writer.writerow(['phi', 'T (K)'] + mix2.species_names)
    for i in range(npoints):
        writer.writerow([phi[i], tad2[i]] + list(xeq2[:, i]))
print('Output written to {0}'.format(csv_file2))

if '--plot' in sys.argv:
    import matplotlib.pyplot as plt
    plt.plot(phi, tad, label='adiabatic falme temperature')
    plt.plot(phi, tad2, label='theoretical', linestyle='--')  # tad2のプロットを追加
    plt.xlabel('Equivalence ratio')
    plt.ylabel('Adiabatic flame temperature [K]')
    plt.legend()  # 凡例を表示
    plt.show()

上記のプログラムを実行すると、理論断熱燃焼温度と平衡断熱燃焼温度が当量比ごとに計算されその結果がcsvで出力される。実行の際に、--plotを引数として与えると、グラフのプロットも同時に行われる。以下に、計算結果のプロットを示す。

破線のプロットが理論断熱燃焼温度であるが、熱解離を考慮した平衡断熱燃焼温度の方が温度は低くなっている。
熱解離が起こる場合は、エネルギーを吸収して安定化するために理論断熱燃焼温度よりも低下するためである。