学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

jh170040 - NAH

反応・相変化を伴う多分散系混相流の大規模シミュレーション
松尾亜紀子（慶應義塾大学）
相変化・反応によって多分散性を有する分散相と、連続相の二相からなる混相流問題について、3 つの研究課題に取り組んだ。本課題研究では、「A：爆轟燃焼波と液滴との干渉」、「B：粉体の燃焼現象」、「C：水中高速飛翔体周りのキャビテーション現象」の 3点を対象とした昨年度の「A:爆轟燃焼波と液滴との干渉」の解析コード MPI 実装化に引き続き、今年度は「B：粉体の燃焼現象」及び「C:水中高速飛翔体周りのキャビテーション現象」について、MPI の実装化を実現した。さらに、解析コードのベクトル化率及び並列化効率の向上を行うことにより、数値解析の高速化に取り組んだ。その結果、計算コードの３次元化や計算時間の短縮等、大規模スケールで高効率な計算が可能となり、工学的スケー
ルにおける現象解析につなげることができた。

1. 共同研究に関する情報

2. 研究の目的と意義(1) 共同研究を実施した拠点名

1）研究目的東北大学サイバーサイエンスセンター

相変化・反応によって多分散性を有する分(2) 共同研究分野

散相と、連続相の二相からなる混相流問題に■ 超大規模数値計算系応用分野

ついて、大規模スケールでの数値解析を行う□ 超大規模データ処理系応用分野

ための OpenMP 並列／SIMD チューニングを□ 超大容量ネットワーク技術分野行うことを目指す。これらの現象は混相流の□ 超大規模情報システム関連研究分野

中でも、相反応・反応が生じたことによって(3) 参加研究者の役割分担

分散相と連続相が同時に存在するために非・ 松尾亜紀子（慶應義塾大学理工学部）

常に複雑な流れとなり、現象を工学的に再現研究総括することが非常に難しい。そのため、本研究・ 江川隆輔（東北大学サイバーサイエンス

ではこのような流れ場に対する大規模スケセンター）並列化処理ールでの数値解析手法を構築し、現象解明に・ 小林広明（東北大学サイバーサイエンスつなげることを目的としている。
本課題研究では、「A：爆轟燃焼波と液滴とセンター）並列化処理の干渉」、「B：粉体の燃焼現象」、「C：水中高・ 磯部洋子（東北大学サイバーサイエンス速飛翔体周りのキャビテーション現象」の 3センター）並列化処理
・ 奥野航平（慶應義塾大学理工学研究科）

点について扱う。A は数μm の液滴が流れ場「C：水中高速飛翔体周りのキャビテーシに存在し、燃焼波と液滴が干渉する気液混相

ョン現象」担当者

反応問題であり、B は数μm の粉体が燃焼す

・ 志村啓（慶應義塾大学理工学研究科）

るとともに既燃ガスと相互干渉をする固気「B：粉体の燃焼現象解明」担当者

混相反応問題、C は蒸発によって生じた気泡・ 渡部広吾輝（慶應義塾大学理工学研究科）が周囲液体と干渉する気液混相問題である。

「A：爆轟燃焼波と液滴との干渉」担当者

以下に各条件における状況及び学術的新規性について述べる。

・ 笠原弘貴（慶應義塾大学理工学研究科）(A)「爆轟燃焼波と液滴との干渉」「A：爆轟燃焼波と液滴との干渉」担当者

マイクロメートルオーダーの液滴が浮遊し1

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月ている状態において、伝播する燃焼波が通過端緒とする。

する。液滴は燃焼波通過後の高温気体によっ

2）研究意義て加熱され蒸発するとともに、相変化時の吸多分散系の混相流は数多く存在するが、そ熱によって燃焼波に対して影響を与える。これらに対して実用的な計算コストで解析をのような流れ場は、エンジン燃焼室内部にお行うモデルは存在しない。本課題研究によっいて多く生じているが、本研究においては燃て工学的スケールにおけるこれらの現象を焼波がデトネーションである場合について定量的に再現することのできる手法が開発

