Moment–of–Fluid 法による沸騰熱伝達現象の三次元詳細数値解析
徳島大学大学院

社会産業理工学研究部

機械科学系教授

太田光浩

1．緒言
高温鋼材を水によって冷却する場合，沸騰冷却によって鋼材が冷却される．高温鋼材の冷却開始時は膜沸騰であり，その後，冷却により鋼材の温度が下がり始めると膜沸騰から遷移沸騰を経て核沸騰に変化する．高温鋼材の冷却過程を深く理解するためには沸騰熱伝達の詳細な知見が必要となる．沸騰に関しては，古くから様々な研究が行われており，これまでに多くの実験研究が報告されている．沸騰現象は時空間スケールが極めて小さい中での非定常現象であるため，実験的に熱伝達機構の詳細を把握するには限界がある．したがって，数値解析(CFD)により現象を的確に再現し, 局所熱流動場を詳細に把握した上で，異相間を介しての熱伝達機構を評価することは有効な手法となる．しかし，沸騰現象は短時間で熱移動を伴いながら気泡が生成・成長する複雑熱流動系となるために CFD が非常に難しい．このような背景の下，相変化を伴う気液二相流の挙動を高精度に再現し，沸騰現象による熱伝達現象を精度良く評価できる数値解析方法の確立が重要な課題となっている．
沸騰の CFD は，近年活発化しており，数値解析による研究報告が増えている[1–9]．
これらの研究では，｢界面でのみ相変化が生じる｣，｢界面温度は飽和温度に保たれる｣と言う仮定の下で熱流動の支配方程式が数値的に解かれた．すなわち，相変化を伴う沸騰現象での気液界面において満足するべき質量，運動量，エネルギーの不連続条件を考慮しながら，質量，運動量，エネルギーの保存式が液相および蒸気相に対して数値的に解かれてきた．沸騰では気泡成長を伴うために，界面の時間発展過程を数値的に 追 跡 す る 必 要 が あ り ， 先 の 研 究 で は ， Volume–of–Fluid(VOF) 法 [10] ， Front–Tracking(FT)法[11]，Level–Set(LS)法[12]，二相系格子ボルツマン(LB)法[13]などが主として用いられた．各方法とも長短所を併せ持つが，本研究では Moment–of–Fluid(MOF)法[14, 15]を用いて界面の時間発展を捕獲する．MOF 法は VOF 法の一種であるが，計算格子内の気液界面勾配(区分的界面近似直線)を決定する際の手順に独自性がある．界面近似直線は，計算セル内の流体体積を厳密に保存しつつ，近似直線で構築された流体の重心位置が計算セル内の実際の流体の重心位置とできる限り一致するように決定される．したがって，MOF 法では計算格子 1 セルの物理情報から，界面近似直線を二次精度で決定できる利点を有し，多流体・多相流れにも容易に拡張できる．沸騰のように加熱壁に気液界面が接触しながら時間発展する場合，計算格子 1 セルの物理情報から界面近似直線を決定できることは計算を進める上で大きな長所となる．また，加熱壁が傾斜して配置されているような場合も計算格子内の気液固界面をモデル化しやすい．

－ 105 －

本研究では，MOF 法を用いて加熱壁面から核沸騰による気泡の生成・離脱過程の数値解析を行い，複雑な沸騰現象を三次元数値解析により再現する．
2．数値解析
支配方程式である連続の式，運動方程式およびエネルギー保存式は，一流体近似に基づいて次式で表される．
(1)

u  0 u u  u   p    2  D  F   g t T u  T     kT tC p 



(2)
(3)

ここで，u：速度 [m/s]，p：圧力 [N/m2]，D：変形速度テンソル [1/s]，ρ：密度 [kg/m3]，

μ：粘度 [Pa･s]，σ：界面張力 [N/m], κ：曲率 [1/m]，F：VOF 関数 [–]，g：重力 [m/s2]，
Cp：定圧比熱 [J/(kg⋅K)]，T：温度 [K]，k：熱伝導率 [W/(m⋅K)] である．
また，密度(ρ)，粘度(μ)，比熱(Cp)，熱伝導率(k)は次式で定義される．

   V 1 F    L F   , , C , k p(4)

