SENAC Vol. 52, No. （3 2019. 7）― 1 ―［共同研究成果］

Flamelet approach に基づくガス燃料を対象とした燃焼シミュレーション
赤尾津翔大：東北大学大学院工学研究科
松下 洋介：東北大学大学院工学研究科
青木 秀之：東北大学大学院工学研究科
Weeratunge Malalasekera：Loughborough University

本研究では，ガス燃料を対象とした乱流燃焼シミュレーションのターゲット・フレームとして最も実績のある乱流拡散火炎の 1 つである Sandia Flame D を対象に Flamelet approach の 1 つである Flame/Progress-Variable approach を用いて乱流燃焼の Large Eddy Simulation を実施し，測定結果と比較することで本乱流燃焼シミュレーションの妥当性を検討した．その結果，混合分率，温度およびラジカルを含む化学種の時間平均値とその分散について解析結果は測定結果とほぼ完全に一致した．今後，同手法を拡張することで NO やすすなどの大気汚染物質の予測精度について検討する予定である．

1.

緒言調整パラメータを必要としない詳細化学反応機構を用いた燃焼シミュレーションは理想的である．しかしながら，最も単純な炭化水素であるメタンの燃焼ですら，例えばこのメタンの燃焼を対象に対して最も実績のある詳細化学反応機構の 1 つである GRI-Mech 3.0 [1]には 53 の化学種と353 の化学反応が含まれる．そのため，詳細化学反応機構を用いて三次元の燃焼シミュレーショ
ンを実施すると，解くべき化学種の保存式の数が多く，計算負荷が膨大となってしまう．また，一般に詳細化学反応機構を用いた反応解析の時間刻みは流体解析の時間刻みより数オーダー小さい．そのため，時間スケールの小さい化学反応の時間刻みに合わせて三次元の燃焼シミュレーシ
ョンを実施すると，やはり計算負荷が膨大となってしまう．さらに，詳細化学反応機構に限らず，有限の反応速度を用いて燃焼シミュレーションを実施する場合，化学種の保存式の生成項にアレニウス型の反応速度が現れる．Reynolds-Averaged Navier-Stokes (RANS)あるいは Large EddySimulation (LES)を用いて乱流燃焼の計算を行う場合，化学種の保存式に時間平均あるいは空間平均を施すが，この化学種の保存式の生成項に現れるアレニウス型の反応速度に時間平均あるいは空間平均を施す方法論が確立されていない．これはアレニウス型の反応速度式に温度に関する指数関数が含まれ，非線形性が強いためである．LES を用いた乱流燃焼では化学反応速度の項にス
ケール相似則を適用する Scale Similarity Filtered Reaction Rate (SSFRR) [2]モデルを用いた計算が報告されているものの，その妥当性は十分に確認されているとは言えず，また，解くべき化学種の保存式の数が多くなることと時間刻みが小さくなることにより計算負荷が膨大となってしまう問題は解決されない．
Flamelet approach [3,4]では，拡散燃焼において形成される火炎が一次元的な火炎片(Flamelet)の集合で表現できるとする．また，化学反応の時間スケールが流体の時間スケールと比較して極めて小さいことを利用し，化学反応と流体をスケール分離する．すなわち，流体の時間スケールでは化学反応が定常状態に達しているとし，あらかじめ 2，3 のパラメータに対してデータベースを作成し，三次元の燃焼シミュレーションではこのパラメータを求め，データベースを参照することで興味ある残りのパラメータすべてを決定する．そのため，Flamelet approach を用いることで，比較的低い計算負荷で間接的ではあるものの詳細化学反応機構を考慮した三次元の燃焼シミュレーションを実施することが可能となる．

― 2 ―

SENAC Vol. 52, No. （3 2019. 7）

本研究では，ガス燃料を対象とした乱流燃焼シミュレーションのターゲット・フレームとして最も実績のある乱流拡散火炎の 1 つである Sandia Flame D [5–7]を対象に Flamelet approach の 1 つである Flamelet/Progress-Variable approach [8]を用いて乱流燃焼の Large Eddy Simulation を実施し，測定結果と比較することで本乱流燃焼シミュレーションの妥当性を検討する．

2.

Flamelet approach に基づく燃焼シミュレーション

