91

蒸気発生器２次系における気液二相流の数値解析
Numerical Analysis of Gas-Liquid Two-Phase Flow in Secondary Sideof Steam Generator
村瀬 道雄（Michio Murase）＊
中村 晶（Akira Nakamura）＊ 八木 良憲（Yoshinori Yagi）＊

要約 加圧水型原子炉（PWR）の蒸気発生器（SG）は１次系と２次系の境界を構成する重要な二相流機器である．本研究では，熱流動解析コードPHOENICSの二流体モデルを用いて（財）原子力発電技術機構（NUPEC）で実施されたSG信頼性実証試験の気液二相流解析を実施した．沸騰
開始点を精度よく計算するには，流速分布に影響する流動抵抗と伝熱量分布を詳細に入力する必要がある．一方，流速分布や伝熱量分布は伝熱管群の上部におけるボイド率にはあまり影響しな
い．ボイド率の計算値は，沸騰開始点近傍の低ボイド率領域を除いて，局所平均で４％，断面平均で２％以内で測定値と一致した．
キーワード

加圧水型原子炉，蒸気発生器，気液二相流，数値解析，二流体モデル，伝熱管群

Abstract
The steam generator (SG) in a pressurized water reactor (PWR) is an important two-phaseflow component as the boundary between the primary loop and the secondary loop. In this study, we
performed gas-liquid two-phase flow analyses for SG reliability tests conduced by Nuclear PowerEngineering Corporation (NUPEC) using the two-fluid model of a thermal-hydraulic computer code
PHOENICS. In order to calculate the location of the boiling initiation accurately, detailed inputs wererequired for the friction coefficients affecting the velocity distribution and the heat transfer distribution.
However, the velocity and heat transfer distributions did not greatly affect the void fractions in the upper
region of the heat transfer tubes. The calculated void fractions agreed with the measured values within 4%
as the local average and within 2% as an average in a cross-section, except the region of low void fractions.
Keywords

pressurized water reactor, steam generator, gas-liquid two-phase flow, numerical analysis,
two-fluid model, heat transfer tubes

定とスリップ式を用いた数値解析が行なわれたas．

1． 緒言

Kondoらdは，SG管群内の熱流動解析用に開発された
加圧水型原子炉（PWR）の蒸気発生器（SG）は，

計算コードを使用してNUPECでのSG信頼性実証試験

１次系と２次系の境界を形成する重要な二相流機器の数値解析を実施した．前者は大型プロジェクトとであり，事故時にも１次系の冷却器として使用されして実施され，後者では専用計算コードが使用される．伝熱管群外を流れる２次系流体は，伝熱管群のている．一方，SGの２次系での二相流データが極め管内を流れる高温加圧水で加熱され，管外面で沸騰て少なかったこと，および汎用熱流体解析コードでする．伝熱管群は管支持板や振れ止め金具で支持さは二相流モデルが十分に整備されていないことから，れており，２次系では，これらの流動抵抗による流汎用コードを用いた数値解析例は見あたらない．

速分布や伝熱量分布が温度分布とボイド率分布に影

そこで，本研究では，汎用コードを使用して，基響し，ボイド率分布が流動抵抗に影響する複雑な気礎研究に基づく各種相関式を適切に選択することに

液二相流になる．

より，SGの２次系における気液二相流をどの程度の経済産業省）の委託を受けて（財）

精度で計算できるかについて検討した．NUPECで実原子力発電技術機構（NUPEC）でSG信頼性実証試験施されたSG信頼性実証試験を計算対象とし，汎用コが実施され，実機相当の条件下で２次系二相流の測ードとして英国CHAM社で開発されたPHOENICSを

通商産業省（現＊（株）原子力安全システム研究所

技術システム研究所クノバ

92

使用した．まず，単相流の領域が広い低伝熱量の実添

験ケースで管群流動抵抗の影響と伝熱量分布の影響

c:

計算値

を評価した．この結果に基づいて100％出力相当の高

f :

液相

伝熱量の実験ケースについて計算し，測定値と比較

g:

気相

検討した．

i, j :

座標方向

k:

相の番号

m:

平均値または測定値

r:

相対値

主な記号
Ai :

気液界面の面積濃度 [m2/m3]

CD :抵抗係数 [-]

Cf :摩擦係数 [-]Cp :比熱 [J/kgK]dh :水力等価直径 [m]字2． 計算方法
2.1

計算対象

dt :

伝熱管の外径 [m]Fkf :摩擦による体積力 [N/ m3]

Fkm :界面抗力による体積力 [N/ m3]

中央部５層の伝熱管群を実寸大で模擬しており，外G:

