本田博司*眞喜志治†山城光# Stability of Vapor Film in Subcooled Film Boiling on a Sphere* Osamu MAKISHI * Hiroshi HONDA † and Hikaru YAMASHIRO# Abstract The previously proposed linear stability theory of vapor film in subcooled film boiling on a sphere was generalized to take account of the interaction between the base flow and the perturbed component.    A disturbance of standing wave type was assumed to be superimposed on the base flows of the surrounding liquid and the vapor film.  For the surrounding liquid, the wave equation was applied to the whole region including the boundary layer and the energy equation was solved analytically by introducing a simplifying assumption.  For the vapor film, the basic equations were solved by the integral method.    By use of the compatibility conditions at the liquid-vapor interface, the solutions for the surrounding liquid and the vapor film were combined to yield an algebraic relation among the vapor film thickness, the order of disturbance and the complex amplification factor of the disturbance.  The numerical solutions of the critical vapor film thickness at which the real part of the complex amplification factor was equal to zero were obtained for the disturbances of the orders of 0, 1 and 2.    The numerical results indicated that the vapor film was most unstable for the disturbance of the 0-th order (i.e., uniform disturbance).    The calculated value of the critical vapor film thickness for the uniform disturbance compared well with the average vapor film thickness at the minimum-heat-flux point obtained from the immersion cooling experiments of spheres in water at high liquid subcoolings. Key Words:  Film boiling, Minimum-heat-flux point, Linear stability analysis, Sphere号

© 2006 The Heat Transfer Society of Japan
a
&
i
極小熱流束点
蒸気，蒸気膜

1格音
膜沸騰蒸気膜の崩壊に伴うクエンチと，それに伴う膜沸騰の極小熱流束（Minimum Heat Flux）点の発生機構の解明は蒸気爆発現象や種々の工業プロセスに関連して重要であり，多くの研究がなされている[1].この問題に関する従来の理論は，伝熱面温度が熱力学的な過熱限界温度まで下降し，固液接触が可能になることによって蒸気膜が崩壊するとする熱力学的不安定説[2-4]と，Taylor不安定[5]に起因する蒸気膜崩壊が起こるとする流体力学的不安定説[6-10]に大別される．本田ら[11-13]らは高温水平細線の浸漬急速冷却実験を行い，そのMHF点の挙動が上述のいずれの理論によっても説明できないことを見出した．ついで水平円筒上の自由対流過冷膜沸騰[13]と強制対流過冷膜沸騰[14]，および球上の自由対流過冷膜沸騰[15]について，気液界面の擾乱が印加された場合の蒸気膜と周囲液からなる系の安定性解析を行った．そして，蒸気膜の振動モードに応じた限界蒸気膜厚さる。が定まり,蒸気膜厚さるがる> δ。を満足する場合には蒸気膜は安定であり，δ<δ,の場合には不安定になることを示した．しかし，この解析には膜沸騰の定常流（ベースフロー）と擾乱との干渉の影響が無視されるなどの問題点があった.本報では，球上の自由対流過冷沸騰についてより一般性のある解析を行い， δ 。の計算値をDhir- Purohit[16]，西尾・上村[17]，菊池ら[18]による加熱金属球の水中浸漬冷却実験におけるMHF点の平均蒸気膜厚さおよび前報[15]の結果と比較する.
理論探折
2
球上の自由対流過冷膜沸騰について考える.物理モデルと座標系を図1に示す。半径座標をr，角座標をφ ,球の表面から計った外向き座標をyとする. 蒸気膜厚さδおよび蒸気膜にそう液の境界層厚さは球の直径dに比べて十分小さいものとする．周囲液の周方向および半径方向速度成分をu, , v, ,蒸気膜内のそれらをu, v,とする。周囲液および蒸気膜内のベースフローに鉛直軸に対称な定在波型の気液界面の微小擾乱が印加された場合を考える.ベースフローの諸量を添字。援乱成分を添字'をつけて表す.従ってδ = δ o + δ ' ,  u i = u io + u i等と表される。なお，物性値については添字。を省略する。
2.1周囲液
周囲液は境界層領域とその外部領域に分けられ,外部領域はほぼ静止液の状態にある．微小圧力擾乱が印加されると，境界層内ではベースフローと微小擾乱との干渉が起こるが,外部領域では干渉がない.仮定により境界層厚さは十分小さいので，ここでは境界層領域を無視し，圧力擾乱の分布は外部領域の解によって与えられるものと仮定する．基礎式は
Pi
ここにρ,は液密度,  tは時間.   P,は液圧力の擾乱成分である。断熱変化を仮定すると，ρとβの間には
の関係がある.ここにc ,は液の音速である.式(1)-(4)よりu 1 ,  v 1 ,  p 1を消去すると，波動方程式
が導かれる。境界条件は
(6)
式(5)-(7)の解は(8) P i = A x h x (i sr/c i )P x (cos φ )e s" ,   k = 0, 1, 2...ここにh,は球ハンケル関数,  P,は第一種ルジャンドル関数, s = s, + is, , s,は擾乱の時間増幅係数, s,は擾乱の振動数である.気液界面において，νとδの間には次の関係があタ
式(8)のPおよび式(9)のvを式(3)に代入するとδに関する微分方程式が導かれる。その解は次式で与えられる.
ただし
と定義する． |c t / sṛ|>>1の場合には，式(10)は
と簡単化される．後述の計算結果によれば， 2c, / sd |は10²のオーダーであり，第一近似としては式(12)で十分である.液境界層は気液界面の振動に伴って振動する．そこで座標yの代わりに気液界面を原点とする外向き座標y 'を用いると，周囲液のエネルギー式は擾乱成分に関して次式で表される.
ここにT 10はベースフローの液温度, T ,は液温の擾乱成分，κ,は液の温度伝導率である．対流項中の￠成分を無視すると，式(13)は次式で近似できる.
ここにq x = q s −q 。は蒸発熱流束の擾乱成分, q s = −λ, (∂T, /∂y) s , q c = −λ, (∂T, /∂y*) o , q evo = λ,ΔT,sar / δ o + α r ΔT s w − α r ΔT s wはベースフローの蒸発熱流束,えは液の熱伝導率，えは蒸気の熱伝導率， h,は蒸発潜熱, α, = σ, (T, 4 − T, 4 )/ ΔT, a,は放射熱伝達率, α,は周囲液への熱伝達率，σ,はステファン・ボルツマン定数, ΔT … = T … – T ,は伝熱面過熱度, ΔT seep = T s – T sは液の過冷度，T，は伝熱面温度，T，は飽和温度，T，はバルク液温度である。またTの分布として次式を仮定する.
ここによ= λ / α は温度境界層厚さである．境界条件は
ここにdT, / dP,は飽和蒸気圧曲線の勾配である．式(14)の解は
s >> (κ, / δ, − q, o / ρ,h, δ, )の場合には，式(18)は次式で近似できる.
計算結果によれば, ∫b / a , s /(κ, /δ, −qc,o /ρ,h,δ,)はそれぞれ10³，10²のオーダーであり，第一近似としては式(19)で十分である。式(19)をq。の定義式に代入すると次の解を得る.

