Thermal Science & Engineering Vol.15 No.4 (2007)

離散気泡モデルによる二相流ダイナミクスの解明

網 健行†

梅川 尚嗣†

小澤 守†

庄司 正弘‡

Investigation on Two-phase Flow Dynamics with Discrete Bubble Model
Takeyuki AMI†, Hisashi UMEKAWA†, Mamoru OZAWA† and Masahiro SHOJI‡

Abstract
Conventional modeling including drift-flux model and two-fluid model is based on “continuous flow hypothesis”,
being constructed by time-averaging, and thus both phases are defined in every spatio-temporal space. This makes it
possible to apply to a variety of two-phase flow dynamics, while the intrinsic void fraction fluctuations, typicallyobserved in slug and churn flows, are hardly simulated. In order to break through such a problem caused bytime-averaging, discrete bubble model based on one-dimensional mass conservation equation, i.e. void propagation
equation, has been developed. This model takes into account, as momentum effects, the wake effect induced bypreceding bubbles, the local pressure fluctuation and the compressibility of gas phase together with the phasere-distribution due to geometrical constrains. Thus obtained spatio-temporal fluctuation characteristics of void fraction
well simulated inherent two-phase behavior not only in a steady flow but also in an oscillatory flow.
Key Words: Two-phase flow dynamics, Discrete bubble model, Void wave propagation, Flow pattern,
Probability density function, Oscillatory two-phase flow

記 号
A
: 面積
c
: 定数
CD : 抗力係数
C0 : 分布パラメータDb : 気泡直径
Dp : 管直径f: 管摩擦係数g: 重力加速度GN : メッシュ数
J: 体積流束JT: 総体積流束L: 長さ
Lmx : 空気混合部長さ
Lp: 圧力測定区間
LT: テストセクション全長
mxi : 空気混合部入口座標
mxo : 空気混合部出口座標

[m2]
[-]
[-]
[-]
[m]
[m]
[-]
[m/s2]
[-]
[m/s]
[m/s]
[m]
[m]
[m]
[m]
[m]
[m]

P
q
t
U
u
ur
uw
vgj
Vb
y
z

α
αcr
∆J
∆t
∆y
µ

:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:

圧力
湧き出し
時間
相速度
相対速度
スリップ速度
ウェーク速度
ドリフト速度
気泡体積
距離
距離
ボイド率
限界ボイド率
脈動振幅
タイムステップ
単位長さ
静粘性係数

[Pa]
[1/s]
[s]
[m/s]
[m/s]
[m/s]
[m/s]
[m/s]
[m3]
[m]
[m]
[-]
[-]
[m/s]
[s]
[m]
[Pa・s]

受付日: 2007 年 2 月 2 日, 担当エディター: 門出 政則
関西大学大学院 工学研究科 (〒564-8680 大阪府吹田市山手町 3-3-35)‡
神奈川大学 (〒221-8686 神奈川県横浜市神奈川区六角橋 3-27-1)†- 197 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

Ltt

:
:
:
:

密度
脈動周期
Lockhart-Martinelli パラメータ
二相摩擦乗数

添字
0
G
L
max

:
:
:
:

定常値
気相
液相
最大

tt

１

緒

[kg/m3]
[s]
[-]
[-]

言

気液二相流は気相・液相が混在した流れで，そこで形成される気液界面は時空間にわたって大きく変動しながら変形・合体・分裂しており，その結果形成される幾何学的構造のパターンである各流動様式
において，気相・液相の配置がそれぞれ異なるため，
その熱伝達特性は大きく異なることになる．このような観点に基づき，従来から，原子炉やボイラなどの沸騰関連機器において設計基準として使用されている計算コードでは計算精度を向上させるために，
まず流動様式を同定し，次にその流動様式に対応す
る熱伝達を含む各種相関式を使用する方法がとられている．その流動様式の同定に際しては定常状態あるいは，時間的に平均化された相関式・構成式に基づいているため，気液二相流に本質的に内在する変動特性やそれと同程度の時間スケールの過渡現象には上記計算コードでは対応できないことになる．
特に，近年のボイラの小型化・高性能化に伴い，高性能小型ボイラは低圧・低質量流束条件下にも関わらず，高熱流束条件下で運転される傾向にある．
このような場合，
伝熱面配置や収熱量分布によって，流量の不均一分配が発生したり，気液密度比が大き
いため密度波振動に代表される流動脈動を伴うような不安定流動を発生しやすい状態にある．しかもそれは水管長さが短いことにより，スラグ流に代表されるような流れの変動時間スケールと同オーダーの同期を有する特徴を持つ．
このような不安定流動が発生すれば，必然的に気泡の停滞・逆流が発生する．大気泡が加熱壁面を停滞気味に通過する場合には[1,2]，低クオリティ領域においても気泡周囲の液膜がドライアウトして限界
熱流束（Critical Heat Flux : CHF）に至ることが知られている[3]．このような二相流が元来有している時
間変動に強く依存するクリティカルな問題は，従来から用いられている時間平均化された，すなわち特にスラグ流やチャーン流において時空間にわたって断続するような流れを，ボイド率が 0.3～0.5 の滑ら
かなボイド率分布を有する流れとみなしてしまうよ
うなモデルでは評価・予測しきれないと考えるのが自然である．
本一連の研究は，著者ら[2,4,5]が開発した離散気泡モデルをさらに発展させることで，二相流特有の変動特性が原因となる沸騰二相流の限界熱流束を評
価・予測することを最終目標としたものであり，本報ではその第一段階として，流動脈動を含めた等温二相流の時間変動特性を数値的に実現したので報告する．なおここで用いた離散気泡モデルは，複雑系に対するシミュレーション手法でよく用いられる
CA 法のようにメタファモデルを用いるのではなく，できるだけ物理的なイメージを保持した形でパター
ンダイナミクス的アプローチとして開発されたものである．
２

