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

高速二相流中に発生する 2 次元膨張波に関する研究
原田 敦史‡

中川 勝文†

Two-dimensional Rarefaction Waves in the High-speed Two-phase Flow
Masafumi NAKAGAWA† and Atsushi HARADA‡

Abstract
Two-phase flow nozzles are used in the total flow system for geothermal power plants and in the ejector of
the refrigerant cycle, etc. One of the most important functions of a two-phase flow nozzle is to convert thethermal energy to the kinetic energy of the two-phase flow. The kinetic energy of the two-phase flow exhausted
from a nozzle is available for all applications of this type. There exist the shock waves or rarefaction waves atthe outlet of a supersonic nozzle in the case of non-best fitting expansion conditions when the operation
conditions of the nozzle are widely chosen. The purpose of the present study is to elucidate theoretically thecharacter of the rarefaction waves at the outlet of the supersonic two-phase flow nozzle.
Two-dimensional basic equations for the compressible two-phase flow are introduced considering theinter-phase momentum transfer. Sound velocities are obtained from these equations by using monochromaticwave approximation. Those depend on the relaxation time that determines the momentum transfer. The
two-phase flow with large relaxation times has a frozen sound velocity, and with small one has an equilibriumsound velocity. Rarefaction waves which occurred behind the two-phase flow nozzle are calculated by the CIP
method. Although the frozen Mach number, below one, controls these basic equations, the rarefaction waves
appeared for small relaxation time. The Mach line behind which the expansion starts depends on the inlet
velocity and the relaxation time. Those relationships are shown in this paper. The pressure expansion curves are
only a function of the revolution angle around the corner of the nozzle outlet for the relaxation time less than 0.1.
For the larger relaxation time, the pressure decays because of internal friction caused by inter phase momentum
transfer, and the expansion curves are a function of not only the angle but also the flow direction. The calculated
expansion curves are compared with the experimental ones. The calculations considered only the compressibilityof the gas phase, but those resembled with the curve of the experiment in which the phase changes occurred.
Key Words: Two-phase flow, Nozzle, Supersonic, Rarefaction wave, Mist, Mach number
記 号
a
:
d
:
f
:
k
:
:
Κ
L
:
:
Μ
p
:
r
:
t
:
v
:

音速
液滴径
速度ポテンシャルのθ 関数部
伝播数
定数 α0p0/(1-α0)
入口長さ
マッハ数
圧力
半径座標
時間
速度

[m/s]
[m]
[m/s]
[rad/m]
[Pa]
[m]
[Pa]
[m]
[s]
[m/s]

x

:

x0

:
:
:
:
:
:
:
:
:

α
γ
µ
θ
ρ
τ
φ
ω

[m]
空間座標，ノズル出口からの距離
入口クオリティ
ボイド率
比熱比
粘性係数
[Pa･s]
[rad]
角座標
3
[kg/m
]
密度
[s]
緩和時間
[m2/s]
速度ポテンシャル[rad/s]
角周波数

受付日: 2007 年 3 月 19 日，第 44 回日本伝熱シンポジウムより受付, 担当エディター: 河村 洋
豊橋技術科学大学 機械システム工学系 (〒441-8580 愛知県豊橋市天伯町雲雀ヶ丘 1-1)‡
豊橋技術科学大学 機械システム工学専攻 (〒441-8580 愛知県豊橋市天伯町雲雀ヶ丘 1-1)†

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

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

添字
e
f
g
l
0
1

緒

そこで，本研究では相変化を起こさず等温変化を
:
:
:
:
:

平衡状態
凍結状態
気相
液相
入口

すると仮定した高速ミスト流中の相間の運動量緩和現象に注目して理論解析を行い，運動量の緩和を伴う流れ場に発生する膨張波の流体力学的性質を明ら
かにする．
２

論

