SR8000/mpp による乱流燃焼の直接数値計算
東京工業大学
店橋

大学院理工学研究科

護*，名田譲，斎藤

機械宇宙システム専攻
敏彦，宮内敏雄* mtanahas@mes.titech.ac.jp
１．はじめに流体機械や実用燃焼器の設計・開発において，それらの流れ場（多くの場合乱流状態）を数値シュミレーションを用いて予測することは，設計・開発コストの削減に非常に重要である．乱流場の数値シュミレーションには幾つかの方法があるが，対象となる実機の高精度化や環境対策の強化などから，より精度の高いモデルの開発が必要となってきている．高精度乱流モデルや乱流燃焼モデルを開発するには，乱流中の微細構造や火炎の局所構造などを明らかにする必要があり，近年，直接数値計算の結果を用いたモデルの開発が多く行われている[1][2]．
直接数値計算 (Direct Numerical Simulation: DNS) は，モデルを一切導入せずに支配方程式を高精度離散化法により直接解く方法である．乱流や乱流燃焼などの支配方程式は，非線形偏微分方程式群から構成されており，DNS はこれらの非線形支配方程式の厳密解に限りなく近い解を求めることができる．しかし，乱流場を厳密に解くには，広範囲の時間・空間スケールにわたる速度変動を分解する必要があるため，乱流の DNS には非常に多くの格子点数と膨大な計算時間が必要とされる．
特に乱流燃焼の DNS の場合，乱流スケールだけでなく，火炎における化学種分布の変動も考慮に入れなければならない．例えば，図 1 は水素・空気層流予混合火炎における火炎近傍の化学種濃度分布を示している．化学種ごとにその分布の幅（空間スケール）は大きく異なっており，HO2 のように分布幅の狭い化学種も存在している．詳細化学反応機構[3]を用いた水素燃焼の DNS の場合，12 種類の化学種が考慮に入れられているため，それぞれの化学種の分布に対して十分な空間分解能を確保する必要がある．また，化学反応機構を構成する多くの化学反応の時間スケールは，乱流運動の時間スケールに比べて短い場合が多く，化学反応の時間変動を十分に分解するために非常に小さな時間積分幅（通常はナノ秒オーダー）を用いる必要がある．
このように DNS は膨大な計算時間を必要とするが， DNS により乱流中の微細構造や乱流燃焼における火炎構造などが明らかにされている[4][5]．今回の報告では，著者らにより行われてきた DNSの中で乱流予混合燃焼の DNS についてその計算方法と主な計算結果について報告する．
2.0 10-40.012

1.5 10-4XO

XHO20.016

0.008

1.0 10-40.004
0.000-0.75

5.0 10-5-0.50-0.25

x0.00(a)O のモル分率

0.25

0.0 100-0.75-0.50-0.25

x0.00(b)HO2 のモル分率

図 1 火炎近傍における各化学種の濃度分布0.25表 1 水素・空気乱流予混合火炎の DNS コードに用いた各サブルーチンにおける計算時間の割合
サブルーチン名

User time (%)

FFT51.9

Chemkin32.8

ノード間データ転送サブルーチン

8.35

質量，運動量，化学種及びエネルギ保存式

3.58

境界条件0.92

その他

2.4

図 2 計算領域

２．水素・空気乱流予混合火炎の直接数値計算
本計算では，質量，運動量，エネルギー及び化学種の保存方程式からなる 17 の非線形偏微分方程式群を基礎方程式として用いた．計算は図 2 に示すような三次元一様等方性乱流中を伝播する水素・空気乱流予混合火炎を対象としている．また，513×128×128 の格子点を用いており，総格子点数は 840 万点程度である．境界条件としては，x 方向の流入・流出境界条件には NSCBC を用い，y 方向と z 方向には周期境界条件を用いた．離散化は，x 方向に 4 次精度中心差分法を，y 方向と z方向にはスペクトル法を用いている．本計算の場合，y 方向と z 方向は周期境界条件を用いている
ため，速度ベクトルの各成分，温度及び各化学種濃度などがフーリエ級数に展開されている．化学反応機構としては，H2, O2, H2O, H, O, OH, HO2, H2O2, N, NO2, NO, N2 の 12 種類の化学種と 27 の素反応からなる詳細化学反応機構[3]を用いた．
表 1 に乱流燃焼の DNS コードの各サブルーチンにおける計算時間(user time)の割合を示している．
計算時間の 80%以上は FFT と，物性値，輸送係数及び反応速度を計算する Chemkin[6]サブルーチンによって占められている．この FFT は独自に開発した任意基数型の FFT である．2, 3, 4, 5, 6, 8, 9 及
び 12 基数からなり，N=2n13 n24 n35 n46 n58 n69 n712 n8 個の倍精度浮動小数点データの FFT が可能なものである．支配方程式はフーリエ空間において時間微分項のみを残した常微分方程式となるが，時間微分項は 3 次精度ルンゲクッタ法を用いて積分されており，支配方程式は初期値問題として解かれる．また，Chemkin サブルーチンはもとは Sandia National Laboratory で開発された汎用の燃焼計算コードであるが，著者らはこの Chemkin サブルーチンを水素燃焼専用サブルーチンとして，ベクトル化と並列計算が容易に行えるように変更を行った．水素燃焼の場合，Chemkin サブルーチンの計算時間に占める割合は 32.8%と FFT より小さくなっているが，メタンなどの炭化水素燃焼の場合，化学種及び素反応数が数百にも達するため[7]，その計算時間に占める割合は 70%程度まで増大する．
これらのプログラムは日立系のベクトル型スーパーコンピュータ用にコーディングされている．
日立系のベクトルプロセッサーはリスト・ベクトルに優れていたため，ベクトル化効率を上げるためにリスト・ベクトルが多用されている．バンク競合に関しては，配列の大きさを 2n+1 に設定することにより回避している．また，最も重要な FFT のプログラムは日立系のベクトル・プロセッサーでチューニングされており，インライン強制ベクトル化命令が含まれている．SR8000/mpp のUnburnt