質量流量 [kg/s]

径22.23mmの伝熱管46本がピッチ32.54mmの正方格子

g :重力加速度 [m/s2]で配列されている．逆U字形の伝熱管の直管部長さはj :見かけ流速 [m/s]

約９m，最大曲率半径は1.52mである．２次側の幅はK:

熱通過率 [W/mK]

約3.1m，高さは約11m，奥行きは伝熱管ピッチの５倍L:

長さ [m]（約163mm）であり，壁面の影響を除いて２次元構造

Nt :伝熱管数 [-]になっている．１次系流体には高温・高圧水，２次p:

圧力 [Pa]系流体にはR123が使用されている．２次系流体は，

pt :伝熱管ピッチ [m]下端の両側から水平に流入し，管群を通過して気水q:

熱流束 [W/m]

分離器に流入するが，気水分離器より上方は模擬さ

Re :レイノルズ数 [-]れていない．管群は，流量分配板，７カ所の管支持

Rk :k 相の容積率 [-]板，３カ所の振れ止め金具で支持されている．本計r:

半径 [m]

算では，２次系流体の下端水平流入部から気水分離t:

時間 [s]

器下方の流出部までを対象にした．

u:

流速 [m/s]Vgj :ドリフト速度 [m/s]x :距離 [m]

Y:

水平方向の座標 [m]Z:

垂直方向の座標 [m]α:ボイド率 [-]

β:式¡6で定義される気相の容積流速比 [-]βtsp :

流路面積比 [-]

Γ:

相変化量 [kg/ m3s]ε:

計算誤差 [-]

θ:

角度 [度]

ρ:

密度 [kg/ m3]σ:

表面張力 [N/m]τ:

応力 [Pa]

計算対象の概要を図１に示す．実験装置 aはSGの

気水分離器
振れ止め金具

伝熱管本数： 46×５
伝熱管外径： 22.23mm
伝熱管ピッチ：32.54mm
最小曲率半径：55.7 mm
最大曲率半径： 1.52m

3. 1m
11 m
高温側 低温側

管支持板

Z

流量分配板

Y

２次系流体
１次系流体
図１

計算対象の概要

93

2.2

水力等価直径dhは管群内部，中心部，両端部で異なる．

基礎式

管群上端の半円部では，直交流に対する式sと式d，
本計算では汎用熱流体解析コードPHOENICSを使用した．PHOENICSには曲線座標系，直交座標系，

平行流に対する式sと式fを水平方向と鉛直方向に分解して加算した．
流量分配板，管支持板には多孔板の流動抵抗係数g

円筒座標系があり，３次元の二流体モデルが組み込まれているが，実験装置の構造から直交座標系２次元の二流体モデルを使用した．基礎式は一般に使用を使用した．
（１-βtsp）
βtsp-2
Cf,tsp = 1.4g

されている二流体モデルと同じであるため，運動量βtspは開口率である．振れ止め金具では管群との交差式のみ示す．

により多数の流路孔が形成されるため，多孔板で近似し，式sと式gを水平方向と鉛直方向に分解して∂ρk Rk uki∂ρk Rk uki ukj
+
=
∂t
∂xj

流動抵抗係数を計算した．管群部には空間容積率を入力したが，流速はセル境界で定義されており，空∂p∂Rkτk ρ++ k Rk g
∂xi
∂xi

-Rk

間面積率により多孔板での実流速が計算される．管
群上端の半円部ではセルの境界面が管群と直交せず，+Fkm + Fkf + Γk uk,intaRk [-]は体積分率，τk [Pa]は相内部のせん断応力，Fkm[N/m3]は相間の界面抗力，Fkf

空間面積率を定義することが困難であるため，空間面積率は入力せず，抵抗係数に空間面積率を含めた．

[N/m3]は管群などによる抵抗力，Γ [kg/m3s]は相変化量である．第１相（k=k１）を液相，第２相（k=２）を気相とし，液相にはkε乱流モデルを適用した．第２相の圧力は第１相にs

気液界面抗力

運動量式aにおける界面抗力Fkmは次式で定義した．
hFkm = -0.5CDρc|ur|ur Aiρc [kg/m 3 ]は連続相の密度であり ρc = ρf とした．u r

等しいとして１圧力で計算した．

[m/s]は気液相対速度である．抗力係数CD [-]は，気泡

2.3

相関式

流，スラグ流，環状流に区分し，気液界面の面積濃度Ai [m2/m3]との積で与えた．

相関式として，流動抵抗，伝熱量分布，気液界面
抗力をユーザサブルーチンに組み込んだ．