Flamelet approach では，三次元の燃焼シミュレーションに先立ち，あらかじめ一次元の対向流拡散火炎や予混合火炎などの単純な火炎を対象に詳細化学反応機構を用いて一次元の燃焼シミュレーションを実施し，混合分率 Z や反応の進行を表す Progress Variable (C)と呼ばれるパラメータに対して，温度，化学種の質量分率，密度や粘度などの物性値および C の正味の生成速度をFlamelet table と呼ばれるデータベースに保存する．三次元の燃焼シミュレーションにおいてこの混合分率 Z や C の保存式を解き，あらかじめ作成した Flamelet table を参照することで，燃焼場における温度，化学種の質量分率，密度や粘度などの物性値および正味の C の生成速度などの変数を決定する．なお，Flamelet table から参照する変数のうち，三次元の燃焼シミュレーションに影響を及ぼす変数は密度や粘度などの物性値および C の正味の生成速度である．

2.1 解析対象と解析条件
本研究では，ガス燃料を対象とした乱流燃焼シミュレーションのターゲット・フレームとして最も実績のある Sandia Flame D を対象とする．図 1 に解析対象と解析領域の概念図を示す．Sandia
Flame D は，開放空間において(a)内側のノズルから温度 294 K，体積分率でメタン 0.25 と空気 0.75の混合気を断面平均流速 49.6 m/s で，(b)その外側からメタンを当量比 0.77 で燃焼させた温度 1880
K のパイロット火炎の燃焼ガスを断面平均流速 11.4 m/s で，(c)周囲から co-flow として温度 291 Kの空気を断面平均流速 0.9 m/s で供給することで形成される部分予混合火炎である．なお，燃料ノズルを基準とすると，レイノルズ数は Re = 22,400 である．混合分率は酸化剤側の境界面で 0，燃料側の境界面で 1 となる変数であり，境界条件として，(a)では Z = 1，C = 0，(b)では Z = 0.27，C= 0.21，(c)では Z = 0, C = 0 となる．なお，計算格子には六面体格子のみを用い，総計算格子数を128 万あるいは 300 万分割とした．

図 1 解析対象と解析領域の概念図 [5–7]Flamelet approach に基づくガス燃料を対象とした燃焼シミュレーション― 3 ―2.2 Sandia Flame D に対する Flamelet table の作成
一次元対向流拡散火炎を対象に，詳細化学反応機構にメタンの燃焼に対して多くの実績があるGRI-Mech 3.0 [1]を採用し，FlameMaster V3.3.10 [9]を用いて一次元の燃焼シミュレーションを実施した．熱損失ゼロおよびいわゆる unity Lewis number を仮定し，混合分率空間において種々のスカラー消散率について一次元対向流拡散火炎の定常解を求め，混合分率 Z とスカラー消散率に対して温度，化学種の質量分率，密度や粘度などの物性値および正味の Progress Variable C の生成速度などの変数で表される二次元のデータベースを保存する．その後，混合分率 Z とスカラー消散率に対する変数で表される二次元のデータベースを混合分率と Progress Variable C に対する変数で表される二次元のデータベースに変換する．さらに，次節で述べるとおり，Large Eddy
Simulation を用いて三次元の燃焼シミュレーションを実施するため，混合分率 Z の確率密度関数に関数，Progress Variable C の確率密度関数に関数を仮定し，Favre 平均を施した混合分率，その分散と Progress Variable C に対する変数で表される三次元のデータベースを作成する．なお，
Flamelet table の作成方法の詳細は著者らの総説[10]で解説したとおりである．
2.3 三次元の燃焼シミュレーションの解析方法
Sandia Flame D を対象とした三次元燃焼シミュレーションの基礎式は低マッハ数近似に Favre
フィルタを施した連続の式 (1)，運動量保存式 (2)，混合分率の保存式 (3)と Progress Variable (PV)の保存式 (4)である．
��̅����̅ ��� � � ��� ���������̅�
���� ���� 2 �����̅ �� � ���̅ ��� ��� � ���̅ ����� ��� �� � ���� ��� ������ �� � ������ ��� 3��������� ���������� � ����� � ����� �� ���̅ ��� ���̅ ����� � ���̅ ���������������� �� ������� � ����� � �� ��� �� ���̅ �� � ���̅ �� ��� � ���̅ ����� � �������
����������� �� ���(1)
(2)
(3)
(4)

