蒸気爆発は、熔融金属と水など大きな温度差の液と液が直接接触した際に生じる激しい爆発現象である。この現象は原子力の苛酷事故において発生することがあるため、近年、原子力安全に関連して研究が進められている。実験と並び、近年の蒸気爆発の素過程についての知見の蓄積や多相流シミュレーション技術の向上などを背景として、数値解析も進められている。汎用の多相多次元コードCHAM PAGNEを基礎とした,よりメカニスティックな蒸気爆発シミュレーションコードCHAMP/VEを開発することを目的として研究を行なった。

　蒸気爆発の過程は、複数の物質と複数の相を含む(図1)。したがって、多相の熱流動問題を扱うことになるが、従来の現象に比べ速い相変化と体系の変化を伴うために、非常に複雑な物理現象が現れる。そのような複雑なメカニズムがなかなか解明されなかったため、どのような物理モデルを選択するかということも難しい問題である。

　また、ある素過程については、数値解析コードに適用できる簡単なモデルは未だに存在しない。よって従来は、比較的簡単な粗混合と微粒化過程しか数値計算を行なっていない。

　更に、激しい蒸発のために体系内の変数の空間的・時間的変化が強いことから、数値不安定と数値拡散の問題が生ずる。したがって、蒸気爆発の解析では、基礎式とモデルの他に、数値解法や計算技術などにも特に厳しい要求がなされる。

　この研究では、これらの問題への対策として、様々なモデルや手法の開発を行なった。これらの方法については後に簡単に説明する。

　数値シミュレーションに用いている現象的物理モデルは相間相互作用である。相間相互作用には相の分類が深く関係している。従来の蒸気爆発素過程の計算では、3〜4相を仮定したものであった。すなわち低温液と高温液(または高温の細粒)と蒸気である。しかし、粗混合過程において主要なモードは高温液滴とそれを取り巻く蒸気膜である。相分布は高温液の周囲に集中している。この場合には、蒸気膜の伝熱方式や高温液の抵抗力などが重要な要素になり、体系全体の相分布を支配する。しかし、膜表面から脱走した気泡は、蒸気膜における伝熱・流動方式とは全く異なるものである。よって、同じ蒸気ではあるが、それらは別々のものとして取り扱うことが重要である。従って、本研究では、蒸気膜と気泡を別にして、5〜6相からなる体系を考えた。

　伝熱と流動モデルについては、本研究で主要な相間熱伝導、質量輸送のモデルは他のコードで採用されているものと殆んど同じであるが、抵抗力については蒸発の影響を良く考慮した。これは、高速伝熱の場合には、熱源として働く高温粒子の抵抗力は常温におけるそれと全く異なることを考慮したものである。常温においては粒子と周囲の低温液の速度差が主要な要素と考えられるが、高温蒸発状態では、粒子周辺の不均一な蒸発が大きな抵抗力を生むと考えられるからである。

　従って本研究では、蒸気膜の質量保存式と運動量保存式から理論的に2つの蒸発抵抗力モデルを導出した。これらは高温液滴が自由表面を通る時と低温液体の内部を移動する時の抵抗力に対してのものである。結果では、極高温の場合(Tm>2500K)には、この抵抗力はかなり大きいことが明らかになった。これは有名なMIXA粗混合実験の結論と一致する(図4,5)。

　本研究は微粒化後の過程における伝熱モデルも検討した。ここでは、液体の過熱により蒸気のフラッシュが生じる現象も考慮している。更に、少量分散の高温液滴間の微粒化伝播のシミュレーションモデルも開発した。

　モデルの構築も困難であるが、計算の時に最も困難なのは数値計算上の問題である。蒸気爆発のような変化が激しい現象では、数値発散はよくある問題である。よって本研究では、これらの対策として、3つの数値アルゴリズムの開発を行った。これらは分散流モデルとマルチ領域スキーム、それに自由表面処理と呼ぶものである。

　中小規模の蒸気爆発では通常高温相の体積率が非常に小さいが、相変化に大きな影響を与えるのもこの相である。したがって、この相の体積率は、数値的な小さな振動でも大きな不安定を生じやすい。しかし、多相流数値計算で体積率が小さい相を扱う際にはいくつかの問題がある。その一つは、分散する相の適切な数学的表現ができないことである。多相流は数学的には、分散流についても、連続相として表現する。このため実際には存在しない空間勾配による拡散が分散相に対しても評価されてしまう。特に高温相の擬拡散は、周囲に擬蒸気発生を引き起こしてしまう。これに対して、分散流モデルでは、分散相の基礎保存式の物理拡散項を省略し、体積率が小さいところで相間熱・運動量輸送を無視することとした。

　もう一つの問題は圧力の計算における問題である。CHAMP/VEで使う数値解法はSIMPLE解法から発展させたもので、圧力補正式を用いて速度の修正を繰り返す。通常、圧力補正式は全ての相の和に基づいて得られるが、例えば蒸気膜などの分散相には大きな数値振動があるので、圧力の計算がよく発散を引き起こす。このため、ここでは連続相だけに基づいて圧力補正式を立てる計算スキームとした。

