流体基礎および複雑現象の大規模シミュレーション
航空宇宙技術研究所
溝渕泰寛
１．水素噴流浮き上がり火炎の数値シミュレーション
浮き上がり火炎は部分予混合火炎の代表例であり、基礎的な研究という面から見てもまた実用的な研究の面から見ても興味深い研究対象である。浮き上がり火炎の研究における主たる研究課題は、燃焼場の詳細構造の解明および火炎が吹き飛ばされずに定在的な浮き上がり火炎が保たれるための安定化機構の解明であり現在まで主に拡散火炎の消炎、三重火炎構造という視点から多くの研究がなされている。しかしながら従来の研究は定常あるいは二次元解析が中心であり、三次元性非定常性が強い浮き上がり火炎の構造の解明は未だ十分には行われておらず、また安定化機構についても結論が出ていない。浮き上がり火炎の直接数値シミュレーション(DNS)を行い、そのデータを解析することにより、それらの現象の解明を試みた。
計算対象の火炎は空気中に噴き出された亜音速水素噴流の燃焼によって得られる浮き上がり火炎である。噴き出し口径Dは2mm、噴き出し速度は680m/sec、レイノルズ数は13600である。実験では浮き上がり高さ7Dの浮き上がり火炎が観察されている。
燃焼流計算におけるDNSでは流体の微細な運動とともに局所的な燃焼反応の構造もとらえる必要がある。水素の燃焼モデルとしては詳細な9化学種17反応モデルを用い、各輸送係数および各化学種のエンタルピー等は厳密に温度等の関数として与える。計算格子幅は噴流中心近くでは0.05mmとしている。これは実験で計測された着火点近くのKolmogorovスケールの約2.5倍、また当量比1の一次元層流水素・空気予混合火炎の反応帯幅の約1/10となっており、流体および燃焼反応の微細構造をとらえられることが可能である。火炎中心から離れたところでは計算格子は徐々に粗くなり、おおよそ一辺5cmの立方体の計算領域に用いた格子点数は約2300万点である。計算手法は有限体積法であり対流項は三次精度風上数値流束で評価し、時間積分は二次精度としている。時間積分幅は約6nsecであり燃焼反応をとらえるのに十分小さい積分幅である。
図1はある瞬間における温度1200Kの等値面図であり、数値シミュレーションにより浮き上がり火炎が再現されていることが分かる。浮き上がり高さは時間とともに変動するが時間平均的には約5.5Ｄであり実験結果と比較すると若干浮き上がり高さが短いが問題の複雑さを考慮すれば良好な一致である。これはDNS的なアプローチで浮き上がり火炎をとらえた世界で最初の例である。
図2は水素のモル分率80%の等値面であり表面の色は温度を表している。極めて三次元的な乱流噴流がとらえられている。この乱れは水素噴流の不安定性に起因している。未燃部分の等値面は非常に複雑な形状を有しており、そのフラクタル次元は多くの乱流噴流実験での計測値に近い2.3-2.4となっている。既燃部では温度上昇に伴う輸送係数の増加により乱れが減衰していることも分かる。
図3は水素消費速度104mol/m3/sec の等値面であり表面の色は発熱率を表している。この浮き上がり火炎には三つの火炎要素があることが分かる。すなわち火炎先端部、内側にある核状の火炎、外側にある島状の火炎である。
|
|
|
|
|
図1：温度等値面(1200K)
|
図2：水素モル分率等値面(80%)
|
図3：水素消費速度等値面(104mol/m3/sec)
|
|
火炎構造を調べるために無次元Flame Indexを導入した。Flame Index F.I.は F.I.＝∇YH2･∇YO2
と定義され、値が正の場合は局所的な火炎構造が予混合火炎、負の場合は拡散火炎となる。このF.I.を、正の場合は局所燃料混合分率に対応する一元層流予混合火炎におけるF.I.で、負の場合は一次元対向流拡散火炎の消炎時のF.I.で無次元化する。この無次元F.I.(N.F.I.)は予混合火炎に対しては局所的な燃焼速度の指標となり、また拡散火炎に対しては消炎の指標となる。図4は水素消費速度104mol/m3/secの等値面であるが表面の色は燃焼モードに対応している。モードの分類は局所的なF.I.の符号および燃料混合分率zの値を用いて行い、赤い部分はF.I.＞0,
z＞zst すなわち過濃予混合火炎、青い部分F.I.＞0, z＜zst すなわち希薄予混合火炎、緑の部分はF.I.＜0
すなわち拡散火炎を表している。zstは量論混合分率である。内側の核状の火炎は過濃予混合火炎、外側の島状の火炎は拡散火炎、火炎先端部は過濃/希薄予混合火炎と拡散火炎の複合体であることが分かる。
図5は二次元断面内のN.F.I.の分布である。温度が600Kに満たない領域においてはN.F.I.を0としている。N.F.I.が正値の等値線は0.4から10.0まで0.4刻みで実線で,
負値の等値線は-0.0002から-0.001まで0.0002刻みで破線で描いている。太線はz＝zstの等値線である。先端火炎は三重火炎的な構造を持っていることが分かる。動画1
を観察すると、外側の島状の拡散火炎はz＝zstの線に沿ってゆっくりと下流に流れて行き、内側の過濃予混合火炎は乱れが非常に強く不安定な火炎であることが分かる。一方で先端火炎は安定であり、この安定な先端火炎が全体の火炎の安定性をもたらしている。
先端火炎の希薄予混合火炎部(図5中A)ではN.F.I.がおおよそ1であり局所燃焼速度が層流燃焼速度であることを示している。一方で軸方向の流体の速度の時間平均はおおよそ層流燃焼速度と同程度となっている。すなわち先端火炎は希薄予混合火炎部での流体の速度と火炎の燃焼速度が釣り合うことによって安定化されていると考えられる。
この他に、外側の島状火炎は先端火炎付近の火炎の非定常な動きと局所的な消炎により生成されること、また島状火炎での燃焼は分子拡散に支配されていること、内側の予混合火炎では乱れの影響が火炎内部構造にまで及び従来LES等の計算に用いられてきたLaminar
Flameletが適応出来ない火炎となっていることも明らかになった。
今後はDNSで得られた膨大なデータの解析を進め、火炎構造および火炎安定機構の解明を更に進めるとともに、乱流燃焼のサブグリッドモデル構築に有効な知見を抽出していかなければならない。
|
|
 |