ここで，式 (2)–(4)の右辺第 3 項をそれぞれ式 (5)–(7)で近似する．
������ ���� 2 ����̅ ����� ��� � � �� �
�
� �
�
(5)
� �� � ���� ��� 3 �� ���� ����� � ����� � � ���̅ ���(6)����� ���� ����� � �� ��� � � ���̅ ���(7)����� ���式 (5)–(7)中のt は乱流粘性係数であり，式 (8)に示す Smagorinsky モデル[11]を用いて推算した．
(8)�� � ���� ∆�� �
なお，CS は Smagorinsky モデルの唯一のモデル定数であり，本研究では 0.1 で一定とした．また，ScZ,t と ScC,t,は混合分率と Progress Variable の乱流シュミット数であり，本研究では 0.7 で一定とした．混合分率の分散は式 (9)を用いて推算した．
����
�
���
�� � �� ∆ � �
���

(9)

以上より，求めた混合分率，混合分率の分散と Progress Variable を用いて前述した Flamelet tableを参照し，密度と Progress Variable の正味の生成速度を求めるとともに温度や化学種の質量分率を決定する．
非構造格子の有限体積法に基づき，運動量保存式の対流項は三次風上差分法を 5%ブレンドした二次中心差分法，拡散項二次中心差分法を用いてそれぞれ離散化し，時間進行法には二次のAdams-Bashforth 法を用いた．計算時間を短縮するため時間刻みを可変とし，クーラン数が 0.4 と― 4 ―

SENAC Vol. 52, No. （3 2019. 7）なるように決定した．なお，乱流燃焼場が定常的な挙動を示す際の時間刻みは約 1 s である．圧力の解法には Simplified Marker And Cell (SMAC)法[12]を適用し，SMAC 法における圧力補正値に関するポアソン方程式の解法には Algebraic Multigid Solver [13]を用い，各タイム・ステップにおいて収束解が求められるまで反復計算を実施した．混合分率と Progress Variable の保存式の対流項は流束制限関数に min-mod 関数[14]を適用した Total Variation Diminishing を，拡散項は二次中心差分法を用いて離散化し，時間進行法には一次の陰解法を用い，多項式前処理付き安定化双共役勾配法[15]を用い，各タイム・ステップにおいて収束解が求められるまで反復計算を実施した．
計算開始から時刻 0.2 s まで乱流燃焼場を発達させ，時刻 0.2–0.3 s のデータを平均化することで各変数の時間平均値を求めた．METIS-5.1.0 [16]を用いて解析領域を分割し，Message Passing
Interface (MPI)を用いて領域分割に基づく並列計算を実施した．また，通信の多い領域間から可能な限り同時に通信する独自のアルゴリズム[17]を採用し，通信に要する時間を最小化し，オーバヘッドを最小化することで並列化効率を向上した．なお，東北大学サイバーサイエンスセンター所有の並列コンピュータ LX 406Re-2 を用い，4 ノード 48 コアを用いた並列計算を実施した．

3.

結果と考察

3.1 混合分率，温度および主要な化学種の質量分率分布の瞬時値
図 2 に t = 0.2 s において，z = 0 の断面における混合分率，温度および主要な化学種の質量分率の分布を示す．混合分率は酸化剤側の境界面で 0，燃料側の境界面で 1 となる変数であり，ノズルから供給された燃料の噴流が減衰するとともに周囲のガスと混合することが示されている．温度については，上流において燃料ノズルの周囲から供給された高温のパイロット火炎の燃焼ガスが高温を示し，また，ノズルから供給された燃料がこの燃焼ガスと混合することで加熱されるとともにパイロット火炎の周囲から供給された co-flow と混合することで燃焼が進行し，中流から下流に向かって高温領域が形成されている．また，燃焼の中間生成物である CO は比較的燃料過濃な領域で，OH は比較的燃料希薄な領域で生成している．さらに，燃焼の最終生成物である CO2が中流から下流に向かって生成している．これらの分布は定性的な比較ではあるものの既往の研究[18]と概ね一致している．次節において解析結果を測定結果と詳細に比較する．

図 2 速度，混合分率，温度および主要な化学種の質量分率の瞬時値の分布

