平成２６年度先端的計算科学研究プロジェクト成果報告書

LES/Flamelet approach を用いた産業用実機燃焼器における乱流燃焼場の大規模解析
高橋裕介1，佐々木天詩2，村瀬景信3，大島伸行4（北海道大学），岩井保憲5（株式会社東芝）

1. はじめに	 	  産業用ガスタービンで用いられる実機燃焼器は，発電や輸送機関の核であり，その性能向上（とくに熱効率）は産業界において重要な案件である．一方，環境問題よりNO! 排出には厳しい規制が課せられている．熱効率向上とNO! 低減化は相反する設計パラメータであり，装置の最適化問題となる．近年，産業用燃焼器の分野において希薄予混合燃焼が注目を浴びている．これは，高熱効率を維持した上で，NO! 排出の大幅な低減が期待されており，産業界・学術分野などの関連分野において研究開発が進行している．しかしながら，希薄予混合燃焼は，逆火や燃焼振動などの問題が生じやすい．そのため，燃焼器内部におけるこのような複雑な現象の挙動を明らかにする需要は現在高まっている[1][2]．	  実機燃焼器内部は強い乱流が燃焼を維持することに特徴がある．乱流現象を明らかにする上で，数値解析は有効な技術の１つであることは，この数十年間の研究開発の積み重ねから示されてきた．また，近年の計算機環境の発展に加え，燃焼解析技術も急速に進歩して お り ， 乱 流 燃 焼 場 を 精 度 よ く 解 析 す る 道 が 拓 け て き た と 言 え る ． 例 え ば Flamelet
approach と呼ばれる燃焼モデルは，火炎面を取り扱うことで燃焼場を表現し，数百・数千種類もの化学反応を直接組み込むモデル（詳細反応計算等）に比べて計算コストが小さい．
そのため非定常乱流場を精度良く捉えることができる Large-eddy simulation (LES)と合わせることで，上述のような実機燃焼器における乱流燃焼場を現実的な計算コストで調べることが可能になりつつある[3]．	  LES を 利 用 す る 場 合 ， 乱 流 燃 焼 場 の 解 像 度 を 高 め る こ と は 重 要 で あ る ． と く に
LES/Flame approach において，精度を高めることは計算格子を細分化し微細な渦構造，及びそれによって維持される燃焼構造を直接解像することに帰結する．このようなことから，本研究目的を実機燃焼器における強い乱流燃焼場について，LES/Flamelet approach1

大学院工学研究院，ytakahashi@eng.hokudai.ac.jp
工学部
3 大学院工学院
4 大学院工学研究院
5 電力・社会システム技術開発センター2平成２６年度先端的計算科学研究プロジェクト成果報告書

および大規模解析を用いて調査することに設定した．	  これまでの当研究グループの研究では，FrontFlow/Red と呼ばれる乱流燃焼場解析コードを用いて計算要素数 4500 万(770 万節点数)の計算対象に対して解析が行われた．ここでは自動細密化ルーチンによって細密化した計算格子を用いて燃焼器内部の乱流燃焼場の解析を実施した．およそ 4500 万要素数だったベース格子は，細密化ルーチンを利用することで 3.6 億要素に細密化される．

2. 基礎方程式	 	  本研究では Large-Eddy Simulation および 2-scalar flamelet アプローチによる数値解析手法を用いた．同手法の基礎方程式となる空間平均操作を施した NS 方程式系，およびflamelet 概念に基づいた G 方程式モデル，保存スカラー方程式モデルについて述べる．なお，数値解析には文部科学省 IT プログラム「革新的基盤ソフトウェアの研究開発」成果をもとに北海道大学が改良開発した有限体積法ベースの流体解析コード「FrontFlow/redver3.1」[4]を使用した．

