極低温推進薬の軌道上貯蔵・輸送に向けた減圧沸騰現象の解明と予測手法の開発
谷 和 磨
東京大学大学院工学系研究科

１．はじめに宇宙分野においても極低温流体の大規模輸送・貯蔵技術の利用が期待されている。近年、米国を中心に月・火星有人探査に向けた動きが活発になってきており、そうした宇宙空間での長期間ミッションの実現に際し、軌道上燃料再補給船 (Space depot)[1] や液体水素を燃料とする
核熱推進ロケット (Nuclear Thermal Propulsion Rocket)[2] などの利用も検討されており、極低温燃料を軌道上や低重力環境において長期間管理する技術が重要な役割を占めるようになると予想されている [3,4,5]。
また、近年の環境問題やエネルギー問題の解決に向けた活発な動きに伴い、これまで主に液体ロケット等の燃料として利用されてきた極低温流体の大規模利用の場面が増えつつある。水素社会を実現するために、日本政府は水素基本戦略を 2017 年に策定し、2030 年頃に商用サプライチェーンを構築し、2017 年現在の年間 0.02 万トンから 30 万トンへ水素流通を増やすことを掲げている[6]。この目標を達成するために構築される商用サプライチェーンの一つが The
Hydrogen Energy Supply Chain である[7,8]。この計画では、2030 年に液化水素をオーストラリアで褐炭から生産し、160000 m3 タンクを保有する液化水素運搬船で日本へ輸送することを目指している。
こうした将来的な軌道上燃料再補給船や国際水素サプライチェーンなど、大規模かつ長期間極低温流体を貯蔵し、輸送するプロジェクトの共通する課題の一つとして、減圧沸騰に関する問題がある。具体的には、Fig.1 に示すように、外部からの入熱で発生したボイルオフガスによりタンク圧力は高くなり、タンクの限界圧力を超える危険性がある。これを防ぐために、放圧 (ベント) と呼ばれるタンク内の気体を外に逃すことでタンク内を減圧する。その際、減圧に伴う飽和温度の低下により、液体が過熱状態になり沸騰する現象である「減圧沸騰」が発生する。この現象により、タンク内の圧力制御が難しくなる上、突沸が発生しタンクの限界圧力を超えてしまう可能性や、発生した気泡により液面が押し上げられ液流出に繋がる危険性がある。宇宙においては、地上での課題に加えて、界面張力が重力を卓越するため気液の分離が難しい課題も付随する。
従来の減圧沸騰を対象にした研究は、高々0.1m3 程度の小寸法のタンクでの水・液体窒素を対象にした実験が中心である[9-20]。タンクの安全な運用の面で問題となりうる非定常な圧力回復

Fig. 1 Schematic of liquid state inside liquid hydrogen tanker during shipping.

スーパーコンピューティングニュース- -1 48- -Vol. 25, No.22023

現象(圧力が減少後上昇する)やそれに伴う液面上昇も確認されたが、その発生機構についての研究は少なく、大型タンクでの同様の実験事実が得られるのか十分に明らかにされていない。また、極低温流体を用いた研究も少なく、物性の影響に関しても不明な点が多い。よって、これらの問題を解決するためには、極低温流体の大型タンクにおける減圧沸騰現象に対する数値流体解析を実施する必要がある[21,22]。しかし、相変化現象を取り扱う難しさから大型タンク内部の詳細な熱流体挙動の分析は少なく、その現象予測手法の開発が求められている上、減圧が数秒で完了する小タンクと異なり大型タンクをフルスケールで数値流体解析を行うのは計算コストの面で難しい。よって、大型タンク内部の熱流動現象を低コストで予測可能な手法の開発が急務である。
以上より、本研究の目的は、極低温推進薬の大規模・貯蔵輸送技術に向けて、大型タンク内の減圧時の圧力・液面挙動の予測技術を確立することを目的とする。

２．数値流体解法
本研究では、圧力回復現象の解明を目的に数値流体解析(CFD)を実施した。解析には、Himenoら[23]や Umemura ら[24,25]によって構築された自由表面流数値解析手法(CIP-LSM: CIP-basedlevel set & MARS)を用いた。
均質二相流の支配方程式は、質量・運動量・内部エネルギーの式