添え字は，V：蒸気相，L：液相 を表す．VOF 関数の移流方程式は次式で定義される．

F u  F  0t





(5)

気液界面(Γ )において満足するべき質量，運動量，エネルギーの不連続条件は次式となる．なお，TInt(界面温度) = TSat(飽和温度) とする．
(6)
(7)
(8)
上式において，[ ]Γ は[f ] Γ = f V – f L を表す演算子である． ：沸騰質量流束[kg/(m2⋅s)]，n：界面の法線ベクトル， un n ：界面での法線方向速度の微分，LV：潜熱 [J/kg]である．また，界面において VOF 関数の移流方程式は次式を満足する．

F
 uInt  F  0t(9)uInt は界面速度[m/s]であり，次式で表される．

－ 106 －

(10)以上の支配方程式を二次精度のアルゴリズム[16, 17]を用いて解析する．
沸騰の数値解析を実施する上で固体壁と気泡の接触領域の間に形成されるマイクロレイヤーの取り扱いは非常に重要となる．実験的にマイクロレイヤーの存在は確認され，気泡成長速度に大きな影響を与えることが報告されている[18]．マイクロレイヤ
ーは，発泡点からの距離が増大するにつれて大きくなることが知られている．本研究では，計算時間の短縮のために境界適合細分化格子(AMR)[19]を適用し，気液界面や固気界面付近では計算格子をより細かく動的に設定して解析を行った．マイクロレイ
ヤーも高次の AMR を適用し，計算格子をマイクロレイヤーに割り当てることでマイクロレイヤー内の熱流動場を厳密に解析することを試みたが，マイクロレイヤーの厚さが非常に薄いために現実的には適用は難しかった．そこで，本研究では，Vahab らのマイクロレイヤーモデル[20]を適用して数値解析を行った．
3．結果と考察
3.1 数値解析法の妥当性
本研究で用いる数値解析手法の有効性を検証するために Dhir[21]が行った水の沸騰実験と数値解析との比較を行った．Dhir は，加熱壁温度(TW)を 380 K(過熱度 ΔT = 7 K)にして単一蒸気泡の生成過程を実験的に観察し，LS 法を用いた数値解析を行った．本
研究での数値解析でも Dhir の実験条件に基づいて物性・物理条件を設定した．すなわ
ち，加熱壁と気泡の接触角(θ)を θ= 50°，TW = 380 K，TS(飽和温度) = 373 K とし，t =0 s において D0(初期気泡核直径) = 0.4 mmを設定した．なお，t = 0 s において加熱壁面付近において，ある程度，温度境界層が発達していることとし，温度境界層厚さ(δT)
は Dhir[21]が使用した相関式を用いた．
数値解析における解析系を図 1 に示す．
数値解析は，二次元軸対称解析とし解析系の大きさは R(半径方向長さ) = 5 mm，Z(解析系高さ) = 4R= 20 mm と設定した．また，液相は H(液相高さ) = 2R = 10 mm まで満たし，解析系の上部領域は気相とした．なお，計算格子の最小メッシュサイズ(Δx)は，解析結果のメッシュ依存性を検討した上で，

図1

二次元数値解析解析系

Δx ≈ 0.02 mm とした．
核沸騰過程の数値解析結果を図 2 に示す．解析結果では気泡核を 55 ms 毎に設定し－ 107 －図2

二次元数値解析結果ており，加熱開始の 0 ms から 4 個目の気泡が生成する 190 ms までの気泡形状と温度分布である．温度分布は，青色(373 K)から赤色(380 K)の範囲でカラーコンターで表わされている．図より明らかなように核沸騰による気泡生成の様子が再現されている．
最初の気泡が離脱した後も，気泡が安定して連続的に生成していることが分かる．気泡成長に関して定量的に考察するために，図 3 に気泡直径の時間変化を示す．1 つ目の気泡が離脱するまでの時間(tD)と離脱後の気泡径(D)は，(1)tD = 52 ms，D = 2.93 mmであった．2 つ目から 4 つ目までに関しては，(2)tD = 46 ms，D = 2.97 mm，(3)tD = 46ms，D = 2.98 mm，(4)tD = 47 ms，D = 3.01 mm となった．図 3 から明らかなように，2 つ目以降の気泡について着目すると，気泡径に大きな変化はないが気泡離脱までの時間が短くなっていることが分かる．つまり，気泡成長速度は 1 つ目の気泡よりも 2 つ目以降の気泡の方が大きいことを意味し，2 個目以降の気泡成長速度はほぼ一定であると言える．これは，気泡核周りの熱流動場の発達と関係している．図から明らかなように 2 つ目以降の気泡核は，気泡周囲が 1 つ目の気泡核よりも高温領域とな－ 108 －っているためであり，1 つ目の気泡が生成した後に熱流動場の状態が定常状態に達していると言える．