2.1.乱流モデル：Large	 Eddy	 Simulation	 	  従来の Reynolds 平均量（時間平均量）の支配方程式を解析する手法（RANS モデル：
Reynolds Averaged Navier-Stokes model）と異なり，Large Eddy Simulation（LES）では空間フィルタのみを施した空間的な平均量に対する支配方程式を解く．これにより，乱流場を特徴付ける大規模乱流変動に対して時系列を追った直接的な解析が可能となる．一方，空間フィルタリングによる平均操作によって消失する小スケールの変動成分の影響については，SGS（Sub-Grid Scale）乱流モデルを用いて表現する．LES では時間平均操作によるモデリングを行わないため，乱流モデルが現象の時間的な変化に忠実な解析が可能であり，非定常乱流現象の解析に適している．	  低 Mach 数近似を適用した圧縮性流体の支配方程式に空間フィルタを施した LES の基礎
方程式は𝜕𝜌 𝜕𝜌𝑢!
+=0
𝜕𝑡
𝜕𝑥!
𝜕𝜌𝑢! 𝜕𝜌𝑢! 𝑢!
𝜕𝑢! 𝜕𝑢!
𝜕𝑝𝜕+=−+𝜇+− 𝜏𝑖𝑗!"!
𝜕𝑡
𝜕𝑥!
𝜕𝑥! 𝜕𝑥!
𝜕𝑥! 𝜕𝑥!となる．本解析で対象としているのは多くの化学種とそれらによる化学反応を含んだ系であるため局所密度が空間分布を持つ．したがって支配方程式に空間フィルタ操作を施すと

平成２６年度先端的計算科学研究プロジェクト成果報告書多くの密度相関項が現れてしまう．これを回避するために本解析では Favre 平均（荷重密度平均）を用いる．上式において物理量φ に対して空間フィルタを施した値を𝜙，密度荷重
空間フィルタ（Favre filter）を施した値を𝜙 = 𝜌𝜙 𝜌という記号でそれぞれ表す．また，燃焼反応による発熱で大きな温度変化を伴うため，粘性係数の温度依存性を考慮する必要がある．本解析では次式で示される式を用いて粘性係数を定めた．
𝜇=𝐴𝑇 !
𝑃𝑟𝑇!

ここで，𝐴 = 2.58×10!!   kg (m. s)，参照温度𝑇! = 298  K，𝛾 = 0.7とした．また𝑃𝑟はプラントル数である．	  SGS乱流モデルには，SGS 運動エネルギの生成と散逸がほぼ等しいとする局所平衡状態を仮定して得られる標準Smagorinskyモデル[5]を用いた．同モデルは，以下の式のように表される．
!"!
𝜏!"= −2𝜇 !"! 𝑆!" ,𝜇 !"! = 𝜌 𝐶! ∆ ! 2𝑆!" 𝑆!" ,𝑆!" =1 𝜕𝑢! 𝜕𝑢!
+2 𝜕𝑥! 𝜕𝑥!

𝐶! は Smagorinsky 定数，∆は空間フィルタ幅（解析格子幅）である．本解析では Smagorinsky
定数𝐶! の値は 0.15 とした．	 

2.2.

燃焼モデル	 :	 2-scalar flamelet approach	 	  予混合火炎モデル，非予混合火炎モデルのそれぞれに対する flamelet aprroach について述べる．さらに，本解析で仮定する部分予混合火炎に対して，2 つのモデルをカップリングした 2-scalar flamelet approach について述べる．

2.2.1. 予混合火炎モデル	  予混合火炎面の挙動に関しては，火炎面位置をその等値面(G=G0)として表現する指標ス
カラーG に対する方程式（G 方程式）[6][7][8]を用いて表現する．ここでは Liu らによって提案された火炎厚さを考慮した修正 G 方程式[9]を用いる．Favre フィルタを施した G 方程式を以下に示す．
𝜕𝜌𝐺 𝜕𝜌𝑢! 𝐺𝜕𝐺𝜕 𝜇 !"!
𝜆 𝜕𝐺𝜕+=++2𝜌! 𝑆! 𝐺𝜕𝑡
𝜕𝑥!
𝜕𝑥!
𝜕𝑥! 𝜎!
𝐶! 𝜕𝑥! 𝜕𝑥!
LES における乱流燃焼速度𝑆! は，空間フィルタリングにより粗視化された GS 火炎面の面積に対して，SGS の変形により実際の火炎面積は増大する影響を考慮するものである．そのモデルとして，本解析では Daniele [6]による乱流燃焼速度モデルを用いた．これは 1 気平成２６年度先端的計算科学研究プロジェクト成果報告書

