1.緒言

　二相流や沸騰現象の数値解析は、原子炉の安全にとって非常に重要である。しかし、二相流では界面の形状が複雑であり、しかも急激な変化を伴う。また、大気圧の水のように、液相と気相の間の密度差が大きい場合、界面で数値不安定性が起きやすい。そのため、二相流や沸騰現象の数値解析は現状では困難である。

　MPS-MAFL法は移動境界を持つナビエ・ストークス方程式を解くために開発された数値解析方法である。これは、ラグランジュ計算のための粒子法であるMPS法(Moving Particle Semi-implicit method)と、グリッドレス対流計算法であるMAFL法(Meshless Advection using Flow-directional Local-grid method)で構成され、任意ラグランジュ・オイラー(ALE, Arbitrary Lagrangian-Eulerian)の計算が可能である。

　そこで本研究では(1)MPS-MAFL法を用いて急激な過渡条件での核沸騰の解析を行い、(2)相変化を伴い、密度差が大きい二相流が扱える解析コードを開発し、(3)それを用いて、大気圧での水の膜沸騰の解析を行う。

2.数値解析方法

2.1支配方程式

　連続の式、ナビア・ストークス方程式、エネルギー保存式を式(1)〜(3)に示す。

2.2計算アルゴリズム

　MPS-MAFL法はLagrangian phase, re-configuration phase, convection (Eulerian) phaseで構成される。まずLagrangian phaseでは、圧力項を除いたナビエ・ストークス方程式を陽的に解き、仮の速度と仮の計算点位置を得る。その後,ポアソン方程式を用いて圧力場を陰的に解き、ここで得られた圧力場により、仮の速度と位置が更新され、ラグランジュ計算フェーズでの速度と位置が求められる。また、エネルギー保存式を陽的に解き、計算点の温度が得られる。各支配方程式での微分演算子は、後で説明される粒子間相互作用モデルを通じて計算される。

　なお、re・configu岨tionphaso(計算点再配置フェーズ)では、計算点の再配置が行われる。この時、MPS-MAFL法は非均一粒子配置が可能であるため、計算点を自由に配置できる。本研究では、計算の精度を上げるために、界面の形と温度分布を考慮し、計算点の再配置を行うことにした。計算点が再配置されると、計算点の速度Ucと任意対流速度Uaが決められる。各速度と各位置の関係を図1に示す。

　Convection phase(対流フェーズ)では、再配置された計算点での物理量を決めるための対流計算を行う。本研究では、Yoonらによって開発されたMAFL法を用いた。

2.3 Lagrangian Phase

　前述したようにMPS法では粒子間相互使用モデルを通じて、支配方程式での微分演算子を離散化する。このモデルで、各粒子は重み関数の与える重みで周辺の粒子と相互作用する。本研究では式(5)の重み関数が用いられている。

　重み関数の半径reを有限とすることによって、計算を行う粒子が近接の粒子に限定されるので、大規模な計算の時でも計算時間を短縮することができる。また、MAFLのように計算点を非均一的に配置する手法では、有効半径内に同じ数の計算点が入るように、計算点ごとに異なる有効半径を使う。

　次に各微分演算子の離散化について説明する。まず、ナビエ・ストークス方程式の中の圧力項での勾配演算子は、各粒子の物理量と重み関数によって式(6)で表される。

　ここで、規格化係数は〓で定義される。有効半径はre,ij=(re,i+re,j)/2であり、各計算点の有効半径が異なる場合を考慮したものである。なお、dは空間の次元数である。

　発散演算子も勾配演算子のように離散化される。ベクトルψの発散は、ψj-ψiとrj-riの内積を用いて式(8)のように書ける。

　拡散はある粒子iの物理量φiを近接の粒子jに分配することで表す。

　中心極限定理により、上記の式が成立するためには、時間△t秒後におけるある粒子から近接粒子への物理量の分配の総計が2dv/△tと一致する必要があり、λは〓となる。

2.4 Re-Configuration Phase

　このフェーズでは、気体-液体境界面や境界条件などを考慮して、計算点を再配置する。ここで、二相流の境界面に常に計算点を配置するため、気泡表面の運動を直接追跡できる。また、計算点を非均質に配置できるため、温度境界層内に計算点を集中することにより、精度の高い結果が得られる。計算点の具体的な配置方法は、問題によって異なるため、各解析の条件とともに説明する。

