先端研究施設共用イノベーション創出事業【産業戦略利用】『みんなのスパコン』TSUBAMEによるペタスケールへの飛翔
利用報告書 平成 20年度新規利用拡大採択課題 i08ne

3次元圧縮性流体コードによる高精度混相流シミュレーションコード超並列化手法の開発
Feasibility study of parallel processing in three dimensional multi-phase compressible fluid flow
桑原 諸兄
Moroe Kuwahara
株式会社計算流体力学研究所 混相流解析グループ
Institute of Computational Fluid Dynamics
www.icfd.co.jp
環境に配慮したエンジンの開発では、燃料の噴霧特性について混相流の高精度解析が必要であるが、その場合、流体の圧縮性は無視できない。このため、混相状態にある噴霧流の流動特性を、圧縮性を再現しながら高精度で解析し、かつ超並列化により高速に計算できる 3次元圧縮性混相流計算コードの開発を目指した。2次元、3次元の従来の直交格子系での圧縮性流体シリアル計算コードに並列処理を施したところ、TSUBAME でのベンチマークでは、プロセッサ数に対してほぼ線形の並列化効率が得られた。
さらに、圧縮性流体計算コードの応用可能性を広げるため、一般座標系に同コードが適用可能できるかを検証した。計算対象は、比較可能な計算結果が豊富な NACA0012 型の翼周りの流れである。同コードを用いることにより、高速かつ高精度な計算が一般座標系でも可能なことが示された。

Studying characteristics of fuel flows of environmental friendly engines require high resolution simulations of multiphase flows.
Fluid Compressibility is significant in those simulations. We attempt to develop computer programs, which also consider
compressibility to simulate three-dimensional multiphase flows for super parallel computers with high resolution rapidly.
Current two- and three-dimensional serial computer programs for simulating a compressible fluid in orthogonal grids are
modified to run on multiple processors. Benchmark tests in the TSUBAME system confirm nearly linear speed-up dependingon the number of processors. Then, these programs were further developed in generalized coordinates to enhance the accuracyof the simulation of flow. In order to validate our numerical method, we calculate the flow around a two-dimensional
NACA0012 wing since there are abundant data. The results show excellent agreement with existing data.
Keyword: Compressible Fluid, Upwind Differential, Multi-Directional Differential, General Coordinate
・背景と目的新しい化学製品、薬品等の開発において、より高純度の生成物創出の妨げとなっている気泡混入問題の解析や、より環境に配慮したディーゼルエンジン等の新しい環境負荷低減を実現するシステム開発には、混相流の流動特性の精密な解析が必要である。
本プロジェクトは、これらの問題に対応するため、従来の手法よりもはるかに高精度な計算コードの開発を目的とする。高圧下での混相流の解析には、圧縮性流体の計算が必須だが、高精度の実現には、時間刻みを非常に小さくする必要があり、同スケールの非圧縮性流体を扱う場合と比較し、必要な計算量は数桁程度増加する。

しかし、音速が無限大となる非圧縮性流体の計算では、ある地点の物理量の計算にはシステム全体の情報が必要なのに対し、圧縮性流体は音速が有限なため、ある地点での物理量は周囲の地点の情報にのみ依存する。
そのため、ノード間の通信量をなるべく小さくしたい超並列計算では、並列化効率をあげるために、圧縮性を考慮した方がむしろ有利と考えることもできる。
流体計算において、任意形状の物体の取り扱いは常に問題となる。既存の圧縮性流体計算コードでは、一様格子の格子幅をより細かくし、その物体を格子空間でマスクデータとして与えるアプローチをとってきた。しかし、さらに高い精度が要求される系では、物体の境界に適合する一般座標系を形成し、その座標系で計算を行う。格子幅を細かくするアプローチは、計算コードの超並列化により解決の目途がたっているが、同じ計算コス
トで、特定形状の物体周りの流れをより高精度に解くには、一般座標系に計算コードを拡張した方が有利な場合もあると考えられる。そこで本報告では、既存コードを一般座標系でも計算できるように拡張した結果を示す。

・概要
高精度混相流シミュレーション実現のための課題は様々ある。新たな相の導入により、異なる相の境界の扱いが、混相流を解く上で一番重要な問題であると考えられる。
本報告では、そのような混相流特有な問題に立ち入る前段階として、より基礎部分のシミュレーション過程を構築した。具体的には、3 次元の直交格子系の単相流体計算プログラムの並列化とそのベンチマーク、及び 2次元の一般座標系単相流体並列計算プログラムの開発である。我々は、本プロジェクト実行以前に圧縮性、直交格子系での 2 次元と 3 次元の単相流体シリアル計算プログラム、非圧縮性、及び一般座標系のメッシュ生成プログラム、シリアル計算プログラムを完成させており、新たに開発する圧縮性の計算コード開発にそれらを応用できる[1]。本来、目的としている計算対象はエンジン内の燃料の噴霧を想定しているが(図 1)、本報告では、全システムの連携を図る前段階のシミュレーション結果を示す。