Fig. 2 Single bubble hypothesis.

Fig. 1 Flow model.

離散気泡モデルモデルの詳細は伊藤ら[4]，Ito ら[5]，Ozawa ら[2]に詳しいが，読者の便宜を考えて，ここで簡単に説明しておく．
2.1 計算モデル
本研究で対象とする流れ場は Fig. 1 に示すような垂直上昇等温二相流であり，垂直管下部より液相を体積流束 JL で流入させ，空気混合部より気相を体積流束 JG で供給し，そこで生成された気液二相混合物

Fig. 3 Flow configurations.
Phenomenon (left) and corresponding model (right)- 198 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

はテストセクションを下流へ移流する系になっている．
離散気泡モデルでは流れ場を Fig. 2 のような高さ
y と内径 Dp が等しい円柱の単一セルを最小単位と考え，セル内に存在する気相を体積積分した代表長
さが Db のセルと相似の仮想単一気泡として定義している．これにより，セル内でのボイド挙動やその分布は空間分解能 Dp のスケールで集中系近似していることになる．この仮想単一気泡のボイド率 は，
以下のように単純に表すことができる．

Db
Dp

3

2.2 気液間の相対速度
まず，基本となる流れ場について述べる．気相速
度 UG は総体積流束 JT と気相相対速度 uG の和として与えられる．

(1)

UG

このような仮定に基づいて，各流動様式を模式的に表すと Fig. 3 のように流路軸方向に離散化（気泡流に対しては Dp のスケールで集中化）していることになる．
一般に一次元流れのボイド率挙動は，式(2)で表される気相の一次元質量保存式，いわゆるボイド波伝
播方程式によって記述される．

t

行気泡によるウェークの影響」を考慮した流れ場を定義し，仮想気泡を気相の非圧縮性を仮定した式(2)に基づいて伝播させたのち，（Ⅱ）
「気相の圧縮性の影響」により気泡を膨張させる．そして（Ⅲ）「幾何学的制限による相再配分機構」により，相の再配分を行うことで気泡の成長を表現している．以下において，これらの機構を簡単に説明する．

UG
y

q

uG

JG

JL

uG

(3)

ここで，気相相対速度はボイド率と各相間のスリッ
プ速度 ur の関数として，以下のようになる．

uG

1

ur

1

UG

UL

(4)

ここで，UL は液相速度である．
また, スリップ速度は単一気泡に作用する浮力と抗力の釣り合いを考慮すれば以下のようになる．

(2)

2

ur

ここで，UG は気相の速度，q は気相の湧き出しに相当し，空気混合部の長さを Lmx(= mxi - mxo)とすれば，空気混合部では q=JG/Lmx と表され，空気混合部以外
では q=0 とする．
本離散気泡モデルでは上述のボイド波伝播式を
グローバルルール（全流れ系を通して一義的に適用される）として用いる．気液間の相対速度，気泡の膨張や収縮，合体・分裂，圧力損失等運動量保存に関わる特性を組み込まないと，有意な流れの変動を再現することができないのは当然である．これらは次節以降に述べる気液間の相対速度及び最も基本的と考えられる三つの機構を，上記グローバルルール
に対して，ローカルルールとして組み込む．その計算概要は，例えば静圧低下に伴い気泡が膨張する場
合については，Fig. 4 に示すように，まず，（Ⅰ）
「先Fig. 4 Sequence of simulation.

JT

L

G

CD

gDb

(5)

L

ここで，抗力係数 CD=0.44 とする．
2.3 先行気泡によるウェークの影響
二つの気泡が縦に並んで上昇している場合，後続気泡が先行気泡のウェークの影響を受けて，先行気泡に追いつき，合体することが知られており，この現象は流動様式の遷移を左右する重要な要因となる．
ここではこの先行気泡によるウェークは，先行気泡
が後続気泡よりも大きい場合に発生するとして，以下のように与える．
乱流域における最大ウェーク速度は（CDA /z2)1/3に比例することが知られている[6]．ここで，A は物
体正面の投影面積，z は物体からの距離である．これを本モデルに適用して，二つの気泡間で発生する
最大ウェーク速度は以下のように与えることとする．

u w, max,ic1u r ,i 1 C D Db,i 12u w, max,i

0

y2

13

(
(

i 1

i)

i 1

(6)
)
i

ここで，i は Fig. 5 に示したような i 番目のセルでの値を意味し，係数 c1=0.715 とする．
複数の気泡が連続して配置されている場合には，ウェークの影響を受けた流れ場にある気泡によって，さらにウェークが発生することになるので，ここで- 199 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

は単純に先行気泡によるウェークの影響は上流側に進むにつれて指数関数的に減衰すると仮定し，対象気泡よりも下流側 c3 y 内の気泡が作り出すウェー
クが重ね合わさって，対象気泡に作用するとして，
ウェーク速度 uw,i を以下のように与える．
c3

u w, iu w,max, i

u w,max, i j exp

c2 j 1(7)

j 1

ここで，本シミュレーションでは，係数 c2=0.14，
c3=20 とする．
そして一旦，ウェーク速度が上述のように決定されれば，先行気泡のスリップ速度はウェーク速度を考慮した形で，以下のように与えられる．

u r ,i 1

u w,i c1 C D Db2,i 1

y2

13