着目し、デトネーション波と液滴の干渉現象できれば、混相流現象を解明する一端とするを明らかにすることによって、デトネーショことが出来る。

ンの工業利用に貢献するとともに数値解析

本課題研究では特に 3 テーマに絞って解析手法の精度向上を目指す。

を行う事から、波及効果として以下の 3 点が(B)「粉体の燃焼現象」

期待できる。

マイクロメートルオーダーの粉体が堆積① デトネーションエンジンの実用化に向けしており、その上部を衝撃波が通過することた知見を収集し、推進機関の効率向上にによって粉体が巻き上げられ浮遊する。粉体貢献する。

が可燃性粒子であった場合、粉体が高温気体② 粉塵爆発現象の一端を解明し爆発現象にに晒されることによって燃焼を開始し、粉塵対する安全・予防技術の発展に貢献する。

爆発を生ずる。炭鉱の坑道においては微粉炭③ キャビテーションの有効利用に向けた提の燃焼による被害が繰り返し発生しており、言を行い、輸送機関の高速化に貢献する。

この抑制は炭鉱施設において常に課題として挙げられている。そのため、粉塵爆発現象

3. 当拠点公募型共同研究として実施した意義の詳細を把握することは、安全防護技術を発本研究課題における研究対象は、産業・自展させる上で必要とされている。本件では、然界において広く散見される現象であり、小粒子群の挙動を代表粒子によって表すこと規模スケールでの解析のみでは現象解明にで粉体の動きを精度よく再現する計算手法つながらない。数値解析によって得られた知を確立することが目的の一部である。

見を反映することによって混相流体現象に(C)「水中高速飛翔体周りのキャビテーショ対する研究の進展に貢献するためには、工学ン現象」的スケールにおける現象解析が必要である。

水中を 100 m/s 以上で移動する物体の表面

本課題研究において用いる混相流解析コ付近において生じるキャビテーション現象ードの構築は既に着手しており、それぞれのに焦点を当てる。水中飛翔体の周囲では、一問題に対して二次元での解析能力があるこ

時的な圧力の低下によってキャビテーショとは小規模計算において確認している。しかン気泡が発生し、飛翔体はキャビテーションし、工学的スケールにおける大規模計算を効に包まれることとなる。この際、気泡と水に率的かつスピーディーに行うためには、ベクおける粘性の差異によって飛翔体に加わるトル型コンピューターによる高並列計算を水中抗力は大きく変化する。本件では、この行う必要がある。そこで、ベクトル型コンピような高速飛翔体に加わる流体抵抗に対しューターの並列化性能を発揮させるため、自て、キャビテーション現象が与える影響を調主開発コードの OpenMP 並列／SIMD チュー査するとともに、水中高速輸送機関の実現へニングを行うものである。

向けたキャビテーションの有効利用研究の2

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

図１ 燃料液滴パージ法を用いた作動概要図
図 3 先端形状の後に分散し炭塵が燃焼するという層状炭塵爆発の過程をとる。本研究では、層状炭塵爆発の初期過程を二次元数値解析により再現し、火炎が加速する過程を観測することでその機構解明を試みた。
図 2 に 15.5 ms における反応面付近の流れ場を示す。形成されている火炎（x = 9.5 m）

図 2 15.5 ms における反応面付近の流れ場

は、予め流れ場に存在する空気と炭塵の揮発によって放出された CH4 が接する領域で生

4. 前年度までに得られた研究成果の概要じていることから、拡散火炎であるといえる。

(A)「爆轟燃焼波と液滴との干渉」

火炎後方上部の領域（x = 7.7 – 9.4 m）ではデトネーションを利用したエンジンにCH4 移流速度は 1000 m/s と火炎移流速度PDE(Pulse Detonation Engine)がある。PDE は、（600 m/s）を上回る速度で移流しており、デトネーションを繰り返し燃料器内で発生させ、CH4 移流が拡散火炎の超音速での移流を支推力を得る内燃機関である。PDE の燃料供給持していることが分かった。

法に燃料・酸化剤以外の物質を用いない燃料(C)「水中高速飛翔体周りのキャビテーショ液滴パージ法が提案されている（図 1）。

ン現象」