2.2  蒸気膜
擾乱成分に関する連続の式，運動量式，エネルギー式を蒸気膜厚さにわたって積分するとそれぞれ次のようになる.
ここに
P,は蒸気圧力の擾乱成分, P,は気液界面における液圧力の擾乱成分,T,, = (T,, +T,)/2である.本論文ではδ << dを仮定しているのでP li = A 2 h k (isd/2c i )P k (cos φ )e st (28)と近似できる．さらに式(25)中の表面張力による蒸気膜の内外圧力差を無視するとP i = P iとなる．後述の計算結果によれば，表面張力項はP,の10 – 6 – 10 –5のオーダーであり，この仮定は少なくともd > 1 mmの範囲では十分成立する.式(26)のeの値は伝熱条件によって変化するが，ここではe=1を仮定する.擾乱の複素増幅係数sの絶対値が十分大きい場合には，式(22)，(23)左辺の慣性項，対流項のうちの空間微分項は時間微分項に対して無視することができ,それぞれ次式のように簡単化される.
計算結果によれば，慣性項，対流項中の空間微分項と時間微分項の比は10~²のオーダーであり，第一近似としては式(29), (30)で十分である.  ų,   T,の分布について
を仮定し，式(29)，(30)に代入すると
-86- © 2006 The Heat Transfer Society of Japanを得る.さらにρ , sδ² /12μ, << 1 ,  sδ² /12κ, << 1の場合には，式(33), (34)はそれぞれ
と簡単化される。計算結果によれば，ρ sδ² /12μ , sδ² /12κ,ほ10 −3 −10 −2のオーダーであり，第一近似としては式(35)，(36)で十分である．従って
(38)となる。ついで式(19)のTおよび式(36)のTをq 。の定義式に代入すると
を得る．ここでα çΔT no / ρ ιλι hi ε | s | << 1の場合には，式(39)は
と簡単化される．計算結果によればα 2 ΔT sub / ρ 1 λ 1 h 12 s iは10~5 –10~のオーダーであり，式(40)は十分正確である。さらに式(37)のu ，式(40)のq 。を式(21)に代入すると次式を得る.
式(41)にP ′ ,  δ ′ ,  ρ ′ ,  Τ ′の解を代入するとδ ′ ,  s, λに関する代数式が導かれる.その式にs = is iを代入し,実数部と虚数部に関する式を連立させて解くと,各kに対するδの安定下限界δ cと対応するs,の値が求まる.本論文ではφ＝0, k = 0, 1, 2の場合についてsとδ。の関係を求める．解くべき式は
© 2006 The Heat Transfer Society of Japanまた，φ＝0, k＝0,1,2とφ＝π/2, k＝0の場合についてようの解を求める。解くべき式はφ＝0についト
φ = π / 2 , k = 0について
ここにT , =(T , +T , )/2である.式(42)-(44)中のq , eを定義式より求めるにはα、の値が必要であるが，使用できるα ,の表示式は存在しない.そこで,ここでは9 goの面積平均値q goと球を離脱する蒸気流量との関係から求まる次式(式の導出については付録参照)を代わりに使用する.
ここにα,は飽和沸騰の場合の対流熱伝達率であり,菊池ら[18]の式（放射伝熱を無視し、蒸気膜内に直線温度分布と放物線速度分布を，また液境界層内に放物線速度分布と放物線速度分布を仮定した近似解析解）にΔT sub = 0を代入した値を用いる． 0.88α,は放射伝熱の補正項であり，西尾ら[17]の実験式にならっている．なお，式(45)を式(43)，(44)に適用する際にはδ = δ ςとおく．ついでq ∞の定義式にこのq enoを代入すると
が得られる．また,式(44)中のdδ。/ dφには気液界面でu, = 0を仮定した積分法によるベースフローの数値解を用いた.
３　計算値と実験値の比較
実際の計算は Dhir-Purohit[16], 西尾・上村[17],菊池ら[18]による加熱金属球の水中浸漬冷却実験と同一条件について行った．なお，計算に際して蒸気および液の物性値はそれぞれの膜温度で評価した.引用したデータの概略を表1にまとめて示す. 金属球の支持方法として, Dhir-Purohit[16]は直径3mmの鉛直棒の上に支持する方式をとり，球の下部５ mmの位置にカラーを取り付けている。このカラーは支持棒上の蒸気膜崩壊が球に到達するのを防ぐことを目的としている．一方，西尾・上村[17]，菊池ら[18]は細い支持棒の下に球を吊り下げる方式を取っている。菊池ら[18]はシリコーン系耐熱途料を薄く塗布した球に関する実験も行っている．塗布厚さは5-27 μmである．使用された金属球の表面は，加熱温度が高い場合にはいくぶん酸化していると考えられる。沸騰條相の親察結果によると，実験条件によって蒸気膜が全表面上で瞬時に崩壊する場合と，支持棒と球との継ぎ目付近で固液接触を開始し，伝播的に蒸気膜が崩壊する場合とがある．ここで採用したデータは瞬時崩壊の場合のものである.図2は引用した全実験データについて, MHF点の伝熱面過熱度  ΔT sat,MHF  と液の過冷度ΔT,,,の関係を示す. 図中にはLienhard  [19]の式による熱力学的限界過熱度と，Dhir-Purohit [16]の実験式も併せ示している. Dhir-Purohit [16]のデータはΔT ω 。の増大につれてΔT sai MHFも増大している．西尾・上村［17］のデータも同様の傾向を示すが， ΔT ω 。の影響はDhir- Purohit [16]のデータに比べてかなり小さい．菊池ら[19]のデータのうち，清浄面に関するものはΔT withの
影響がほとんどなく，一点を除いて過熱限界温度よりも低い．一方，耐熱塗料を塗布した面は相対的にΔΤ ,よりも高い．図２より，△T wi,мнғは実験条件によって大幅に変化し，ΔT weの関数としてまとめることは出来ないことがわかる.図 3 は沸騰液体が水でd =10 mm,   ΔT set = 200 K, ΔT sub = 70 K, φ = 0の条件について，式(42)から求まふs, ,  s,のδによる変化をk = 0, 1, 2の場合について示す．なお，式(42)はAT supを陽に含まないが，液の物性値はAT supに依存するので,この解は非常に弱いΔT supの関数である．k = 0, 1, 2に対するs,の曲線は一度だけ水平軸s, = 0と交差する. 従ってこの点のδοからδεである。δο <δεではs,>0, δε <δοではs, < 0であるから, δ 0 < δ 。では蒸気膜は不安定であり，δς <δ 0では安定である．δςの値はkが小さいほど大きい.従って蒸気膜の安定性はｋ＝０の擾乱,すなわち一様擾乱に対するδこによって決定される.この値をδρと表すことにする.ではs, <0であるから，δρ <δρでは蒸気膜は不安定であり，δρ <δρでは安定である。δεの値ははが小さいほど大きい．従って蒸気膜の安定性はk = 0の擾乱,すなわち一様擾乱に対するこによって決定される。この値をδρと表すことにする.図4は西尾・上村[17]の実験条件について本解析と前報[15]の解析によるる。とΔT sanの関係を比較したものである．前報[15]では表面張力による気液の圧力差を考慮しているが，計算結果によればこの圧力差は無視できる.従って表面張力項を無視すると,式(41)に相当する前報[15]の式は
となる．図中には¢＝0，k＝0,1,2と¢＝π/2，k＝0に対する本解析の解と，k = 0,1,2に対する前報[15]の解を示している．両者を比較すると，δεの勾配は前報[15]のほうが大きく,両曲線はΔT = 100 Kのあたりで交差する。 δςの値はΔT , α > 100 Kでは本解析のほうが小さく，両者の差はΔT, , ,の増大につれて増大する．しかし，その差はΔT sm = 350 K でも4 %以下であり，ベースフローの影響は小さい.図５は図４中のk＝0に対する解る。の一部を拡大して示す。 φ＝0ではΔT …の影響は無視できるほどである．φ＝π/2radではΔT …の増大につれてδ。も増大するが，φ＝0の解との差は±0.24%以下である．図6, 7, 8, 9それぞれ西尾・上村[17], Dhir- Purohit[16],菊池ら[18]の清浄面および菊池ら[18]の耐熱塗料塗布面の実験データから求めたδ MHFの値とδρの計算値をΔT …の関数としてプロットして比較したものである。なお，ΔT,が大きい場合には膜沸騰におよぼす放射伝熱の影響が無視できないので，δмнғの計算に際して考慮している．また，図中には菊池ら[18]の膜沸騰熱伝達率α,の式から求めた平均蒸気膜厚さるのΔT ,による変化をΔT , ,をパラメータとして示している．なお，菊池ら[18]の解析では放射伝熱が無視されているので，放射伝熱を考慮したるの値を
によって推測した．この式は放射伝熱を考慮した膜沸騰熱伝達率と蒸気膜厚さの関係式
より導かれる．計算結果によれば，式(48)右辺の第２項は実質的に無視できる.図6の西尾・上村[17]の実験データとの比較についてみると，δмнғはΔT satの増大につれて減少し， ΔT sab = 70–73Kに対するδmHFはδc,oとよく一致している．また，δもΔT, , ,の増大につれて減少し，同時にδмнгとの一致がよくなっている。 ΔT wi-mhгの値は各ΔT , ,ごとに一定の範囲に分布しており,その値はΔT  , ,の増大につれて増大する傾向がある．δの曲線の勾配はふ。の曲線の勾配より大きく，両曲線は一点で交差する．従って，理論的にはこの交点よりΔΤ ωの大きい領域では蒸気膜は安定であり，この交点において蒸気膜が不安定となって崩壊する。しかし実際には，ΔT…が小さい場合には発生蒸気の一部は蒸気泡を形成して球から離脱し，その際に発生する気液界面擾乱によって早めの蒸気膜崩壊が起こると考えられる.  ΔT, が増大すると蒸気泡の離脱はなくなり，擾乱の振幅も減少する．そして， ΔT sub 70–73Kでは本理論で想定した状況に近い微小擾乱になり，その結果δмнғとδρがほぼ一致したものと考えられる.なお, ΔT ser,мнғの値に幅があるのは,気液界面擾乱の強さにばらつきがあるためと考えられる。また，蒸気膜崩壊後の固液接触部では瞬間的に非常に高い蒸気発生があり，その際の反発力によって液が持ち上げられて蒸気膜が再び形成される可能性がある．従って，実際のMHF点は球表面の満れ性にも関係すると考えられる.図７のDhir-Purohit[16]のデータも基本的には西尾・上村[17]のデータと同様の傾向を示している.しかし実験範囲がΔT s ≤ 50 K に限られているので,全データについてδмнғ > δσ 0となっている．彼らはステンレススチール，銅および銀製の球に関する実験を行っているが，銀球に関するるNHFの値は他の球に比べて低い.この原因は不明である.図8は菊池ら[18]の清浄面に関する実験データとの比較を示す。この場合にもδμμはΔT smbの増大につれて減少し，ΔT sat = 80−85 Kではδ MHFがδτοよりも小さくなっている．その原因は明確でないが，彼らの試験球の上端部は一体加工された直径約8 mm,高さ10 mmの鉛直円筒部に接続しており，試験球は鉛直円筒部の中間まで浸漬されるので，伝熱面形状の違いが伝熱特性に影響している可能性がある．またらもΔT subの増大につれて減少し，ΔT sub = 82 Kではδ<δαρとなっている．菊池ら[18]の式は近似解析解であるが，その詳細が不明であるので，この問題についてここで議論することは出来ない.図9は菊池ら[18]の耐熱途料途布面に関する実験データとの比較を示す.この場合, δmHFに対応するΔT windeの値は清浄面に比べて大幅に高くなっている。実験データの傾向は清浄面の場合と同様であるが, ΔT ω = 80 Kに対するδмнғとδ ωの差はΔT ωの増大につれて減少し，高過熱度では両者が良好に一致している.図8と9の清浄面と耐熱塗料塗布面のAT sel,mhf の範囲をΔT … = 80–85 Kの場合について比較すると,清浄面は77–241 K，耐熱途料途布面は174–533 Kである．従って，ΔT ω = 533 Kでは清浄面についても気液界面の不安定に起因する局所的な固液接触がすでに始まっていると考えられるが，熱流束の上昇を伴う大規模な固液接触には至っていない．固液接触直後の伝熱面過熱度ΔŤ sat,мнғは非定常熱伝導の解より
で与えられる．ここにＣは比熱である．この式に△T sat,мнғの測定値を代入すると，銀球についてΔŤ sat MHF =73.2–229.2 Kとなる．また，耐熱途料の物性値をシリコーン樹脂の値で近似すると，耐熱塗料塗布球についてΔŤ sat , MHF  = 50.4 – 158.6 K  となる.したがって△Ť sat , MHFの範囲は後者のほうが小さめである.そして，ΔŤ sat,мнғの最大値は水の限界過熱度217.3 Kに近い. この結果は，ΔŤ sat , MHFが液体の限界過熱度以下になることが瞬時的固液接触開始の条件であるとする西尾・上村[17]の説を支持するもののようである．一方，Dhir- Purohit[16]のデータについは， ΔΤ sai,mhfの最大値は銀球で256.8 K，銅球で337.2 K,ステンレススチール球で450.9 Kであり,限界過熱度を超えている．しかし，彼らのデータのΔT sql MHF対ΔT subの関係は，西尾・上村[17]が蒸気膜の伝播的崩壊を観察した実験条件である球の浸漬深さが3.0dの場合と類似しているので，必ずしも瞬時的崩壊ではなかった可能性がある．この点については今後の研究に待たねばならない.
4　結　自
球まわりの自由対流過冷膜沸騰について，伝熱と流動の複合した蒸気膜の安定性解析を行い，安定な蒸気膜厚さの下限界が存在することを示した．ついで，この理論による限界蒸気膜厚さδ。。の計算値をDhir-Purohit[16]，西尾・上村[17]，菊池ら[18]による加熱球の水中浸漬冷却実験データから求めた，瞬時的蒸気膜崩壊が起こる場合のMHF点の平均蒸気膜厚さらмнгと比較した。また，ベースフローの影響を無視した前報[15]の安定性解析の結果とも比較した.前報[15]と本解析のδ,の値を比較すると,  ΔT,,, > 100 Kでは本解析のほうが小さい.しかし両者の差は小さく，ベースフローの影響は実質的に無視できる程度である．また，δρに及ぼすΔT …。やφの影響も小さい。 δ MHFの値はΔT , ,の増大につれて減少し, ΔT sws > 60 Kではδ oとほぼ良好な一致を示した．また，西尾・上村[17]，菊池ら[18]のデータによれば, ΔŤ sat,MHFの値は限界過熱度を大幅には上まわらなかった．以上の結果より，本論文に述べた蒸気膜の不安定性による固液接触開始は高過冷度膜沸騰のMHF 点条件の重要な要素であり， AT sat , MHFが限界過熱度以下になることが全面的な固液接触開始の条件であると考えられる.
I evoの表示式の導出
村家
蒸発熱流束の面積平均値q 。。と球を離脱する蒸気流量M の間には
の関係がある．M yはφ = π / 2における周方向蒸気速度u oを用いて
と表すことが出来る。ここにaは比例定数である.飽和沸騰では全伝熱量が蒸発に使われるので,
が得られる．このaを式(A2)に代入し，次いで得られたM,を式(A1)に代入すると，可∞の一般的な表示式として
が得られる.