過したボイド率 i=( I - cr)を下流側へ再配分し，それに伴って押し出された余剰の液ホールドアップ
L,i を上流側で均一に再配分することでこの問題を解決している．各セルに均一に再配分される液ホ
ールドアップ L,i は， i に相当する液ホールドア
ップの一部が上流方向に再配分されるとして，その割合を cr の関数であると仮定して，以下のように与
える．
L ,i

L ,i

i

i exp

c 4 cri(11)

ここで，本モデルでは係数 c4=5 とし，限界ボイド率は環状流の液膜厚さを考慮して， cr=0.95 としている．

1

( u w, i

0)
(8)

2

u r ,i 1

G gDb ,i 1LCD

( u w, i0)L

2.4 気相の圧縮性の影響
気泡は移流前後の静圧の比に基づき膨張しながら流動しており，いま対象としている系では等温膨張しているとして，以下の式が成り立つ．

Vb 2

Vb1 P1 P2(9)

ここで Vb は気泡の体積，P は静圧であり，1，2 の添
え字はそれぞれ移流前，
移流後を示している．
また，
静圧は，等温二相流で加速項による寄与は非常に小さく，以下のように重力項と摩擦損失により与えら
れる．
GN

Pi

PABS

i

G ,i

1

i

L g

y

i
GN

fi
i

Dp

Db ,i 1
Dp
2

Fig. 5 Simulation model.

(10)
LU L ,i U L ,i

ここで，PABS は系の代表圧力で，ここで対象とする系では大気圧 0.1013 MPa と設定しており，f は管摩擦係数である．管摩擦係数は層流域ではHagen-Poiseuille の式，乱流域では Blasius の式をそれぞれ使用している．なお i，GN については Fig. 5
参照のこと．
2.5 幾何学的制限による相再配分機構
気相の圧縮性に基づき気泡を膨張させる際，管直径を超えて気泡が膨張することが計算上起こる．そこで本モデルでは，ある限界ボイド率 cr を定義し，対象気泡がその限界ボイド率を超えた場合，その超

この液相再配分は，例えばスラグ流中の大気泡を観察すれば分かりやすい．大気泡は下流に存在する液相を押しのけながら移流し，押しのけられた液相
は大気泡周りを流下液膜のように上流に下降していることを考慮すれば，この再配分機構は実現象にもある程度整合したものである．
2.6 計算条件
数値計算はボイド伝播方程式を有限差分近似化した差分方程式を使用しており，その差分モデルをFig. 5 に示す．計算メッシュにはスタガード・メッ
シュを使用し，時間項は前進差分，対流項は風上差分とした．計算条件は表 1，2 に示すとおりで，JG0
は時間平均気相体積流束，JL0 は時間平均液相体積流束である．計算対象の幾何学的条件は後述の実験結- 200 -C 2007 The Heat Transfer Society of Japan
○

3

1

3

2

2

3

3

2

4

3

2

1
0.8
0.6
0.4
0.2
0

2

PDF

Axial coordinate y m

Thermal Science & Engineering Vol.15 No.4 (2007)

1

1

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

1

1
0.8
0.6
0.4
0.2
0

1
JG0=0.02m/s
JL0=0.08m/s

0 0.2 0.4 0.6 0.8 13
JG0=5.00m/s
JL0=0.08m/s

1
0.8
0.6
0.4
0.2
0

0 0.2 0.4 0.6 0.8 12
JG0=0.80m/s
JL0=0.08m/s

0 0.2 0.4 0.6 0.8 1

4
JG0=20.0m/s
JL0=0.08m/s

0 0.2 0.4 0.6 0.8 1

Void fraction

Void fraction(1) Snap shot

(2) PDF
Fig. 6 Characteristics of void fraction.
Table 3 System configurations.
LT [mm]3000
mxi～mxo [mm]150～300
Lp [mm]
1110～2100
GN [-]
200
Dp [mm]15
15
y [mm]
0.05
t [ms]

果と比較するため表 3 に示すとおりとした．
ここで，
LT はテストセクション全長，Lp は圧力測定区間，GNはテストセクションのメッシュ数，計算メッシュ長さ y は管内径 Dp と等しく，タイムステップ t はクーラン条件を十分に満たすものである．また， JLと は既報で述べた気液二相流動特性実験[7]における液相体積流束の振幅と周期を示している．
Table 1 Simulation conditions.
Steady flow
Air‐Water
System pressure [MPa]0.1
JG0 [m/s]
0.020～20.0
JL0 [m/s]
0.030～3.00(Probability Density Function: PDF)を示している．
計 算 結 果 を 順 に 見 て い くと ， 気 相 体 積 流 束 をJG0=0.02 m/s とした Fig. 6(1)-1 の場合，小さな三角形状のボイド率分布（ボイド波）がテストセクション
軸方向にまばらに存在していることがわかる．このとき，Fig. 6(2)-1 を見ると，低ボイド率側に単一ピークが存在する．
次に，気相体積流束を JG0=0.8 m/s に増加させた

Table 2 Simulation conditions.
Oscillatory flow
Air‐Water
System pressure [MPa]0.1
JG0 [m/s]
0.001～20.00
JL0 [m/s]
0.010～0.090
0.327
0.164 0.109
JL [m/s]2
4
6s]Calculation results
BubblySlug
Churn
AnnularFlow pattern transitions
Mishima et al.
Taitel et al.

10

Air-Water 0.1MPa

計算結果及び考察

3.1 定常等温二相流における流動特性(1)ボイド率変動挙動および流動様式線図
Fig. 6(1)は初期状態より十分に時間が経った，
ある瞬間におけるテストセクション軸方向のボイド率分布を示しており，図中に示す上流部の 2 本の水平線
は空気供給部に相当する．計算条件は図中に示すとおりで，JL0=0.08 m/s に対し JG0 を 0.02, 0.8, 5.0, 20 m/sと変化させている．一方，Fig. 6 (2)は Fig. 6 (1)の各条件における y=1.65 m(LPDF)の位置におけるボイド率の時間変動データから計算した確率密度関数- 201 -