燃焼器形状が一次元的であり、デトネーシ

水中高速飛翔体周りに発生するキャビテョンの現象も一次元的に捉える事が出来るーション現象の数値解析を行った。水中高速事より、一次元数値解析を実施した。酸化剤

飛翔体周りにスーパーキャビテーションが供給全圧に対して燃焼器平均圧が高くなる

現れる場合、気泡は先端部分から発生し、先圧力ゲインな燃焼器実現のためには、燃焼端形状により気泡半径や流体抵抗値が異な器を短く、酸化剤供給全圧を低く、点火位置ることが先行実験により明らかになっていから酸化剤供給側に燃料・酸化剤混合気がる。そこで本研究では、気泡半径及び流体抵配置されている距離を大きく、当量比を高抗値を評価指標とし、水中高速飛翔体形状のくすることが有効であることが分かった。

最適化について検討した。最適化の検討には、(B)「粉体の燃焼現象」

他の値との優越を比較することで最適値を炭塵爆発では炭塵は重力沈降によって堆決定するパレート最適手法を用いた．

積していることが多いため、最初に燃焼波と本研究では、先端形状として図 3 に示すよの干渉により堆積した炭塵が巻き上がり、そうな尖頭形状（Ogival nose）、平頭形状（Flat3

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

図 5 気相デトネーション計算条件
図 4 Cd 値と気泡半径nose）及びスパイク形状（Spike nose）を用いた。
Flat nose は、Ogival nose の先端部を平らに切断したもので、切断面の半径（平頭半径）rf をパラメータとた。 Spike nose は Ogival nose の肩部の外周を直角に切り欠いた形状をしている。切(a) 𝜑𝜑 = 0.9, 𝑌𝑌𝐻𝐻2𝑂𝑂 = 0

り欠いた後の残った先端部をスパイク部と(b) 𝜑𝜑 = 0.9, 𝑌𝑌𝐻𝐻2𝑂𝑂 = 0.0701呼ぶ。Spike nose では平頭部の半径 rf 及びスパ

図 6 最大圧力履歴

イク部の半径 rs をパラメータとした。
図４より、Spike nose の平頭半径やスパイク径の値を調整することで、パレート最適解上で任意の Cd 値と気泡半径を実現する水中高速飛翔体形状を設計することが可能である。

5. 今年度の研究成果の詳細(A)「爆轟燃焼波と液滴との干渉」
水液滴が噴霧された C2H4-air 混合気中を伝播する気相デトネーションに関する二次元数値解析を実施した．本研究の目的は水液滴の挙動と水液滴が気相デトネーションの伝播挙動とセル構造に与える影響を明らかにする事である．計算対象と境界条件を図５に示す．
水液滴が存在する未燃混合気中を伝播する場合でも気相デトネーションは安定に伝播し、液滴が存在しない場合の CJ 速度と比

図 7 流れ場（𝜑𝜑 = 0.9, 𝑌𝑌𝐻𝐻2𝑂𝑂 = 0.0701）

較して 3.2%までの速度欠損が確認された．デトネーションと水液滴の干渉により主に衝する気相デトネーションに速度欠損を生じ撃波面と音速線までの領域で気相から液滴させたと考えられる．

へとエネルギーが輸送される．エネルギー欠

図 6 にセル構造を示す最大圧力履歴を示す．

損は水液滴が存在する未燃混合気中を伝播

水液滴の存在によってセル構造が変化した．
4

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

デトネーションと水液滴との干渉によって弱い三重点は減衰した．速度欠損によって反応誘導距離が増加し、セル幅がわずかに増加したと考えられる．
図 7 にデトネーション波が水液滴が存在する混合気中をおよそ 100 mm 伝播した際の瞬

図8

粉体の燃焼現象の計算対象間の流れ場を示す．衝撃波後方約 10 mm の地点で主に水液滴は蒸発し、蒸発している領域はデトネーション波と同じ速度で伝播した．
デトネーション後方の速度、渦度および温度の変動は水液滴の蒸発によって抑制された．
(B)「粉体の燃焼現象」
二次元 CFD-DEM シミュレーションを実施し、堆積炭塵による炭塵爆発現象に対する解析

