気液二相流のパターンダイナミクスモデリング

日本マイクログラビティ応用学会誌 Vol. 29 No. 2 2012 (84–91)
特集：理論化・モデリング
（解説）
その 2
気液二相流のパターンダイナミクスモデリング
小澤
守 1・網 健行 2
Pattern Dynamics Approach to Two-Phase Flow Dynamics
Mamoru OZAWA 1 and Takeyuki AMI2
Abstract
Conventional modeling of two-phase flow, including homogeneous, drift-flux and two-fluid models, is based on timeaveraged properties, even in the case of simulation of transient behavior. Then the inherent fluctuations of two-phase flow
are neglected so as to ensure smoothed and continuous properties of parameters, such as void fraction and velocity. Gasliquid two-phase flow is a typical complex system, and thus to simulate inherent fluctuation, e.g. of slug and churn flow, is
rather hard with above-mentioned conventional modeling. This paper describes one of the approaches to such simulation
aiming at pattern formation and/or evolution of void fraction fluctuation of two-phase flow. The present simulation is
conducted based on the void propagation equation and additional limited number of momentum effects, so that the
complex properties are successfully realized on a computer.
1.
基本的な考え方
気液二相流はボイラや復水器など蒸発・凝縮を伴う流
れを代表する総称であり，そこで形成される気液界面は
変形・合体・分裂を伴いながら時空間にわたって大きく
変動する．その結果形成される幾何学的構造のパターン
である流動様式に応じて気相・液相の配置並びに体積率
分布（ボイド率）が様々に変化し，熱伝達特性や圧力損
失が大きく影響される．
気液二相流の重要性は水管ボイラの開発当初 1,2) から
水循環問題に関連して認識されていたが，適切な計測技
術や解析方法がない時代にあって，当初は Münzinger の
“Dampfkraft”3) でも明らかなように，相対速度の存在
は認識されながらも理論的扱いは均質流モデル即ち気液
間の相対速度がないとしたモデルが用いられていた．当
然ながら実際とは一致しなかったが，その当時のボイラ
の燃焼室熱負荷が現在に比してかなり低く，熱計算にお
いてはガス側の放射伝熱や対流伝熱，水管壁における熱
伝導を把握しておれば十分通用したのだろう．つまり管
壁を通しての熱移動はほとんどガス側に律速されており，
管内沸騰熱伝達は熱抵抗としてほとんど有意な大きさを
持たなかったのである．ところが第 2 次世界大戦後，米
国において軽水炉開発が進展するに伴って，ボイラとは
格段に発熱密度の大きい原子炉炉心冷却問題が顕在化し，
沸騰研究と相まって二相流研究が大きく進展した 4,5)．
わが国においていち早く気液二相流研究を開始した石
谷 6) はその当時の研究を「単一境界面理論」（分離流モ
デルが代表例）とし，今後の展開の方向として「複合境
界面理論」構築に向けて進むことを予測した．当時の二
相流の主流はスラグ流などのような間欠流的な変動を無
視した時間平均モデルが大前提であり，さらに石谷は二
相流の本質について「しかるに従来の研究では，圧力や
流量の脈動成分は殺してしまい，平均値だけを問題とす
るのが例であった．．．圧力脈動には二相流の本質に関す
る重要な情報が含まれている．．」6) との見解を提示した．
つまり二相流の本質はボイド率や圧力の変動であると述
べたのである．一方においてボイラの水循環設計などに
供するには詳細な評価方法でなくても時間平均値として
のボイド率相関式・理論式があればよく，変動特性まで
も組み込んだモデリングに対する要求はそれほど高くな
かったのだろう．
ではそもそも石谷の指摘した二相流に内在する変動特
性は単一境界面理論を用いて予測可能なのだろうか．否
である．つまり境界面が一つあるいはたとえ複数であっ
ても固定的であればボイド率の時空間変動など起こりえ
ないのであって，問題解決には，少なくとも複数の界面
があり，それらが流動のダイナミクスに従って相互に干
渉すること（これが「複合」の意味）が必要であろう．
1
関西大学 社会安全学部 〒569-1098 高槻市白梅町 7-1
Faculty of Safety Science, Kansai University, 7-1 Hakubai-cho, Takatsuki, Osaka569-1098, Japan
2 関西大学 システム理工学部 〒564-8680 吹田市山手町 3-3-35
Faculty of Engineering Science, Kansai University, 3-3-35 Yamate-cho, Suita, Osaka 564-8680, Japan
(E-mail: [email protected])
− 84 −
気液二相流のパターンダイナミクスモデリング
界面の定義方法を決めてやれば二相流固有の変動，すな
わち気泡の合体や複雑な相分布など詳細な三次元流(3D)
の数値シミュレーションを行うことは可能であるが，一
般には膨大な計算をこなさなければならない．それに対
して本稿で説明する離散気泡モデルは上記のような数値
計算とは汎用性や計算から得られる情報量あるいは分解
能はかなり限定されるが，石谷のいう「複合」境界面の
考え方を簡便な方法で実現すべく開発された一次元流
(1D)モデルで，気液二相スラグ流やチャーン流などにお
けるようなボイド率変動を数値的に実現するものとして
構想された 7)．
ところで気液二相流の流動様式は気液界面の形状の分
類はある瞬間の，いわば凍結した，幾何学的な形状を目
視やカメラなどを用いて同定するのが通例であった．一
方，Jones & Zuber 8) は，ボイド率変動の確率密度関数
（PDF）が流動様式に応じて特徴的な分布を示すことを
実験によって明らかにした．管内流であれば，流動様式
固有の空間的なボイド率分布を有する流れを，ある特定
の測定箇所で観察すればその流動様式に応じて固有の
PDF を示すのは当然ともいえる．つまり瞬間写真的な流
動様式同定も，Jones & Zuber 流の同定も本質的には同
じものと考えてよい．言い換えればボイド率変動やその
確率密度関数がどのようなパターンを形成するかを求め
ることは，とりもなおさず流動様式を同定することに対
応し，モデリングと実現象の物理的な対応関係が明確で
あれば，実現象のシミュレーションにも利用できるだろ
うというのが著者らの基本的なスタンスである．
さてパターン形成という視点から考えれば，典型的に
は周辺の要素間の干渉ルールに基づいて時間発展経過を
見るセル・オートマトン(CA)法や Coupled Map Lattice
(CML)などがよく知られているが，これらは実現象への
見通しが悪いという問題がある．一方，粒子系でよく用
いられている離散要素モデル(Discrete Element Method)
では粒子間の動的な干渉をバネとダッシュポットで表現
して粒子挙動の時間発展過程を追うという手法が取られ
ている．この手法では例えば粒子が壁に及ぼす力など物
理的に意味のある特性量を求めることが可能である．本
稿で説明する離散気泡モデルは物理的な意味を持たせな
がら，ボイド率変動から見たパターン形成過程を数値的
にトレースするとともに，圧力損失やボイド率変動，さ
らにはより大規模な不安定流動などのシミュレーション
にも利用できるようにしたいと考えており，その意味で
は粒子系の離散要素モデルと似たアプローチであるとも
いえる．
従来型の Drift-flux モデルや二流体モデルはたとえ非
定常状態を扱うにしても基本は時間平均化モデルで，ボ
イド率や流速を時空間すべての領域で定義できる必要が
あり，それによって流体力学的解析手法が適用可能とな
っている．本モデルは対象とする管内流を多数のセルに
分割し，それぞれに気泡（ボイド率）を定義することに
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
よって，離散化し，長大な気泡を気泡列で代表させ，
個々の気泡の挙動をトレースすることによって，気泡の
離合集散を可能にした．これが本モデルを離散気泡モデ
ルと名付けた所以でもある．この気泡の離散化によって
二相流に内在するボイド率の変動特性，ひいては流れ場
の局所瞬間的な変動を再現しようとするもので，時間平
均モデルとは立場が異なっている．また離散要素モデル
のように個々の粒子そのものをトレースするのではなく，
ここではボイド波の伝播をトレースするという点におい
て物質的（ラグランジュ的）ではなく，波動的である．
気液界面は複雑に変形し，その意味で多自由度系であ
り，決定論的に界面形状や界面の空間分布を求めること
は困難である．そこでアプローチの方法としては CA 法
と同様に，複雑な機構を前提とせず，少数の基本となる
機構のみを組み込んだ数値シミュレーションを行うこと
によって，つまり制約を減らしフレキシビリティを高く
することによって，複雑な現象をコンピュータ上で実現
しようというわけである．
2.
モデリング
二相流の基本的な機構をどのように与えるかには問題
なしとはしないが，ここでは以下のような考えでモデル
化を行った．まず質量保存は系全体として厳に守らなけ
ればならない要件（これを global rule と呼ぶ）とし，さ
らに二相流の本質として気液間の相対速度，先行する気
泡のウェーク効果，局所的な圧力変動に伴う気泡の膨
張・圧縮を考慮する．即ち，運動量保存については局所
的な関係式，すなわち local rule として計算に組み込む．
時間発展を追うことから質量保存式（実際にはガス側の
保存式で，ボイド波伝播式）は 1D の偏微分方程式の形
で与えるが，それ以外の運動量や力などの local rule に
は準定常関係式を用いる．
ここで用いるモデルは 1D モデルである．すなわち，
流れ方向に離散化した各セルにそこでのパラメータを積
分した結果を配置して定式化行うものである．これは通
常の水管ボイラや原子炉蒸気発生器など多くは 1D モデ
ルで十分であるという背景に準拠している．断面方向に
パラメータの積分を行うので空間分解能は当然，管直径
のオーダーになる．したがって比較的長い管路，すなわ
ち L/D が大きい管路を考える（例えば数 10 から 100 以
上）であれば，現象の不確実さから考えて，大きな問題
はない．
2.1
計算モデル
ここで代表として説明する流れ場は Fig. 1 に示すよ
うな垂直上昇等温二相流であり，垂直管下部より液相を
体積流束 JL で流入させ，空気混合部より気相を体積流束
JG で供給し，そこで生成された気液二相混合物はテスト
セクションを下流へ移流する系になっている．本モデル
− 85 −
小澤
守，網
Fig.1 Flow model and single bubble hypothesis.
では流れ場を Fig. 1 のような高さ y と内径 Dp が等しい
円柱の単一セルを最小単位と考え，セル内に存在する気
相を体積積分した代表長さが Db のセルと相似の仮想単一
気泡として定義している．これによりセル内でのボイド
率挙動やその分布は空間分解能 Dp のスケールで集中系近
似していることになる．この仮想単一気泡のボイド率は，
以下のように単純に表すことができる．
  ( Db / Dp )3