1

JL0 m/s

３

0.1

0.01
0.01

Bubbly

Slug
0.1

Churn
1

Annular10

JG0 m/sFig. 7 Flow pattern map.
C 2007 The Heat Transfer Society of Japan○100

Thermal Science & Engineering Vol.15 No.4 (2007)

スラグ流：PDF において高ボイド率・低ボイド率側の双峰ピーク，かつ低ボイド側のピークが顕著
スラグ流－チャーン流遷移：PDF のボイド率の最小値が 0.1 以上であれば，チャーン流と判定する．
チャーン流：PDF において高ボイド率・低ボイド率側の双峰ピーク，かつ高ボイド側のピークが顕著
チャーン流－環状流遷移：PDF のボイド率の最小値が 0.5 より大きければ，環状流と判定する．
環状流：PDF において高ボイド率側に単一ピーク
このような条件のもとに各流動様式を同定し，JG0－JL0 グラフ上にプロットしたものが Fig. 7 で，同時に Mishima ら[9]（実線），Taitel ら[10]（破線）による流動様式遷移線を示す．これらの遷移境界はこ
こでのシミュレーションで同定された流動様式とよく一致していることがわかる．液相体積流束が大きい領域で Mishima らによるとスラグ流となる領域において，本シミュレーションではチャーン流となる
結果が得られた．この領域は Taitel らによればチャーン流に相当する．よって，この領域では二つの理論線で流動様式が一致しないこと，実際の実験においてもこの領域の流動様式を明確には判断できないことを考慮したとしても，本モデルは妥当な結果を与える．
ここで，本シミュレーションに用いた係数や定数について若干コメントしておく．まず，抗力係数 CDについては TRAC コードなどで粒子レイノルズ数が大きい場合の値 CD=0.44 を用いた．ウェーク速度に関係する係数 c1～c3 は気泡流，スラグ流，チャーン流のパターンに主として影響し，相の再配分に関わる係数 c4 は気泡流以外の流動様式に影響を与える．

Axial coordinate y m

Fig. 6(1)-2 の場合，限界ボイド率に達する，つまりテイラー気泡に相当すると考えられる短矩形状のボ
イド波が形成されている．その一方で，ボイド率が0 となる液相のみの部分もある．
このとき，Fig. 6(2)-2を見ると，低ボイド率側と高ボイド率側の双峰にピークが認められ，低ボイド率側のピークが顕著である．
さらに気相体積流束を JG0=5.0 m/s に増加させた
Fig. 6(1)-3 の場合，さきほどよりも長い区間，限界ボイド率を保持する長矩形状の波形がテストセクシ
ョン軸方向にほぼ規則的に存在しており，その長矩
形波の前後には低ボイド率へとなだらかに遷移する区間が存在する．このとき，Fig. 6(2)-3 を見ると，低ボイド率側と高ボイド率側の双峰にピークがあり，高ボイド率側のピークが顕著である．
最後に，気相体積流束が JG0=20 m/s と最も大きい
Fig. 6(1)-4 の場合，上述の条件で認められた低ボイド率となる区間は存在せず，限界ボイド率で保持される一様なボイド率分布のみが認められる．このときは，Fig. 6(2)-4 に示すように，高ボイド率側のみに単一ピークが存在する．
これらの結果より，各条件における流動様式を判別するが，本研究では，Jones-Zuber[8]などによって示された流動様式ごとの PDF の特徴に加え,ボイド
率変動波形の最大値，最小値にしきい値を設定し，以下の条件を用いて判別した．
気泡流：PDF において低ボイド率側に単一ピーク
気泡流－スラグ流遷移：ボイド率変動波形のボイド率の最大値が 0.6 を超えたとき，スラグ流と判定する．

3

Recording period 0.2s

2

1
Air
injection

00 0.5 0 0.5 0 0.5 0 0.5 0 0.5

Void fraction
Fig. 8 Characteristics between void fraction and pressure drop
(Air-Water 0.1MPa: JG0=0.3m/s, JL0=0.08m/s).
- 202 -

Fig.9 Void fraction fluctuations.

C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

本計算で用いた各係数の値は，主として実験及び既存の流動様式遷移条件で与えられる気泡流から環状流までの流動様式をできるだけバランスよく再現できる値として選んだものであるが[11]，その適否の判定は多分に定性的である．したがって，ここで用いた数値及びその組み合わせが真に最適かどうかについてはコメントできない．ここで用いたのはチューニングによって決定した数値ではあるが，全ての流動様式に普遍的に用いることができ，また本論文
におけるような空気‐水等温系のみならず，蒸気‐水等温系[2]に対しても妥当な結果を与える点は留意する必要がある．

8

6

JL0=0.03m/s
JL0=0.08m/s
JL0=0.30m/s
JL0=0.80m/s
JL0=3.00m/s+20%

uB m/s

-20%
4

2

LG m