高速二相流ノズルは，二相流の運動エネルギを利用する装置において重要な役割を果たす．
たとえば，
冷凍サイクルに適応される二相流エジェクタ[1]では，冷媒の膨張時に放出されるエンタルピ落差を二相流の運動エネルギとして効率よく変換させなければならない．二相流の音速は，相変化を起すとき数m になり[2]，相変化を起さない場合でも数十 m のオーダーであり，非常に低い．このような二相流は少が音速を大きく変えるため，圧縮性流体の特性に大きな影響を与えることが明らかにされてきた．本研
究では，付加的な複雑さを避けるため，連続相が気相で分散相が液相(液滴)であるミスト流中に存在す
る運動量緩和現象のみを考慮し，2 次元流れの膨張現象について調べた．
気相と液相のそれぞれの質量保存を表す連続の式
は以下のようになる.

く知ることが重要となる．著者ら[3]はこれまでに冷t

凍サイクルのエジェクタを開発してきたが，ノズル

流膨張波は超音速の性質と亜音速の性質の両方を持gtて，ノズルから二相流が噴出するときの挙動を詳しの特性をまず実験[4]で明らかにしてきた．この二相式

これまでの 1 次元二相流の研究[5]で運動量緩和現象

二相流の運動エネルギを利用する機器の設計におい背後に発生する膨張波の特性を知る必要があり，こ礎

相変化はエネルギ変換機器で非常に重要となるが，し加速されるとすぐに超音速流のような振る舞いをする．二相流エジェクタや二相流タービンのような

基

1

v

xi

g gi

l

xi

(1)

0

1

v

(2)

0

l li

ここで，二相流中の気相の体積分率を表すボイド率を α，密度を ρ，流速(vx,vy)とおく．ただし，添え字
g および l はそれぞれ気相と液相を，i は 2 次元の各方向を表す．

ち備えていることを実験で示した．
また，二相流ノズル内の流れ特性は，これまで 1次元流れを用いて理論的に研究[5,6]されてきた．高速二相流では，
特に相間の運動量輸送が重要であり，各相の加速はその緩和時間に大きく依存することが示された．また，衝撃波管を用いた 1 次元気泡流の膨張波の研究[7]なども行われてきた．一方，2 次元圧縮性二相流の研究では，ミストの量が非常に少な
く，流れが圧縮性気体に支配されている領域に関するもので，タービンブレード[8]，膨張波[9,10]に関しての研究があるが，液体の量が少ないとして液体

熱水タービンや冷凍機のエジェクタの二相流のク
オリティは 0.1 のオーダーでミスト流の領域にある．
液滴径が十分小さい場合，液滴にはストークス抵抗による相間の摩擦力が働く．この力を単位体積当たりに働く液相の外力として書き表すと式(4)の右辺の第 1 項の形になる．ただし，緩和時間 τ は
τ=ρld2/(18μg)で液滴径と物性値で表される．気相には逆の力が働き，式(3)の右辺の第 1 項となる．したがって，気相と液相のそれぞれの運動量保存を表す運動量式は次式のように与えられる．
v giの運動量を無視するか小さい等の近似を行っている．

g

また，気泡流領域での 2 次元二相流では翼回りの流

v

g gj

t

れ[11]があるが，基礎式で気体の運動量等が無視さ

xj

v gi
vgi

1

れている．本研究の対象はクオリティが 0.1 から 0.5程度の低クオリティのミスト流で，これらの近似が1成立しない液体と気体の両方の性質が表れる二相流であり，この領域の 2 次元圧縮性二相流の研究はほ

l

v li
t

1
1

とんどない．

vli

p
xi

l

v

l

l lj

xj

vgi

vli

(3)

v li
1

p
xi

(4)

ここで，p は圧力であり，これらの基礎式は代表的- 224 -C 2007 The Heat Transfer Society of Japan○

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

な流体力学の数値解析の本[12]と同様のものである．
ae

二相流のエネルギ式は各相の温度変化を求めるた
めに必要となるが，本研究で対象とするクオリティ

p g

x

1 x

x

g

l

が 0.1 程度のミスト流では液滴が大きな熱容量をもp1