(1)
このような仮定に基づいて，各流動様式を模式的に表す
と Fig. 2 のように流路軸方向に離散化（気泡流に対して
は Dp のスケールで集中化）していることになる．なお
1D モデルであるにもかかわらず 3D 形状の気泡を考慮す
るところに矛盾するように思われるが，ボイド率の時空
間変化だけでは物理的イメージをつかみにくいこと，の
ちに述べる先行気泡のウェークや気泡の体積膨張・圧縮
効果を組み込むためには 3D の形状が必要であることな
どによる．
一般に 1D のボイド率挙動は，式（2）で表される気相
の質量保存式，いわゆるボイド波伝播方程式によって記
述される．
 U G 

q
t
y
健行
の移動量を表すことになる．その場合には当然ながら質
量流束の減少する液相，あるいは気相の質量流束の変化
を考慮しなければならない．
二相流の重要な特徴の一つである気液間の相対速度は
特に垂直流の場合には重要で，流体中の終端速度と同様
な関係で近似的に与えることができる．これは Nicklin
ら 9) や Zuber ら 10) のドリフト速度の考え方と同様であ
る．さらに気泡挙動特に気泡合体過程から，先行気泡の
ウェークの効果，つまり気泡がある程度先行気泡に接近
したら後続する気泡が加速・合体する観察結果を組み込
みたい．さらに時間平均的な挙動から見た二相流ではほ
とんど問題にならないが，局所的には気泡はその場の静
圧を受けてその体積を変化させているはずである．した
がって管長に沿う静圧分布を求めるための圧力降下計算
過程とそれによる気泡の体積変化を組み込む．これらは
気泡の運動に関する因子としてそれほど無理なく理解が
得られるであろう．
本離散気泡モデルを用いた計算は，当然ながら数値的
に時間発展を追うことになり，その計算過程において，
組み込んだ各種因子を同時的に解くのではなく，収束の
ための繰り返し計算過程などは使用しない．そこで例え
ばあるセルにおける離散気泡が体積変化によって膨張し，
セルの寸法を超えることもありうる．体積に関する整合
性を確保するために，そこで本モデルでは超過した体積
（ボイド率ならびに液相ホールドアップ）を対象とする
セルの前後に再配分する機構を組み込むことにする．こ
れらは上記ボイド伝播方程式が global rule であるのに対
して，局所的に整合させるので local rule と呼んでいる．
言い換えれば，質量保存については流れ場全体として扱
い，流れのダイナミクスに関与する運動量保存あるいは
相変化の場合ではエネルギー保存については局所的な作
用として準定常的関係で組み込むことになる．この点で
質量，運動量，エネルギーの保存式を連立して解いたも
のとは異なるが，基本と考えられる機構は含まれている．
本モデルにおける時間発展の計算手順は，まず（i）
先行気泡によるウェークの影響を考慮した流れ場を定義
し，仮想気泡を気相の非圧縮性を仮定した式（2）に基づ
いて伝播させたのち，（ii）気相の圧縮性の影響により気
泡を膨張・圧縮させ，そして（iii）流れ場の幾何学的寸
法とのマッチングを図る．
(2)
ここで，UG は気相の速度，q は気相の湧き出しに相当し，
気液混合部の長さを Lmx (  mx i  mx 0 ) とすれば，気液
混合部では q
 J G / Lmx と表され，気液混合部以外では