Average period TSlug s(2)スラグ流域における流動特性
本モデルでは，管径オーダーの気泡を最小単位と
0
0
2
4
6
8
考えていることから，大気泡が流れ場の主要を占めるスラグ流を基本としたモデルといえる．そこで，
uNicklin m/s
本モデルにおけるスラグ流域の流動特性を検討する．
Fig. 10 Comparison of slug velocity to
Nicklin’s correlation.
まず，スラグ流域特有の静圧脈動現象について検討する．Fig.8 に JG0=0.3 m/s, JL0=0.08 m/s におけるテ
1.5
Simulation
ストセクション差圧およびボイド率の時間変動波形
JL0=0.03 m/s
JL0=0.08 m/sを示す．Fig.8-(a)は差圧測定区間 LP=0.99 m(LP1-LP2
JL0=0.30 m/s
JL0=0.80 m/s
JL0=3.00 m/s
間)のテストセクション差圧波形であり，Fig. 8-(b),Correlation
Akagawa(c), (d)はそれぞれ差圧測定区間入口 LP1(y=1.11 m)，
1
差圧測定区間のほぼ中央の位置 LPDF(y=1.65 m)，そして同出口 LP2(y=2.10 m)の位置におけるボイド率
JL0=0.03m/s
変動波形である．また，図中の一点差線は気泡の上昇速度に相当する．この差圧変動波形をボイド率変
JL0=0.08m/s
0.5
動波形と共に見ていくと，まず，①差圧測定部に大JL0=0.3m/s
気泡が流入してくると，液体スラグが測定部から流JL0=0.8m/s
出するため，相対的にテストセクション差圧は減少
JL0=3.0m/s
する．②測定部に大気泡が流入し終わるとテストセ
0
クション差圧の減少は止まる．③測定部に液体スラ0123グが流入すると，大気泡がこれに押されて測定部かJG0 m/s
ら流出するためテストセクション差圧は増加する．
Fig. 11 Characteristics of average period of slug bubble.
④大気泡が流出し終わると，測定部上下流部では液1相の流入出があるだけの期間がしばらく続くため，Simulation
テストセクション差圧はほとんど変化しない．このJL0=0.03m/s
JL0=0.08m/s
JL0=0.03m/s
ように流速が比較的低く，摩擦による圧力損失の影
JL0=0.30m/s
JL0=0.80m/s
0.8
J=3.00m/s
L0
響が小さい場合の差圧脈動は，対応するボイド率変Correlation
Street動から説明することができる．このような差圧脈動0.6に関する定性的な特徴は赤川ら[12]による実験におJL0=0.08m/s
いても確認されている．
JL0=0.3m/s
JL0=0.8m/s
次に，大気泡の上昇速度について検討する．Fig. 90.4に初期状態より十分に時間が経った，ある瞬間から
0.2 s ごとのテストセクション軸方向のボイド率分0.2
布を示す．計算条件は JG0=0.8 m/s, JL0=0.08 m/s であJL0=3.0m/s
る．このようなボイド率分布を用いることで，テス
トセクションに存在する大気泡（図中破線囲み）の
0
0
0.5
1
1.5
2
2.5
軌跡（図中一点鎖線）から伝播速度の評価が可能とJG0 m/s
なる．本モデルにおける結果を Fig. 10 にプロットし，次式で与えられる Nicklin[13]の式と比較する．
Fig. 12 Characteristics of length of slug bubble.

- 203 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

102

30
Simulation
BubblySlug
Churn
Annular

Simulation
BubblySlug
Churn
Annular
Correlation
Chisholm

20
JL0=0.8m/s

JG0/

101

Ltt

m/s

JL0=0.03m/s
JL0=3.0m/s

Annular flow

10

Slug flow00100

Drift-flux model
Slug flow
Annular flow

10

20

10-2

30

10-1

100

101

JT m/s0.35 gD p

0.11 0.09 J G
J L 0.2 J G(12)

離れていくとは言うものの，特別に式(13)や(14)の関係を用いなくてもおおよそのスラグ挙動を表現できたことになる．
(3)時間平均特性
ボイド率の時間平均特性として，Zuber ら[14]によるドリフトフラックスモデルに基づく整理をしたものが Fig. 13 で，総体積流束に対する気相の速度 JG/の関係を示している．計算結果は流動様式ごとに記号を変えてプロットし，図中実線および破線はそれ
ぞれスラグ流および環状流における相関式，式(15)を示している．

JG
C0 J T v gj

(13)

Slug flow

液相体積流束が大きい場合は赤川らの相関式と良く一致しており，エラーバーの幅も小さいのに対し，液相体積流束が小さい場合は相関式を下回る値を取り，
エラーバーの幅はかなり大きな値をとっている．
最後に先ほど評価した周期と上昇速度より大気泡の長さの評価を行う．Fig. 12 に気相体積流束に対する平均周期を示し，Street ら[13]により提案されている以下の気泡長さの相関式を実線で示す．

C0

LStreet0.29 J G
J L 0.12

104

Fig. 14 Characteristics of average frictional pressure drop.

計算結果と Nicklin の式とを比較すると，低液相体積
流束域で相関式よりも少し大きく，高液相体積流束域で少し小さくなっているが，全体としては±20%程度で再現できている．
次に，Fig. 8-(c)に示したボイド率の時間変動波形
より，液相スラグ，大気泡の存在時間を求め，従来の式と比較してみる．Fig. 11 に気相体積流束に対する平均周期の値をプロットし，赤川ら[11]により提案されている以下の相関式を実線で示す．なお，平
均値はプロットで，最大値・最小値の幅をエラーバ
ーで表現している．

T Akagawa

103

tt

Fig. 13 Characteristics of average void fraction.

u Nicklin 1.2 J T

102

(14)

気相体積流束が大きくなるほど，また，液相体積流束が小さくなるほど気泡長さが長くなるという定性的な特徴については本モデルで再現できた．液相体
積流束が小さくなるにつれて計算結果は相関式から

1.2 ， v gj

(15)

0.35 gD pL

GL

Annular flowC0

1.0 ， v gj

23