つことから気相はほぼ等温変化をするとみなせるたg

め，温度一定の状態方程式を用いた．液体の圧縮性(6)lは気相のそれに対して小さく無視できるので，液相

ここで，音速を a，添え字 e は平衡状態であることの密度 ρl は一定とした．また，緩和時間 τ は一般的を表す．

には流れ場の関数となるが，本解析では大きく変化しないとして定数として扱った．

また，式(5)において，τ→∞とすると凍結流の音速が得られる．

これらの基礎方程式は運動量緩和時間 τ によって
af

特徴づけられる．すなわち，τ が小さくなると相間

p g

x

1 x

x

2
g

2
l

の運動量輸送が活発になるため速度差が小さくなり，
τ→0 となる極限では，相間の運動量輸送が平衡に達し速度差が無くなる平衡流となる．逆に，τ が大きくなると相間の運動量輸送が疎になり各相間の影響

p

p1

g

l

(7)

力が少なくなり，τ→∞となる極限では，相間の運動

ここで，添え字 f は凍結状態であることを表す．

量輸送が無くなり各相が独立に流れる凍結流となる．

ρl>>ρg であるので，凍結音速は気体の音速に等しく凍結という用語は運動量輸送にふさわしくないが，無次元音速 a’=a/(p/ρg)1/2 はほぼ 1 になる．

緩和現象を扱うときに良く用いられるので，ここで式(5)から得られる微小擾乱の伝播速度，すなわち音速を無次元数 ωτ を横軸にとって Fig. 1 に示す．式も用いた．
基礎式(1)～(4)の微分項は，τ→∞として得られる凍(5)は気体の等温音速の 2 乗 p/ρg で無次元化すれば，結流の基礎式とまったく同じである．すなわち，こ

ボイド率 α と密度比 ρl/ρg がパラメータとなるが，以れらの方程式を支配する無次元数は，凍結流を支配前われわれが行った熱水の実験ではクオリティがする無次元数である凍結マッハ数となる．

0.1 で，密度比が 103 のオーダーであったのでこれから x/(1-x)=αρg/{(1-α)ρl}の関係を用いてこのパラメー

３

音

速

タを定めた．また，この図で縦軸は気体の等温音速

非定常基礎式(1)～(4)を線形化し，角周波数 ω の1

単色波がこれら式を満足する条件から微小擾乱の伝播速度が求まる．これを実行すると，この単色波のaf

複素伝播数ベクトル(kx,ky)が満足しなければならな

0.8

a'

い次の関係式が得られる．
i

2

k

2
x

p

1

p
g

l

k

i

g

g

1

0.6

(5)

l
2
y

x=0.1
ae

l

土方ら[13]は 1 次元の二相流の音波の伝播速度を

l

/

g

1000

0.4

求めているが，本論文では，2 次元の問題を扱うため，kx と ky が表れる．また，そこでは温度や相変化

0.1

1

10

100

の緩和現象も扱われているが，ここでは運動量の緩和現象のみが考慮された．
式(5)において，τ→0 とすると平衡流の音速が得られる．

- 225 -

Fig. 1

Sound velocity of two-phase flow with
momentum relaxation.

C 2007 The Heat Transfer Society of Japan○

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

で無次元化されている．ωτ が小さいときは平衡音速

入口のマッハ数 M0 は 1.27 となる．
超音速流の場合，に，大きいときは凍結音速に漸近する．この結果か下流の影響が上流に伝わらないのでノズル内の流線

ら，音速は ωτ が 1 より小さいときはほぼ平衡音速は変化しない．超音速流れ場のマッハ角 sin-1(1/M0)

で，20 より大きいときはほぼ凍結音速で表せる．

は 52°になり，このマッハ角の膨張波の下流で減圧する．Fig. 3 には，このときの中心線に沿う圧力分

４布を実線でノズル出口からの距離 x に対して示した．

気体力学の膨張波との比較