q  0 とする．なおこの式は沸騰二相流や凝縮二相流の
ように相変化を伴う系にも適用できる．その場合には湧
き出し項 q は液相から気相へ，あるいは気相から液相へ
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
− 86 −
Fig.2 Discrete bubble model (left: image of flow,
right: corresponding model).
気液二相流のパターンダイナミクスモデリング

uw,max,i  c1ur ,i1 CD Db,i1 y 2
2

13
uw,max,i  0
(  i 1
 i )
(  i 1
 i )
(6)
ここで，i は Fig. 3 に示すように，i 番目のセルでの値を
意味し，係数 c1=0.715 とする．複数の気泡が連続して配
置されている場合には，ウェークの影響を受けた流れ場
にある気泡によって，さらにウェークが発生することに
なるので，ここでは単純に先行気泡によるウェークの影
響は上流側に進むにつれて指数関数的に減衰すると仮定
し，対象気泡よりも下流側 c 3 y 内の気泡が作り出すウ
ェークが重ね合わさって，対象気泡に作用するとして，
ウェーク速度 uw,i を以下のように与える．
c3
u w,i  u w,max,i   u w,max,i j exp  c2  j  1
(7)
j 1
本シミュレーションにおいては感度解析結果などに基
づいて c2=0.14，c3=20 としている．一旦，ウェーク速度
が上述のように決定されれば，先行気泡のスリップ速度
はウェーク速度を考慮した形で，以下のように与えられ
る．
Fig. 3 Simulation model.