LJL
G Dp

L

G
L

ここで，C0 は分布パラメータ，vgj はドリフト速度でZuber ら[14]によって，上述のように与えられている．
計算結果は低総体積流束側で多少，相関式よりも高い値となっているが，全体的にはうまく再現できている．
次に，圧力損失の時間平均特性として，
Lockhart-Martinelli 相関[15]に基づいて，摩擦損失の計算結果の整理を試みたのが Fig. 14 である．計算結- 204 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

果は流動様式ごとにプロットされ，Chisholm[16]による相関式は実線で示されている．ここで， tt と Ltt
はそれぞれ以下のように定義している．

Ltt

2

2
tt

Pf

Pf

y

y

Pf

Pf

y

(16)
L

y

L

(17)
G

また，Chisholm による相関式は液相および気相が乱流であるとした以下の式を使用した．
2

21

1

Ltt

1

(18)

2
tt

tt

チャーン流および環状流領域で計算結果の分散が若干大きいものの，全体的には Lockhart-Martinelli 相関と良く一致していることが分かる．
以上，時間平均値あるいは統計的特性について既存の相関式等との比較を通じて，本モデルの適用性を実証することができた．次節では時間変動挙動そのものについて実験と本モデルの比較を行う．
3.2 流動脈動等温二相流における流動特性
著者らは既に流動脈動下の流動様式遷移に着目した実験を行い，脈動に伴う流れの停滞あるいは逆流時に大気泡が過渡的に形成され，定常流を基本と

3

0

1/5

2/5

3/5

4/5

(1)ボイド率変動挙動
液相体積流束 JL を以下に示すように入口流量に脈動流を印加した流動脈動下のボイド率変動挙動についてシミュレーションを行った．

JL

J L0

Bubbly flow
Phase
6/5

7/5

8/5

9/5

J L sin

2

t

(19)

ここで，JL0 は定常値であり， JL は振幅である．
計算結果は Fig. 6 のような定常流で見られたボイド率分布に加えて，入口流量が逆流もしくはかなり流量が小さい範囲（1/5 ～6/5 ）と，入口流量が増加して，最大となり，減少に転じる範囲（7/5 ～2 ）のそれぞれの範囲において異なる流動様式が確認できるようなボイド率分布変動が得られた．そこで，その例として Fig. 15，16 に初期状態より十分に時間が経った，ある瞬間から 0.2 s ごとのテストセクショ
Annular flow

2

Axial coordinate y m

Axial coordinate y m

Slug flow

した流動様式とは，特に脈動周期の短い場合には大きく異なること，また加熱系であれば，CHF を引き起こすなどの問題の可能性を指摘している[7]．ここでは既報[7]の実験に対応した条件で本モデルを用いたシミュレーションを行い，ボイド率変動挙動や流動様式の変化について実験結果との比較を行う．
なお，実験装置には内径 15 mm の透明円管を用いた
空気‐水垂直上昇流であり，空気流量は一定，水流量は定常流にシリンダ－ピストンによる強制脈動を印加した脈動流とした．実験装置・方法や結果の評価については既報[7]を参照されたい．

2

3

0

1/5

2/5

3/5

4/5

Churn flow
Phase
6/5

7/5

8/5

9/5

2

2

1

1

Air
injectionAir
injection

00

0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 1

0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 0 0.5 1

JL m/s

0.2
0
-0.2
0

2

1

2

Time s3

4

Air-Water 0.1MPa
JG0=0.06 m/s
JL0=0.06 m/s
JL=0.327 m/s=2 s

Fig. 15 Void fraction fluctuation in oscillatory
two-phase flow (Air-Water: JG0=0.06
m/s, JL0=0.03 m/s, JL=0.327 m/s, =2 s).

0

0.2
0
-0.2
0

Void fracrion2

JL m/s0

Void fracrion

12Time s3

4

Air-Water 0.1MPa
JG0=6.00 m/s
JL0=0.06 m/s
JL=0.327 m/s=2 s

Fig. 16 Void fraction fluctuation in oscillatory
two-phase flow (Air-Water: JG0=6.00
m/s, JL0=0.03 m/s, JL=0.327 m/s, =2 s).

- 205 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)(a) Calculation results(JG0=0.12 m/s).

(b) Experimental results(JG0=0.10 m/s).

Fig.17 Void fraction fluctuation in oscillatory two-phase flow(Air-Water: JL0=0.03m/s, JL=0.327m/s, =2s).
ン軸方向のボイド率分布（実線）を入口液相体積流束の変動波形（破線）とともに示す．
Fig. 15 では，入口流量が増加する範囲で 7/5 に示
すように，テストセクション全体に小気泡がまばらに存在する気泡流が確認できる．そして，入口流量が減少し，逆流する範囲では 3/5 に示すように，テストセクション上流部では小気泡がまばらに存在する気泡流であるものの，入口流量の逆流により，テストセクション下流部において気泡合体により大気泡が確認できるスラグ流に遷移している．
また Fig. 16 では，入口流量が増加する範囲で 7/5に示すように，テストセクション全体に限界ボイド率に達する長矩形状のボイド率分布がほぼ規則的に形成され，その間を低ボイド率側になだらかに遷移するチャーン流が確認できる．そして，入口流量が減少し 1/5 に示すように，空気混合部で限界ボイド率を一様に保持するボイド率分布が形成されると，入口流量の逆流に伴い，3/5 に示すようにそれがテストセクション全体を覆うような環状流へと遷移している．このように入口流量変動の 1 周期内を通して，Fig. 15 では気泡流からスラグ流，Fig. 16 ではチ
ャーン流から環状流へと遷移しており，ここではこれらの流動状態をそれぞれ，遷移スラグ流，遷移環状流と定義する．