Burnt(a) 一様等方性乱流中の微細渦(b) 燃焼場における微細渦と密度の等値面

図 3 乱流中の微細渦と密度の等値面

Unburnt

Unburnt

Burnt

Burnt

(a) O 濃度(b) HO2 濃度
図 4 各化学種の濃度

1node あたりの実効演算速度は約 3.4GFlops である．
３．水素・空気乱流予混合火炎の計算結果
図 3 は一様等方性乱流中の微細渦(a)と燃焼場における微細渦と密度の等値面(b)を示している．乱流場は図 3(a)に示すような三次元管状の微細な構造によって構成されており，この構造は乱流場の種類やレイノルズ数に依存しない普遍的な構造であることが明らかにされている[5]．また，乱流エネルギやエネルギ散逸のほとんどを担っており，乱流中のスカラー輸送に対しても重要な役割を果たしている．図 3(b)の乱流燃焼の場合，燃焼による温度上昇にともない粘性が増加するため，渦は高温側で急激に減衰する．また，火炎自体も乱流中の微細渦の影響を受けており，密度の等値面は大きく湾曲している．
図 4 は O と HO2 の濃度分布を示している．図中右下に乱流運動の影響を受けていない層流の場合の濃度分布を比較のために示す．図 1 に示したように化学種によりその分布の幅は大きく異なっているが，乱流運動によっても化学種の分布の幅は大きく変化している．このように乱流運動により化学種の分布の幅やその生成速度も大きく変化するため，DNS ではその変化も考慮に入れて時間・空間分解能を設定しなければならない．

case1
case2
case3

310u'rms /SL

210

distributed
reaction zones

well-stirredreactor

110

corrugated flamelets

1wrinkled laminar flame

1110210l/δL

310410

図 5 Combustion Diagram

図 5 は乱流スケール l と火炎厚さδL の比(l/δL)と乱流強度 u’rms と層流燃焼速度 SL の比(u’rms/SL)による火炎構造の分類を示している．図 5 に示すような火炎構造の分類は Combustion Diagram と呼ばれており，異なった領域に分類される火炎では化学種の分布などに大きな差違が現れると考えられている[8]．著者らは，これらの各領域の火炎構造を明確にするために図中の点で示した条件においてDNS を行っている．図 3 と図 4 に示した計算結果は図中の case 1 の計算結果であり，case 1 から case
3 のそれぞれの条件の DNS に必要な CPU time は SR8000/mpp の 1CPU 換算で 2.4, 9.7 及び 39 万時間である．ところが，実用燃焼器などの火炎構造は点線で示した領域に分類されており，乱流スケールは火炎厚さに対して非常に大きい．この領域の条件において DNS を行う場合，乱流スケールの増加のために計算領域を拡大する必要があるが，空間分解能自体は火炎のスケールに依存する．このため，計算領域中の格子点の数が増大し，点線の領域の計算時間は今まで行った DNS の計算時間の 102 から 106 倍に達する．著者らはより実機の条件に近い乱流スケールの大きな条件の DNS を計画中であるが，さらなる DNS コードの高速化と 128node 運用などの有効利用により計算時間の短縮を計りたいと考えている．
４．まとめ
ここでは乱流燃焼の DNS について，コードの最適化と最近行った計算の結果を簡単に説明した．
一連のコードの最適化に際して，東京大学情報基盤センターの多大なご協力を頂いた．
参考文献
[1]

A. Trouve, T. Poinsot: J. Fluid Mech., 278 (2000), p.1[2]

C. Meneveau, T. Poinsot: Combust. Flame, 86 (1991), p.311[3]

E. Gutheil, G. Balakrishnan and F. A. Williams, Reduced Kinetic Mechanism for Applications inCombustion Systems, Springer-Verlag (1993), p. 177[4] M. Tanahashi, M. Fujimura and T. Miyauchi: Proc. Combust. Inst., 28 (2000), p.529[5]

M. Tanahashi, S. Iwase, T. Miyauchi, Advances in Turbulence,VIII (2000), p.655[6] R. J. Kee, G. Dixon-Lewis, J. Warnatz, M. E. Coltrin and J. A. Miller, Sandia Report,
SAND89-8009B(1998)[7]

http://euler.Berkley.EDU/gri_mech/[8]

N. Peters, Lect. Notes Eng., 40 (1989), p.242