ur ,i1  uw,i c1 CD Db2,i1 y 2
2.2
流れ場
気相速度 UG は総体積流束 JT と気相相対速度 uG の和
として与えられる．
U G  JT  uG  J G  J L  uG


1 3 1
( u w,i
 0)
ur  [2 L  G G Db /(CD  L )]1/ 2 ( u w,i  0 )
(8)
2.4
(3)
気相の圧縮性
ここで，気相相対速度はボイド率と各相間のスリップ
速度 ur の関数として，以下のようになる．
気泡は移流前後の静圧の変化に応じて等温変化に従う
とすれば，以下の式が成り立つ．
(4)
(9)
uG  1   ur  1   U G  U L 
ここで，UL は液相速度である．
また, スリップ速度は単一気泡に作用する浮力と抗力
の釣り合いを考慮すれば以下のようになる．
ur  [2 L  G G Db /(CD  L )]1/ 2
(5)
Vb 2  Vb1 P1 P2 
Vb は気泡の体積，P は静圧であり，1，2 の添え字はそれ
ぞれ移流前，移流後の値であることを示している．なお
静圧の計算にあたっては，等温二相流では加速項による
寄与は非常に小さく，以下のように重力項と摩擦損失に
より与えられる．
ここで，抗力係数 CD=0.44 とする．
2.3
Pi  PABS    i  G ,i  1   i  L Gy 
GN
気泡によるウェーク
i
気泡の合体は流動様式の遷移を左右する重要な要因で
ある．先行気泡によるウェークは，先行気泡が後続気泡
よりも大きい場合に発生するとして，以下のように与え
る．
乱流域における最大ウェーク速度は (CD A / z
2 1/ 3
)
に
比例することが知られている
．ここで，A は物体正
面の投影面積， z は物体からの距離である．これを本モ
デルに適用して，二つの気泡間で発生する最大ウェーク
速度は以下のように与えることとする．
11)
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
D
GN
f
i
i
p
 Db,i  1
Dp
2
(10)
 LU L ,i U L ,i