マッハ角で示される斜め膨張波が到達した点から減
本解析では，ノズル出口に発生する膨張波を理論圧し始める．一方，亜音速流の場合，ノズルから流的に調べるため，Fig. 2 に示すような形状で解析を出した後に速度が低下するので圧力は上昇する．Fig.

行った．上部はノズルの流れ方向の中心線に対応し3 には同じ体系で計算された非圧縮性流体のポテンているとして対称対象境界条件を与えた．2 次元ノ

シャル流から得られる中心線に沿う圧力分布を破線

ズル出口高さの半分を代表長さ L にとり，ノズル出で示す．ノズルから出た後に圧力が低下し流体が加口の手前は長さ(1/2)L だけノズル内部とし，ノズル速することは，超音速流の特徴である．

を出た後の背圧室は 5L×(9/2)L の大きさで計算が行次に，二相流の場合の計算結果を Fig. 4～Fig. 6 にわれた．下部および右側は自由流出の境界条件を与示す．入口クオリティ x0 を 0.1，密度比 ρl/ρg0 を 103，

えた．

気体と液体の入口が同じであるとして入口無次元速

距離を L で，圧力を入口圧力 p0 で，速度を入口の

度 u’を 0.5 とした計算結果である．無次元緩和時間

で無次元化して計算が行われた．

τ’=τ/{L/(p0/ρg0)1/2}は，それぞれ 0.01，0.1，1.0 とした．

計算はメッシュ間隔を 0.1L とし，CIP 法[14]を用い

図中，青線で気体の流線を，赤線で液体の流線を示た非定常 2 次元計算が行われた．十分時間を経て波す．これらの流線は前述のように各相の質量流量を形が変わらなくなったところで，定常解として結果示すものである．前節で述べたように本解析で用いを得た．

た基礎式を支配する音速は凍結音速で，ほぼ気体の等温音速(p0/ρg0)

1/2

Fig. 2 には，この体系で計算された比熱比 が 1.4等温音速に等しい．入口速度がこの音速の半分であである理想気体の膨張波の流線がプロットされているにもかかわらず，超音速流れが達成されている．

る．ただし，本論文では曲線間の質量流量が一定でこのことは，緩和時間が小さい場合，流れは平衡流ある曲線を流線と定義した．したがって，定常流でに近いので平衡音速のマッハ数が 1 より大きければ，あるが圧縮性流体であるので接線は流速の方向を正支配方程式から見れば亜音速であるにもかかわらず

確には示していない．等温音速で無次元化された入超音速の振る舞いをすると解釈できる．本解析でも1/2

口速度 u0’=u0/(p0/ρg0) は 1.5 である．この条件では2outlet

1.6

nozzle

1.2

5L

p/p0

L

L
2

u'0=1.5

0.8

1.4
u0 ' 1.5

0.4

9
L
2
Fig. 2

Potential flow
Supersonic flow

0

2

4

x/L

Configurations of analyses. And Stream-line of
rarefaction waves of an ideal gas.

- 226 -

Fig. 3Pressure profiles of a supersonic andincompressible potential flow along the center
axis of the nozzle.
C 2007 The Heat Transfer Society of Japan○

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

平衡音速より小さい入口速度で膨張する解は得られ

なわち，緩和時間の増加によってマッハ数が減少すなかった．したがって，二相流膨張波の特徴は，凍るような形態を示す．緩和時間への依存性について

結音速のマッハ数が亜音速でも超音速の形態を示すは次節で詳しく述べる．また，緩和時間はミスト径ことである．

の 2 乗に比例しているので，緩和時間が大きいこと

Fig. 7 は，Fig. 4～6 のノズル中心線上の圧力分布は液滴径が大きいことを意味している．Fig. 4 からをノズル出口からの距離 x をとって示したものであFig. 6 になるに従い，ミストの慣性力が増加するたる．ノズル出口の角から出る膨張波が始めてノズル

め，液体の流線の曲率半径は大きくなる．二相流膨の中心線まで到達する位置は，入口流速が同じであ