図3気泡成長速度

図 4 に Dhir の実験・数値解析との比較を示す．Dhir は 3 回の実験を行っており，

図4

Dhir の実験・数値解析との比較

図中の黒シンボルは 3 回の実験結果に対応している．実線は数値解析結果であり，紫の実線が Dhir の数値解析結果，赤(1 つ目の気泡)と青(4 つ目の気泡)の実線が本研究での数値解析結果である．本研究の数値解析結果は，Dhir の実験結果および数値解析結果と良好な一致を示していることが分かる．特に 4 つ目の気泡の成長速度は，時間初－ 109 －期を除いて実験結果と良く一致する．数値解析結果と実験結果の時間初期における相違は，実際のマイクロレイヤーの大きさとその時間的発達に起因していると考えられるが，全体として Vahab らのマイクロレイヤーモデル[20]を適用した MOF 法による数値解析方法は核沸騰問題の計算に有効であることを確認できる．
3.2 三次元数値解析
解析方法の有効性が確認できたため，三次元数値解析により複雑な沸騰現象の再現を試みる．解析系は直交座標(x, y, z)系とし，X(x 方向長さ) = Y(y 方向長さ) = 10 mm，Z(z 方向長さ) = 20 mm とし，液相を H(液相高さ) = 10 mm まで満たした．前節の二次元数値解析と同様に θ = 50°，TW = 380 K，TS = 373 K とし，D0 = 0.6 mm の気泡
核を 2 ms 毎に加熱面上にランダムに配置して解析を行った．計算格子の最小メッシュ
サイズは，計算時間を考慮して Δx ≈ 0.156 mm とした．
図 5 に三次元数値解析から得られた t = 175 ms までの気泡挙動の鳥瞰図を示す．図において加熱壁面は赤色で示されている．解析結果から明らかなように，非常に短い時間の間に複数個の気泡が成長しながら合一し，合一後は大きく変形しながら上昇しているのが分かる．上昇した気泡は上部の気相中へと抜けるため液層の気液界面もダイナミックな挙動を示す．

図5

三次元数値解析結果(気泡挙動)－ 110 －

図 6 に x–z 中央断面での温度分布を示す．温度分布に関しては t = 100 ~ 175 ms での結果を示す．気泡の上昇運動に伴い気泡背後の水は上方に運ばれ，その流れの影響を受けて高温度域が上方に引き上げられている領域を確認できる．しかし，ここで再現された解析結果は，短い現象時間での沸騰現象であるために，系内の温度場は十分に発達していない．より長時間の沸騰現象を再現すれば，より発達した温度場が観察できると考えられる．引き続き，より長い時間の沸騰現象を再現し，複雑な熱流動場の形成過程について考察する予定である．

図6温度分布(x–z 中央断面)

4. 結論
MOF 法を用いた数値解析により加熱壁面から核沸騰による気泡の生成・離脱過程を再現可能である．今後，より長い時間の複雑な沸騰現象を三次元数値解析により再現し，沸騰熱伝達の詳細メカニズムを検討する予定である．
謝辞
本研究は，公益財団法人 JFE21 世紀財団の 2017 年度技術研究助成による支援を受けて実施した．ここに深く感謝の意を表する．
参考文献[1] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and Heat Transfer AssociatedWith a Single Bubble During Nucleate Boiling on a Horizontal Surface, J. HeatTransfer, 121, 623(1999).
[2] G. Son, N. Ramanujapu, V.K. Dhir, Numerical Simulation of Bubble Merger
Process on a Single Nucleation Site During Pool Nucleate Boiling, J. Heat－ 111 －Transfer, 124, 51(2002).
[3] A. Mukherjee, V.K. Dhir, Study of Lateral Merger of Vapor Bubbles During
Nucleate Pool Boiling, J. Heat Transfer, 126, 1023(2004).
[4] D. Li, V.K. Dhir, Numerical Study of Single Bubble Dynamics During Flow
Boiling, J. Heat Transfer, 129, 864(2007).
[5] 梅村悠, 姫野武洋, 渡辺紀徳, 加熱面上における気泡成長・離脱過程の数値解析,