ここで， PABS は系の代表圧力で管路の出口圧力と考
えてよい． f は管摩擦係数である．管摩擦係数は層流域
では Hagen-Poiseuille の式，乱流域では Blasius の式を
それぞれ使用している．セルの番号 i ， GN については
Fig. 5 を参照のこと．なお管路の入口，出口においては
連続流入，流出条件を置いている．
− 87 −
3
1
3
2
2
3
3
2
3
2
守，網
4
1
1
2
1
1
Air
injection
0
0
0.5
1
0
0
0.5
1
0
0
0.5
1
0
0
0.5
健行
1
0.8
0.6
0.4
0.2
0
PDF
Axial coordinate y m
小澤
1
JG0=0.02m/s
JL0=0.08m/s
1
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
1
0.8
0.6
0.4
0.2
0
3
JG0=5.00m/s
JL0=0.08m/s
0 0.2 0.4 0.6 0.8 1
1
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
1
Void fraction 
本モデルではある限界ボイド率  cr を定義し，対象気
泡がその限界ボイド率を超えた場合，その超過したセル i
  i   cr を下流側へ再配分し，それ
に伴って押し出されることになる余剰の液ホールドアッ
プ  L,i の一部が上流側のセル群に一様に再配分する方
法をとっている．つまり各セルへの配分は
 L,i  { i exp( c4 cr )} / i