Flamelet approach に基づくガス燃料を対象とした燃焼シミュレーション― 5 ―3.2 中心軸上における各変数の時間平均値の比較
図 3 に中心軸上における混合分率，その Root-Mean-Square (RMS)，温度および化学種(CH4, CO,CO2, OH, H2O, H2)の質量分率の時間平均値を測定結果とともに示す．解析結果および測定結果はいずれも上流から下流に向かって 1 から徐々に減衰している．これはノズルから供給された燃料が拡散しながら周囲のガスと混合するためである．混合分率の時間平均値の解析結果は測定結果とほぼ完全に一致していることから，本解析はノズルから供給された燃料の噴流が乱流燃焼場において減衰し，周囲のガスと混合する過程を良好に表現していることを意味する．
混合分率の RMS については解析結果および測定結果はいずれも上流から下流に向かって上昇し，その後減少している．解析結果は測定結果を概ね再現しているものの，上流において解析結果は測定結果より若干低い値を示している．これは，解析では流入境界面において空間的な速度分布を考慮しているものの，その時間的な変動を考慮していないためであると考えられる．理想的には本解析の流入境界面より上流の燃料を供給するノズルを解析対象に含めることで，本解析の流入境界面における時間的な速度の変動も再現すべきである．しかしながら，LES において壁せん断が支配となるノズル内の乱流流れの解析を実施するためにはさらに多くの計算格子を必要とし，計算負荷が大きくなってしまうため現実的ではない．ノズルの下流に位置するある断面における速度などの物理量をノズルの上流に位置する流入境界面に複製するいわゆる Recycle
boundary を用いることで解析すべきパイプの長さを短くすることはできるものの，それでもなお多くの計算格子を必要とする．一方，Pierce and Moin [19]は事前にノズルのみを対象とした解析を実施し，実際の解析の流入境界面に相当する断面における速度などの空間的な分布とその時間的な変動を含む物理量をすべて保存し，実際の解析でこの保存したデータを流入境界条件として用いることを提案している．この方法では，流入境界面において妥当な空間的な分布と時間的な変動が求められることが期待されるが，保存すべきデータ量が膨大となるだけでなく，異なる解析対象間に生じるメッシュの不整合を処理しなければならない．また，パイプ内の流動が下流の影響を受けないことも確認しなければならないため，用いるためには解決すべき課題も多い．実用的な方法として，流入境界面において人工的に時間的な変動を生成する手法[20]が数多く提案されている．しかしながら，本解析対象と同じ Sandia Flame D を対象に流入境界面に人工的に時間的な変動を生成する手法を適用した解析結果が報告されているものの，解析対象ごとに調整すべきパラメータが存在し，汎用性に欠ける．そのため，LES において妥当な流入境界条件を与えるのは極めて困難であり，今後の研究課題であると言える．
温度の時間平均値については解析結果および測定結果はいずれも上流から下流に向かって上昇し，最高温度を示した後に緩やかに減少している．これは，ノズルから供給された燃料がその周囲から供給されたパイロット火炎の高温の燃焼ガスと混合するとともにパイロット火炎のさらに周囲の co-flow と混合し，燃焼するためである．解析結果は測定結果とほぼ完全に一致している
ものの，最高温度を示した後に解析結果は測定結果と比較して若干過大評価している．これは，解析結果と測定結果に差が生じる領域では，燃焼ガスの主成分であり吸収性ガスである CO2 とH2O の質量分率が高く，ふく射伝熱により熱損失を考慮していないためであると考える．実際，
本解析対象と同じ Sandia Flame D を対象とした研究においても同様の傾向が報告されており[18]，今回用いた燃焼モデルの範疇ではむしろ妥当な解析結果が求められていると考える．しかしながら，NO など温度に敏感な化学種の質量分率を正確に予測するためには温度をさらに正確に予測する必要があり，ふく射伝熱を連成し，Enthalpy defect [21]と呼ばれる熱損失を考慮可能な Flameletapproach を実施する必要があると考える．

4.

結言

本研究では，ガス燃料を対象とした乱流燃焼シミュレーションのターゲット・フレームとして最も実績のある乱流拡散火炎の 1 つである Sandia Flame D を対象に Flamelet approach の 1 つであ― 6 ―

SENAC Vol. 52, No. （3 2019. 7）

る Flame/Progress-Variable approach を用いて乱流燃焼の Large Eddy Simulation を実施し，測定結果と比較することで本乱流燃焼シミュレーションの妥当性を検討した．その結果，混合分率，温度およびラジカルを含む化学種の時間平均値とその分散について解析結果は測定結果とほぼ完全に一致した．今後，同手法を拡張することで NO やすすなどの大気汚染物質の予測精度について検討する予定である．