図4：火炎構造 赤:過濃予混合火炎、
青:希薄予混合火炎、緑:拡散火炎
|
図5：二次元断面内N.F.I.分布
動画表示
|
|
２．RANSによる燃焼流の数値シミュレーション
燃焼流は乱流に大きく影響を受ける。燃料と酸化剤の混合は促進され、また燃焼速度は増加する。燃焼の実用計算においてはこれら乱流の影響をモデル化し計算に組み込む必要がある。ここでは時間平均的な解を求めることを目的としてRANSによる乱流燃焼計算コードを開発した。乱流輸送係数については二方程式の乱流モデルから渦粘性係数を求め、乱流プラントル数、シュミット数一定の仮定のもと乱流拡散係数を求めた。また反応速度を平均温度と温度変動の関数とすることにより乱流が化学反応に及ぼす影響についても考慮した。その際温度変動は、温度変動およびその散逸率に関する二方程式の乱流モデルを用いて求めた。
図6はスクラムジェットエンジン内部のシミュレーションに適用した例であり水素の燃焼によって生成された水の分布を示している。NWT20PEを用い、計算格子点数は約150万点である。水素・空気の反応モデルとしては9化学種17反応モデルを用いている。図7は壁面圧力分布である。Φは燃料当量比で噴出し流量に対応するが、噴出し流量が大きく再循環域が保炎用のステップを越えてくる場合でも実験と良い一致を示している。
RANSで乱流燃焼のシミュレーションを行う場合、乱流輸送係数および反応速度をどのように決めるかが課題である。乱流プラントル数、シュミット数の値をいくつにすべきか分かっていないし、それらが一定であるという仮定の密度変化が大きな場での妥当性は極めて低い。また温度変動についての方程式についても議論の余地は限りなく大きい。今後DNSや実験のデータを基にこれらの値について調べていく必要がある。
|
|
|
| 図6：スクラムジェットエンジン内のH2O分布 |
図7：壁面圧力分布 |
|
３．衝撃波反射形態に関する数値シミュレーション
衝撃波の反射形態には正常（Regular）反射とマッハ（Mach）反射がある。入射角が小さい場合は正常反射となるが角度を大きくするとマッハ反射に遷移する。また逆に角度を小さくしていくとマッハ反射から正常反射に遷移する。二次元定常を仮定した理論によれば前者はdetachment
条件から決まる入射角、後者はvon Neumann 条件から決まる入射角で遷移する。すなわち反射形態の遷移にはヒステリシスが存在し、どちらの形態も可能なdual-solution
領域と呼ばれる条件が存在する。実験における正常反射からマッハ反射への遷移はdetachment条件に至る前にこのdual-solution 領域で起こる。その原因としては入射角変化による流れ場の非定常性、三次元性、擾乱の影響が考えられている。ここではピッチング運動をする三次元楔周りに生じる衝撃波の反射形態の遷移を観察することにより、非定常性および三次元性が遷移に及ぼす影響について調べた。
計算は非粘性計算であり用いた計算格子点数は約100万点である。図8はマッハ数4の一様流内におかれた二つの楔間に生じたマッハ反射の様子を示している。左図、右上図および右下図はそれぞれ楔間の水平面内、楔直後の面内および対称面内のマッハ数分布を示している。下側の楔角は22.6deg.に固定されており、上側の楔角を26deg.から30deg.まで徐々に大きくすることによって正常反射からマッハ反射へと遷移した後の結果である。楔は長さc=120mm、幅w=330mmであり楔間の距離は213.8mmである。遷移の様子は動画2
で観察することが出来る。図9は楔のアスペクト比と遷移時の上側の楔角の関係を示している。これらの結果は、非定常性は遷移に影響しないこと、また流れ場の三次元性は正常反射からマッハ反射への遷移を遅らすことを示している。
この研究では非定常性及び三次元性は正常反射からマッハ反射への遷移がdual-solution領域で起こる原因にはならないことが示された。今後は風洞内の微小粒子が衝撃波を通過する際の遷移への影響などを調べていく必要がある。 |
 |
 |
図8：マッハ反射時のマッハ数分布
動画表示
|
図9：楔のアスペクト比と遷移角の関係
|
|
参考文献
- 水素噴流浮き上がり火炎の構造について、日本燃焼学会誌 第44巻
128号 pp.73-79, 2002.
- A Numerical Analysis on Structure of Turbulent Hydrogen Jet Lifted FLame,
29th Proceedings of Combustion Institute, to appear, 2002.
- Numerical Analysis of Mach/Regular Reflection
FLow Around Three-dimensional Wedges, Proceedings of 22nd International Symposium
on Shock Waves, pp.1273-1278, 1999.
|