(11)
ここで c4=0.5 とした．
液相と気相を置き換える操作において，超過した液ホ
ールドアップがそっくりそのまま上流のセルに負わせる
のではなく，多くは下流に押し出され，残りの一部を全
体で少しずつ負担させている．このことは比較的長いス
ラグ気泡の観察結果とも整合している．上流側に位置す
る管路の入口でポンプなどによって流速が固定されてい
るのが通例で，自由に流出する下流側に向かって気泡の
膨張があるのは当然の成り行きであろう．なおここでい
う限界ボイド率は環状流の液膜厚さのオーダーを考慮し
て，  cr =0.95 としている．
3. 定常等温二相流における流動特性
以上説明した離散気泡モデルを用いて，様々な条件の
もとで数値シミュレーションを行った．ここでは空気－
水二相流の 1 例を示しておこう．
Fig. 4(1) は，計算開始から十分に時間が経ったある時
刻におけるテストセクション軸方向のボイド率分布を示
しており，図中に示す上流部の 2 本の水平線は空気供給
部に相当する．計算条件は図中に示すとおりで，
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
0 0.2 0.4 0.6 0.8 1
(2) PDF
JL0=0.08m/s に対し JG0 を 0.02, 0.8, 5.0, 20m/s と変化
管の幾何学的寸法とのマッチング
のボイド率  i
4
JG0=20.0m/s
JL0=0.08m/s
Void fraction 
(1) Snap shot
Fig. 4 Void fraction fluctuation.
2.5
2
JG0=0.80m/s
JL0=0.08m/s
させている．一方，Fig. 4(2) は Fig. 4(1)の各条件におけ
る y=1.65m(LPDF)の位置におけるボイド率の時間変動デ
ー タ か ら 計 算 し た 確 率 密 度関数（Probability Density
Function: PDF）を示している．
計算結果を順に見ていくと，気相体積流束を JG0=0.02
m/s とした Fig. 4(1)-1 の場合，小さな三角形状のボイド
率分布（ボイド波）がテストセクション軸方向にまばら
に存在していることがわかる．このとき，Fig. 4(2)-1 を
見ると，低ボイド率側に単一ピークが見られる．
次に，気相体積流束を JG0=0.8m/s に増加させた Fig.
4(1)-2 の場合，限界ボイド率に達する，つまりテイラー
気泡に相当すると考えられる短矩形状のボイド波が形成
されている．その一方で，ボイド率が 0 となる液相のみ
の部分も見られる．このとき，Fig. 4(2)-2 を見ると，低
ボイド率側と高ボイド率側の双峰にピークが見られ，低
ボイド率側のピークが顕著である．
さらに気相体積流束を JG0=5.0m/s に増加させた Fig. 4(1)-3
の場合，さきほどよりも長い区間，限界ボイド率を保持
する長矩形状の波形がテストセクション軸方向にほぼ規
則的に存在しており，その長矩形波の前後には低ボイド
率へとなだらかに遷移する区間が存在する．このとき，
PDF は低ボイド率側と高ボイド率側の双峰にピークが見
られ，高ボイド率側のピークが顕著である．
最後に，気相体積流束が JG0=20m/s と最も大きい
Fig. 4(1)-4 の場合，上述の条件で認められた低ボイド率
となる区間は存在せず，限界ボイド率で保持される一様
なボイド率分布のみが認められる．このとき，Fig. 4(2)-4 を
見ると，高ボイド率側のみに単一ピークが見られる．
各条件における流動様式の判別に際しては，本研究で
は，Jones-Zuber 8) による各流動様式に対応するボイド
率の PDF の特徴を条件とするのに加えてボイド率変動
− 88 −
気液二相流のパターンダイナミクスモデリング
Calculation results
Bubbly
Slug
Churn
Annular
Flow pattern transitions
Mishima et al.
Taitel et al.
10
Air-Water 0.1MPa
JL0 m/s
1
0.1
Bubbly
0.01
0.01
Slug
Churn
0.1
1
Annular
10
100
JG0 m/s
Fig.5 Flow pattern map.
波形の最大値，最小値にしきい値を設定した．
気泡流：低ボイド率側に単一ピーク
気泡流－スラグ流遷移：ボイド率変動波形のボイド率の
最大値が 0.6 を超えたとき，スラグ流と判定する．
スラグ流：高ボイド率・低ボイド率側の双峰ピーク，か
つ低ボイド側のピークが顕著であること．
スラグ流－チャーン流遷移：ボイド率の最小値が 0.1 以
上であれば，チャーン流と判定する．
チャーン流：高ボイド率・低ボイド率側の双峰ピーク，
かつ高ボイド側のピークが顕著であること．
チャーン流－環状流遷移：ボイド率の最小値が 0.5 より
大きければ，環状流と判定する．
環状流：高ボイド率側に単一ピークを有すること．
上記の判定条件を用いて流動様式を同定し，JG0－JL0
グラフ上にプロットした結果を Fig. 5 に示し，同時に
Mishima & Ishii 12)（実線），Taitel ら 13)（破線）による
30
流動様式遷移線を示す．Fig. 5 に示した流動様式判定結
果は，既存の遷移境界モデルの結果とよく整合している
というより，実際にはいくつかの local rule の定数・係
数を調整して Mishima & Ishii の遷移境界にほぼ一致す
るように設定したという方が正確であろう．なお液相体
積流束が大きい領域では Mishima & Ishii によればスラ
グ流となる領域において，本シミュレーションではチャ
ーン流と判定している．定数や係数を調整して整合性を
高めることは可能ではあるが，それによって他の領域で
予測性能が低下する問題があり，ここで採用している定
数・係数は最適とは言わないまでも許容できる程度に調
整した結果である．なおこの領域は Taitel らによればチ
ャーン流に相当し，二つの予測方法で流動様式が一致し
ないこと，実際の観察においてもこの領域では流動様式
の同定が難しいこと，そもそもチャーン流とスラグ流の
区別に任意性が残っていることなどを考えれば，シミュ
レーション結果と Mishima & Ishii 線とのずれについて
議論する必要もないだろう．
以上で説明した離散気泡モデルの妥当性は実験観察結
果との比較によって検証される必要があるが，汎用性の
非常に高いと定評のある Mishima & Ishii の流動様式遷
移境界とほとんど一致させていることからして，本モデ
ルの有用性はかなり高いと考えている．なおこのモデル
を水平流に拡張し，CO2 高圧二相流の流動様式との比較
検討によってその妥当性を確認しているので，詳細につ
いては参考文献 14-19) を参照されたい．
本モデルはボイド率の時空間変化におけるパターン形
成時間発展をトレースするのが本来の趣旨であるが，計
算過程でパラメータの時間平均を求めれば，既存の相関
式などと比較検討可能である．ここではボイド率と圧力
損失の時間平均について述べる．ボイド率の時間平均特
性として，総体積流束 J T と気相の流速 J G /  の関係を
Fig. 6 に示す．離散気泡モデルによる結果は流動様式ご
とに記号を変えてプロットし，図中実線および破線はそ
れぞれスラグ流および環状流に対する相関式 10, 20) を示
Simulation
Bubbly
Slug
Churn
Annular
102
Simulation
Bubbly
Slug
Churn
Annular
Correlation
Chisholm
JG0/ m/s
20
JL0=0.03m/s
JL0=0.8m/s
Ltt
JL0=3.0m/s
Annular flow
10
Slug flow
0
0
Drift-flux model
Slug flow
Annular flow
10
20
100
30
10-2
JT m/s
10-1
100
101
102
103
104
tt
Fig. 6 Time-averaged void fraction.
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
101
Fig. 7 Time-averaged pressure drop (friction).
− 89 −
小澤
守，網
している．離散気泡モデルは低体積流束側で多少ではあ
るいが相関式よりも高い値となっていることを除けば，
全体的にはかなりうまく再現できている．
次 に ， 圧 力 損 失 の 時 間 平 均 特 性 と して，LockhartMartinelli 相関（L-M 相関）21) の方法に準拠して，計算
結果から得られる摩擦損失の時間平均値を整理した結果
を Fig. 7 に示す．離散気泡モデルによるシミュレーショ
ン結果は流動様式ごとにプロットされ，Chisholm 22) に
よる相関式は実線で示されている．ここで Ltt と  tt は
それぞれ以下のように定義している．
 P   P 
 Ltt 2   f   f 
 y   y  L
 Pf 