張波のもう 1 つの特徴は緩和時間によって流れの状るのに緩和時間が増加すると上流側に移動する．す態が大きく変わることである．
５

膨張開始点の緩和時間への依存性
膨張波が発生し始める条件は前述のように緩和

時間に依存し，このような高速二相流流れの形態を大きく左右する．ノズルの中心線に沿う減圧特性として現れる膨張波の開始点はマッハ数が大きいとき

u0 ' 0.5' 0.01

気体力学では変曲点となって現れるが，二相流膨張
波では滑らかに変化する．たとえば Fig. 8 の実線は入口速度 u0’=0.9，緩和時間 τ’=0.1 の場合のノズル中Fig. 4

Stream-lines of rarefaction waves of two-phaseflow with relaxation time τ’=0.01 relaxation.

心線上の無次元圧力 p/p0 をノズル出口からの距離
x/L の関数として表したものである．本研究では Fig.
8 のように，50％膨張した点で膨張曲線を近似する 2次曲線（図中に点線で表す）が p/p0=1 の線と交わる点を膨張開始点と定めた．厳密な意味での膨張開始点ではなく，膨張曲線の大きな変化が仮想的に起こり始める点とした．
ノズル出口の角から出る二相流の膨張波も，気体

u0 ' 0.5' 0.1

力学の膨張波から類推できるように，膨張波に流入
1

Stream-lines of rarefaction waves of two-phaseflow with relaxation time τ’=0.1 relaxation.

'
0.01
0.10
1.00

0.8

p/p0

Fig. 50.6

u0'=0.5

0.4

0.2

u0 ' 0.5
' 1.0

0

1

2

3

4

x/L
Fig. 6

Stream-lines of rarefaction waves of two-phaseflow with relaxation time τ’=1.0 relaxation.

- 227 -

Fig. 7

Pressure profiles along center axis
for the various relaxation time.

C 2007 The Heat Transfer Society of Japan○

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

する速度，すなわち流入速度の膨張波の波面に垂直に大きく依存し緩和時間が小さいとき，平衡音速にな成分が音速となる性質を持っていると考えられる．

近く，大きくなるにしたがい，徐々に凍結音速へとそうすれば，上で求めた膨張開始点とノズルの角に変化していく．入口流速の依存性も，入口速度が大引いた直線（膨張波）が流れ方向となす角 θ はマッ

きいほうが音速が速くなることを表している．

ハ角（sinθ=1/M0）となるので，入口流速を用いると緩和時間の小さい場合は，ほぼ 3 節で求めた音速開始点から流れの音速を計算することができる．Fig.

と一致するが，無次元緩和時間が 2 以上で，この整9 は膨張開始点から計算される音速 a を Fig. 1 と同理方法との一致が悪くなる．これは，緩和時間が大様に ωτ を横軸にとって示したものである．
ここで，きい場合，液体の慣性力のため液体は曲がらずほと波の角周波数 ω として，代表長さであるノズル幅 L

んど膨張しないため，液相の流れから出たときに膨と代表速度となる入口速度 u0 から 2πu0/L が用いられ張することと，緩和時間の大きな流れ場では気液のた．図中のそれぞれの記号は無次元緩和時間 τ’をパ

速度差が大きく，気液間の速度差による内部損失とラメータにとってある．求められる音速は緩和時間

圧力が減少するので，圧力波形から直接音速を予測

p/p0

することは難しいためであると考えられる．
1.1
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3

６

二相流膨張波の減圧特性
気体力学において，半無限流体中にある角を回る

流れはプラントル・マイヤーの膨張波として知られている．長さの次元が消えるので，膨張時の圧力はx/L

角の回り角だけの関数として表される．
本解析でも，
膨張特性を明らかにするため，膨張時の圧力とノズ
ルの出口の角からの角度の関係を調べた．

θ
sin

比較のため，気体の圧縮性のみを考慮した平衡二1M0