圧から 20 気圧までの幅広い圧力領域に対して，乱流燃焼速度が検証されており，本解析対象のような高圧での燃焼に対しても適用することが可能である．Daniele ら[10]の定式を以下に示す．
𝑆!
𝑢!
= 𝑚𝑎𝑥 𝑎
𝑆!
𝑆!

!.!"

𝐿 ! !!.!" 𝑝 !.!" 𝑇! !!.!", 𝛼!"# ,𝜎!
𝑇!
𝑝!

ここで，𝐿!
𝑝 !.!!
=𝑏, 𝑎 = 337.45, 𝑏 = 8.3, 𝑝! = 0.1  MPa,   𝑇! = 1  K, 𝛼!"# = 173.4exp   1.3424𝑝
𝜎!
𝑝!

2.2.2. 非予混合火炎モデル	  非予混合燃焼（拡散燃焼）の取り扱いには，化学種およびエネルギの輸送を保存スカラーによって表現する保存スカラーアプローチを用いる．すべての化学種の拡散係数が等しい，Lewis 数 Le = Sc/Pr = 1，低 Mach 数近似の仮定の下，壁面での伝熱および輻射の影響を無視すると，元素の質量分率𝑍の保存式とエンタルピℎの保存式が相似形となり，他の境界条件についても相似であれば 1 つのスカラー量𝜉によって両者を表現することが出来る．
このときスカラー𝜉は化学反応による生成消滅を伴わない保存スカラーである．火炎の構造はこの保存スカラー𝜉のみによって表現され，その形状は化学量論混合比の等値面𝜉!" で表される．その輸送方程式に Favre フィルタを施すと以下のように記述できる．
𝜕𝜌𝜉 𝜕𝜌𝑢! 𝜉𝜕 𝜇 𝜕𝜉
𝜕+=−𝜌 𝑢! 𝜉 − 𝑢! 𝜉𝜕𝑡
𝜕𝑥!
𝜕𝑥! 𝑆! 𝜕𝑥!
𝜕𝑥!
保存スカラー𝜉 は酸化剤と燃料の混合を表現しているため以降「混合分率」と呼ぶ．上式右辺第２項の SGS スカラー流束項にはＧ方程式と同様に勾配拡散モデルを用いた．
𝜌 𝑢! 𝜉 − 𝑢! 𝜉 = −𝜇 !"! 𝜕𝜉𝑆!!"! 𝜕𝑥!

ここで𝑆!!"! は乱流 Schmidt 数である．混合分率  𝜉と密度・温度などの各スカラー変数との関係は，後述の flamelet data を参照することにより与えられる．

2.2.3. 部分予混合火炎モデル	 	  前項までで述べた予混合燃焼モデルと非予混合燃焼モデル，すなわち G 方程式と𝜉方程式による解析モデルを組み合わせることにより，両燃焼形態を含んだ系に対する解析モデルが提案されている．このモデルでは，スカラーG を未燃，既燃の判別パラメータとして用平成２６年度先端的計算科学研究プロジェクト成果報告書

い，一方，未燃，既燃の両状態での化学組成は混合分率によって酸化剤と燃料の混合状態を決定することにより定める．これらを結合することによって，局所温度，局所密度，局所燃焼速度を決定するモデルが得られる．すなわち，局所における化学種𝛼の質量分率は𝑌! = 1 − 𝜃 𝐺

