1.はじめに
加圧水型原子炉で想定される事故，例えば小破断冷却材喪失事故の条件下では，炉心で発生した蒸気が一次系配管や蒸気発生器（SG）の伝熱管に流入し，放熱により生じる凝縮水が蒸気の流れに対向して重力落下する気液対向流状態になり，蒸気速度が大きいと凝縮水の落下が制限されるフラッディングが生じる.事故解析における配管系の計算では一次元二相流モデルが使用されるが,フラッディングによる落下水流量の計算は不確かさが大きいため,蒸気の見かけ速度♫と落下水の見かけ速度♫との関係を与える気液対向流制限（counter-current flow limitation: CCFL)の相関式が使用されている.著者らの研究グループでは,神戸大学で取得したデータと従来データおよび解析結果を活用し,ホットレグ,加圧器サージ管, SG伝熱管を対象にしてCCFL相関式の信頼性向上と不確かさの定量化に取り組み(1)~(3), CCFL相関式の改良は一段落した.気液対向流ではCCFLに次いでボイド率αの評価が重要であり，現在はαを含む流動特性について評価している.鉛直管での気液対向流について1980年ころまで多くの研究が実施された(4)(5).しかし,気液対向流でのαの測定例は少なく, Imuraら(6)は自由落下液膜厚さの相関式を用いてαを評価した. Bharathanら(7)は空気・水系で液膜厚さδと圧力勾配dP/dzを測定して壁面摩擦係数f,と界面摩擦係数fを求めたが，fが負になる場合が生じた.そこで, Bharathanら(8)およびBharathan and Wallis (9)はは,下端でのフラッディングによる粗い液膜(rough film: RF）を対象にf, = 0と近似してdP/dzデータからf,を求めて相関式を提案した. Abeら(10)は, Bharathanら(7)(8)によるJς, JιとdP/dzのデータを用いてα, f,とf,を求めf,は大きくなることを示したが，f,の予測値に対する信頼性は不明である. このような背景から,神戸大学の富山研究室では,空気・水系でJcに対するJr, dP/dzとαを測定し,環状流モデルを用いてf,とf,を求めた(11)(12).著者らの研究グループでは,神戸大学で取得した基礎データを活用して，フラッディング状態における流動特性に関する相関式を整備することを目的とし,相関式の適用範囲を大口径管や高温高圧条件に拡大することを目標にしている. 従来研究ではαデータが少ないことから, Takakiら(13)はCCFLデータ（JςとJ,の関係）とWallisのエンベロープ法(4)を用いてaとfを求める方法およびdP/dzデータとfの相関式を用いてαとf,を求める方法について検討し,前者は不確かさが大きく後者は評価精度がよいことを示した.また, Bharathanら(7)による空気・水系での) dP/dzデータおよびIlyukhinら(14)による蒸気・水系でのdP/dzデータを用いて上端フラッディングにおけるαを求めてふの相関式を提案した(13).下端フラッディングについて，佐野ら(15)はBharathanら(8) 9)による空気・水系でのdP/dzデータおよびIlyukhinら(14)による蒸気・水系でのdP/dzデータを用いてf,を求めて相関式を提案した.本報では，神戸大学での空気・水実験(11)(12)および上端フラッディングでのδの相関式(13)と下端フラッディングでのf,の相関式(15)について紹介する.なお,下端フラッディングに関する実験(11)は神戸大学での単独実施であり，上端フラッディングに関する実験(12)は神戸大学との共同研究の一環として寒施した.
2　実験装置と環状流モデル
2.1 実験装置
鉛直管での空気・水フラッディング実験装置(12)を図1に示す. 鉛直管の上端が曲面で下端が直角の下端フラッディング実験(11)と上端が直角で下端が曲面の上端フラッディング実験(12)で装置の構成は同じであり,曲面の曲率Rは鉛直管の直径Dの1/2(R = D/2)である.下端フラッディング実験(11)ではD = 20, 40 mm,上端フラッディング実験(12)ではD = 20 mmであり，本報では主にD = 20 mmでの測定値について紹介する.実験装置は,下部タンクへの空気の供給系,上部タンクへの水の供給系,鉛直管で構成され，鉛直管の直径と高さおよび上下端の形状を変更することができる.実験では,空気流量と水流量を一定値にして,鉛直管に流入しない水をドレンして上部タンク水位を一定値（本実験では0.1 m）に保つ.落下水の流量と見かけ速度∑は下部タンクの水位上昇率から測定する. αは弁締切法を用いて60回測定した平均値を使用している. 弁の締切時間は1/30秒以下であり，α測定の不確かさは（1-α）＝0.05～0.13に対して0.0025 ～ 0.0040である．弁体は厚さ2 mmのスライド板であり，締切弁の設置がJ,に影響しないことを確認している.
2.2 環状流モデル
鉛直管内における気液対向流では環状流になり.流動状態の計算には環状流モデルが使用される. 気相流路および流路全体に対する圧力バランスは次式で表せる(11)(13)(15) dP/dz + ρς g + 2f ; ρς[Jς /α - Jι/(1-α)] ²/(D・α 1/2 )= 0  （1） dP/dz +[(ρι. -ρς)(1- α)+ρς]g -2ƒ,. ρι[Jι/(1-α)]2/D = 0 （2) gは重力加速度, Pは圧力, zは鉛直方向座標, pは密度である. zの上向きを正にすると圧力勾配dP/dzは負になる. Jςに対するJι, dP/dz, αを測定すると,式(1)から界面摩擦係数f,を,式(2)から壁面摩擦係数f,を求めることができる.
2.3 実験条件
本報で用いるデータの実験条件を表1に示す.  aとf,の評価が目的であり， CCFL測定のみの実験は対象外である.  Bharathanら(7)は,  D = 50.8 mm一定で，上端と下端に対して曲面（rounded: R）と直角（square: S）の全ての組み合わせで実験している. Bharathanら(8)は,上端と下端とも直角で直径を変更して実験し, RF（粗い液膜）でのdP/dzのみ報告している.  Ilyukhinら(14)の蒸気・水系を除いて空気・水系実験である. 環状流モデルでは, αとらとは相互に変換できる. (3)
S-W: steam-water
3.実験結果
本章では, D = 20 mmでの下端フラッディンケ(11)と上端フラッディング(12)のデータを用いて,フラッディング状態での流動特性について紹介する.
3.1 流動様式
高速度カメラを用いて流動様式を撮影（350 fps)し、得られた画像のうち管軸上の輝度値のみを抽出して時間軸方向に並べた画像（time-strip画像(12)を図2に示す.  D = 20 mm, 上端は直角で下端は曲面である.低Jででは上端で液相の落下が制限されて気液界面が滑らかな液膜（smooth film: SF)が形成される. 液膜が厚い黒い部分が時間経過により下方に移動している. Jcを増加すると下端でも液相の落下が制限されて援乱波が生じ,上昇空気により擾乱波が上昇するが途中で落下する遷移状態（(transition: TR)が生じる. TRでは液相は上端と下端で落下が制限されている.高Jででは下端で発生した援乱波が上部タンクに流出し，RFが形成される.  RFでは液相は下端で落下が制限される. time-strip画像を用いると流動構造を鮮明に可視化できる.
3.2 圧力勾配とボイド率
表1に示した実験条件のD = 20 mmでR/S (11)とS/R (12)における無次元圧力勾配（dP/dz)*と液相体積率（1-α）の測定結果を図３に示す. (dP/dz)*は重力（ρι-ρς)・gで無次元化している. 式(2)から明らかなように, (1-α)と(dP/dz)*の差が壁面摩擦項であり, Bharathanら(8)およびBharathan and Wallis (9)のようにƒ, = 0と近似すると(1-α) ≒ −(dP/dz)*になる.壁面摩擦項はRFで比較的小さいが，SFでは大きく-(dP/dz)*が小さい.流動様式の変化は,図３に示すように−(dP/dz)*の変化で判別できるが，流動様式の境界が不明確な場合には図２に示すように可視化で確認する必要がある.
3.3  CCFL特性
J,を式(2)で計算できるが精度がよくないため,原子力プラントの事故解析では, WallisによるCCFL相関式(4)が広く使用されている.
Hçは無次元速度, J,は見かけ速度, wは代表長さであり，定数と勾配mはCCFL測定値から決定される。一方，Bankoffら(16)は代表長さｗを次式で表している.
Lはラプラス長さ, ♂は表面張力である.式(6)の代表長さｗによりDがCCFL特性に及ぼす影響が決まるため、CCFL相関式の導出では式(6)のβの選定が重要である.式(4)の無次元速度Ｈ,*は,β = 0でWallisパラメータJ๊に, β = 1でKutateladzeパラメータK๊になる。 J 。とK 。は無次元直径D を用いて相互に変換できる. J k = K */D*¹ ² (k = G or L), D* = D/L
事故解析では, 例えば計算モデルでCCFLが生じる位置に定数Cと勾配mおよびβを入力してCCFL相関式を適用する. D = 20 mmでR/S (11)とS/R (12)におけるCCFL特性を図4に示す.  Kusunokiら(17)は下篇フラッディングに対して次式を提案した.
(9)式(9)の不確かさ± 0.05にはIlyukhinら(14)のD = 20 mm,  P = 0.6 ～ 4.1 MPaでの測定値を含む. Wallis相関式(4)はJ *1/2とJ *1/2との線形式であるが,式(9)のようにJc*1/2をJr*1/2の三次関数で表す方がよい場合がある.  Takakiら(13)は烏村ら(12)による上鶏フラッディングでの測定値に対し式(10)を導出し, Ilyukhinら(14)のD = 20 mm,  P = 0.6 ~ 4.1 MPaでの測定値にも適用できることを示した. J G *1/2 + 1.02 J *1/2 = 0.83±0.05 ( D = 20 mm) (10)式(10)はD = 20 mmのみ適用でき,  D = 30 mmでははK パラメータ支配になる（Muraseら(18)).
3.4 壁面摩擦
壁面摩擦を直接測定するのは困難であるため,一般的にはJςに対するJς, dP/dz, αを測定して式(2)から広を求める.気液対向流に対する広を評価した例は極めて少ないため，近似的に単相流に対するんの相関式を用いる場合がある（数土(19) ; Takaki 5 (13)). f , = max (16/Re , 0.079/Re , 25 ), Re , = J , D/v ,    (11) Re ,は液相レイノルズ数, ν,は動粘度である.  Abeら(10)はBharathanら(7)(8)によるdP/dzデータを用いてか,を評価して次式を提案した.
D = 20 mmでR/S (11)とS/R (12)におけるƒ,を図５に示す. 上端フラッディングのSFではf,を式(11)で近似できる.下端フラッディングではƒ,が式(11)より大きくなるため, Godaら(11)は居流に対して次式を提案した.
式(13)はRFを対象としており測定値とよく一致する.
3.5界面摩擦
Bharathanら(8)およびBharathan and Wallis (9)は，下篇でのフラッディングによるRFを対象にf, = 0と近似してD = 6.4～152 mmでのdP/dzデータ
からf,を求めて次の相関式を提案した.
fi = 0.005+A (δ/L) β, log 10A = -0.56+9.07/D*, B = 1.63+4.74/D*
数土(19)は,式(14)をベースにして, Richter (20)によるD = 19～140 mmでのCCFLデータとWallisのエンベロープ法(4)を用いて, CCFL計算値が測定値と一致するf,を求めて次の相関式を提案した. f = 0.008 { 1 + m (δ/D)"}, m = 41.31D*(x+0.25)109.07/D*, n = 1.63 + 4.74/D*      (15) D = 20 mmでR/S (11)とS/R (12)におけるfを図6に示す.  SFよりRFの方が流路内の水量が多く液膜が厚いためf,が大きくなる. 式(14)と(15)はは，δの関数であるためR/SとS/Rで計算値が異なり, S/Rを対象に計算した.式(14)はRFを対象に提案されたがSFの測定値と比較的よく一致している。式(15)の導出ではエンベロープ法を使用してJr測定値に一致するf,を求めたため過大なf,をあたえる。 RFでの測定値とよく一致する従来相関式はなく, Godaら(11)はRFに対して次式を提案した. f i = 0.11J c * d (1-a), a = -1.96～1.98, b = 0.42～0.43
式(16)はRFを対象とし、RFでの測定値とよく一致する.
4.ボイド率と界面摩擦の評価方法
気液対向流条件での鉛直管内における流動特性は落下水が制限される位置により異なり，流路内での水量(1-α)がδ, ƒ…, ƒに影響する.著者らの研究グループの目的は，神戸大学で取得した基礎データを活用して，フラッディング状態における流動状態に関する相関式を整備することである。流動状態に関する相関式の適用範囲を大口径管や高温高圧条件に拡大することを目標にしている.しかし,表1に示したように従来研究ではαやらの測定例が極めて少ない．そこで，本章ではCCFLデータやdP/dzデータからαやfを評価する方法について紹介する。 αやf,の評価方法を検討するために，ここではD = 20 mm,  S/R (12)での測定値を使用する（詳細はTakakiら(13)を参照).
4.1    CCFLデータとエンベロープ法
式(1)と(2)からdP/dzを消去して無次元化すると,次のJςとJςの関係式が求まる.
J k はWallisパラメータである.  Wallisは, f iとf ,の相関式を用いαを変数としてJgとJ,の関係を計算し、そのエンベロープ（図７を参照）がCCFL曲線になるとした(4) .しかし, Bharathan and Wallis (9)は,エンベロープ法によるJ,の計算値は過大であり, J,の理論的な最大値を与えるとしている.エンベロープ法は次の最大流量条件と等価である. ∂J¸*/∂α = 0 or ∂J¸*/∂α = 0
α計算の検証を目的とするため，αの測定値があるD = 20 mmでS/R (12)を対象にエンベロープ法を適用した。 fには図５より式(11)を使用し、fには図６のS/RでSFでのデータを用いて最小二乗法により次式を導出して使用した.
各計算の番号と使用した式を表２に示す. Wallisのエンベロープ法(4)に従って，式(17)を用いαを変数にしてJcとJcの関係を計算した表2のNo. 1の結果を図７に示す. JcとJcの関係のエンベロープは，式(19)の最大流量条件を用いたNo. 2の計算結果と一致する。式(19)でαJ,*/∂α = 0を用いるとJ,* = 0で計算が不安定になるため，著者らは∂Jς*/∂α = 0を用いた.数土(19)は，式(17)と(19)および式(11)と(14)を用いてJ cとJの関係を計算したが，Richter (20)によるD = 19～140 mmでのCCFLデータと計算値が一致するfを求め式(14)を式(15)に変更した.数土(19)と同様に，図４に示したS/RでSFでのCCFLデータに一致するfを求めた. f i = 0.008 + 2000[(1-α 1/2 )/2] 3   (D = 20 mm)      (21) No. 2の式(20)を式(21)に変更したNo. 3の計算結果はCCFL測定値と一致している.エンベロープ法を用いた（1-α）の計算値と測定値(12)の比較を図８に示す.エンベロープ法を用いると（1-α）を過大に計算する.
B. eq., a  eq., fi eq. and f., eq. denote equation number for basic equation, α , f i , and f ,, respectively.式(2)でdP/dz, J cとJ l , f vを与えればαが求まる。つまり，JとJの関係をCCFLデータで与え,広の相関式とdP/dzデータを用いれば式(1)からαが求まる. αが求まれば式(1)からf,が求まる. f,の影響を評価するため、表２に示すようにf, = 0, Bharathan and Wallis (9)が推奨するf, = 0.005, および式(11)について検討した. dP/dzデータを用いた（1-α）の計算値を図9に示す. ƒ = 0.005を用いたNo.５の計算値は，Wallis (4)が推奨するSFのよに対する次の相関式とよく一致する. δ/D = 0.063/ , * 2/3 (22) dP/dzデータと式(11)を用いたNo. 6の計算値は(1-α)の測定値とよく一致する.計算値が中にで(1-α)を過小に評価するのは,図５で式(11)が層流と乱流との遷移域でf,を過小に評価することによる.
5.上端フラッディングでの液膜厚さ
亀井ら(21)は気液対向流条件でのδ測定を初めて試みたとしている. 亀井ら(21)は天秤で試験部の重量を測定して試験部内の液重量を測定しているが，ガスを流動させると天秤の精度が悪く定量的に言及できないとしている.  Feind (22)やHewitt  and Wallis (23)は気液対向流条件でふを測定したが、前者はフラッディングが生じる前，後者は唇流条件で測定値は自由落下液膜厚さとほぼ等しい. Imuraら(6)はSFとRFを区分せず，Feind (22)による相関式を用いて気液対向流条件でのαを評価した.自由落下液膜厚さに対して多くの相関式が提案されているが，比較用としてNusselt (24)とFeind (22)による式を組み合わせて用いた. δ/(ν ι ² /g) 173 = max  {(3Re ι /4) 173 , 0.266Re 172 } (23) Bharathanら(7)によるる測定以後，神戸大学でのα測定(11)(12)まで気液対向流条件でのαもしくはδの測定は見当たらない. Takakiら(13)は，図9に示したようにdP/dzデータと式(11)を用いるとSFでのαを精度よく評価できることを示し、表１に示したBharathanら(7)の空気・水系データとIlyukhinら(14)の蒸気・水系データを用いてSFでのαを評価した(Bharathanら(8)はD = 6.4 mm ～ 152 mmでdP/dzを測定しているがRFでの測定値しか報告していない） αの評価値から求めたるを図10に示す. 品村ら(12)はα測定値から,その他(7)(14)はdP/dzデータから求めたらである.  Takakiら(13)は，図10に示した評価値から乱流に対して次式を提案した. δ/(v L 2 /g) 173 = 0.091 Re L 0.64 (Re L ≥ 1800)
上端フラッディングについては,図５に示した層流から乱流への遷移域でのf,に対する相関式の導出，図10に示した層流から乱流への遷移域でのらに対する相関式の導出，fの相関式の導出などの課題がある. これらの課題については,  D = 40 mmで追加実験を行い，追加データを含めて検討する予定である.
6.下端フラッディングでの界面摩擦係数
佐野ら(15)は,表１に示した従来研究によるdP/dzデータを用いて、式(2)からαを、式(1)からfiを求めた。 ∑の相関式には，式(11)の層流項をGodaら(11)による式(13)に置き換えて使用した. f , = max (2.86×10 4 /Re , 1.96 , 0.079/Re , 0.25 ) (25)多くの実験データを用いて相関式を導出して適用範囲を広げると不確かさが大きくなる.そこで, SG伝熱管寸法であるD = 20 mmでのf。空気・水系でのf,直径と流体物性値の影響を含むf,に区分して相関式を導出した. D = 20 mmでのf,を図11(a)に示す.圧力条件ごとに最小二乗法を用いてf,の指数関数を求め,係数と指数を気液密度比ρς/ριもしくは気液粘性比μς/μιで整理して最小二乗法を用いて累乗関数を求めた. R²値が１に近いμς/μιを選択し、次の相関式を導出した.
式(26)による計算値fc。と測定値fc 。の比較を図11(b)に示す.測定点数の95％を含む式(26)の不確かさは,空気･水データ(11)に対して± 4.2 %,蒸気･水デー γ (14)に対して± 15 %である. D = 6.4 ～ 152 mmでのfを図12(a)に示す. 図11 (a)に示したJç*を用いるとKç* = D*¹/2 Jç*でf,を良好に整理できた.直径ごとに最小二乗法を用いてf,の指数関数を求め,係数と指数をDで整理して最小二乗法を用いて累乗関数を求めた. Kçの係数と指数をD°の累乗関数にした次のf,相関式を導出した. f (= 0.30D* exp(-1.90K c *), 6.4 mm ≤ D ≤ 152 mm（27）式(27)は式(15)や(16)より簡潔である.式(27)による計算値f cと測定値f cの比較を図12(b)に示す.測定点数の95%を含む式(27)の不確かさは44 %である.式(27)でD = 20 mmとし,式(26)の導出と同様に,係数と指数をμς/μιで整理して最小二乗法を用いて累乗関数を求めた.  D*とμς/μιを含むf,相関式として次式を導出した. f i = 0.157( μς/μι) + 166D* exp {-1.14( μς/μι) + 133K ς*} (28)空気・水系に対しては，式(28)は式(27)と一致するため式(28)による計算値と測定値の関係は図12(b)と同じである.式(28)による計算値f cとD = 20 mmでの測定値fi, wの比較を図13に示す.測定点数の95％を含む式(28)の不確かさは，空気・水データ(7)(8)(11)に対して±51 %,蒸気･水データ(14)に対して±10 %である.下篇フラッディングについては, αもしくはδの相関式化は未検討である。図10と式(24)で示した上篇フラッディングではδをRe」の関数で整理できたが，下端フラッディングでのαもしくはδについては実験ごとの相達が大きく主要な影響因子を特定できていない.これまでに上端が直角の上端フラッディングに対する実験（D = 40 mm）をほほ終了し、データ評価を実施中である.現在は上端と下端が曲面での管内フラッディングに対する実験(築)を実施している.
7.まとめ
本報では，鉛直管でのフラッディングを対象として，神戸大学での空気・水実験(11)(12)および上篇フラッディングにおける液膜厚さるの相関式(13)と下端フラッディングにおける気液界面摩擦fの相関式(15)について紹介した.鉛直管でのフラッディングについてWallis研究室で多くの実験が行われたが，信頼できるポイド率αの測定値がなく壁面摩擦係数を広≒ 0 ((1-α)が無次元差圧-(dP/dz)*とほぼ等しい(1-α) ≒−(dP/dz)*)と近似して界面摩擦係数f,が求められた.つまり,フラッディング条件におけるα, f ,, ƒ,の評価の不確かさが不明であった.そこで,神戸大学の富山研究室では弁締切法を用いてαを測定し，環状流モデルを用いてƒ„とƒ„を求め，これらを定量化した.著者らの研究グループでは,神戸大学で取得されたデータを用いて従来研究による測定値を有効に活用する方法について検討し、αとf iの相関式の適用範囲の拡大を目指している. Wallisのエンベロープ法を用いると(1-α)を過
ギリシャ文字
添字
略号