2.5 Convection Phase

　流体の物理量は、流線に沿って対流するので、流線に沿った計算格子を作ると、多次元の対流問題も一次元計算になる。MAFL法では、計算点ごとに局所的な1次元格子を、その計算点の持つ速度ベクトルの方向に生成する。そして、仮の計算点をその格子の上に等間隔に配置し、そこでの変数値を周囲の粒子の持っている変数値より内挿する。なお、その仮の計算点での値を用いて、1次元差分スキームを使って移流計算を行う。本研究では式(11)のupwind schemeを使った。

　ここでq=|ua|△t/△rであり、△rは仮の計算点間の距離である。MAFL法の概念図を図2に示す。

2.6非圧縮性モデル

　連続の式に圧縮性の項を入れることにより、擬似的に圧縮性の効果も扱うことが可能である。ナビエ・ストークス方程式の陰的な圧力勾配項に、等エントロピー変化の式と、連続の式を代入すると、圧力のポアソン方程式が得られる。本研究では、液相では非圧縮性流体の計算式(12)を、気相では圧縮性流体の計算式(13)を使用する。

2.7表面張力モデル

　水-蒸気境界面にある粒子にかかる体積力Fsは以下の式で与えられる。

　ここで、σは表面張力係数,κは曲率,nは単位垂直ベクトルである。表面の曲率半径は三つの点を過ぎる円の半径であり、単位垂直ベクトルは表面から中心へ向かうのベクトルの単位ベクトルである。

3.過渡状態での核沸騰の解析

3.1反応度事故と過渡沸騰

　軽水炉の炉心で制御棒が何らかの理由で引き出されると、反応度が投入され、炉心出力が上がり、燃料集合体の温度も上昇する。これを反応度事故という。このような急激な出力上昇時には、燃料の温度上昇によるドップラー効果と、炉心冷却水のボイド発生により、負の反応度が投入され、炉心出力は低下する。しかし、従来の相関式では過渡沸騰時のボイド率の予測が困難であり、高熱流束および高サブクール条件での過渡沸騰現象についても解明されていない。そのため現在の原子炉の安全解析では、ボイドの効果が無視されている。そこで、MPS-MAFL法を用いて、高熱流束および高サブクール条件での過渡沸騰時における気泡の発生と成長過程の2次元数値解析を行うことにした。

3.2解析条件

　本解析での解析条件を表1に示す。この解析条件は低温反応度事故を仮定した山田らの実験を元にして決めたものである。また、初期状態での気泡の形は円であり、気泡の表面と加熱面の間の接触角は45度である。薄い温度境界層を精度高く解析するため、水・蒸気境界面の形状と温度境界層と考慮して計算体系を四つの領域に分け、気泡周辺に計算点を多数配置する。

　気泡の中の物理量は均一であると仮定し、蒸気中には計算点が配置されてない。また、蒸発と凝縮による蒸気量の変化と界面の運動による体積の変化により、蒸気の圧力が更新される。蒸気の温度は飽和温度である。なお、界面上の計算点の温度と圧力は蒸気と同じである。ただし、圧力では表面張力の影響は考慮する。熱伝導計算の後、界面上の計算点の温度が飽和温度と異なる場合には、その分のエネルギーは全て水の蒸発または蒸気の凝縮に使われることにする。ただし、温度境界層がまだ発達していない初期状態では、界面の運動と凝縮量は計算しない。底面と側面ではnon-slip境界条件を使い、計算体系の下にある加熱面では一定熱流束になるように温度分布を与える。

3.3解析結果

　初期半径が50x10-6mである場合、気泡の形と半径の変化を図3に示す。気泡は0.0162秒から成長を始め、0.1363秒に加熱面を離脱する。その後、加熱面からの熱流入の減少と高いサブクール度による熱流出の増加のため、急激に凝縮する。

　気泡体積の変化(図4)を見ると、初期気泡半径が小さいほど早く成長が始まることが分かる。ただし、気泡成長開始時間は初期半径に対してそれほど大きな依存性はない。定常沸騰では、初期気泡半径が大きいほど低い熱流束で成長が始まることが知られており、これとは傾向が異なっている。この理由は、過渡沸騰の場合では、加熱面付近の温度が局所的寺こ高くなるため、小さい気泡核が先に成長開始の条件を満たすためである。また、気泡離脱時の体積は、初期気泡半径に依存せず、ほぼ同じである。