日本マイクログラビティ応用学会誌, 29, 2(2012).
[6] G. Tryggvason, J. Lu, Direct Numerical Simulations of Flows with Phase
Change, Procedia IUTAM, 15, 2(2015).
[7] A. Urbano, S. Tanguy, G. Huber, C. Colin, Direct Numerical Simulation ofNucleate Boiling in Micro-Layer Regime, J. Heat Mass Transfer, 123,1128(2018).
[8] Y. Sato, B. Niceno, Pool Boiling Simulation Using an Interface Tracking
Method: From Nucleate Boiling to Film Boiling Regime Through Critical Heat
Flux, Int. J. Heat Mass Transfer, 125, 876(2018).
[9] T. Sun, A Numerical Study on Dynamics Behaviors of Multi Bubbles Merger
during Nucleate Boiling by Lattice Boltzmann Method, Int. J. Multiphase

Flow, 118, 128(2019).
[10] C.W. Hirt, B.D. Nichols, Volume of fluid(VOF)method for the dynamics of freeboundaries, J. Comput. Phys., 39, 201(1981).
[11] S.O.

Unverdi,G.Tryggvason,

A

Front-Tracking

MethodforViscous,Incompressible, Multi-Fluid Flows, J. Comput. Phys., 100, 25(1992).
[12] M. Sussman, P. Smereka, S. Osher, A Level Set Approach for Computing
Solutions to Incompressible Two-Phase Flow, J. Comput. Phys., 114, 146(1994).
[13] Z. Dong, W. Li, Y. Song, A Numerical Investigation of Bubble Growth on andDeparture from a Superheated Wall by Lattice Boltzmann Method, Int. J. Heat

Mass Transfer, 53, 4908(2010).
[14] H. Ahn, M. Shashkov, M. Christon, The Moment-of-Fluid Method in Action,Commun. Numer. Methods Eng., 25, 1009(2009).
[15] M. Jemison, E. Loch, M. Sussman, M. Shashkov, M. Arienti, M. Ohta, Y. Wang,
A Coupled Level Set-Moment of Fluid Method for Incompressible Two-Phase
Flows, J. Sci. Comput., 54, 454(2013).
[16] G. Li, Y. Lian, Y. Guo, M. Jemison, M. Sussman, T. Helms, M. Arienti,
Incompressible Multiphase Flow and Encapsulation Simulations Using theMoment-of-Fluid Method, Int. J. Numer. Methods Fluid., 79, 456(2015).
[17] Y. Lian, D. Liao, H. Qiu, M. Sussman, M. Vahab, Y. Hussaini, M. Arienti,－ 112 －Experimental and Numerical Investigation of Icing Process of a Liquid Droplet,
9th AIAA Atmospheric and Space Environments Conference, 4481(2017).
[18] M.G. Cooper, A.J.P. Lloyd, The Microlayer in Nucleate Pool Boiling, , Int. J.

Heat Mass Transfer, 12, 895(1969).
[19] M. Sussman, A. Almgren, J. Bell, P. Colella, L. Howell, M. Welcome, An
Adaptive Level Set Approach for Incompressible Two-Phase Flows, J. Comput.

Phys., 148, 81(1999).
[20] M. Vahab, K. Shoele, M. Sussman, Interaction of an Oscillating Flexible Plate
and Nucleate Pool Boiling Vapor Bubble: Fluid-Structure Interaction in a
Multimaterial Multiphase System, AIAA 2018 Fluid Dynamics Conference,
AIAA 2018–3718.
[21] V.K. Dhir, Numerical Simulations of Pool-Boiling Heat Transfer, AIChE J., 47,813(2001).

－ 113 －