相流の膨張波の関係式を示す．気体の速度と液体の速度が等しく，気体の密度のみが圧力に比例すると仮定して導いた．速度が速度ポテンシャル φ=rf(θ)

Fig. 8Calculation method of sound velocity
using Mach angle.

で表せるとき，r-θ 座標の二相流の連続の式と運動量
の式から
df
d

1.2
Theoretical solution
'
0.01
0.02
0.05
0.08
0.10
0.20
0.50
0.80
1.00
2.00

1

a'

0.8

0.6

のであり，膨張波の波面に垂直な 方向の速度が音速になるということを意味している．その導出方法
は，たとえば[15]に示されている．ポテンシャル流
れのとき，平衡二相流の運動量方程式は積分でき，

1

10

100

u
2 0
L
Fig. 9(8)αρg+(1-α)ρl である．式(8)は，文献[4]の式(5)と同じも

xj
0.1

dp
d

の関係が得られる．ここで，ρ は気体の平均密度

0.4

0.2
0.01

2

Sound velocities obtained from Machangle for the various relaxation time.

vi2
2

p K log p
K g 0 p0l0(9)の形になり括弧の中が保存される．ここで，二相流の連続の式(1)，(2)と状態方程式から
p
vj
0
xj 1

(10)

が得られ，括弧の中が流線に沿って保存されること- 228 -C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No.4 (2007)が用いられた．αp/(1-α) は定数となるのでこれを K

等温変化のものである．

とした．(8)式の右辺の音速の 2 乗は
dp
d

p K
K K g 0 p0ノズル出口の角からの距離を r，ノズル壁の水平

2

(11)
l

線からの角度を θ として，無次元緩和時間 τ’が 0.1，
無次元入口流速 u0’が 0.96 の解析結果の圧力を θ の関数として Fig. 10 にプロットする．各記号は図中にと書かれ(p+K)の関数となる．(9)式は
f 2 ( p K ) 2 2 K p K log p2
K g 0 p0l

記されたそれぞれの半径である．点線は上で求めた
const. (12)となり p が f だけの関数で表せるので微分方程式(8)平衡二相流のプラントル・マイヤーの膨張曲線である．両者はよく一致し，無次元緩和時間が 0.1 より小さいときは平衡理論の膨張曲線で表される．

は解くことができ，圧力が角度 θ の関数として求め

Fig. 11 に無次元緩和時間が 10 の場合の結果を示られる．式(9)，(10)，(11)の導出は付録に示されていす．
圧力の減衰が大きく，平衡理論から離れている．

る．[4]では，相変化がある場合の関係式が示されてまた，角からの距離 r が大きいほど圧力が下がっているが，ここで得られた関係は相変化は起こらない

おり，角度 θ の関数で表すことができない．無次元
緩和時間が大きいとき相間の速度差が大きくなるの

1

で気体と液体間の摩擦が大きくなり，結果として二p/p0

相流の内部摩擦が増加し，このために圧力が減衰し0.8

Equilibrium

0.6

r/L = 1r/L = 2r/L = 3r/L = 4ていると考えられる．
７

これまでに，熱水の二相流ノズル出口に発生する

0.4

膨張波を調べる実験[4]を行ってきた．本解析は相変' = 0.1化を考慮していないので直接比較することはできな

0.2
0
90

いが，二相流膨張波の 2 次元的特徴を比べることができる．Fig. 12 は横軸にノズル出口からの距離をと
60

30

0

-30

-60

-90

[deg]

Fig. 10

沸騰水の膨張波の実験結果との比較って，ノズルの中心軸に沿う圧力分布を示したものである．ノズル入口温度は 130 ℃で，各記号は質量

Expansion curves versus the rotational
angle for a small relaxation time.

流量を変えたときの圧力分布であるが，図中にはノズルのど部の流路断面積で割った質量フラックスの形で数値が示されている．図中の各線は気体と液体の温度と速度が等しい平衡状態で膨張するとしたと

1