・結果および考察
既存の直交格子系の単相流体シリアル計算プログラムを並列化し、その並列化効率を検証した。この計算プ
ログラムの特徴は、圧縮性流体計算における高精度維持に一般的な、TVD と 3 次精度の風上差分（MUSCLE内挿）に加え、非圧縮性流体計算で実績のある多方向差分を組み合わせた。これにより、圧縮性のある流体の亜音速から超音速にまたがる非定常問題を高い精度で解けるようになった。
2 次元と 3 次元のシリアル計算プログラムに対し、MPIにより並列化を施した。簡略化のため、全計算領域において格子点の最も多い軸（X 軸）に沿ってのみ領域分割し、そのメッシュ上のデータを各並列プロセスに与える。プロセス間で通信されるべきデータは、プロセスの個別領域の境界にある 1 点、あるいは端から 2 点の部分のみであることから、通信負荷を最小限に低減できる。
3 次元の直交格子によるダクト内流れ解析について、カラーシェードよる圧力分布の例を図 2 に示す。

図 2 3 次元による計算結果例

図1 エンジン内の噴霧流の例（出典:Kawahara、 et.al. 2007、 The 19th internal
combustion engine symposium)
一般座標系単相流体並列プログラムでは、従来航空分野で多くの実験、数値シミュレーション（特に定常等エンタルピー流れの計算）が行なわれてきた翼周りの流れを解析することで、新たな計算コードの計算結果の妥当性を検証する。特に高精度シミュレーションの必要性が高い遷音速スケールの領域に着目した。

圧縮性単相流れについて並列化したコードを用い、256×128×128 の格子サイズでベンチマークを行った。この格子サイズは、実用的な計算上最低限のサイズと考えられるが、メモリの利用を分散できる並列計算コードにより、初めて計算可能になったものである。図 3、図 4に、並列化による計算速度向上の度合いを示す。図 4では、プロセス数の増加による計算速度がほぼ線形に向上しており、並列化により、非常に効率よく計算速度が向上している。ただし、1 ノードに多くのプロセスを割り当てるよりも、多くのノードにわたって少数のプロセスを割り当てた方が実行速度が速い。同じプロセス数であっても、ノードの利用方法により、実行速度に 1.5 倍程度の差が生じている。しかし、後者の利用方法は、現状の利用時間制限枠内において多くの利用時間を消費してしまい、計算資源の無駄遣いとも言える。実際の計算におけるプロセスの割り振りには、問題に応じた対処が必要である。

800
700speed(arb. unit)

600500400300200
1000081624324048566472プロセッサ数

図 3 プロセッサ数とノード数を変えた 1,000step の計算にかかる時間。
横軸は（プロセッサ数）/（ノード数）を示す。

1000時間（ｓ）

900800

図 5 一般座標系における NACA0012 型の翼周りの遷音速流れ。
圧力（赤：高、青：低）のコンターを示している。

700
600
500
400
300
200100・まとめ、今後の課題

P64(32noP6de4(s)16node64s
P(8n )
6 4 odesP(4n )od32e
P(3 2 s)n
3 2 od es
P(4n )
3 2 odesP(2n )
od
16es
P()
16n
1 6 od es)
P(2n1 6 odes
P(1n )
odes)
8p(8no8 P de s)(1nodes)

0

プロセッサ数

図 4 プロセッサ数を変えた場合の速度向上
次に、2 次元並列計算において、一般座標系でも計算できるように計算コードを拡張した。メッシュの生成は、従来の非圧縮性の計算コード用のものを適用している。
ここでは NACA0012 型の翼の周りに O 型の 360（翼の周囲）×60 のメッシュを生成し、マッハ数 0.8、迎え角1.25 度の流れ場を計算した。よく知られているように、この流れは翼の上部で衝撃波を形成する。図 5 は計算結果の圧力のコンターを示す。時間刻みは 0.0001 秒で14.8 秒後の結果である。衝撃波とともに、翼の表面で目立った渦の生成が捉えられた。衝撃波の生成される位置は、先行研究と比較して妥当な結果であった。

本報告の範囲では、混相流そのものの計算については報告しないが、一般に圧縮性流れを高解像度で計算できることを示した。3 次元では、並列化により実用的な計算が可能になったとともに、TSUBAME によるベンチマークでは、当初の期待通り、プロセッサ数に対してほぼ線形に計算速度が向上した。また計算コードは一般座標系での計算にも適用でき、翼周りの流れを高解像度で計算できた。
今後も継続して、より良い計算コードの開発を進める。
[1]参考文献
1. Mi Young LEE, Tetuya KAWAMURA, Kunio
KUWAHARA and Junichi OOIDA, Simulation ofUnsteady Flows at High Angle of Attack, AIAA 20042136.