図9

衝撃波と燃焼波の伝播速度履歴塵粒子との干渉によって減衰され衝撃波に

図 10

燃焼波付近の流れ場（16.0 ms）影響を与えなかった。

は、燃焼波に対して右側から左側へ流れていを行った。計算条件および境界条件を図 8 に示す。本条件は、先行研究における実験を模擬して小スケールで再現したものである。解析領域は 18 m x 0.05 m であり、底面から 2mm まで炭塵が空隙率 0.47 で配置されている。解析領域の左端から 2 m は CH4-Air 予混合気が充填されており、以降の 18 m は空気で満たされている。
図 9 に衝撃波と燃焼波の伝播速度履歴を示す。燃焼波の伝播速度は 350 – 500 m/s であり、大規模スケールでの実験結果と一致した。一方で、衝撃波の伝播速度は 500 m/s となり、数値解析による先行研究よりも低い値を示した。これは、炭塵粒子の巻き上げ挙動が異なっているためである。本計算結果では炭塵粒子が上壁に達するまで巻き上がったため、燃焼波から生じた圧縮波は浮遊する炭図 10 に燃焼波付近の流れ場を示す。燃焼

くのみであり、拡散火炎の構造とは一致しな波面に対して、流体は前後から流入している。

い。炭塵粒子は燃焼波面上で揮発を開始する

このとき、燃焼波の左側からは揮発によってのではなく、燃焼波の通過後に揮発を開始す

生じた CH4 が流入し、右側からは予め存在しる。このため、揮発によって生じた CH4 が燃ている空気が流入する。これにより燃焼波は焼波の左側から燃焼波へ流入する形となる。

拡散火炎の構造を持って伝播する。炭塵粒子

炭塵粒子は燃焼波に対して 400 m/s の速度で5

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

流入するため、Dust-Air 予混合気を伝播する

場合は急激に増加した。

層流火炎とは異なり、燃焼波の前面に予熱帯

水中移動体速度をパラメータとした気相体積分率の結果を図 12 に、抗力、揚力、転を形成しない。

倒モーメントの値を表 2 に示す。抗力は水中(C)「水中高速飛翔体周りのキャビテーショ移動体速度に依存するが、キャビテーションン現象」

発生場所は速度にはほとんど影響を受けな

これまでは２次元解析によって得られたいことが分かる。このことから、スーパーキ流体力を基に剛体運動方程式を立て、飛しょャビテーションを利用した水中移動体の設う軌跡及び姿勢を求めた。しかし、迎角を持計においては、使用速度ごとに形状を変更すった軸対象物体の飛しょう安定性を検討するという必要は少ないと考えられる。

るためには３次元解析による流体力の算出

6. 今年度の進捗状況と今後の展望が必要である。そこで、本研究ではスーパーキャビテーションを伴う水中移動体の飛し昨年度は「A:爆轟燃焼波と液滴との干渉」

ょう安定性について、3 次元計算コードによの解析コードについて MPI の実装を実施した。

り水中移動体の 3 次元的な姿勢及び速度の影

今年度は「B：粉体の燃焼現象」及び「C:水中響を解析した。

高速飛翔体周りのキャビテーション現象」に

迎角をパラメータとした気相体積分率の結果を図 11 に、抗力、揚力、転倒モーメン
トの値を表 1 に示す。迎角にかかわらず抗力はほぼ同じ値となった。揚力及び転倒モーメ
ントは水中移動体が気泡内にある場合は迎角に比例して増加するが、水からはみ出した

図 12
表１

迎角をパラメータとした気相体積分率6

揚力、抗力、転倒モーメント（各迎角）

AoA

Drag [N]Lift [N]Moment [Nm]

0˚
5˚
10 ˚
45 ˚

26,17325,781
26,28923,58705,321
10,040132,138

01,9253,64128,394

表2

図 11

水中移動体速度をパラメータとした
気相体積分率

揚力、抗力、転倒モーメント（各速度）

Velocity

Drag [N]Lift [N]Moment [Nm]450 m/s
600 m/s
1000 m/s