きの理論膨張曲線である．その理論計算ではノズル
内も平衡状態にあるとしてノズル出口の圧力が実験

p/p0

0.8
0.6

Equilibrium

値と同じ値になるまで膨張したと仮定している．実

r/L = 1r/L = 2r/L = 3r/L = 4

験値と平衡理論の膨張波の理論の大きな差は，理論値にはノズルの角から出る膨張波がはじめて中心線まで達する間フラットな圧力部が存在することであ

0.4

り，実験値にはこれが無い．また，理論膨張曲線は'=10

常に曲率が正の凹状をしているが，実験値の膨張曲

0.2

線は始め曲率が負の凸状で後半が凹状をしている．
Fig. 13 に，Fig. 12 と同じ条件である入口速度無次

0
90

60

30

0

-30

-60

-90

[deg]

Fig. 11

Expansion curves versus the rotational
angle for a large relaxation time.

元緩和時間 τ’を 10 としたときの中心線に沿う圧力分布をノズル出口からの距離 x をとって実線で示す．
この高速二相流の液滴径を 10 m のオーダーとすれ
ば緩和時間は 1 ms のオーダーとなり，無次元緩和時- 229 -C 2007 The Heat Transfer Society of Japan○

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

間は 10 のオーダーと考えられる．
図中の点線は前節８

結論で示した平衡理論から得られる膨張波による理論膨張曲線である．Fig. 11 では角度 θ で表したがこれをノズル出口に発生する二相流膨張波を理論的に調x で展開した曲線である．平衡理論では膨張波が到べるため，相間の運動量輸送を考慮した圧縮性二相達するまでにフラットな圧力部が存在するが，本解流の基礎式を導き，これを数値解析することにより

析結果にはそれが無く，Fig. 12 で示す実験値に近い．

二相流膨張波の特性を求め，以下の結論を得た．

緩和時間の大きな二相流では気液の速度差が大きく，

1．非定常基礎式から微小擾乱の伝播速度が求めらその内部損失によって圧力は減少し，フラットな部れた．この音速は緩和時間と波の角周波数の積に依分は消滅したことを本解析が示している．また，解存し，これが小さいとき運動量輸送が平衡状態にあ析結果の膨張曲線の曲率の変化も実験値と同じ傾向

る平衡音速に，大きいときは凍結状態にある凍結音を示している．

速に近づく．
2.

音速が 1 よりも小さい入口条件でも，緩和時間が小p [Pa]

20×10415×10

二相流膨張波は基礎式の支配音速である凍結

さいとき存在しうることを明らかにした．
Equlibrium
1.03×104[kg/(m2s)]
1.13×104[kg/(m2s)]
1.22×104[kg/(m2s)]
Experimental result
1.24 ×104[kg/(m2s)]
1.80×104[kg/(m2s)]
2.79 ×104[kg/(m2s)]

43.

高速二相流に現れる最初の膨張波のマッハ角は緩和時間と波の角周波数の積 ωτ の関数で与えられる音速から求められる．ただし，波の角周波数には二相流の入口速度 u0 とノズル出口の半値幅 L から得られる 2πu0/L をとる必要がある．
4.

緩和時間が小さいとき，圧力変化はほぼノズ

ル出口の角からの角度で表され，平衡理論に近い．
無次元緩和時間が 1 より大きくなると相間の運動量交換による内部摩擦が増加し，圧力は流れ方向に減衰するため，この角度だけでは表されない．

10×104
0

1

2

3

4

5

x/L
Fig. 125.

実験で得られる熱水の二相流ノズル出口の二相流の無次元緩和時間はほぼ 10 程度であるが，この条件で得られた本解析結果の膨張曲線の形状は実験

Experimental expansion curves along thecenter axis of the nozzle.

値の膨張曲線と類似している．
参考文献

1

p/p0

0.8
Calculated solution
Equilibrium

0.6

0.4

0.2
0

1

2

3

4

x/L
Fig. 13

Calculated expansion curves along
the center axis of the nozzle.