𝜌𝜌𝐷𝐷𝐷𝐷= −𝜌𝜌∇ ∙ 𝑢𝑢)⃗𝐷𝐷𝐷𝐷𝐷𝐷𝑢𝑢)⃗= −∇𝑝𝑝 + ∇(𝚻𝚻𝝂𝝂 + 𝚻𝚻𝝈𝝈 ) + 𝜌𝜌𝑔𝑔⃗
𝐷𝐷𝐷𝐷
𝐷𝐷𝐷𝐷= −𝑝𝑝∇ ∙ 𝑢𝑢𝜌𝜌)⃗ + Θ̇
𝐷𝐷𝐷𝐷

(1)

2𝜇𝜇
+⃗)𝐈𝐈 + 𝜇𝜇(𝐃𝐃 + 𝐃𝐃" ),𝚻𝚻𝛎𝛎 = − (∇ ⋅ 𝑢𝑢
3𝐃𝐃 = ∇u+⃗,++++⃗𝑛𝑛𝚻𝚻𝝈𝝈 = 𝜎𝜎𝛿𝛿$ (𝐈𝐈 − 𝑛𝑛$ ++++⃗),$

Θ̇ = (𝚻𝚻𝝂𝝂 + 𝚻𝚻𝝈𝝈 )∇ ⋅ 𝑢𝑢+⃗ − ∇ ⋅ 𝑞𝑞⃗.

を解く流体解法の部分と、気液を識別する関数𝐻𝐻# の移流方程式
𝐷𝐷𝐻𝐻$ 𝜕𝜕𝐻𝐻$=+ (𝑢𝑢+⃗ ⋅ ∇)𝐻𝐻& = 0
𝐷𝐷𝐷𝐷
𝜕𝜕𝜕𝜕(2)

𝐻𝐻& = 0.5
𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿𝐿
𝐻𝐻& = 0
𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼
𝐻𝐻& = −0.5
𝐺𝐺𝐺𝐺𝐺𝐺(3)を解く界面追跡法の部分で構成される。

CIP-LSM で は 、 流 体 解 法 と し て CIP-CUP 法[26] の 一 種 で あ り ， 温度を 独 立 変 数 と す る
TCUP(Thermo CIP−CUP)法[23]を，界面追跡法として PLIC-VOF 法の一種である MARS 法[27]を採用している。また、CLSVOF(Coupled Level Set Volume-of-Fluid)法[28]と同様に、界面再構築には密度関数法の一種である Level-Set 法[29]を用いて、数値拡散を防ぎ界面を高精度で追従する仕様となっている。まず相変化を考慮した TCUP 法について説明する。局所熱平衡を仮定して相変化項を加えると、式(1)は基本量𝑄𝑄 = (𝒖𝒖, 𝑇𝑇, 𝑝𝑝)$ を状態量とする次式に変換することができる。

スーパーコンピューティングニュース- 2- - 49

Vol. 25, No.22023𝜌𝜌$⃗𝐷𝐷𝑢𝑢= −∇𝑝𝑝 + ∇(𝚻𝚻𝝂𝝂 + 𝚻𝚻𝝈𝝈 ) + 𝜌𝜌𝑔𝑔⃗ + 0𝑢𝑢$⃗#$% − 𝑢𝑢$⃗&'( 1𝑚𝑚̇δ)𝐷𝐷𝐷𝐷
𝜌𝜌+ 𝐷𝐷𝐷𝐷𝐷𝐷𝐷𝐷̇𝜌𝜌𝐶𝐶*= −𝑇𝑇+Θ𝜌𝜌 𝐷𝐷𝐷𝐷𝐷𝐷𝐷𝐷𝜌𝜌+ Θ̇
111 𝐷𝐷𝐷𝐷= −∇ ∙ 𝑢𝑢$⃗ −−:−; 𝑚𝑚̇𝛿𝛿%𝜌𝜌#$% 𝜌𝜌&'(𝜌𝜌𝐶𝐶%, 𝐷𝐷𝐷𝐷𝜌𝜌 𝜌𝜌𝐶𝐶*(4)Θ̇ = (𝚻𝚻𝝂𝝂 + 𝚻𝚻𝝈𝝈 )∇ ⋅ 𝑢𝑢+⃗ − ∇ ⋅ 𝑞𝑞⃗ − Jℎ𝑔𝑔𝑔𝑔𝑔𝑔 − ℎ𝑙𝑙𝑙𝑙𝑙𝑙 K𝑚𝑚̇ δs .

このうち、相変化を考慮することで生成された項
0𝑢𝑢$⃗#$% − 𝑢𝑢$⃗&'( 1𝑚𝑚̇δ)−>ℎ#$% − ℎ&'( ?𝑚𝑚̇δ)
11−:−; 𝑚𝑚̇𝛿𝛿%𝜌𝜌#$% 𝜌𝜌&'((5)

は順に、運動量損失と生成、潜熱、気液密度差による体積湧き出しを示す。
TCUP 法では、支配方程式(4)の左辺の流体微分項を推定する移流項と右辺の非移流項に大別される。移流項では、矢部ら[26]が提案した双曲型偏微分方程式を解く高次精度差分法である CIP法(Constrained Interpolation Profile Scheme)により、各状態量の値を更新する。非移流項では粘性項と伝熱項を陰的に評価し、速度と温度をそれぞれ修正する拡散段階と、それに続いて圧力項を陰的に評価して圧力を修正する音響段階に分けて解く。音響段階でポアソン方程式を解く
際には、前処理として SA-AMG 法[30]を用い、非定常反復法の安定化双共役勾配法(Bi-CGSTAB 法)で解いている。また、大規模問題に対応するために、MPI を用いた領域分割法による解法の並列化も行っている。また、界面追跡法の部分に関しても、相変化による界面からの体積湧き出しを考慮すると次式を導くことができる。
𝜕𝜕𝐻𝐻$
𝑚𝑚̇+ (𝑢𝑢𝑛𝑛+⃗ ⋅ ∇𝐻𝐻$ = 0+⃗$ ⋅ ∇)𝐻𝐻& −𝜕𝜕𝜕𝜕𝜌𝜌'($ $𝜕𝜕𝐻𝐻$
𝑚𝑚̇+ (𝑢𝑢𝑛𝑛+⃗ ⋅ ∇𝐻𝐻$ = 0+⃗$ ⋅ ∇)𝐻𝐻& −𝜕𝜕𝜕𝜕𝜌𝜌)*+ $(6)相変化がない場合は𝑚𝑚̇ = 0で、かつ各相の速度と界面の速度が一致する(𝑢𝑢)⃗%&' = 𝑢𝑢)⃗()* )ので、)⃗# = 𝑢𝑢