4,4147,949
26,173000000学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

ついて、MPI の実装を行う予定としていたと
ころ、予定通り両方の解析コードへの MPI の実装に成功した。
また、解析コードのベクトル化率及び並列化効率の向上を行うことにより、数値解析の高速化に取り組んだ。「B：粉体の燃焼現象」
図 13 サブルーチンの初期コードを例に、その概要及び結果を示す。
（１）ベクトル化率の向上「B：粉体の燃焼現象」の解析コードの中には、ベクトル化率が低く、ベクトル長が短いサブルーチンが存在していた。当該サブルーチンの内部を詳細に解析したところ、ループ内に依存関係があり、最内ループがベクトル化されていないことが分かった(図 13）。また、最内ループは外側のループ毎にループ長が異なり、かつそのループ長は短い。外側ループのループ長が長いため、外側ループを最内ループに移動することで性能改善が期待できる。

図 14 ベクトル化後のコード

そこで、ループの一部を分割し、最内ループのループ長のうち最大になるものを検索する。次に、最内ループの終端を検索した最大ループ長に変更し、ループ間の依存関係を解消する。最大ループ長に満たないループはマスクにより処理を行わないようにしている。
依存関係が解消したため、ループ交換を行い、外側のループを最内側ループに移動することでループ長を長く確保し、ベクトル化率を高めた(図 14)。さらに、ベクトルプロセッサに搭載されている ADB を有効に活用できるように最内ループのブロック化を行う。 SXACE は１命令で 256 要素を同時に演算できるため、ブロック長は 256 とした(図 15）。
高速化後の性能を表 3 に示す。高速化によりベクトル化率が 12.23%から 98.91%に改善し、ベクトル長が 12.3 から 151.3 に伸びること
で、サブルーチン単体で約 6 倍の性能改善を
図 15 高速化後のコード得ることができた。

7

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

や、ノード内の並列化に OpenMP を用いた

表 3 サブルーチンの高速化効果

Hybrid MPI などの検討が必要である。

7.

研究成果リスト(1) 学術論文•奥野 航平, 松尾 亜紀子, 水中移動体の先端形状が飛しょう安定性に与える影響に関する数値解析, 日本機械学会論文集, Vol.83, No.856,
2017(2) 国際会議プロシーディングス該当なし(3) 国際会議発表
図 16 並列化性能

・Kei Shimura, Akiko Matsuo, NumericalInvestigation on the Effect of Particle（２）並列化率向上

Collisions in the Process of Dust

並列化にはノード内のコア毎にプロセスをLifting, The Ninth JSME-KSME Thermal

処理する MPI 並列を行う Flat MPI を用いた。

and Fluids Engineering Conference

本コードの計算領域は 1,901×161 の 2 次(TFEC9), Okinawa, Japan, Oct. 27 - 30,

元空間で、この空間中に粒子が存在する。こ2017の 2 次元空間中の熱や圧力等の計算後、粒子・Hiroaki Watanabe, Akiko Matsu, Numericalの相互作用の計算を行う。並列化をするにあinvestigation

ontheたり、計算領域を分割し分割後の各領域にbehaviorofhighpressureMPI のプロセスを割り当て、各領域の熱や圧releasedfrom力等と存在する粒子に対して相互作用の計International算を行うように並列化を行った。

Safety (ICHS 2017), Hamburg, Germany,

図 16 に並列化後の性能を示す。2 次元分割

Sept. 11 - 13, 2017

を行うと計算の初期で粒子が存在しない領・

theself-ignition

tube,

ConferenceonhydrogenThe

7th

HydrogenKen Matsuoka, Haruna Taki, Jiro Kasahara,

域が多数存在し、プロセス間のインバランス

Hiroaki Watanabe, Akiko Matsuo, Takuma

が大きくなるため、計算領域を x 方向(1,901)

Endo, Study on a High-Frequency Pulse

での 1 次元分割とする。高速化後のシングルDetonation Operation, 1st International

コア実行に比較し、64 プロセス実行で約 9Constant倍の性能改善を得るにとどまった。実行後、Combustion Workshop, Poitiers, France,

時間が経過すると、各領域に存在する粒子数

June 13 – 16, 2017に偏りが発生し、プロセス間のインバランス

・VolumeandDetonation

Hiroaki Watanabe, Akiko Matsuo, Haruna

により、性能向上が抑制されていると考えら

Taki,える。

Numerical Investigation on the Effect of

今後、さらなる高速化のため、多数の粒子がFuel Supply in High Frequency Pulse

存在する領域をさらに細かく分割すること

Detonation Operation, 1st International
8KenMatsuoka,

JiroKasahara,

学際大規模情報基盤共同利用・共同研究拠点 平成 29 年度共同研究 最終報告書 2018 年 5 月

Constant

VolumeandScitech 2018, Florida, USA, Jan 8-12,

Detonation

2018

Combustion Workshop, Poitiers, France・June 13 – 16, 2017・Kohei

Okuno,

AkikoMatsuo,‘The

Haruna Taki, Kazuyuki Takao, Naritoshi

numerical analysis of the nose shape

Hirota, Ken Matsuoka, Jiro Kasahara,

effect on the flight stability of theHiroaki Watanabe, Akiko Matsuo, Takuma

underwater

Endo, Investigation of High-Frequency

International Conference on Numerical

Pulse Detonation Cycle with Fuel Phase

Methods

Transition, 31st International Symposium

Japan, June 26-29, 2017

inprojectile’,
Multiphase

The

Flows,

3rdTokyo,on Shock Wave, Nagoya, Japan, July 9 –14,・2017

HarunaTaki,Naritoshi

Hirota,(4) 国内会議発表

Ken・Matsuoka, Akira Kawasaki, Jiro Kasahara,Hiroaki Watanabe, Akiko Matsuo, Takuma

ンによる衝撃波によって誘起される堆積炭塵

Endo,の燃焼に関する解析, 第 31 回数値流体力学シ

Pulse

Detonation

Operation

atンポジウム, 京都, 2017 年 12 月

Kilohertz Frequency, 26th International
・Colloquium on the Dynamics of Explosions

・笠原次郎，低背圧環境下で作動するパルスデ

30 – August 4, 2017

トネーション燃焼器の推力性能に関する数

Kei Shimura, Akiko Matsuo, Numerical

値解析，第 49 回流体力学講演会/第 35 回航

Investigation of Dust Lifting Induced by

空宇宙数値シミュレーション技術シンポジ
ウム，東京，2017 年 6 月
・Symposium on Shock Wave, Nagoya, Japan,広吾輝，松尾亜紀子，遠藤琢磨，パルスデKei Shimura, Akiko Matsuo, Numerical

トネーション燃焼器の高周波数作動に関する実験的研究，第 49 回流体力学講演会/第35 回航空宇宙数値シミュレーション技術シ

of Layered Coal Dust Combustion, 26thンポジウム，東京，2017 年 6 月

International Colloquium on the Dynamics
ofExplosions

andReactive

・System,

奥野航平, 松尾亜紀子, 水中移動体後端部

Boston, USA, July 30 – August 4, 2017が受ける水圧力による飛しょう安定性に関Hirotaka

する数値解析, 日本混相流学会 混相流シン

Kasahara,Akiko

Matsuo,ポジウム 2017、 東京、 2017 年 8 月

Numerical Investigation on the EffectsofAir

Dissociation

upon

Hypersonic

Projectile in Standard Atmospheric Air,(5) その他（特許，プレス発表，著書等）

31st International Symposium on Shock

該当なし

Wave, Nagoya, Japan, July 09 – 14, 2017
・瀧春菜，松岡健，川崎央，笠原次郎，渡部July 09 – 14, 2017
Investigation on the Initial Development

・渡部広吾輝，松尾亜紀子，瀧春菜，松岡健，

and Reactive Systems, Boston, USA, July

Vertical Shock Wave, 31st International

・志村啓, 松尾亜紀子, CFD-DEM シミュレーショ

Hirotaka Kasahara, Akiko Matsuo, The
Effect of Hypersonic Projectile Shape on
AerodynamicandThermal

Performance

Launched by Ground-based Railgun, AIAA
9