Cantera：感度解析

※Canteraバージョン：3.1.0

Canteraで感度解析を行う。ここでは、メタンの燃焼における各素反応式について、温度に対する感度を計算する。

# conditions
gas = ct.Solution('gri30.yaml')
temp = 1500.0
pres = ct.one_atm
phi = 1.0

gas.TP = temp, pres
gas.set_equivalence_ratio(phi, 'CH4', 'O2:1.0, N2:3.76')
r = ct.IdealGasConstPressureReactor(gas, name='R')
sim = ct.ReactorNet([r])

GRI-Mech3.0の反応メカニズムを使う。条件は、温度1500K、圧力1atm、当量比1.0を設定する。ゼロ次元定圧反応器モデル（IdealGasConstPressureReactor）を使う。

# all reactions with respect to the temperature sensitivity
nreac = 325
for i in range(nreac):
    r.add_sensitivity_reaction(i)

反応メカニズムに含まれる325の素反応を感度解析にセットする。

# set the tolerances for the solution and for the sensitivity coefficients
sim.rtol = 1.0e-6
sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6

計算用のトレランスを設定する。

# time loop
smax = [0]*nreac
for t in np.arange(0, 2e-3, 5e-6):
    sim.advance(t)

    # sensitivities for each reaction
    for i in range(nreac):
        s = sim.sensitivity('temperature', i) # sensitivity of temperature to reaction
        if abs(s) > abs(smax[i]): # maximum sensitivity
            smax[i] = s

    print('time: {:12.7f}'.format(t))

時間積分の計算を行う。ここでは、0から0.002秒まで5e-6秒刻みで計算する。

2つ目のfor文で、各反応式の感度を計算している。sim.sensitivityで、i番目の反応式の温度に対する感度を取得する。smaxには時間を通して、絶対値が最大となるときの感度が格納される。

print('\n---------------- Sensitivity for all reactions --------------- \n') 
for i in range(nreac):
    print('{:3d} {:<35s} {:14.6f}'.format(i + 1, sim.sensitivity_parameter_name(i), smax[i]))

すべての反応式の感度を出力する。

---------------- Sensitivity for all reactions ---------------

---------------- Sensitivity for all reactions ---------------

  1 R: 2 O + M <=> O2 + M                     0.000161
  2 R: H + O + M <=> OH + M                   0.000609
  3 R: H2 + O <=> H + OH                      0.140857
  4 R: HO2 + O <=> O2 + OH                   -0.106612
  5 R: H2O2 + O <=> HO2 + OH                 -0.000479
...
321 R: C3H7 + H <=> C2H5 + CH3                0.000013
322 R: C3H7 + OH <=> C2H5 + CH2OH             0.000022
323 R: C3H7 + HO2 <=> C3H8 + O2               0.000014
324 R: C3H7 + HO2 => C2H5 + CH2O + OH         0.000018
325 R: C3H7 + CH3 <=> 2 C2H5                 -0.000113

# top 10
ntop = 10
smax_sort = sorted(smax, key=abs, reverse=True) 
smax_indx = [smax.index(x) for x in smax_sort]

rlabel = []
sp = []

print('\n---------------- Sensitivity for Top 10 reactions --------------- \n') 
for i in range(ntop):
    rlabel.append(sim.sensitivity_parameter_name(smax_indx[i]))
    sp.append(smax_sort[i])
    print('{:3d} {:5d} {:<35s} {:14.6f}'.format(i, smax_indx[i] + 1, rlabel[i], sp[i]))

ここでは、感度の大きな反応式を10個取り出して表示してみる。

まず、感度が大きい順にソートする。smax_sortは感度の絶対値でソートしたデータ、smax_indxはその時の式番号のデータが格納される。

次のfor文で、感度の大きな10個のデータを取得する。rlabelには反応式が入る。spには感度が入る。print文で結果を出力する。

---------------- Sensitivity for Top 10 reactions ---------------

  0   155 R: CH3 + O2 <=> CH3O + O                 16.836198
  1   158 R: 2 CH3 (+M) <=> C2H6 (+M)             -13.609158
  2    38 R: H + O2 <=> O + OH                     13.094800
  3    53 R: CH4 + H <=> CH3 + H2                 -10.150386
  4   156 R: CH3 + O2 <=> CH2O + OH                 6.775377
  5   119 R: CH3 + HO2 <=> CH3O + OH                6.266314
  6   159 R: 2 CH3 <=> C2H5 + H                     4.342360
  7   161 R: CH2O + CH3 <=> CH4 + HCO               4.026277
  8    32 R: CH2O + O2 <=> HCO + HO2                4.023569
  9   118 R: CH3 + HO2 <=> CH4 + O2                 3.434750

# plot
plt.title('Sensitivity for Temperature')
plt.barh(rlabel, sp)
plt.xlabel('Sensitivity')
plt.gca().invert_yaxis()
plt.gca().set_axisbelow(True)
plt.grid(axis='x')
plt.tight_layout()
plt.show()

最後にグラフにプロットする。$\mathrm{CH_3}+\mathrm{O_2} \leftrightarrow \mathrm{CH_3O}+\mathrm{O}$の反応が温度に対する感度が最も高い。次の$ 2 \mathrm{ CH_3}+ \mathrm{M} \leftrightarrow \mathrm{ C_2H_6}+ \mathrm{M}$は負の感度が高く、温度に対してマイナスに寄与することがわかる。

温度だけでなく化学種など各変数に対する感度も同様に計算できる。