CD Ai = min｛（CD Ai）b, （CD Ai）s ｝（α<0.8）=（CD Ai）a （α> 0.8）j

気泡流にはMinatoらの式kを使用し，ドリフト速度Vg jaにはIshiiの式lを使用した．

流動抵抗

運動量式aにおける抵抗力Fkfは次式で定義した．
Fkf = -0.5Cfρk|uki| uki /Δ xis

Δ xiは流れ方向の長さである．管群に直交する流れの（１-α）g/Vg j2（CD Ai）b = ２αVg j = 1.41［ gσ（ρf -ρg）/ρf2］1/4

kl

スラグ流にはIshii and Mishimaの式¡0を使用した．

流動抵抗係数にはJakobの式fを使用した．
（CD Ai）（１-α）3s =9.80.08（pt /dt）
Cf,t =４Nt 0.044+
Re-0.15（0.43+1.13dt /pt）p｛（ t /dt）-１｝dαgs（１α4.5（α-αgs） ６
- ）
+
d（１
-αgs） dsm（１-αgs）
h¡0

dh [m]は流路の水力等価直径，αgs [-]はスラグ後流の小気泡のボイド率である．小気泡の平均直径d sm [m]式sと式dの流速に関して，機械工学便覧gでは管群にはTaylor不安定で形成される気泡径を参照して次式への近寄り流速（最小流速）で定義されているのにを使用した．

対し，伝熱工学資料hでは管群間の最大流速で定義さ

dsm =２［σ/g（ρf -ρg）］1/2¡1

れているため，最小流速と最大流速を用いて計算し，環状流にはWallisの式¡1を使用した．

圧力差を測定値と比較した．

（CD Ai）a = 0.02｛１+75（１-α）｝α0.5/dh ¡2管群に平行する流れの流動抵抗係数は，水力等価
直径を用いて管内流に対するBlasiusの式jを使用した．
Cf,p = 0.361 Re-0.25（Z /dh）f

実験圧力0.73MPaでの飽和温度における流体R123の物
性値を用いると式¡1は d sm=1.57mmになる．管群内での水力等価直径dh=32.4mmで代表し，式kと式¡0がα94=0.33近傍で繋がるようにαgs =0.02として，抗力係数m）での飽和温度373K，出口温度370Kで直線近似しを計算して使用した．

て伝熱量分布を数値計算した（１次系流体の入口温（１-α）,（CD Ai）= min｛1257α

度は376K）．この場合，総伝熱量の計算値は測定値と3（64α+75）｝
9.8（１-α）（α<0.8）

約３％で一致した．

0.5 （= 0.617｛１+75（１-α）｝α

α>0.8）

式¡4に よる近似計算をCase-A，数値計算をCase-B¡3として伝熱量分布を図３に示す．Case-A，Case-Bと式¡3に よる抗力係数を図２に示す．界面抗力の式hも高温側と低温側でそれぞれ水平方向に平均した値と抗力係数の相関式¡3は 鉛直方向だけでなく水平方を使用した．管群上端の半円部では高温側と低温側向にも適用した．

でそれぞれ平均値を使用した．Case-Aでは熱通過率

1000と２次系流体温度を一定値で近似しているため，１100ている．Case-Bの沸騰領域（高さ約５m以上）では，
上部ほど２次系流体温度が低くなるため，上部で伝気泡流
Minato（８）熱量が高くなっている．
Case-Ａ
Case-Ｂ1210
環状流Wallis（１１）

10

1
0

0.2

0.4
0.6
ボイド率 α（−）

図２

気液界面抗力係数0.81

高さ位置（m）

抗力係数 CD Ai（1/m）

次系流体温度の低下に伴って伝熱量が単調に減少し
スラグ流
Ishii/Mishima（１０）

低温側

8
6

高温側

4
2

d

伝熱量分布

0
0

熱通過率と２次系流体温度を一定と近似すると伝0.5

熱量の２次元分布を解析的に解くことができる．

11.5
伝熱量の比（−）

2

図３ 伝熱量分布（矢印：１次系流体の流れ方向）
K L
qexp（-Km x/CpG）= m m
qm CpGm｛１-exp（-Km Lm /CpGm）｝¡4

G/Gm =（Lm /L）0.5

2.4

計算条件

¡5
qm [W/m]は伝熱管単位長さ当たりの平均熱流束，Km気相の見かけ流速比 β は管群上端半円部の最外周伝[W/mK]は平均熱通過率，Lm [m]は伝熱管の平均長さ，熱管部分で定義されているa．

G m [kg/s]は伝熱管１本当たりの１次系流体の平均流計算対象とした２次系流体の条件を表１に示す．