一般的な識別関数の移流方程式となる。この式の解法には、二値関数の有限体積的な移流に優れた MARS 法を採用している。

相変化現象は、その複雑な機構から支配方程式で現れる𝑚𝑚̇を正確に計算するためには格子解像度を高くして計算する必要がある。しかし、相変化計算をする際には格子解像度を高くすると計算コストが増大するため、格子解像度が低くても相変化量を見積もれるモデルを導入する必要がある。そのため、Fig.2 のような 3 段階の相変化モデルを採用する。まず CFD で捉えられる相変化と、CFD では捕らえられないほど小さい相変化現象に分けて考え、前者は CFD により直接計算し、後者はサブグリッドスケールの相変化モデルで計算する。サブグリッドスケールの相変化モデルは 2 段階のモデルとなっており、壁面での気泡生成などを壁面沸騰モデルで、液相中での相変化を温度回復法ベースのモデルで計算する。これにより計算コストを抑えつつ、壁面セルでの気泡生成・成長から液相中での気泡成長をサブグリッドスケールモデルでも扱うことができる。
第一段階の壁面沸騰モデルでは、核沸騰と膜沸騰の両方を扱うことができるが、本研究では高い熱流束が必要な膜沸騰は無視する。核沸騰フェーズでは気泡生成の活性化条件Δ𝑇𝑇 K が既存
気泡核における活性化条件
Δ𝑇𝑇 =!
!
/$.
"#$% "&'( )*+ 0+,-