𝑌!,! + 𝜃 𝐺 𝑌!,!となる．本解析では，未燃，既燃の両状態の結合を最も単純なスカラーG による線形結合として仮定し，
𝛼 𝐺 =𝐺とした．また，本解析では解析の高速化のため，温度についても，未燃，既燃の両状態での値をスカラーG によって線形結合して算出する仮定をおき，密度についてはスカラーGによる調和平均によって算出されるものとした．
𝑇 = 𝐺𝑇! 𝜉 + 1 − 𝐺 𝑇! 𝜉 ,𝜌=𝜌! 𝜉 𝜌! 𝜉𝐺𝜌! 𝜉 + 1 − 𝐺 𝜌! 𝜉2.3.火炎データ	 
flamelet approach では化学反応を解く代わりに混合分率𝜉 の輸送方程式，および，火炎面レベルセット方程式（Ｇ方程式）を解くことによって燃焼状態を表現することを述べた．
ここでは flamelet approach によって燃焼を表現する際に必要な混合分率に対する温度や密度などの熱物性，および，層流火炎速度といった変数を関連づける火炎データ（flameletdata）について説明する．これらの火炎データを求めるために，予め化学反応計算ソフト
ウェア CHEMKIN ver4.1.1[11][12]を用いて計算を行った．この計算に用いた物性データ，化 学 反 応 デ ー タ は ア メ リ カ ガ ス 協 会 (GRI) に よ っ て 公 開 さ れ て い る デ ー タ ベ ー ス
GRI-MECH3.0 [13]を使用した．これはメタンからプロパンまでの炭化水素系燃料を用いた
詳細化学反応解析を行うための化学反応データセットであり，53 化学種，325 反応式が含まれている．
混合分率に対する温度及び密度のデータについては，等エントロピ，等圧条件における無次元化学平衡計算を行い，化学平衡状態の密度，温度を求めた．混合気組成については，混合分率を区間 0～1 において 0.001 刻みで変化させた初期条件を用いた．そして，得られた温度，密度のデータを元に gnuplot による重み付け最小二乗法によるフィッティング機
能を利用し多項式近似を行い，混合分率の多項式関数として近似した．
また，予混合気の層流燃焼速度については，予混合気中の自由伝播火炎に対する一次元詳細化学反応解析を行い，その伝播速度と混合分率との関係を得た．ただし数値解析上の平成２６年度先端的計算科学研究プロジェクト成果報告書

簡便のため，データテーブルは層流燃焼速度𝑆! と未燃密度𝜌! との積，すなわち質量流束の

図 1.1. 層流燃焼速度

図 1.3(a) 温度(ケース A)	 図 1.2. 密度（比体積）

図 1.3(b) 温度(ケース B)

形で構築した．これは，予混合火炎面の単位面積辺りの燃料消費速度と解釈する事もできる．解析では，この質量流束と混合分率の関係について，温度，密度と同様の手順で混合分率の多項式として近似した．本解析で用いたテーブルを図 1.1-1.3 に示す．

2.4.NO! 予測モデル	 	  Flamelet approach でNO! を評価する場合，NO生成速度は流れ場特性時間に比べても遅く，Fast chemistry が適用できないと考えられる．ただしNO生成量は他の化学種に比べて小さく，燃焼場に大きく影響は与えないものと仮定し，パッシブスカラ量として以下のNO輸送方程式（化学種保存則）を解いて評価する．
𝜕𝜌𝑌!" 𝜕𝜌𝑌!" 𝑢!
𝜕 𝜇 𝜕𝑌!"+=+ 𝜔!" 	 𝜕𝑡
𝜕𝑥!
𝜕𝑥! 𝑆! 𝜕𝑥!
なお，Fuel NO は無視し，Prompt NO と Thermal NO の生成量のみを計算する．	  それぞれの生成速度を計算するにあたって，N2, O, O2, Fuel(CH4)の質量分率𝑌を評価する必要がある．これはフレイムレットテーブルと，局所的な無次元温度(G)と混合分率(𝜉)より求める．

平成２６年度先端的計算科学研究プロジェクト成果報告書

2.4.1. Thermal NO	 
本研究においては，以下の３つの反応モデル（拡大 Zeldvich 機構）を考慮することによって Thermal NO を評価する：
N! + O ↔ NO + N
N + O! ↔ NO + O
N + OH ↔ NO + H
生成速度は次の式で求める，
𝜔!"#$%&'_!" = 2𝑘!" 𝑌!! 𝑌! ,      𝑤ℎ𝑒𝑟𝑒  𝑘!" = 1.84×10! exp −38370 𝑇

2.4.2. Prompt NO	 
Prompt NO に関しては，生成速度を次の式[14]で求める：
𝜔!"#$!%_!" = 𝑘!" 𝑌!!! 𝑌!! 𝑌!"#$ ,      𝑤ℎ𝑒𝑟𝑒  𝑘!" = 1.2×10! exp −30200 𝑇
ここで，酸素濃度指標𝑏は以下で与える，
𝑏 = 1,𝑖𝑓  𝑌!! ≤ 4.1×10!!