β= jg （/ jf + jg）¡6量，x [m]は伝熱管の中心軸に沿った入口からの長さ

β= 0.7における総伝熱量は約２MW，β= 0.9におけるである．

総伝熱量は約５MWであり，β= 0.9が100%出力運転

単相流領域と沸騰領域で２次側の熱伝達率が異なり位置により熱通過率が異なる．また，２次系流体相当の条件である．
水平方向には伝熱管中心をセル境界として100メッの温度は単相流領域では伝熱により下流ほど高くな

シュ，鉛直方向には105メッシュとし，流量分配板とるが，沸騰領域では下流ほど圧力が低く飽和温度が管支持板の間を複数のセルに分割し，管群上端半円低くなる．そこで，ボイド率の測定値から沸騰開始

部では伝熱管３ピッチを１セルにした．したがって，点を求め，単相流熱伝達率にDittus-Boelterの式h，プ

ΔZ /ΔY =３〜４の縦長のセルになっている．時間幅ール核沸騰熱伝達率にKutateladzeの式hを使用し，２0.05秒で2000ステップの過渡計算を行い，定常状態に次系流体の温度を入口温度371K，沸騰開始点（Z≒５

達していることを確認して測定値と比較した．

95

計算パラメータと計算ケースを表２に示す．単相

響は小さい．各ケース間の圧力の相違は沸騰開始点

流領域での伝熱流動計算の妥当性と沸騰開始点の計とボイド率の相違により，沸騰開始点が下流になる

算精度を検討するために，低伝熱量である β = 0.7にほど静水頭が大きくなり，圧力の計算値が高くなっ

ついて，管群に直交する流れの流動抵抗の影響ている．計算では沸騰領域の圧力勾配を約20％過小（Case-１と２）と伝熱量分布の影響（Case-２と３）評価した．

について検討した．これらの結果に基づいて100％出力相当である β = 0.9（Case-４）の計算を実施した．

50

Case-３とCse-４で伝熱量分布の計算方法は同一であ分布はCase-３と異なる．
表１

計算対象の２次系流体の条件a

２次系流体

R123

気相の見かけ流速比β[-]

0.7

0.9

入口流量[kg/s]

240

174

入口温度 [K]371

371

出口圧力[MPa]

0.73

0.73

出口温度[K]

370

370

30
20

水平部

鉛直部

0
0

0.5

1

1.5

2

入口からの距離（m）
図４

２次系流体入口近傍での圧力差

計算条件

表１の
β

式s〜d
の流速

図３の伝熱量分布

Case-10.7

最小流速

Case-A

Case-20.7

最大流速

Case-A

Case-30.7

最大流速

Case-B

Case-40.9

最大流速（Case-B）*

12
高温側

10

Data（1）Case−1Case−2Case−38
64＊計算方法が Case-B と同じで伝熱量分布は異なる。

2

3． 計算結果

0

3.1

Data（1）Case−1Case−2Case−310

高さ位置（m）

表２

高温側

40

圧力差（kPa）

るが，Cse-４では伝熱量が大きく沸騰開始点と伝熱量

700

750

管群流動抵抗の影響

図５

管群に直交する流れの流動抵抗の影響が大きい２
点が（Y,Z）=（-1.5m, 0.08m）と（-0.79m, 0.81m）であ

最小流速を用いたCase-１では水平方向の圧力差を１
桁過小評価したが，最大流速を用いたCase-２とCase３では圧力差の計算値が測定値とよく一致した．したがって，管群に直交する流れの流動抵抗の計算式
sと式dでは最大流速を用いるのが妥当である．
管群に直交する流れの流動抵抗は，図４に示した

圧力分布

8
高温側

6

Data（1）Case−1Case−2Case−3

4
2
0
0

0.2

ように管群下部の水平方向圧力分布に大きく影響するが，図５に示すように鉛直方向の圧力分布への影

900

10

高さ位置（m）

を結ぶ直線上の計算値を示した．流動抵抗の計算に

850

12

次系流体入口近傍の圧力差を図４に示す．圧力測定
るため，これらの測定点と中間点（-0.79m, 0.08m）

800

圧力（kPa）

0.4

0.6

ボイド率（−）
図６

鉛直方向ボイド率分布

0.8

96

圧力差
（Pa）

温度
（℃）

高温側

流速
（m/s）

ボイド
率（-）

低温側

11m
3.1m

２次系流体
１次系流体

a液相の流速分布3.2s液相の圧力分布
d液相の温度分布
図７ Case-３における２次元分布

fボイド率分布

伝熱量分布の影響