(a) Z とその RMS(b) 温度(c) CH4(d) CO(e) CO2(f) OH(g) H2(h) H2O
図 3 中心軸上における混合分率，その RMS，温度および化学種の質量分率 [5–7]謝辞
本研究は，東北大学サイバーサイエンスセンターの並列コンピュータ LX 406Re-2 を利用することで実現することができた．また，本研究の一部は特別研究員奨励費(18J11135)，JSPS 科研費(JP18K03964)，東燃ゼネラル石油研究奨励・奨学財団および東北大学若手リーダー研究者海外派
遣プログラムの助成を受けたものであり，ここに謝意を表す．
参考文献[1] Smith, G. P. et al., GRI-Mech 3.0, http://www.me.berkeley.edu/gri_mech/[2] DesJardin, P. E. and Frankel, S. H., Large eddy simulation of a nonpremixed reacting jet: application
and assessment of subgrid-scale combustion model, Phys. Fluids A, 10, 2298–2314 (1998)[3] Peters, N., Laminar diffusion flamelet models in non-premixed turbulent combustion, Prog. Energy
Combust. Sci., 10(3), 319–339 (1984)

Flamelet approach に基づくガス燃料を対象とした燃焼シミュレーション― 7 ―[4] Peters, N., Laminar flamelet concepts in turbulent combustion, Symp. (Int.) on Combust., 21, 1231–1250 (1998)[5] Barlow, R. S. and Frank, J. H., Effects of turbulence on species mass fractions in methane/air jetflames, Proc. Combust. Inst. 27, 1087–1095 (1998)[6] Barlow, R. S. et al., Piloted methane/air jet flames: Scalar structure and transport effects, Combust.
Flame, 143, 433–449 (2005)[7] Schneider, Ch. et al., Flow field measurements of stable and locally extinguishinghydrocarbon-fuelled jet flames, Combust. Flame, 135, 185–190 (2003)[8] Pierce, C. and Moin, P., Progress-variable approach for large-eddy simulation of non-premixedturbulent combustion, J. Fluid Mech., 504 73–97 (2004)[9] Pitsch, H., and Bollig, M., FlameMaster; A Computer Code for Homogeneous Combustion andOne-Dimensional Laminar Flame Calculations, Institut für Technische Mechanik, RWTH, Aachen.
[10] 松下洋介ら, Flamlelet Model に基づく乱流燃焼シミュレーション, 金属, 85(11), 915–921(2015)[11] Smagorinsky, J., General circulation experiments with the primitive equations, Mon. Weather Rev.,
91(3), 99–164 (1963)[12] Harlow, F. H. and Welch, J. E., Numerical calculation of time-dependent viscous incompressible flowof fluid with a free surface, Phys. Fluids, 8 2182–2189 (1965)[13] 藤井昭宏ら, 領域分割による並列 AMG アルゴリズム, 情報処理学会論文誌コンピューティングシステム（ACS）, 44, SIG06(ACS1), 9–17 (2003)[14] Roe, P. L., Characteristic-based schemes for the Euler equations, Annu. Rev. Fluid Mech., 18, 337–365 (1986)[15] Van der Vorst, H. A., Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for theSolution of Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comput., 13(2), 631–644 (1992)[16] Karypis, G. and Kumar, V., A Fast and Highly Quality Multilevel Scheme for Partitioning IrregularGraphs, SIAM J. Sci. Comput., 20(1), 359–392, (1999)[17] Matsushita, Y. et al., Efficient communication strategy in parallel computation based on domain
partitioning, J. Chem. Eng. Jpn., 51(1), 79–82 (2018)[18] Raman, V. and Pitsch, H, A consistent LES/filtered-density function formulation for the simulation ofturbulent flames with detailed chemistry, Proc. Combust. Inst., 31(2), 1711–1719 (2007)[19] Pierce, C. and Moin, P., Method for generating equilibrium swirling inflow conditions, AIAA J., 36(7),1325–1327 (1998)[20] Klein, M., A digital filter based generation of inflow data for spatially developing direct numerical orlarge eddy simulations, J. Comput. Phys., 186(2), 652–665 (2003)[21] Hossain, M. et al., Modelling of a Bluff-Body Nonpremixed Flame using a Coupled
Radiation/Flamelet Combustion Model, Flow, Turbul. Combust., 67(3), 217–234 (2001)