(,*-

1.

,(7)

により計算され，液体の過熱度がこれより高い場合に気泡が生成される。

スーパーコンピューティングニュース- -3 50- -Vol. 25, No.22023そして、壁からの入熱により成長する気泡(Fig.3)の成長速度は下式により計算する。
𝑚𝑚̇(𝜑𝜑) =2&'( ($/*00 .$)*+ )(,*- ×1(678 9.678 :)(8),+B2&'( ($/*00 .$)*+ )1
0CD678 9𝑉𝑉̇ = ; ∮=>??@A 𝑚𝑚̇(𝜑𝜑)𝑑𝑑𝑑𝑑(𝜑𝜑) =log C .
;(#$%#$% ,*-(9)

気泡がある程度成長すると、第 2 段階の温度回復法ベースの沸騰モデルに移行し、
𝑚𝑚̇ =;.100 E.100 F- GHI (J$&'( .$)*+ K,M$)2- )(,*- × NO,(10)を適用し沸騰量を計算する。温度回復法を、格子内でのボイド率が 0.1%から 50%の場合の格子かつ液温が飽和温度より高いものに適用することで、壁面沸騰モデルから CFD の相変化計算に移行できるようにした。
こうして気泡が生成・成長して、一格子あたりに占める気体体積の割合が 50%を超えると、CFDの相変化の計算に移行する。Fig.4 のように CFD での相変化量の計算では、式(11)を用いて格子点上での温度と界面での飽和温度の温度勾配から沸騰量を計算する。𝜙𝜙は Level Set 法の概念に従った気液界面からの距離関数である。

𝑚𝑚̇ =JPQ⃗&'( .PQ⃗#$% K⋅TQ⃗

𝑞𝑞⃗ = −κ∇T~ −(,*-,

κ(𝑇𝑇 − 𝑇𝑇#UO )
∇𝜙𝜙,
𝜙𝜙

(11)

(12)

Fig. 2 Phase change model.

Fig. 3 Definition of a coordinate axis.

Fig. 4 Temperature distribution around liquid-vapor interface.

スーパーコンピューティングニュース- 4- - 51

Vol. 25, No.22023

3．解析条件
計算対象は、Fig.5 に示すように、JAXA 能代ロケット実験場の 30m3 液化水素タンク(直径 2.3m・高さ 7.622 m)[21,22]の断面形状を模した 2 次元計算格子と軸対称 1/4 領域の 3 次元計算格子とした。2 次元格子は，全部で 33 の領域から成り立っており、総格子点数は 551,520 点である。3次元格子は、129 領域で、総格子点数は 1,968,131 点である。
設定した初期条件は、液位は 6.248m で、系内の温度は実験値から線形補間して与え、初期圧力は 321.8kPaG で液体に対しては静水圧も課している。境界条件について、2 次元格子に関してのみ、下から 3 個目のドメインから 12 個目の中心に等温の粘着壁を設定し，仮想的に温度計プローブが存在するとした。同様に、タンク側壁にも等温の粘着壁条件を課している。3 次元格子に関しては、タンク壁に等温の粘着壁条件、半径方向壁面は滑り壁を与えている。流出境界条件
は、2 次元と 3 次元格子に関して共に Fig.6 のように実験で得られた質量流量をタンク上部の管の先端に与えている[22]。計算は、全て東京大学情報基盤センターが所有するスーパーコンピュータ Oakbridge-CX を用いて、二次元計算は 256 並列で３次元計算は 1024 並列で実施した。
4．解析結果・考察
実験[21,22]及び CFD で得られた圧力履歴とそれに対応させた CFD 可視化結果を Fig.7 に示す。
圧力履歴の結果より、定性的には数値計算により圧力回復の傾向をとらえることに成功を収めた。
2 次元計算の可視化結果と圧力履歴を比較しつつ考察していくと、圧力が落ちていく 0 から 6秒付近までは沸騰が起こっているものの、流出流量による減圧が支配的であるため、圧力は減少していく。液位に関しては 6 秒で最大値をとるとともに、圧力も極小値をとる。そして、6 秒後(a) 2-D(b) 3-D

Fig. 5 Computational domain and initial condition (H: liquid occupation ratio. H=0.5: Liquid, H=-0.5: Gas)

Fig. 6 Time variation of vent gas mass flow rate

スーパーコンピューティングニュース- -5 52- -Vol. 25, No.22023

以降から圧力が上昇に転じ、CFD においても圧力回復現象が確認できる。Fig.7(a)の可視化画像の下段の右端に示すように、7 秒後以降では気泡により液体がよく混ざり圧力がゆっくり減少していくことがわかる。つまり、圧力回復現象は、気液界面面積が増加するまでに時間がかかるため、ベント量が沸騰量に上回ってしまうことが原因で発生する。そして、液相上部が気泡で攪拌され気液界面面積が十分に大きくなれば、系全体が平衡状態に到ることが考察される。3 次元計算においても、Fig.7 (b)に示すように同様の圧力回復の傾向が見られ、2 次元計算で考察されたことが 3 次元計算でも成立することを示す。
ただし、3 次元 CFD 計算において実験より圧力回復の幅が大きいことが見て取れる。これは、今回の CFD 計算ではサブグリッドスケールの第 2 段目の相変化モデルでは凝縮の影響を考慮していないため、沸騰により系の圧力が高くなっても気相部での凝縮が進まず、沸騰の影響を緩和できていないためと考察される。よって、今後の研究では第 2 段目に温度回復法ベースの凝縮モデルを追加して計算を実施する必要があると言える。また、3 次元 CFD 計算の可視化結果より、液面の上昇が起こっていないことがわかる。そのため、今回の減圧速度の範囲では、液面が上昇して外部に液化水素が流出することなどは起こりにくいことが示唆される。

5．結論
本研究では、大規模液化水素タンク内の高液位条件下での減圧沸騰現象を明らかにする CFD解析を実施した。圧力回復現象の２次元および３次元 CFD 計算から、減圧開始時には壁からの(a) 2-D results(b) 3-D results
Fig.7 CFD and experimental result of time variation of pressure and visualized results.

スーパーコンピューティングニュース- 6- - 53

Vol. 25, No.22023

沸騰がないため、圧力低下が支配的であることがわかった。その後、沸騰による気液界面積の増加に伴い圧力が上昇し、最終的に平衡状態に達することが示された。また、今回の減圧速度の範囲では、圧力回復時に液界面の急激な上昇が発生せず、液の外部流出の危険性や気相部が押し潰されて圧力が急上昇することは起こり得ないことがわかった。
また、サブグリッドスケールの沸騰モデルを二段階に分け、一段目に壁面沸騰モデル、二段目に温度回復法に基づく沸騰モデルを導入することで、減圧沸騰試験での圧力回復現象を 3 次元
CFD 解析により再現した。ただし、定量的な観点では実験値との解離があるため、凝縮モデル等に関する考察を行う必要があると言える。

参1.

考文

献

Goff, J., Zegler, F., Bienhoff, D., Marchetta, J., Kutter, B., & Chandler, F..
Realistic near-term propellant depots: implementation of a critical spacefaring
capability. In AIAA Space 2009 Conference & Exposition (p. 6756) (2009).

2.

Tanner, C.L., Young, J.J., Thompson, R.W., & Wilhite, A.W. IAC-06-D 1 . 1 . 01 OnOrbit Propellant Resupply Options for Mars Exploration Architectures (2006).

3.

Muratov, C. B., Osipov, V. V., & Smelyanskiy, V. N. Issues of long-term cryogenic
propellant storage in microgravity (No. ARC-E-DAA-TN4295). (2011).

4.

Chato, D. J. Cryogenic fluid transfer for exploration. Cryogenics, 48(5-6), 206-209.

5.

Meyer, M. L., Chato, D. J., Plachta, D. W., Zimmerli, G. A., Barsi, S. J., Van(2008).

Dresar, N. T., & Moder, J. P. Mastering cryogenic propellants. Journal of Aerospace
Engineering, 26(2), 343-351 (2013).
6.

Hydrogen and Fuel Cells Strategy Office, Advanced Energy Systems and Structure
Division, Energy Conservation and Renewable Energy Department, Agency for Natural
Resources and Energy Japan. Basic Hydrogen Strategy. 1–37 (2017).

7.

Kamiya, S., Nishimura, M. & Harada, E. Study on introduction of CO2 free energy toJapan with liquid hydrogen. in 25th International Cryogenic Engineering Conference

and the International Cryogenic Materials Conference in 2014, ICEC 25–ICMC 2014 vol.
67 11–19 (Elsevier B.V., 2015).
8.

Takaoka, Y., Saeed, A., Nishifuji, K. & Kanbe, K. Design and Operation of the First
LH2 Carrier. in Gastech 2019 Exhibition & Conference 1–10 (2019).

9.

Miyatake, O., Murakami, K., Kawata, Y. & Fujii, T. Fundamental Experiments of Flash
Evaporation. Bull. Soc. Sea Water Sci. Japan 26, 189–198 (1972).

10. Miyatake, O., Fujii, T., Tanaka, T. & Nakaoka, T. Flash evaporation phenomena of poolwater. Kagaku Kogaku Ronbunshu 1, 393–398 (1975).
11. Saury, D., Harmand, S. & Siroux, M. Experimental study of flash evaporation of awater film. Int. J. Heat Mass Transf. 45, 3447–3457 (2002).
12. Saury, D., Harmand, S. & Siroux, M. Flash evaporation from a water pool: Influence ofthe liquid height and of the depressurization rate. Int. J. Therm. Sci. 44, 953–965(2005).

スーパーコンピューティングニュース- -7 54- -Vol. 25, No.2

2023

13. Cheng, H. P. & Lin, C. T. The morphological visualization of the water in vacuum
cooling and freezing process. J. Food Eng. 78, 569–576 (2007).
14. Zhang, D., Chong, D., Yan, J. & Zhang, Y. Study on steam-carrying effect in staticflash evaporation. Int. J. Heat Mass Transf. 55, 4487–4497 (2012).
15. Wang, C., Xu, R., Chen, X., Jiang, P. & Liu, B. Study on water flash evaporation
under reduced pressure. Int. J. Heat Mass Transf. 131, 31–40 (2019).
16. Hewitt, H. C. & Parker, J. D. Bubble Growth and Collapse in Liquid Nitrogen. J. HeatTransfer 90, 22–26 (1968).
17. Shingo, Y., Yutaka, H. & Ikuo, T. An Experimental Study on Flashing Phenomenon ofLiquid Nitorgen. Trans. Japan Soc. Mech. Eng. Ser. B 58, 1498–1503 (1992).
18. Toshiaki, W., Yutaka, H. & Ikuo, T. Flashing Phenomena of Liquid Nitrogen in aPressure Vessel. Trans. Japan Soc. Mech. Eng. Ser. B 61, 1849–1854 (1995).
19. Toshiaki, W., Yutaka, H. & Ikuo, T. Flashing Phenomena of DepressuTized Liquid
Nitrogen in a Pressure Vessel. Trans. Japan Soc. Mech. Eng. Ser. B 63, 2148–2153(1997).
20. Takeda, M., Usui, T. & Maekawa, K. Study on boiling behavior of pressurized liquid
nitrogen under rapid depressurization. IOP Conf. Ser. Mater. Sci. Eng. 502, (2019).
21. Tani, K. et al. Prediction of Pressure Reduction Rate in 30m3 Liquid Hydrogen Tank
Based on Experimental And Numerical Analysis. in International Conference on Hydrogen

Safety 2019 1–12 (2019).
22. Tani, K. et al. Pressure Undershoot During Pressure Reduction Experiment with LargeScale Liquid Hydrogen Tank. 23rd World Hydrog. Energy Conf. (2020).
23. HIMENO, T. et al. Numerical Analysis of Free-Surface Flows under Various Conditionsin Acceleration : Improvement of CIP-LSM : CIP-Based Level Set & MARS(The Forefront ofMulti-Physics CFD/EFD). Trans. Japan Soc. Mech. Eng. Ser. B 76, 778–788 (2010).
24. Umemura, Y. et al. Liquid nitrogen chill-down process prediction by direct interface
tracking approach. 53rd AIAA/SAE/ASEE Jt. Propuls. Conf. 2017 1–12 (2017).
25. Umemura, Y., Himeno, T. & Watanabe, T. Numerical analysis of boiling flow in surfacetension-dominated environment. 49th AIAA/ASME/SAE/ASEE Jt. Propuls. Conf. 1 PartF, 1–14 (2013).
26. Yabe, T. & Takei, E. A New Higher-Order Godunov Method for General HyperbolicEquations. J. Phys. Soc. Japan 57, 2598–2601 (1988).
27. Kunugi, T. MARS for Multiphase Calculation. Comput. Fluid Dyn. J. 9, 563–571 (2001).
28. Sussman, M. & Puckett, E. G. A Coupled Level Set and Volume-of-Fluid Method for
Computing 3D and Axisymmetric Incompressible Two-Phase Flows. J. Comput. Phys. 162,301–337 (2000).
29. Osher, S. & Sethian, J. A. Fronts propagating with curvature-dependent speed:
Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 79, 12–49 (1988).

30. Vaněk, P., Mandel, J. & Brezina, M. Algebraic multigrid by smoothed aggregation forsecond and fourth order elliptic problems. Computing 56, 179–196 (1996).

スーパーコンピューティングニュース- 8- - 55

Vol. 25, No.22023