　なお、山田らの実験ではボイド率の上昇は0.08〜0.1secに観測されており、本解析における気泡の成長のタイミングとよく一致している(図5)。発泡点密度を4x105/m2にすると定量的にボイド率が一致するが、現状では過渡沸騰時の発泡点密度を予測することが難しい。

4.大気圧での水の膜沸騰の解析

4.1膜沸騰現象

　原子炉では、事故時に炉心出力の増加や炉心冷却の悪化により、炉心内で膜沸騰が生じる可能性がある。また、膜沸騰現象を理解し、さまざまな条件における膜沸騰時の熱流束などを予測することは、原子炉の安全にとって重要である。しかし、膜沸騰では核沸騰と異なり、蒸気膜内の伝熱と流動を解析する必要がある。そこで、密度差が大きい二相流が扱えるMPS-MAFLコードを開発して、大気圧での水の膜沸騰の数値シミュレーションを行うことにした。

4.2解析条件

　膜沸騰における数値シミュレーションの計算条件を表2に示す。本計算では、体系最上部の圧力を1気圧で維持する。また、流体の下にある加熱面の温度は200〜400℃であり、境界面の初期位置は式(15)とする。

ここで、ycは境界面平均位置、εは変位である。

4.3シミュレーション結果

　計算結果のうち、加熱面温度の500度の場合の蒸気膜内の温度変化を図6に示す。図6(a)は流体の初期状態で、上から水、蒸気膜、加熱面となっている。加熱初期では、流体の対流がないため熱は伝導によって伝達され、加熱面付近から熱くなる。その後、蒸発により蒸気量が増加すると蒸気膜の平均厚さも増加するが、不安定性により蒸気膜の薄い部分が生じる。そして、この蒸気膜の薄い部分での蒸発量が増加し、蒸気の対流が生じる。また、蒸気の対流が発達すると、蒸気膜上部にあった低温の蒸気が加熱面に供給され、この蒸気の循環によって加熱面付近の温度が低くなるため、熱伝達が増加する。

　加熱面での平均熱伝達率の変化を図7に示す。また、各条件での平均熱流束を表3にまとめた。MPS-MAFL計算結果とBerensonの式での熱流束を比べて見ると、本計算の値が低く、Berensonの式の約40〜60%であることがわかる。また、平均蒸気膜厚さと最小計算点間距離が小さい場合の熱流束が高かった。本計算では、蒸気膜がBerensonの式で仮定されているよりも厚いため、熱流束が低いと思われる。

4.緒言

　任意ラグランジュ・オイラ計算が可能なMPS-MAFL法を用いて、過渡沸騰における気泡成長の数値解析をおこなった。高熱流束および高サブクール条件では、小さい気泡核が大きい気泡核よりも先に成長するのが分かった。これは、過渡沸騰の実験において、小さな気泡が大量に発生することが観察されていることを説明できる。また、気泡成長開始時間は山田らの実験結果と良く一致した。

　次に、MPS-MAFL法を用いて、密度差が大きい二相流が扱える膜沸騰解析コードを開発した。さらにこれを用いて、大気圧での水の膜沸騰の解析を行った。膜沸騰における蒸気膜内の流動を計算することができ、膜沸騰における熱伝達の過程を再現することができた。しかし、計算された熱流束を実験式であるBerenson式と比較したが、定量的な一致はまだ不十分だった。

表1核沸騰解析条件

表2膜沸騰解析条件

表3膜沸騰解析での熱流束

図1各位置と速度の関係

図2MAFL法

図3核沸騰解析結果(rinit=50μm)

図4初期気泡半径の影響

図5山田らの実験結果との比較

図6膜沸騰解析結果

図7熱伝達率

　本論文はメッシュレス手法であるMPS-MAFL法による核沸騰と膜沸騰の数値シミュレーションについて記述したもので論文は5章より構成されている。