2

 tt  
 y  L
 Pf

 y


G
,
(12)
チャーン流および環状流領域において計算結果の分散
が若干大きいものの，全体的には L-M 相関と良く一致し
ている．ボイド率にしても圧力損失にしても本モデルに
は Drift-flux モデルや L-M 相関を組み込んでいるわけで
はないにも関わらず，既存の代表的な相関式との整合性
が高いということは，本モデルのポテンシャルの高さの
証左でもある．
さらに本モデルを用いた時空間発展の計算例と実験に
おける脈動流から得られた結果の定性的な比較を Fig. 8
に示しておく．左の計算結果は，右図の実験結果に合わ
せて得られたボイド率分布から模式的に気泡分布に書き
直したもので，空間座標は実験に併せてある．両者比較
してみると本離散気泡モデルは実現象と非常に良く対応
している．つまり本離散気泡モデルはボイド率変動のパ
ターン形成にとどまらず，二相流の挙動そのものを定常
状態も含めてかなり上手く再現できることがわかる．
また Chisholm による相関式は液相および気相が乱流で
あるとした以下の式を使用した．
Ltt 2  1 
21
1

 tt  tt 2
(13)
健行
4.
離散気泡モデルの展開
本稿では垂直管内等温二相流を例として離散気泡モデ
ルについて説明した．このモデルは垂直流であるため重
力項が支配的であったため問題にならなかったが，水平
Fig. 8 Void fraction fluctuation in oscillatory two-phase flow (left: simulation, right: experiment)
( J L0
 0.03 m/s, J L  0.33 m/s, J G 0  0.12 m/s , oscillation period=2.0 s).
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
− 90 −
気液二相流のパターンダイナミクスモデリング
流では気泡運動に伴う各セル間の圧力のジャンプ条件，
成層化に伴う静圧変化，さらに水平管内における気泡の
ドリフト速度 23), 24) を導入すれば容易に水平流に展開で
きる 16, 17)．また 0.5mm 程度のミニチャンネルへの展開
も基本的には気液間のスリップ相関の変更によって可能
である 18,19)．モデルの詳細についてはそれら文献を参照
されたい．さらには本離散気泡モデルを用いて水平沸騰
流の CHF 挙動，水－ケロシンによる液－液系の流動様
式にもかなりうまく適用できることがわかっている．つ
まり本モデルは完全とは言わないが，ボイド率の時空間
変動のパターン形成を定量的にも定性的にもうまく再現
できるモデルとなっている．1D モデルであり，運動量保
存やエネルギー保存関係を質量保存式と同時に解いてい
るわけではないが，ここで選定した少数の基本メカニズ
ムが二相流に内在する（ボイド率の）変動特性を再現す
るうえでかなりの確度で有効であったともいえる．この
ことは多分野の多くの複雑系に対するアプローチに対応
するが，しかし「ではなぜそのようなパターンになるの
か」という疑問には依然として明確な回答を持ち合わせ
ていない．
参考文献
1)
2)
3)
4)
5)
6)
7)
8)
9)
The Babcock & Wilcox Co.: Steam - Its Generation and
Use, 22nd ed., New York, 1890.
The Babcock & Wilcox Co.: Steam - Its Generation and
Use, 8th British ed., London, 1913.
F. Münzinger, Dampfkraft-Wasserrohrkessel und Dampfkraftanlagen, Verlag von Julius Springer, Berlin (1933).
S. W. Gouse, Jr.: An Index to the Two-Phase Gas-Liquid
Flow Literature, MIT Report No.9, The MIT Press,
Massachusetts, 1966.
赤川浩爾: 気液二相流研究史と関連技術―気液二相流研究
史の私観―，学術出版印刷, 2007．
石谷清幹：気液二相流研究の現状，日本機械学会誌，66528 (1963) 58-66.
伊藤浩二，井上満，小澤守，庄司正弘：気液二相流の流動
様 式 の 簡 易 遷 移 モ デ ル ， 日 本 機 械 学 会 論 文 集 ， 70-689
(2004) 151-158.
O. C. Jones Jr. and N. Zuber: “The Interrelation between
Void Fraction Fluctuations and Flow Patterns in TwoPhase Flow”, Int. J. Multiphase Flow, 2 (1975) 273-306.
D. J. Nicklin, J. O. Wilke and J. F. Davidson: “Two-Phase
Flow in Vertical Tubes”, Trans. Inst. Chem. Engr., 40-1
(1962) 61-68.
J. Jpn. Soc. Microgravity Appl. Vol. 29 No. 2 2012
10) N. Zuber, F. W. Staub, G. Bijwaard and P. G. Kroeger:
“Steady State and Transient Void Fraction in Two Phase
Flow Systems”, GEAP-5417, 1967.
11) H. Schlichting: Boundary Layer Theory, 6th ed.,
McGraw-Hill, New York, 1968, 685-695.
12) K. Mishima and M. Ishii: “Flow Regime Transition
Criteria for Upward Two-Phase Flow in Vertical Tubes”,
Int. J. Heat Mass Transfer, 27-5(1984) 723-737.
13) Y. Taitel, D. Barnea and A. E. Dukler: “Modeling Flow
Pattern Transitions for Steady Upward Gas-Liquid Flow
in Vertical Tubes”, AIChE J., 26-3(1980) 345-354.
14) M. Ozawa, T. Ami, H. Umekawa and M. Shoji: “Pattern
Dynamics Simulation of Void Wave Propagation”,
Multiphase Science and Technology, 19-4(2007) 343-361.
15) 網健行，梅川尚嗣，小澤守，庄司正弘：“離散気泡モデル
による二相流ダイナミクスの解明”, Thermal Science and
Engineering, 15-4(2007) 197-209.
16) T. Ami, H. Umekawa, M. Ozawa and M. Shoji: “Traveling
Void Wave in Horizontal Two-Phase Flow”, Int. J. Heat
Mass Transfer, 52 (2009), 5682-5690.
17) M. Ozawa, T. Ami, I. Ishihara, H. Umekawa, R.
Matsumoto, Y. Tanaka, T. Yamamoto and Y. Ueda: “Flow
Pattern and Boiling Heat Transfer of CO2 in Horizontal
Small-Bore Tubes”, Int. J. Multiphase Flow, 35-8(2009)
699-709
18) M. Ozawa, H. Umekawa, T. Ami, R. Matsumoto and T.
Hara: “Flow Boiling in Horizontal Minichannels: Flow
Pattern of CO2 at High Pressure”, Multiphase Science
and Technology, 22-2 (2010) 115–132.
19) M. Ozawa, T. Ami, H. Umekawa, R. Matsumoto and T.
Hara: “Forced Flow Boiling of Carbon Dioxide in
Horizontal Mini-Channel”, Int. J. Thermal Sciences, 50
(2011) 296-308.
20) 日本機械学会編，改訂 気液二相流技術ハンドブック, コロ
ナ社，2006.
21) R. W. Lockhart and R. C. Martinelli: “Proposed
Correlation of Data for Isothermal Two-Phase, TwoComponent Flow in Pipes”, Chem. Eng. Prog., 45-1(1949)
39-48.
22) D. Chisholm: “A Theoretical Basis for the LockhartMartinelli Correlation for Two-Phase Flow”, Int. J. Heat
Mass Transfer, 10 (1967) 1767-1778.
23) M. E. Weber: “Drift in Intermittent Two-Phase Flow in
Horizontal Pipe”, Can. J. Chem. Eng., 59 (1981) 398-399.
24) T. Sakaguchi, M. Ozawa, H. Hamaguchi and T.
Fukunaga: “Behavior of a Large Bubble in a Horizontal
Channel, 2nd Report: Large Bubble Penetrating into
Running Liquid”, Trans JSME, Ser. B, 56-527 (1990)
1891-1898.
− 91 −
(2012 年 2 月 29 日受理)