先ほどの遷移スラグ流域において，気液二相流動特性実験[7]で指摘しているような様々な問題を引き起こす可能性のある過渡的に形成された大気泡挙動を確認することができたため，以下でそのボイド率変動挙動を詳しく検討する．
Fig. 17-(a)は本モデルでの計算結果を示しており，Fig. 3 で示したのと同様な方法で二次元離散気泡状に表示している．ただし，限界ボイド率に達したものが連続した場合，それは一つの大気泡とみなし，気泡間は連結して表記している．計算条件はJG0=0.12 m/s，JL0=0.03 m/s， JL=0.327 m/s， =2 s である．また，Fig. 17-(b)は実験[7]より得られた流動可視化画像を 0.2 s ごとに並べたもので，実験条件はJG0=0.1 m/s，JL0=0.03 m/s， JL=0.327 m/s， =2 s であり，それぞれの寸法をおよそ合わせて表記しており，各図にスケールを示す．また，局所ボイド率変動（図中に示す位置 z は空気混合部からの位置である）および入口流量変動も同時に示す．
これらを共に見ていくと，①入口流量（Fig. 17(b)ではオリフィス差圧トランスデューサの出力電圧が表示されている）の減少に伴い空気混合部上部で気泡が停滞し始め，②入口流量が最小となるあたりで気泡は合体・成長して大気泡になり，③これが入口流量の増加と伴に上流から下流に向かって伝播して- 206 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

JG0 m/sCalculation results
BubblyBubbly-Slug
Churn
Churn-Annular

Slug
Annular

0.1

Flow pattern transitions
Mishima et al.
Taitel et al.

JL0 m/s

JL0 m/s

0.1

0.01

0.01

Bubbly
0.001

0.01

Slug
0.1

Churn
1

Annular10

JG0 m/s

Bubbly

100

0.001

0.01

0.1

Churn
1

Annular
10100

JG0 m/s(a) Steady flow

JL0 m/s

0.1

JL0 m/s

0.1

Slug

0.01

0.01
Bubbly
0.001

0.01

Slug
0.1

Churn
1

JG0 m/sAnnular10

Bubbly

100

0.001

0.01

(b) Oscillatory flow: JL=0.327 m/s, =2s

Churn
1

Annular
10100

JG0 m/s

JL0 m/s

0.1

JL0 m/s

0.1

Slug
0.1

0.01

0.01
Bubbly
0.001

0.01

Slug
0.1

Churn
1

JG0 m/sAnnular10

Bubbly

100

0.001

0.01

(c) Oscillatory flow: JL=0.164 m/s, =4s

0.1

Slug
0.1

Churn
1

Annular
10100

JG0 m/s

JL0 m/s

JL0 m/s

0.1

0.01

0.01
Bubbly
0.001

0.01

Slug
0.1

Churn
1

Annular10

Bubbly

100

0.001

0.01

JG0 m/s(d) Oscillatory flow: JL=0.109 m/s, =6s

Slug
0.1

Churn
1

Annular
10100

JG0 m/sFig. 18 Flow pattern map(Left: Simulation, Right: Experiment).

- 207 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

いる．計算結果では気泡の成長から下流への伝播までが実験結果と比べて若干早いが，これは実験においては，計算で与えているような正弦波を作れないためであると考えられる．局所ボイド率変動については，入口流量がおよそ最大を迎えるあたりで局所
ボイド率変動の周期的なピークが発生しており，これは実験結果および計算結果においても確認することができる．これは成長した大気泡が入口流量の増加と伴に下流に周期的に伝播するためである．以上のように，本計算結果は流動脈動下においても実験結果をよく再現できている．
このような流量の大規模な変動に伴うボイド波伝播の特性は，例えば既報[17]におけるように本報でもグローバルルールとして用いたボイド波伝播式とドリフトフラックスモデルによるボイド率相関式を用いて再現できるし，二流体モデルによっても可能である．ただしその際には過渡的な流れの状態において流動様式が変化し，それに伴ってボイド率相関式を切り替えるなどの所作が必要となるばかりでなく，その切り替えに伴ってボイド波伝播に不自然な挙動が出現しないように遷移域を設けるなど追加的な操作が必要となる．それに対して今回のモデルでは細部において若干乱暴な取扱いもあるが，そのような操作なしに，いわば自然発生的に過渡的な大気泡の発生を再現できている点に本質的な差異がある．
(2)流動様式線図
先で得たような流動脈動下のボイド率分布を基に，流動様式を脈動周期ごとに JG0－JL0 平面にプロッ ト し， 同時 に定 常流 につ い ても 同様 に示 すと
Fig.18-(a)～(d)のようになる．ここで，Fig. 18 内左図
は本モデルでの計算結果を示しており，右図は実験結果[7]を示している．まず，定常流については，流動様式の遷移を各理論式と実験結果を良く再現できている．次に，流動脈動流における流動様式遷移については，定常流での流動様式遷移と比べて遷移ス
ラグ流，遷移環状流がそれぞれ低気相体積流束側に発達しているのに対し，スラグ流－チャーン流の遷移境界は変化していない．脈動周期が短くなり，脈動振幅が大きくなるほど，遷移スラグ流が低気相体
積流束側に発達している傾向があり，計算結果は実験結果をよく再現できている．
４

結論

既報[2,4,5]の離散気泡モデルに改良を施し，具体的な実験結果と比較することによって，以下のことが得られた．
定常等温二相流下におけるボイド率の時空間変
動挙動を実現することができた．
これより，気泡流，
スラグ流，チャーン流，環状流の各流動様式を確認することができ，PDF を用いた流動様式判別を利用することで，Mishima ら，Taitel らによる流動様式線