　第1章は序論で、沸騰様式として核沸騰および膜沸騰があり、原子炉の除熱などにおいて工学的に幅広く利用されている伝熱形態であると述べられている。こうした沸騰を伴う気液二相流の解析には、従来より集中定数モデルが用いられてきたが、最近は実験相関式を必要とせず基礎方程式のみに基づいた機構論的な解析が試みられている。そして、気液界面の大変形に適用できるメッシュレス法がその有力な手法であるとしている。論文提出者が所属する研究室でこれまで研究されてきたメッシュレス法としてMPS法とMPS-MAFL法を紹介し、本研究の目的はMPS-MAFL法を用いて沸騰の機構論的な解析をおこなうことであるとしている。

　第2章では本研究で用いられている数値解析手法についてまとめている。任意ラグランジュ・オイラー記述の熱流動の支配方程式、MPS法で用いる粒子計算モデル、および粒子の不均一配置に拡張したモデル、計算メッシュを用いない対流項離散化スキームであるMAFL法、その他沸騰の解析に必要な表面張力、相変化、密度比が大きい場合の圧力計算のアルゴリズム等である。このうち、膜沸騰の解析コードは論文提出者が新たに開発しており、その中で用いられている相変化モデルは独自に考案したとしている。

　第3章は原子炉の反応度事故(RIA)などに見られる過渡時の核沸騰のMPS-MAFL法での解析について述べている。RIAでは急速な出力上昇による沸騰により負の反応度が投入されるが、従来の定常沸騰での相関式では評価することができず、高熱流束および高サブクール条件での過渡沸騰現象については解明されていないとしている。なお、核沸騰における単一気泡の成長と離脱は、YoonがMPS-MAFL法を用いて既に定量的解析に成功していることを述べている。この解析では蒸気内の温度分布を一様とし、蒸気流動の計算も必要ない。この成果を踏まえて、実験の条件での過渡沸騰の解析をおこなっている。熱流束が限界熱流束を超えているにもかかわらず、過渡的に気泡の成長と離脱が生じることが計算され、気泡の成長にかかる時間が模擬実験の結果と一致していると述べている。さらに、初期気泡核半径の影響を調べ、気泡の成長にかかる時間は初期気泡核半径にほとんど依存しないという結果を得ている。これは、過渡沸騰実験では多数の細かい気泡が発生することをよく説明するとしている。

　第4章は大気圧下における水の膜沸騰を解析した成果である。従来の研究として、座標変換を用いたSon-Dhirの研究、および界面追跡法を用いたJuric-Tryggvasonの研究を紹介し、大気圧下の水の膜沸騰は、気液の密度比が大きいだけでなく蒸気膜が非常に薄くなるため非常に難しい問題であり、まだシミュレーションされたことが無いとしている。膜沸騰では、蒸気膜の中の熱流動により沸騰熱伝達が支配されているため、新たに、蒸気の熱流動も解析する必要があり、新しい計算コードが開発された。大気圧下の水と蒸気では密度比が大きいため、それぞれの相で別々に圧力場を解くアルゴリズムが導入されている。初期の蒸気膜に擾乱を与え、計算をおこなったところ、安定した蒸気膜が維持されることが示された。熱伝達量をBerensonの実験相関式と比較したところ、40-60%程度であり、過小評価となっている。

　第5章は結論であり、メッシュレス法であるMPS-MAFL法を用いて、核沸騰と膜沸騰の機構論的な数値シミュレーションをおこなったとしている。核沸騰ではRIAの条件で気泡の成長時間が定量的に予測できるとともに、小さな気泡が多数見られる理由を説明できたことが成果であると述べられている。膜沸騰では、大気圧下の水-蒸気系で安定した蒸気膜を計算することが初めてできたが、定量的な伝熱量については過小評価であることが述べられている。

　以上を要するに本論文はMPS-MAFL法を用いて沸騰の数値シミュレーションをおこない、過渡状態の核沸騰解析では気泡成長開始時間について実験と一致する結果を得るとともに、膜沸騰解析では実際の大気圧下の水の膜沸騰を初めて解析している。この成果はシステム量子工学、特に原子力熱流動工学に進展をもたらすのみならず、工学諸分野における沸騰伝熱の数値シミュレーションの進歩に貢献するところが少なくない。よって本論文は博士(工学)の学位請求論文として合格と認められる。