𝑏 = −3.95 − 0.9ln𝑌!! ,𝑖𝑓  4.1×10!! < 𝑌!! ≤ 1.11×10!!

𝑏 = −0.35 − 0.1ln𝑌!! ,𝑖𝑓  1.11×10!! < 𝑌!! ≤ 3.00×10!!

𝑏 = 0,𝑖𝑓  3.00×10!! < 𝑌!!

平成２６年度先端的計算科学研究プロジェクト成果報告書

3. 解析モデル	 
3.1.解析手法	 
解析条件を表 3.1 にまとめる．	 
表 3.1	  解析条件	 

流体の性質	 

低マッハ数近似非圧縮性流体	 

乱流モデル	 

標準スマゴリンスキーLES	 (Cs=0.15)	 

燃料の種類	 

都市ガス 13A	 

メッシュの形状	 

四面体(774 万節点・4548 万要素)	 

時間積分法	 

クランク・ニコルソン法	 

空間補間	 

速度場(対流項)：2 次精度中心差分(95	 %),	  	 1 次精度風上差分(5	 %)	 
速度場(粘性項)：2 次精度中心差分	 
スカラー場：2 次精度中心差分	 

時間刻み	 

1.0e-6[s]:	 非燃焼計算	 
5.0e-7[s]:	 燃焼計算(base	 grid）	 
2.5e-7[s]:	 燃焼計算(refine	 grid)	 

計算ステップ数	 

200,000 ステップ	 	 
3.2.境界条件・解析格子	 	  対象とした希薄予混合燃焼器の解析領域・外観を図 3.1 に，流入部名称を図 3.2 に示す．
境界は流入部・流出部・壁面部で構成される．さらに流入部はメイン予混合バーナー

図 3.1. 解析対象の外観・境界名

図 3.2. 各流入部の名称

平成２６年度先端的計算科学研究プロジェクト成果報告書(Fm-Am)，パイロット予混合バーナー(Fpp-App)，パイロット拡散バーナー（燃料流入部:Fpd・空気流入部:Apd），冷却空気スリットで構成される．流入部において燃料流量と空気(酸化剤)流量およびそれぞれの温度が与えられる．流出部において静圧が与えられており，また逆流なしの流出条件を課している．壁面では速度場にたいしてすべりなし，圧力に対して壁面法線方向の勾配 0，壁面温度は断熱，すなわちスカラー方程式に対して壁面法線方向の勾配 0 である．	 

9	  図 3.3.	  計算格子(base	 grid）の概観	 
計算格子を図 3.3 に示している．各バーナー部が燃焼室内部に接続する領域の着火領域や，燃焼室内の保炎領域における燃焼・乱流の挙動を精度よく捉えるために，そのような領域において計算格子を密にしている．相対的に，出口ダクト領域における計算格子は粗い．また，実機設計における実用性を重視してすべての格子にテトラ格子を採用している．
これによって比較的容易な格子生成が可能になっている．節点数・要素数はそれぞれ 774万節点，4548 万要素となっている．	 

平成２６年度先端的計算科学研究プロジェクト成果報告書

3.3.計算条件	 	  解析は燃料・空気流入量や内部圧力を変更した２つのケースに対して行われる．表 3.2にそれぞれのケースの条件を示している．
表 3.2. 計算条件
CaseA

B

Inlet fuel/air ratio, kg/kg

0.01920.0281

Fuel temperature, K

300.0

300.0

Air temperature, K

606.9

659.1

Total pressure at exit,
MPa

0.9782

1.171

なお，今回の Refine 格子解析を行うにあたって，九州大学の先端的計算科学研究プロジェクトによる計算リソース(FUJITSU PRIMERGY CX-400, 1024 並列, ３ヶ月間専有利用)を利用した．

3.4.ガス温度計測位置	 	  実験において出口ダクト部においてガス温度およびNO! が計測されている．ガス温度の測定には熱電対が用いられている．図に出口ダクト部におけるその計測位置を表している．
解析においても，毎ステップそれらの座標のガス温度，NO! 質量分率を取得している．

図 3.4.ガス温度計測位置	 

平成２６年度先端的計算科学研究プロジェクト成果報告書

4. Refiner 機能による細密格子を用いた乱流燃焼場解析	 	  非燃焼解析による速度場，温度，燃料混合分率の計算結果(流入開始後 0.01 秒)を初期条件として，燃焼解析を継続して行った．非燃焼場における乱流場の挙動は昨年度報告としたものと定性的に同じであるのでここでは割愛する．なお，着火過程の模擬として，まず非燃焼解析結果(0-10,000 step)を初期条件として，メイン予混合バーナーおよびパイロット予混合バーナーの流入部より既燃(G=1)ガスを流入させる．燃焼器内部に燃焼場が十分に形成された後(10,000-30,000 step)に，流入部既燃ガス(G=1)を未燃ガス(G=0)に切り替えて燃焼解析をスタートさせた．ただし，切り替え後しばらくはバーナーダクト部に既燃ガスが残ることから，各流入ダクト内部が未燃ガスで満たされるまで解析を継続する(30,000-50,000 step)．したがって統計量の取得はその後からである．この方法は多少計算時間を要するが，保炎される運転条件であればほぼ安定的に保炎領域が形成されるため，燃焼場を解析する上で有用である．

4.1.乱流燃焼場	 
FrontFlowRed のパッケージに同梱される計算格子 Refiner 機能を利用して，細密格子を作成し，乱流燃焼場解析を行った．ただし，この Refiner 機能は現在開発中のもので，公開版
FFR v3.1 には含まれていない．計算格子の Refine は１本のエッジを２つに分割することで行われる．したがって，１つのセルは８つ分割されることとなり，計算格子数（セル数，ノード数）は 2^3=8 倍となる．本解析においては，Base grid (45 million elements)で乱流燃焼場を十分に発達させた後に，Refiner 機能で細密にした格子(360 million elements)で計算を開始する．安定的な計算を進める上で，クーラン数一定を維持するならば，Refineした際に時間刻み幅をおよそ半分にする必要がある．ここでは時間刻み幅は 2.5e-7[s]とした（Base grid では 5.0e-7 [s]）．ケース B に関しては Step 200,000 まで Base grid で解析し，Step 200,000-279,000 の間で Refine grid で解析を行った．ケース A に関しては，同様に Step 200,000 まで Base grid で解析し，それから Step226000 まで Refine grid で解析した．そのため，ケース A に関しては，十分な議論を行う上では統計量が少ないと判断し，今回は省略した．

平成２６年度先端的計算科学研究プロジェクト成果報告書(a) Refine grid(b) Base grid

図 4.1. レベルセット関数 G 分布（瞬時値）(a) Refine grid(b) Base grid

図 4.2. 混合分率分布（瞬時値）	  図 4.1 と 4.2 に瞬時レベルセット関数 G の分布と瞬時混合分率の比較を示している．(a)が Refine grid の結果，(b)が Base grid の結果で，運転条件はいずれもケース B である．特徴的な点として，格子を細密にすることによって，メイン予混合バーナー近傍のレベルセット関数 G が増加している（バーナー部２，３段目の着火が捉えられた）ことが示されており，この条件において，格子細分化は G 方程式に対して大きなインパクトがあると言える．一方で混合分率の分布については大きな差異は見られなかった．なお，出口ダクト部でのレベルセット関数 G は Refine grid，Base grid いずれもほぼ G=1 であることを示しており，格子を細分化しても既燃状態から変化することはない．

平成２６年度先端的計算科学研究プロジェクト成果報告書	  図 4.3 と 4.4 にはそれぞれ Refine grid(a)，Base grid(b)における瞬時ガス温度の分布と瞬時速度の絶対値の分布を示している．また，図 4.5 と 4.6 には同様にそれぞれ瞬時層流燃焼速度の分布と瞬時乱流燃焼速度の絶対値の分布を示している．これらの諸量に関しては瞬時値については格子の細密化に伴う大きな差異が見られない．

(a) Refine grid(b) Base grid

図 4.3. ガス温度分布（瞬時値）(a) Refine grid(b) Base grid

図 4.4. 速度絶対値の分布（瞬時値）

平成２６年度先端的計算科学研究プロジェクト成果報告書(a) Refine grid(b) Base grid

図 4.5. 層流燃焼速度の分布（瞬時値）(a) Refine grid(b) Base grid

図 4.6. 乱流燃焼速度の分布（瞬時値）	  図 4.7 は出口ダクト部におけるガス温度の Base grid 結果，Refine grid 結果および，実験結果の比較を示している． Refine grid と Base grid の結果で，混合分率分布に大きな差異が現れなかったことと，出口ダクト部においていずれの格子の場合でもほとんど既燃（G=1）であったことから，ガス温度も大きく変わってはおらず，実験結果を若干上回る値となった．ただ，ダクトコア部よりも，ダクト側壁部での反応進行度が進んでおり（G が高い値をとる），Refine grid の方が Base grid よりも大きめなガス温度をとることが示されている．ガス温度を大きく予測する方が実験結果よりも離れてしまうが，細密格子にしたことで新たに誤差を生じるということではなく，これは本来壁面部での非断熱効果が現れているものと考えるのが合理的である．本解析モデルにおいて非断熱効果を表現するモデルが導入されていないので，その導入は将来的に必要な事項になる可能性も示唆する結果平成２６年度先端的計算科学研究プロジェクト成果報告書となった．

図 4.7 出口ダクト部の温度比較（解析：時間平均値）

4.2.燃焼器内部の NO 挙動	 	  Refine grid, Base grid それぞれでの Thermal NO と Prompt NO の結果を図 4.8, 4.9 に示している．運転条件は Case B であり，いずれも時間平均値である．なお，Refine gridに対しては 79000step 分の統計，Base grid については 30000step 分の統計を取得した結果である．Refine grid では Thermal NO の質量分率を大きく予測する傾向にあることが示された．それは保炎領域において顕著である．これはパイロット部近傍の局所的なガス温度が格子を Refine することで増加した結果，Thermal NO の生成量が増加したためと考えられる．一方で，Prompt NO は若干ではあるが Refine grid の方が減少した．	  表 4.1 は出口ダクト部における NO 質量分率 (Thermal NO+Prompt NO) の解析結果(Refine grid, Base grid)，および実験結果との比較である．ただし CFD の結果は実験値で規格化されている．また，解析結果は時間平均値となっている．運転条件は Case B である．
Refine grid の結果は，Base grid の結果と比べて，10%弱ほど大きく予測することが示されている．いずれの解析の場合も実験結果を過大に予測している．ただし，オーダ評価は十分に可能なものと思われる．

平成２６年度先端的計算科学研究プロジェクト成果報告書(a) Refine grid(b) Base grid

図 4.8. Thermal NO の分布（時間平均値）(a) Refine grid(b) Base grid

図 4.9. Prompt NO の分布（時間平均値）

表 4.1. 出口ダクト部における NO の比較（実験値との比）

Average

Experiment,

CFD Refine grid,

CFD Base grid,

ppm/ppm

ppm/ppm

ppm/ppm

1.000

1.708

1.589

平成２６年度先端的計算科学研究プロジェクト成果報告書

4.3.並列計算の効果	 	  本解析モデルの並列計算パフォーマンステストとして，今回利用した九州大学情報基盤研究開発センター高性能演算サーバ(FUJITSU PRIMERGY CX-400: tatara)の他に，同センター高性能アプリケーションサーバ(HA8000-tc/HT210: Hakozaki)および北海道大学情報基盤センタースーパーコンピュータシステム SR16000 を用いた．計算格子として Basegrid（4548 万要素数・774 万節点数）と，それをベースに自動細分化ルーチンを利用して生成した Refine grid（3.639 億要素数・6194 万節点数）の 2 つである．	  図 4.10 に並列数を変更した際の計算実時間を示している．この計算実時間は基準ステップ(1000 step)および基準格子節点数(100 万節点数)で規格化したものである．同じ格子を用いた場合，tatara と hakozaki, SR16000 と hakozaki で計算実行時間に大きな違いは現れなていない．本解析モデルでは流体の運動量保存方程式，燃焼場を表現するためのスカラー方程式の他に，ポアソン方程式型である圧力補正方程式を毎ステップ厳密に収束させる必要がある．実際に，計算時間の大部分が圧力補正方程式の収束に割かれている．その収束手法にあたるマトリクスソルバーには，共役勾配法（ICCG 法）が用いられている．この手法は収束性が良好であるものの，計算格子数の増加に伴って収束までの反復回数が増加することが特徴である．そのため基準節点数で規格化しているにもかかわらず，格子サイズの大きなケースほど計算実時間は増加する．圧力補正方程式の収束に必要な反復回数はBase grid で約 1000 回，Refine grid で約 2200 回となった．

図 4.7 計算実行時間（基準ステップ・基準格子節点数で規格化）

平成２６年度先端的計算科学研究プロジェクト成果報告書

5. まとめ	 	  ガスタービン燃焼器における乱流燃焼場の数値解析を行った．乱流モデルには標準スマゴリンスキーモデル，燃焼モデルには 2-scalar flamelet approach を用いた．Refiner 機能による細密格子（360 million elements）を用いて，ケース B における燃焼器内部の乱流燃焼場解析を実行した．格子を細密にすることで，ケース B におけるメイン予混合バーナー近傍の関数 G がより大きな値を持つようになった．すなわちその領域における燃焼特性の変化が見られた．ただし，出口ダクト部におけるガス温度分布は Base grid/Refine grid で大きな違いは現れない．また，同様に出口ダクト部での NO 質量分率についても格子間での大きな違いは現れなかった．

参考文献	 
[1]

Y. Iwai, F. Maeda, M. Yamada, and M. Ota. “Study on Low Emission Combustorfor Gas Turbine (Optimization of Combustion Characteristics for MinimumEmission)”. Journal of Gas Turbine Society of Japan, 37(2):97–103, 2009.

[2]

F. Maeda, Y. Iwai, and M. Yamada. “Combustion Stability of Dry Low NOxCombustor with Axially Fuel-Staged Combustion”. Journal of Gas Turbine Societyof Japan, 32(2):102–109, 2004.

[3]

Y. Takahashi, N. Oshima, and Y. Iwai. “Large- Eddy Simulation of Transient
Behavior in a Combustion Field for Gas-Turbine Engine”. In 11th World Congress
on Computational Mechanics, 5th European Conference on Computational
Mechanics, 6th European Conference on Computational Fluid Dynamics (WCCM
XI-ECCM V-ECFD VI), Barcelona, Spain, 20-25 July 2014.

[4]

N.

Oshima,“FrontFlow/reddownloadpage”,

https://www.eng.hokudai.ac.jp/labo/fluid/ download/download.htm.
[5]

J. Smagorinsky. “General Circulation Experiments with the Primitive Equations, I.
The Basic Experiment”. Monthly Weather Review, 91(3):99–164, 1963.

[6]

N. Peters. “Local Quenching due to Flame Stretch and Non-Premixed Turbulent
Combustion”. Combustion Science and Technology, 30(1-6):1–17, 1983.

[7]

N. Peters. “Laminar Diffusion Flamelet Models in Non-Premixed Turbulent
Combustion”. Progress in Energy and Combustion Science, 10(3):319–339, 1984.

[8]

Williams,F.A., Combustion Theory(2nd ed.),Addison-Wesley, 1985.

[9]

Y. Liu and N. Oshima. “A new level set approach for a premixed flame based on a平成２６年度先端的計算科学研究プロジェクト成果報告書

new concept of flame speed”. Journal of thermal science and Technology, 6(1):140–153, 2011.
[10] S. Daniele and P. Jansohn. “Correlations for turbulent flame speed of different
syngas mixtures at high pressure and temperature”. In Proceedings of ASME
Turbo Expo 2012, GT2012-69611, Copenhagen, Denmark, June 11-15 2012[11] R.J. Kee, F.M. Rupley, E. Meeks, and J.A. Miller. “CHEMKIN-III: A FortranChemical Kinetics Package for the Analysis of Gas-phase Chemical and PlasmaKinetics”. UC-405 SAND96-8216, 1996.
[12] R.J.

Kee.

“CHEMKIN,Reaction

Design,Inc,SanDiego”.

http://www.reactiondesign.com/ products/chemkin/.
[13] M. Frenklach, T. Bowman, and G. Smith. http: //www.me.berkeley.edu/gri-mech/.
[14] G. G. De Soete, “Overall Reaction Rates of NO and N2 Formation from Fuel
Nitrogen”, 15th Symp,(int.)Comb. pp.1093-1102, 1975.