図と良く一致した．特にスラグ流においては静圧脈
動現象を再現することができ，大気泡の上昇速度を相関式との良い一致が得られたことを確認し，大気泡の存在周期，大気泡の長さにおいては定性的に有意な結果が得られた．さらに時間平均特性としてボ
イド率をドリフトフラックスモデル，摩擦圧力損失
を Lockhart-Martinelli 相関による整理を行うことで定量的に既存の相関式とよく一致する結果が得られた．
流動脈動等温二相流下におけるボイド率の時空間変動挙動を実現することができ，気液二相流動特性実験の結果を良く再現することができた．また，脈動流下における流動様式についても実験結果と良く一致した．
以上の結果から，本離散気泡モデルは流動脈動を含む等温二相流の時空間変動挙動を実現することができ，従来から使われているような時間平均化されたモデルとは異なり，二相流ダイナミクスが原因となるような過渡現象について対応可能なモデルを構築できたと考えている．
ここで用いたモデルは，スラグ流をベースにおいて展開したモデルで，当然のことながら１次元流が対象である．そのようなことから，このモデルはスリップに関する関係式に表面張力効果を組み入れる
ことにより，比較的簡単に細管内二相流のシミュレーションが可能になる．水平管においては細管への適用もさることながら管内で相分離が起こるような場合に興味があり，そのような場合も含めて目下検討中である．また既報[17]にも示しているが，本報における混合部の取扱い方と同様に各セルに気相の湧き出し項をおくことによって沸騰流のモデルが可能になる．ただし本報の緒言で述べたような CHF
問題に適用するには，液膜消失の組込み方や，管内面での伝熱特性と壁面の熱容量の動的干渉を如何に組込むかなどが非常に重要になる[1, 18]．現在，これらの問題を解決すべく鋭意検討中である．
本研究の一部は 2005-2006 年度関西大学学術研究
助成（共同研究）によったことを記し，謝意を表す
る．
参考文献[1] Ozawa, M., “Critical heat flux induced by flow
instability in boiling channels”, ECI Int. Conf. on
Boiling Heat Transf., Spoleto, Paper No.KL-2(2006).
[2] Ozawa, M., Ami, T., Umekawa, H., Shoji, M.,“Pattern dynamics simulation of void propagation”,
4th Japanese-European Two-Phase Flow Group
Meeting (2006), to be published in Multiphase
Science and Technology, 19-4 (2007).
[3] Collier, J. G., Thome, J. R., Convective boiling and- 208 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)

condensation (1994), 375-429, Oxford University
Press.
[4] 伊藤浩二，井上満，小澤守，庄司正弘，“気液二相流の流動様式の簡易遷移モデル”，日本機械学会論文集，70-689, B (2004), 151-158.
[5] Ito, K., Ozawa, M., Shoji, M., “Pattern-dynamics
approach to two-phase flow regime transition”,
Japan-US Seminar on Two-Phase Flow Dynamics,
1 (2004), 103-112.
[6] Schlichting, H., Boundary layer theory, 6th ed.
(1968), 685-695, McGraw-Hill.
[7] Inoue, Y., Okano, M., Umekawa, H., Ozawa, M.,“Flow pattern transition and void fraction behavior
in adiabatic oscillatory two-phase flow”, Proc. the
4th JSME-KSME Thermal Engineering Conference,
1 (2000), 747-752.
[8] Jones Jr., O. C., Zuber, N., “The interrelation
between void fraction fluctuations and flow
patterns in two-phase flow”, Int. J. Multiphase
Flow, 2 (1975), 273-306.
[9] Mishima, K., Ishii, M., “Flow regime transition
criteria for upward two-phase flow in verticaltubes”, Int. J. Heat Mass Transfer, 27-5 (1984),
723-737.
[10] Taitel, Y., Bornea, D., Dukler, A. E., “Modelingflow pattern transitions for steady upwardgas-liquid flow in vertical tubes”, AIChE J., 26-3
(1980), 345-354.

[11] 伊藤浩二，気液二相流の流動様式と遷移モデル，
東京大学学位論文 (2002)．
[12] 赤川浩爾，浜口八朗，坂口忠司，碇哲二，“気液二相スラグ流の差圧脈動に関する研究（第１報，差圧脈動の実験結果）”，日本機械学会論文
集，36-289, B (1970), 1520-1527．
[13] 日本機械学会，改訂 気液二相流技術ハンドブ
ック (2006), 191-248, コロナ社．
[14] Zuber, N., Staub, F. W., Bijwaard, G., Kroeger, P. G.,“Steady state and transient void fraction in twophase flow systems”, GEAP-5417 (1967).
[15] Lockhart, R. W., Martinelli R. C., “Proposedcorrelation of data for isothermal two-phase,
two-component flow in pipes”, Chem. Eng. Prog.,45-1 (1949), 39-48.
[16] Chisholm, D., “A theoretical basis for the
Lockhart-Martinelli correlation for two-phase flow”,
Int. J. Heat Mass Transfer, 10 (1967), 1767-1778.
[17] 平山美緒，梅川尚嗣，岡野充浩，小澤守，“気液二相脈動流の流動様式とボイド率変動特性”，
混相流，17-3 (2003), 285-295.
[18] Ozawa, M., Umekawa, H., Mishima, K., Hibiki, T.
and Saito, Y., “CHF in Oscillatory Flow Boiling
Channels”, Chemical Engineering Research and
Design, 79-A4 (2001), 389-401.

- 209 -C 2007 The Heat Transfer Society of Japan○