伝熱量分布の影響が顕著である鉛直方向のボイド率分布を図６に示す．Case-１とCase-２で伝熱量分布は同一であるが，Case-２ではCase-１より水平方向の流動抵抗が大きく入口側の流速（流量）が大きくな
温度（℃）

るため同一伝熱量でも温度上昇率が小さくなり沸騰
開始点が上方になる．Case-３では図３に示したよう

ボイド
率（-）

に単相流領域で伝熱量が低下しているためCase-２より温度上昇率が小さくなり沸騰開始点が上方になる．
この結果，Case-３ではボイド率分布が測定値と最もよく一致している．管群上部では，各ケースで圧力・温度・総伝熱量に顕著な相違がなくなることから，ボイド率がほぼ等しくなっている．

3.3

高温側

低温側

諸量の２次元分布

伝熱管群に直交する流れの流動抵抗を最大流速で定義し，伝熱量分布に詳細計算（図３のCase-B）を使用したCase-３の２次元分布として，液相の流速ベ
クトル，液相の圧力分布（出口圧力基準），液相の温度分布およびボイド率分布を図７に示す．伝熱管群に直交する流れの抵抗が大きいため，２次系流体は

a液相の温度分布sボイド率分布
図８ Case-４における２次元分布

0.8
0.7
0.6
0.5
Z=8.328m0.4
測定値（1）
0.3
計算値
0.2
0.1
0
-1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6
水平位置（m）

0.8
0.7
0.6
0.5
0.4
Z=5.758 m
0.3
測定値（1）
0.2
計算値
0.1
0
-1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6
水平位置（m）

ボイド率（−）0.5
0.4

ボイド率（−）

0.6
r=1.504m0.4

測定値（1）
計算値

0.2
0
0

30

60

90
120
角 度（度）

150

180

1
0.8
0.6
r=1.056m0.4

測定値（1）
計算値

0.2
0
0

30

60

90
120
角 度（度）

150

180

0.8
0.6

（1）

測定値
計算値

r=0.593m0.4

測定値（1）
計算値

0.2
0

1.2

1.6

0

30

60

90
120
角 度（度）

150

180

1
Z=3.192m
測定値（1）
計算値

0.3

0.8
0.6
r=0.267m

0.4

0.2
0.1

測定値（1）
計算値

0.2

0

0

-1.6 -1.2 -0.8 -0.4 0 0.4 0.8
水平位置（m）
図９0.8

1
Z=4.475m-1.6 -1.2 -0.8 -0.4 0 0.4 0.8
水平位置（m）
0.6

1

ボイド率（−）

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

1.52m

ボイド率（−）

0.8
0.7
0.6
0.5
Z=7.04m0.4
測定値（1）
0.3
計算値
0.2
0.1
0
-1.6 -1.2 -0.8 -0.4 0 0.4 0.8 1.2 1.6
水平位置（m）

振れ止め金具

ボイド率（−）

ボイド率（−）

ボイド率（−）

ボイド率（−）

ボイド率（−）

97

1.2

Case-４における水平方向のボイド率分布（Z：高さ位置）

1.6

0
図10

30

60

90
120
角 度（度）

150

180

Case-４における上端の周方向ボイド率分布（r：曲率半径）98

下端の両側から中心方向上向きに流入し，下端中央

らの長さxが外側の伝熱管の下部ほど長くなり熱流束部では鉛直方向流速が低くなっているが，上部になqが低下するが，計算ではこの影響を考慮していない．

るにつれて流速分布が均一化している．下端では伝したがって，管群下部における伝熱量・温度・ボイ

熱管群に直交する流れにより水平方向に圧力勾配がド率を正確に計算するには，１次系流体と２次系流生じているが，この部分を除いてほぼ１次元的な圧体の温度差を正確に計算して伝熱量分布を求める必力分布になっている．単相流領域では，液相の温度要がある．Z = 4.475mでは，高温側のボイド率分布測は伝熱により下流（上部）ほど高くなり，流速ベク

定値がほぼ一様になっているが，低温側のボイド率

トルにほぼ直交する等温線になっている．水平方向

測定値は外側で低くなっており，計算では再現できの流動抵抗が大きいため，鉛直方向より水平方向のていない．しかし，Z = 5.758mでは高温側と低温側で温度勾配が大きくなっており，中心部に細長い高温

ボイド率が異なるものの，各々ほぼ一様なボイド率領域が生じている．沸騰を開始すると飽和温度になになっている．Z=7.04mとZ = 8.328mでは高温側と低るため，圧力分布と同様にほぼ１次元的な温度分布