[1] 中川，”二相流エジェクタ”，冷凍，79-925 (2004)，
856.
[2] L.D. Landau and E.M. Lifshitz, Fluid Mechanics(1959), 422, Pergamon Press.
[3] 中川，他 3 名，“管内高速ミスト噴流の混合性能に関する研究”，日本機械学会論文集，60-575，
B (1994), 2552.
[4] 中川，他 2 名，
“超音速二相流ノズル出口に発生する膨張波に関する研究”，Therm. Sci. Eng.，14-4(2006), 67.
[5] 中川，杉浦，”速度緩和現象によって発生する高速ミスト流中の衝撃波”， 日本機械学会論文集，66-645，B (2000), 1495.
[6] 杉浦，中川，“ノズル内の亜音速ミスト流中に発生する速度緩和および温度緩和を伴う衝撃波”，
日本機械学会論文集，67-622，B (2001), 2462.
[7] 中川，間宮，”蒸発を伴う気ほう流中の膨張波”，- 230 -C 2007 The Heat Transfer Society of Japan○

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

日本機械学会論文集，52-474，B (1986), 650.
[8] J.B.Young，”Two-Dimensional, Nonequilibrium,
Wet-Steam Calculations for Nozzles and Turbine
Cascades,” Trans. ASME J. Turbomach., 114 (1992),
569-579.
[9] Marble, F.E.，”Dynamics of a gas containing small
solid particles ,”Proc.5th AGARD Combustion(1963), 175-215．
[10] 森岡，他 7 名，混相流体の力学 (1991), 133，朝倉書店.
[11] 松本，西川，大橋，“二相流における翼列計算（第１報，モデル化と数値計算）”，日本機械学会論文集，86-1412，B (1987), 3210.
[12] 保原・大宮，数値流体力学 (1992), 326，東京大学出版.
[13] 土方・他 2 名，日本機械学会論文集，37-295, B(1971), 556.
[14] 矢部・観山・椛島，パソコンによるシミュレー
ション物理 (1992), 138，朝倉出版[15] 森，土方，「熱と流れの工学 II」(1977), 40，共立出版．

式(10)は，αp/(1-α) が流体粒子に沿って保存されることを表している．さらに，プラントル・マイヤーの流れのように入口条件が一様な場合，流れ場全体で保存されることになり，これを K とした．
p
(A4)
K const.
1
平衡二相流の平均密度はこの定数 K を用いて圧力だけの関数として次のように表すことが出来る．
p
K
g0 p
(A5)
l
p K p0p K
これを用いると式(A3)の右辺は次のようになる．
1 p( p k) ppxi
K g 0 / p0
x
l
i
1
(A6)
p K log p
K g0 / p0xil
また，二相流全体の運動量の式(A3)の右辺は，ポテ
ンシャル流れの場合には rot(v)=0 なるので
vj

vi
xj

v 2j
xi

(A7)

2

したがって，式(A3)は積分できる形になり，付録 平衡流の関係式(9)，(10)，(11)の導出
定常の平衡流の基礎式は，基礎式(1)，(2)，(3)+(4)において，非定常項を消し vg = vl = v とすることに
よって得られ，次のようになる．
v(A1)
pv j v j
p
p j 0
xj
xj
xj
xj

vj

(1

)v j

vi
xj

1 p
xi

vj

xj

(1

)

vj
xj

xj

p
K

K log p
g 0 p0l(9)0

式(9)が得られる．この式は単相流のベルヌーイの式に対応し，かっこ内が流れ場全体で一定となる．
音速の 2 乗も式(A5)の関係から
pp

0 (A2)

p
p K

(A3)

K

g0
l

p0

2

p K
K {K g 0 / p 0

式(A1)，(A2)から∂vj /∂xj を消去することによって式(10)が得られる．
p
vj
xj 1

vi2
2

l}

(11)

となり，式(11)が得られる．
0(10)- 231 -C 2007 The Heat Transfer Society of Japan○