　更に、体積率の計算における問題もある。体積率は通常、まず各相の質量保存式から求めて、全体の体積率の和が1になるという条件で個々の相あるいはある相だけの修正をする。しかしこの修正は、体積率が小さい相については誤差が大きくなり易い。これに対して、分散流モデルでは、体積率の修正には高温液滴の誤差が起きないように改良した。

　マルチ領域スキームは分布の不連続界面の数値拡散によって数値上の擬拡散を起こさないために開発された解析手法である。蒸気爆発過程の体系は相分布が位置的に集中しているため、混合領域と水だけの単相の領域の2つの領域に分けられる。従来の方法では、系に含まれるどの成分も体系全体にわたって微小量でも存在するような分布を想定し、全領域で計算していた。そうすると、数値拡散は上で述べた2つの質量の境界で大きくなる。これに対して、マルチ領域スキームでは、この2つの領域を分離して別々に計算する。混合領域では全部の相を計算するが、単相領域では水の単相のみについて計算を行う。この2つの領域の界面は時間とともに移動していくので、時間ステップを進めていく際に界面の位置を知る必要がある。そこで、領域フロントのトラッキングという手法を開発した。この手法は各時間ステップで混合領域の端に位置する擬似粒子の軌跡をトラッキングしていくものである(図2)。この方法では、移動する境界格子系を定義し、そこで、質量・運動量・エネルギーを保存することが必要である(図2)。

　自由表面の問題も非常に顕著である。数値拡散は当然な問題であるが、表面での圧力差が大きく、数値不安定を生じ易いというのはますます大きな問題となる。この問題は、液体内部の蒸気膨張が大きく、内部の圧力が大きくなり、表面で圧力勾配が大きくなるためのものである。液体速度は、圧力勾配に支配されるので、大き過ぎる圧力勾配は大きな速度を引き起こす。そして、計算上表面のセルでは質量の流失が生じることになる。その結果、表面内部の圧力が大幅に下降する。こうして数値的な振動が発生するが、この振動が発散してしまうことが頻繁に生じる。

　更に、蒸発による抵抗力は表面において最も大きいので、表面の位置が確定されなければならない。そのために、表面の処理手法を開発した。まず、表面の位置をトラッキングしていくものである。或は、表面セルを定義し、表面セルの大きさを計算するものである(図3)。このような取り扱いにより、数値拡散の防止が可能である。更に、大きい圧力勾配により生じた数値不安定問題の解決のためには、表面セルの計算を行なう際に低い緩和係数を用いたり、速度を制限するなどの手法を用いた。

　これらのアルゴリズムをコードに組み込み、2つの数値的な問題が有効に抑えられることを確認した(図6,7)。

　前述のように、連続な相における圧力計算や、体積率の小さい相の体積率修正誤差の防止は、計算技術の改良に含まれる。他にも、緩和係数の値の選択も、これに含まれる。また、全ての重要な相と重要な変数の関係式を、異なる領域で別々に取り扱うようにした。更に、計算中の収束傾向により自動的または人工的に改変することができるようにした。その結果、緩和係数は低い値であることが重要であることが示せた。

　コードのその他の改良点は、入力メニューの作成と出力データの図示の自動化と、円柱座標の導入である。円柱座標により、対象的な3次元問題もシミュレーションできるようになった。これはもとのCHAMPAGNEコードで用いている直角座標よりも応用範囲を拡大するものである。

　計算は、粗混合過程を中心に行なった。使用したデータはMIXAの実験データである。計算結果と、この実験の結果との比較を行なうことができた。この他には、微粒化の計算を例として実行した。微粒化については2つのケースを計算した。すなわち、全ての高温相を連続分裂するものとして計算した場合と、逐次分裂するものとして計算した場合の2つである。図8と図9にその計算例を示した。

　蒸気爆発素過程の数値シミュレーションのために、多相多次元コードCHAMP/VEを開発した。物理モデルの不足と数値計算上の問題から、素過程中最も簡単な粗混合と微粒化の2つの過程だけのシミュレーションを行なった。

　物理モデルについては蒸発抵抗力モデルを開発した。

　数値計算手法については分散流モデル・マルチ領域スキーム・自由表面処理、以上3種類の手法を開発した。

　更に、計算技術にもいくつかの改良をした。

　これらのモデルや手法を用いて粗混合と2種類の微粒化のケースをシミュレーションし、その結果をMIXA実験の結果と比較することができた。

Fig.1.A sketch of a vapor explosion system

Fig.2.Mixture region boundary tracking technique

Fig.3.Free surface cell treatment

Fig.4.Comparison of evaporation drag model and a general drag model for a partile penetrating free surface

Fig.5.Comparison of evaporation drag model and a general drag model for a particle moving inside liquid

Fig.6.Void fraction of melt by Multi-Region Method(MRS)and Single Regiou Method(SRM)

Fig.7.Void fraction of melt by dispersed distribution method(DPM)and Continuous Fluid Model(CFM)

Fig.8.Evolution of the volume fraction of liquid coolant in a mixing process

Fig.9.Evolution of the volume fraction of melt drops in a fragmentation process