温側のボイド率の差は小さく，計算値は測定値とよになり，上部ほど温度が低下している．下部の中心く一致している．したがって，本計算のように伝熱部における細長い高温領域に対応して沸騰開始境界

量分布を高温側と低温側で１次元近似しても，管群が形成されているが，ボイド率は上部ほど水平方向

下部での詳細な挙動を除いてSGの２次側における熱に均一化されている．中心部の高温領域と沸騰開始

流動特性を評価できると考えられる．

境界は中心線より高温側にずれている．これは，図図10の曲率半径 r = 0.267mとr = 0.593mでは測定範囲５の圧力勾配から推定されるように，二相流によるが狭いこともあり，測定値，計算値ともに高温側と摩擦損失より静水頭の方が大きいため，静水頭が大低温側のボイド率に有意な相違はない．r =1.056mときい低温側から高温側に向かう流れが生じることに

r =1.504mでは高温側のボイド率が低温側より若干高よる．

くなっており，流体が管群部を通過する距離が長く

100％出力相当のCase-４における２次元分布としなる中央部のボイド率が高くなっている．一方，両て，液相の温度分布およびボイド率分布を図８に示

端では管群部の通過長さが短いため，ボイド率の測す．計算条件は表２に示したように，伝熱管群に直定値は低くなっているが，計算ではボイド率を過大交する流れの流動抵抗を最大流速で定義し，伝熱量評価している．このように，ボイド率の計算値は局分布には図３のCase-Bと同様の詳細計算を使用した．

所的には測定値と異なるが，全体の傾向をよく再現

流速分布と圧力分布は図７に示したCase-３と顕著なしている．なお，r =1.504mでのボイド率の計算値に相違はない．Case-４では伝熱量が大きいため，Case-

おける局部的な低下は振れ止め金具の部分における３より下部から水平方向の温度が均一化し，下部で変化である．管群上部半円部におけるボイド率測定

沸騰を開始しており，上部でのボイド率が高くなっ

値62点に対する計算値の平均誤差は約４％であった．

ている．
1

3.4

ボイド率分布0.8

図10に示す．図９は直管部での水平方向分布であり，
図10は管群上部半円部での周方向分布である．
高さ位置Z=3.192mでは高温側と低温側でボイド率

計算値

100％出力相当のCase-４でのボイド率分布を図９と
0.6

0.4

に大きな相違があり，計算値はこの傾向を再現しているが相違を過小評価している．特に低温側で計算
値と測定値の相違が大きい．この原因として，単相0.2流領域では図８のaに示すように水平方向に大きな温度差があるが，計算では２次系流体の水平方向温00度分布を伝熱量分布に反映していないことが考えられる．また，低温側では式¡4に おいて伝熱管入口か

0.2

0.4

0.6

0.8

測定値
図11

Case-４での断面平均ボイド率1

99

局所的なボイド率の計算値と測定値との差は気液

実験では気液界面速度が測定されているが，測定界面抗力係数だけでなく，流速分布や伝熱量分布に値と計算値の対応関係の検討が必要であることからも依存し，特に下部での差が大きい．一方，各測定比較評価していない．流速は伝熱管に及ぼす流体力

断面での平均値は流速分布の影響と水平方向伝熱量の評価に重要であり，今後の課題である．

分布の影響が緩和される．断面平均ボイド率の計算
値と測定値の比較を図11に示す．測定断面は直管部

4． 結言５ヵ所（図９），上部半円部４ヵ所（図10）である．
計算値は，約９％過大評価したZ=3.192mを除いて，測定値と約２％以内で一致した．

NUPECで実施されたSG信頼性実証試験を対象に，
熱流体解析コードPHOENICSの二流体モデルを使用
して２次側の沸騰二相流を計算し，以下の結果を得

3.5

考察

た．
a 管群に直交する流れの流動抵抗に管群への近寄り

本計算ではSGに固有な相関式は使用せず管群や管速度である最小流速を用いると圧力差を１桁過小

内流に基づく相関式を使用し，管群の流動抵抗計算

評価し，管群の隙間での最大流速を用いると圧力で使用する流速の定義と，伝熱量分布の近似の影響差の計算値が測定値と一致したことから，圧力損について検討した．この結果，管群に直交する流れ

失を最大流速で定義するのが妥当である．

による流動抵抗の計算では最大流速を用いるのが妥

s 熱通過率と２次系流体温度を一定と近似して求め当であることを示した（図４）．また，液単相流伝熱

た伝熱量分布を使用すると液単相流領域での伝熱と沸騰伝熱の相違は伝熱量分布に大きく影響し，管量を過大評価し，測定値より上流で沸騰を開始し群下部での温度分布，沸騰開始位置とボイド率分布

た．単相流伝熱と沸騰伝熱を区分し２次系流体のに大きく影響するが，上部のボイド率にはほとんど

鉛直方向温度分布を考慮した伝熱量分布を使用す影響しないことを明らかにした（図６）．物性値を除ると鉛直方向ボイド率分布の計算値が測定値と比

いて流体R123に固有な相関式は使用していないため，

較的よく一致した．伝熱量分布は管群下部での熱実流体である高温高圧水にも適用できると考えられる．

流動には大きく影響するが，上部でのボイド率へ

気液界面抗力係数については，管内流に基づく相の影響は小さく，上部での二相流評価には近似に関式を使用したが，管群外での流れで管内流に近似よる伝熱量分布を使用できる．一方，下部での熱できること，液単相流領域でも平均流速が約0.6m/s

流動を精度よく評価するには水平方向の伝熱量分であることから，伝熱量分布の影響が大きい下部を布も考慮する必要がある．

除いてボイド率の計算値は測定値とよく一致した

d 気液界面抗力係数には管内流に基づくMinatoらの（図９〜図11）．式¡0に おけるスラグ後流での小気泡式（気泡流），Ishii and Mishimaの式（スラグ流），のボイド率の値 αgs =0.02には根拠がないが，式からWallisの式（環状流）を組合せて使用したが，伝明らかなように大幅に変更しない限り気液界面抗力

熱量分布の影響が大きい下部を除き，測定断面で係数とボイド率への影響は小さいと考えられる．

の平均ボイド率の計算値は約２％以内で測定値と実験では管群上部での二相流挙動が多く測定され一致した．管群上部での局所ボイド率の測定値62ており，Takaiらsによる数値解析，Kondoらdによる点に対する平均計算誤差は約４％であった．した

数値解析とも管群上部での測定値と比較評価されてがって，SGに固有の相関式を使用しなくても２次いる．伝熱量分布，鉛直方向ボイド率分布，管群下部

系での気液二相流を適切に計算できると考えられでのボイド率分布は記載されていないため，詳細なる．

比較は困難であるが，図10に示した管群上部半円部
でのボイド率の計算結果は相互に顕著な相違はない．

謝辞

本計算で使用したPHOENICSコードでは連続相と分散相を区分していないが，１圧力で計算した．し（財）原子力発電技術機構の高井睦夫氏（現たがって，連続相と分散相を区分した汎用熱流体解力（株）），岩瀬敏彦氏，宇和川誠一氏（現

関西電

三菱重工

析コードでも適切な気液界面抗力を使用すればSGの業（株）），名倉三郎氏，嶋村健吾氏，および関西電力

２次側の沸騰二相流を計算できると考えられる．

（株）の高松洋氏には有益な討論と助言をいただいた．

100（株）CRCソリューションズの松永史彦氏には計算入

定値とよく一致したが，式¡3で は流動様式が変化す力の作成とコードの使用法について協力していただ

る際に抗力係数が連続的に変化するように便宜的ないた．また，大阪大学の片岡勲教授と大川富雄講師

数値を使用した．そこで，より物理現象に対応したには，付録で述べるように気液界面抗力係数の相関

抗力係数を導出することを目的として，大阪大学大式を提供していただいた．

学院工学研究科の片岡勲教授，大川富雄講師のグループと共同研究を実施した，抗力係数の導出過程は文献複雑であるため．ここでは概要について述べる．

a（財）原子力発電技術機構，蒸気発生器信頼性実証

a

試験に関する報告書（平成12年３月）.

気液界面抗力係数の相関式
伝熱管群の外側の２次系流体は複数の伝熱管に囲s Takai, M., et al., Thermal Hydraulic Test and

まれた流路を流動することから，管内流で近似する

Verification of Thermal Hydraulic Computer Code for

ことができる．そこで，片岡教授，大川講師のグルー

Two-Phase Flow in U-bent Tube, Proc. ８th Int. Conf.

プでは，管内流でのドリフトフラックスモデルに基on Nucl. Eng.（ICONE-８）, Baltimore, No.8653

づいて気液界面抗力の相関式を導出した．導出され（2000）.

た相関式と計算コードPHOENICSで使用している気d Kondo, Y., et al., Development of Thermal-HydraulicComputer Code（Two-Fluid Model）for Steam

液界面抗力の定義式より，抗力係数は次式で表せる．
３（ρ -ρ ）
2 ｝
ρ
（１-α）
（CD Ai）=｛２αfg g/（ f Vg j ）Generator, Proc. ４th JSME-KSME Thermal Eng.

×|C1 ug - C0 uf （C| 1 ug - C0 uf ）/（|ur| ur）（A1）

Conf., B212, Kobe（2000）.

C1 =（１-αC0 ）/（１-α）（A2）

f Jakob, M., Trans. of ASME, Vol. 60, 384 （1938）.

αはボイド率，urは気液相対速度であり，Vg jとC0 は

g 日本機械学会編，機械工学便覧，丸善（1990）.

ドリフト速度と分布定数である．Vg jとC0 は，流動様

h 日本機械学会，伝熱工学資料，改訂第４版

式に依存するが，気泡流（低ボイド率）からチャー（1986）.
j Blasius, H., Forschg. Arb. Ing.-Wes., No. 131, Berlin（1913）.
k Minato, A., et al., An extended Two-fluid Model forInterface Behavior in Gas-liquid Two-phase Flow, 8thン流（中ボイド率）に対するIshiiの式lを使用すると，．
式（A1）とC0 は次のようになる（Vg jは式lを参照）
３/L ｝
（１-α）
（CD Ai）=｛αL |C1 ug - C0 u f |（A3）
×（C1 ug - C0 uf ）/（|ur| ur）
C0 =｛1.2 - 0.2（ρg /ρf ）｝｛１- exp（-18α）｝（A4）Int. Conf. on Nucl. Eng. （ICONE-８）, Baltimore,
No.8045 （2000）.

LL

=［σ｛
/ （ρf -ρg ）g ｝］0.5（A5）l Ishii, M., One-dimensional drift-flux model andL Lはラプラス係数であり，物性値のみの関数であconstitutive equations for relative motion between

る． 式（A1）は高ボイド率の環状噴霧流にも適用できる

phases in various two-phase flow regimes, ANL-77-47

が，環状噴霧流では液膜と液滴を区分して評価する（1977）.

必要があり，Vg jとC0 は液滴流量比を含む複雑な式に¡0 Ishii, M. and Mishima, K., Two-Fluid Model andなる．また，環状流に対するVg jを使用すると気液界

Hydrodynamic Constitutive Relations, Nucl. Eng. and

面抗力が小さくなり，数値計算が不安定化する．そDesign, 82（1984）.

こで，高ボイド率に対しては噴霧流に対するVg jを使¡1 Wallis, G. B., One-dimensional Two-phase Flow,
McGraw Hill, New York（1969）.

用した．
（ρg /ρf）|C1 ug - C０uf |（CD Ai）=｛α（１-α）/LL ｝
×（C1 ug - C0 uf ）/（|ur| ur）

付録

（A3'）

式（A3）と（A3'）は流速を含み複雑であるが，流路内で一様なボイド率を仮定すると分布定数はC 0 =1.0となり，気液界面抗力係数
本計算では，気液界面抗力を式hで定義し，抗力係数に相関式¡3を 使用した．ボイド率の計算値は測次のように簡略化される．
３/L
（１-α）
（CD Ai）==αL（A6）（CD Ai）==α（１-α）/LL｝（ρg /ρf）（A6'）

SG信頼性実証試験 aでの実験条件に対する２次系流101

体R123の物性値を用いると，式（A6）と（A6'）は次のようにくなったため，図示していない．流速を含まない相なる．

関式¡3と 相関式（A7）では計算が安定であり，相関式¡3と相関式（A7）でボイド率の計算値は顕著な相違がなく，３, 44.7α（１-α）］（CD Ai）= max［1277α（１-α）（A7）

測定値とよく一致している．管支持板部でボイド率の計算値が急激に変化しているが，大きな圧力損失sと圧力変動により計算が局部的に不安定化したこと

計算結果

β=0.9（100％出力）での高温側における鉛直方向によると考えられる．

ボイド率分布を図A１に，上端での周方向ボイド率分

図A１と図A２に示したように計算結果には顕著な布を図A２に示す．気液界面抗力係数に流速を含む相相違はないが，相関式（A7）では相関式¡3よ り物理モデ

関式（A3）と（A3'）を使用すると計算が不安系化し，質量保ルが明確であり，気泡流とスラグ流を区分せず広い

存とエネルギー保存が維持されず，計算誤差が大き

流動様式に適用できる利点がある．

1

12
高温側

測定値（1）
相関式（13）相関式（A7）80.8

ボイド率（−）

高さ位置（m）

10

6
4

0.6
曲率半径:1.056 m

測定値（1）
式（13）
式（A7）

0.4
0.2

2
0

0
0

0.2

0.4

0.6

0.8

1

0

30

ボイド率 （％）
図A１

鉛直方向ボイド率分布

60

90

120

150

角 度 （度）
図A２

上端の周方向ボイド率分布180