Paper No 7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Local Flow Measurement and Drift-Flux Modelling of Subcooled Boiling Flow in

Sub-Channel of Vertical 3 x 3 Rod Bundle

J. Enrique Julia*, Takashi Hibikit, Byong Jo Yun" and Goon Cherl Park: and Mamoru Ishiit

*Departamento de Ingenieria Mechnica y Construcci6n, Universitat Jaume I.

Campus de Riu Sec, 12071 Castellon, Spain

School of Nuclear Engineering, Purdue University.

400 Central Dr., West Lafayette, IN 47907-2017, USA

Korea Atomic Energy Research Institute.

Taejon, 305 600, Korea

Nuclear Engineering Department, Seoul National University.

Seoul, 151 742, Korea

bolivar@emc.uji.es, hibiki@purdue.edu, bi uiil kacii ic ki parkgc@snu.ac.kr, ishii@purdue.edu

Keywords: Drift-flux model, rod bundle, sub-channel, subcooled boiling flow, void fraction

Abstract

In this paper, the interfacial flow structure of subcooled water boiling flow in a sub-channel of 3x3 rod bundles is presented.

The 9 rods are positioned in a quadrangular assembly with a rod diameter of 8.2 mm and a pitch distance of 16.6 mm. Local void

fraction, interfacial area concentration, interfacial velocity, Sauter mean diameter and liquid velocity have been measured by a

conductivity probe and a Pitot tube in 20 locations inside one of the sub-channels. A total of 53 flow conditions have been

considered in the experimental dataset at atmospheric pressure conditions with a mass flow rate, heat flux, inlet temperature and

subcooled temperature ranges of 250-522 kg/m2s, 25-185 kW/m2, 96.6-104.9 C and 2-11 K, respectively. The dataset has been

used to analyze the flow characteristics of subcooled boiling flow and to study the effect of the heat flux and mass flow rate on

the local flow parameters. In addition, the distribution parameter and drift velocity constitutive equations have been obtained for

subcooled boiling flow in a sub-channel of rod bundle geometry. In the case of the distribution parameter the constitutive

equation for subcooled boiling flow in a sub-channel is obtained from the bubble layer thickness model. In this derivation an

existing constitutive equation for subcooled boiling flow in a round pipe is modified by taking account of the difference in the

flow channel geometry between the sub-channel and round pipe. In the case of the drift velocity, the constitutive equation is

proposed based on an existing correlation and considering the rod wall and sub-channel geometry effects. The prediction

accuracy of the newly developed correlations has been checked against the area-averaged data integrated over the whole

sub-channel in the 3 x3 rod bundle, obtaining better predicting errors than the existing correlations most used in literature.

Introduction

There is an increasing interest in both academia and industry

in a variety of engineering systems concerning two-phase

flows for their optimum design and safe operations. Among

them, the modelling of the two-phase flow in a sub-channel is

of great importance to the safety analysis of nuclear power

plants and verification of thermal-hydraulic design codes.

This fact is especially important in boiling water nuclear

reactors (BWRs) since the two-phase flow is involved in its

standard operational conditions.

Nowadays, drift-flux model is widely used to predict the

two-phase flow behaviour in multiple scenarios (Zuber and

Findlay, 1965; Ishii, 1977; Hibiki and Ishii, 2006). Its

importance relies on its simplicity and applicability to a wide

range of two-phase flow problems of practical interest. In

the drift-flux model, two constitutive relations are needed to

close the mathematical formulation and solve the problem.

In this regard, distribution parameter and drift velocity need

to be obtained by appropriate constitutive equations.

Recently, several works have been developed to obtain sound

and accurate distribution parameter and drift velocity

constitutive equations based on extensive dataset and

mechanistic models in different flow conditions (Hibiki and

Ishii, 2002; 2003; 2006). A complete review of the

available constitutive equations can be found in (Hibiki and

Ishii, 2006). However, most of the studies mentioned above

were concentrated on flows in round tube geometry because,

in spite of its simplicity, it is involved in many practical

applications. Though, in many of the nuclear systems, more

complex geometries like separators, fuel bundles and steam

generators are present. Often, studies performed on simple

geometries are not sufficient for the modelling needs in these

complex geometries. This fact is especially important in the

case of the drift-flux modelling in rod bundle sub-channel

Paper No

geometries, since many of the actual computational

thermohydraulic code calculations are based on the drift-flux

model.

In the existing works, modeled distribution parameter and

drift velocity were evaluated by one-dimensional data.

There are very limited attempts to develop mechanistic

drift-flux model considering phase distribution and channel

geometry effect and to validate them with local flow

parameters such as local void fraction and gas and liquid

velocities. In this paper, new constitutive equations for the

drift-flux model developed for subcooled boiling bubbly

flow in a sub-channel of rod bundle geometry are presented.

It will be shown that the constitutive equation of the

distribution parameter for subcooled boiling flow in a

sub-channel is obtained from the bubble layer thickness

model proposed by Hibiki et al. (2003). In this derivation an

existing constitutive equation for subcooled boiling flow in a

round pipe (Ishii, 1977) is modified by taking account of the

difference in the flow channel geometry between the

sub-channel and round pipe. In the case of the drift velocity,

the correlation proposed by Ishii (Ishii, 1977) for round pipes

will be modified to work in subcooled boiling bubbly flow in

a rod bundle sub-channel. The accuracy of the newly

developed correlations has been confirmed with the data

obtained in a 3x3 rod bundle sub-channel as well as with

other published correlations.

Nomenclature

A Coefficient

Ac flow channel area

Aw bubble layer area

Bsf bubble size factor

Co distribution parameter

CO,Ishil distribution parameter given by Ishii's equation

for boiling flow in a round tube

CoQ asymptotic value of distribution parameter

Do rod diameter

Db bubble diameter

DH hydraulic diameter

Dsm Sauter mean diameter

G mass flow rate

j mixture superficial velocity

n Exponent

p Pressure

peit critical pressure

Po pitch distance

Q heat flux

Reb bubble Reynolds number

Re~ bubble Reynolds number in single-bubble system

Rp radius of round tube

Ro rod radius

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

r radial coordinate measured from rod surface

T Temperature

vg interfacial velocity

vf liquid velocity

v, drift velocity

vr relative velocity

vr. bubble terminal velocity

xP, bubble layer thickness

Greek symbols

a void fraction

awp void fraction at assumed square void peak

ATsub

A

liquid subcooling

modification factor

v kinematic viscosity

p density

0o surface tension

Subscripts

g gas phase

f liquid phase

in inlet

sub subcooled

Mathematical symbols

< > averaged value

<<>> fraction weighted mean value

Acronyms

BWR Boling Water Reactor

LSTF Light Water High Conversion Reactor

THTF Thermal Hydraulic Test Facility

TPTF Two-Phase Flow Test Facility

Literature Survey

The basic form of the one-dimensional drift-flux model is

expressed by

(r)

(vs)) (a)

C' (j) + ((V),

where v, j, j a, C, j and v, are, respectively, the

gas velocity, superficial gas velocity, void fraction,

distribution parameter, mixture volumetric flux and drift

velocity, and indicate the void fraction weighted

cross-sectional area-averaged and area-averaged quantities.

The distribution parameter and void-fraction weighted

averaged drift velocity are defined as

Co and

((va) ( (1 a>)

(a) (a)

Paper No

where v, is the relative velocity between phases.

Available distribution parameter and drift velocity

constitutive equations developed for rod bundle geometries

are empirical correlations that have been derived from rod

bundle data (mainly differential pressure measurements) and

have been tested against different rod bundle databases [8] in

a wide range of void fraction conditions.

Bestion (1990) developed a correlation to be used in the

thermalhydraulic code CATHARE. The correlation is based

on experimental data in rod bundles with hydraulic diameters

of 12 mm and 24 mm and also the visual observations given

by Venkateswararao et al. work (1982). The Bestion's work

(1990) only provides the drift velocity correlation, thus a

constant distribution parameter of unity has been assumed.

Chexal et al. (1992) developed a distribution parameter

correlation based on the Zuber-Findlay's drift-flux model

(1965) and the main modifications are focused on the

distribution parameter. In their work, Chexal et al.

developed correlations for both upward and downward flows

in vertical, inclined and horizontal pipes with different fluid

types. In the case of upward water-steam flow in a vertical

pipe, the distribution parameter and drift velocity correlations

are those known as the EPRI correlation that were developed

by the same authors in 1986 (Chexal and Lellouche, 1986).

In addition, the critical pressure parameter is introduced in

the distribution parameter equation. Later, Inoue et al. (1993)

developed a distribution parameter correlation based on the

Zuber-Findlay's drift-flux model that was derived from void

fraction data in an 8x8 BWR facility (Morooka et al., 1991).

Thus, the distribution parameter and drift velocity

correlations were obtained by experimental data fitting and

taking the inlet pressure and mass flux as working parameters.

Finally, Maier and Coddington (1997) extended the

correlation obtained by Inoue et al. to a wider range of

experimental conditions and flow configurations.

Drift-flux Modelling

Co Formulation Assuming Power Law Profile

In 1977, Ishii (1977) proposed a simple model for the

distribution parameter for bubbly, slug and chur-turbulent

flow in adiabatic flow given by

CO =co,0 -(Co, i -1) ,

in this equation p. and p, represent the gas and liquid

phase densities, respectively, and Co,. the asymptotic

distribution parameter value for high void fraction.

Ishii extended the use of Eq. (4) to boiling flow by the

addition of a weighting factor that takes into account the wall

bubble nucleation and makes C -> 0 when a) -> 0 In

this case the distribution parameter is given by

C = -(co, -1) g (1-e ')),

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

where the coefficient A and Co,, are recommended to be

-18 and 1.2 for round pipes, respectively. It has been shown

that the use of Eqs. (4) and (5) can be extended to other flow

channel geometries if the parameters Co, and A are

properly modified.

In the present work, the distribution parameter for a

sub-channel of rod bundle geometry has been derived.

First the Co, value in a subchannel can be obtained

analytically by the use of Eq. (2) and the appropriatej- and

a distributions at high void fraction condition. In order

to obtain the complete Co distribution, the Ishii's equation

for boiling flow in a round pipe, Eq. (5), has been modified

in order to take into account the flow channel geometry

difference. The modification factor has been analytically

obtained by the use of the bubble layer thickness model

(Hibiki et al., 2003).

As mentioned above, the asymptotic value of the

distribution parameter can be approximated by Eq. (2),

using mixture volumetric flux and void fraction profiles at

high void fraction condition. For this purpose the

analytical forms of thej- and a distributions are needed.

The modeled sub-channel, including the coordinate system,

is given in Fig. la).

(i!

b)

Figure 1. a) Modeled sub-channel including coordinate and

b) Physical meaning of assumedj.

(ii ~

Paper No

The mixture volumetric flux, j, is assumed by

j(r,0)= (0) 1- 1 2,r'o I

11 -R, 0

r < R/4 Ro

0 < /4 (6)

where ro is the radial distance measured from the rod surface

andj,(0) is the mixture volumetric flux at a point on the line B.

Rc is defined as the distance between the origin and the point

of the intersection of the line A with line B, Po is the distance

between two rods (pitch) and R is defined as

R = 2R, R. (7)

It should be noted here that R ranges from P Ro to

2P, Ro.

Equation (6) indicates that the mixture volumetric flux along

the line A is assumed to be a power low profile with its

maximum value at (r,0) = (R ,0), namely at a point of the

intersection of the line A with the line B and zero at

(r,) = (Ro,O) namely at a point on the rod surface.

Substituting R in Eq.(6) yields

j(r,0)= j (0) 1-1

Here, jo is assumed by

j, (0) = jo I1- I1

2rcos0

P, 2Rocos0

P, 2Rcos0

(2P- 2R,) coso

where jo is the mixture volumetric flux at the center of the

sub-channel as indicated by open square in Fig. lb), namely

the mixture volumetric flux at (r,0) = (Pb/2, T/4) Eq. (9)

indicates that the maximum mixture volumetric flux on the

line A, namely the mixture volumetric flux at a point of the

intersection of the line A with the line B is assumed to be the

same as the mixture volumetric flux at (r,0) = (R, T/4) .

In other words, the mixture volumetric flux at a point

indicated by solid triangle (or open triangle) in Fig. lb) is

assumed to be the same as that at a point indicated by open

triangle (or open triangle).

Finally, substituting j,, in Eq.(9) yields

rP 2RocosO

j(r,0)= jo 1- 1- ( p_2Ro) cos0

o 2rcos0 I

P, 2Rocos0

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

In a similar way, the void fraction distribution for the

adiabatic flow in the sub-channel can be defined as

S\ 1Po 1 2RocosS

a (r,0) =a0 1- 1 2Ro)S (11)

( P,- 2R,) cos

Po 2RcosO

Thus, the asymptotic value of the distribution parameter can

be calculated by using the Eq. (2) and Eqs. (10) and Eq. (11).

The modeled mixture volumetric flux shows a good

accuracy with the experimental data (Yun et al., 2008).

The calculated asymptotic values of distribution parameter

in sub-channel as parameter of exponent, n, is shown in Fig.

2a). In the calculation, the exponent for void fraction

profile is assumed to be the same as that for mixture

volumetric flux for simplicity and the non-dimensional rod

diameter, defined as Do/Py (where Do is the rod

diameter), is 0.5. As shown in Fig. 2a), as the exponent

increases, or the mixture volumetric flux and void fraction

profiles become flatter, the distribution parameter

approaches 1.0. For n=2, the distribution parameter

reaches almost to 1.2, which is a typical value of the

distribution parameter in a round pipe. However, in real

two-phase flow in a sub-channel, the exponent may be

around 7. Unlike the case for a round pipe, the distribution

parameter in sub-channel may be around 1.04, as pointed

out by the experimental data (Yun et al., 2008). Since

30% change of n only causes a 1.5% deviation from

the value of Co, calculated by using n = 7, a slight

change of n may not affect Co, significantly. Figure 2b)

shows the dependence of the asymptotic value of the

distribution parameter with the non-dimensional rod

diameter, using n = 7 for the calculations. It is possible

to observe that Co, increases with Do/Po The

dependence is weak for all the Do/Po range, but specially

for Do/Po values below 0.7. If a value of 0.5 for Do/Po

is chosen, that corresponds to Co, of 1.04, a maximum

error of 3% is obtained for a rod bundle sub-channels with

Do/Po values between 0.2 and 0.8 calculated by using a

Do/Po value of 0.5.

1.4

1.2

1.0

" 1.0

-o

where ro ranges from 0 to Po/2cosO Ro and jo, can be

obtained by integrating j (r, 0) over the flow channel.

C -104forn-7

------ --------------- -

2 4 6 8 10

Exponent, n [-]

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Heater Rod Surface Outer Sub-channel Surface

SAll/

1.2

| 1.1 *C=104forD

$ 30

0o

1.0-------

--------------------

0.9 1 -

0.2 0.3 0.4 0.5 0.6 0.7 0.8

Non-D Rod Diameter, D/P [-]

b)

Figure 2. Dependence of distribution parameter on a)

exponent n inj-distribution and b) non-dimensional rod

diameter D0/P0.

It should be noted here that the C0,~ values obtained may be

also valid in slug and chur-turbulent flow regimes.

Co Modeling in Subcooled Boiling Flow using

Bubble-Layer Thickness Model

In order to obtain the complete distribution parameter

correlation in the rod bundle sub-channel, the bubble layer

thickness model can be introduced to obtain the modified

C0, and A parameters. The bubble layer thickness model

was successfully introduced by Hibiki et al. (2003) in order to

obtain the distribution parameter in an internally heated

annulus. In this model, see Fig. 3, the subcooled flow path

near a heated rod is divided into two regions, namely (i)

boiling two-phase (bubble layer) region where the void

fraction is assumed to be uniform, and (ii) liquid single-phase

region where the void fraction is assumed to be zero. In Fig.

3, a x, R,, a, x, and R are the local void fraction,

the radial coordinate measured from the center of the heater

rod surface, the radius of the heater rod, the void fraction at

the assumed void peak, the bubble-layer thickness, and the

coordinate of the outer part of the considered sub-channel,

respectively. Consequently, the void fraction distribution

can be assumed as

a = aw for 0r < xwp

(12)

a = 0 for xwp r R Ro.

where the r coordinate is considered from the rod surface.

a Heater Rod Surface Outer Sub-channel Surface

,,,,,,

R, Ro+x,p R

Modeled Subcooled Boling Flow

Figure 3. Basic concept of bubble-layer thickness model.

The distribution parameter for subcooled boiling flow in the

sub-channel can be calculated using Eq. (2) and numerical

integration of Eqs. (10) and (12). Unfortunately, no

analytical solution can be obtained for the distribution

parameter, so only numerical solutions will be provided in

this work. Hibiki et al. (2003) showed that the difference

in the dependence of CO on (a) between the round tube

and other sub-channel geometries may mainly be attributed

to the difference in the channel geometry. This assumption

is valid for A ,,/A values lower than 0.3, where Ap

and Ac are the bubble and channel areas respectively.

Since the product of A,,/Ac and ap is equal to (a),

A,,plA may correlate closely with (a). As a result, the

distribution parameter for subcooled boiling flow in the rod

bundle sub-channel can be obtained from Ishii's equation,

Eq. (5), taking into account of the channel geometry effect

on the distribution parameter as

C = A(1.2 -0.2 /pp )( e 18(),

where A is the modification factor defined by the ratio of

the distribution parameter for the sub-channel to that for the

round tube for the same Awp /A value given by

A W_ (x +D for rod bundle

Ac P2 n(DO/2)2 (14)

sub-channel

-= 1- 1 for round tube (15)

where Rp is the radius of the round tube.

In order to calculate the modification factor, the distribution

parameter of boiling flow for a round tube is needed.

Therefore, the void fraction and j-distributions are defined

as (Hibiki et al., 2003)

a = a, forRp x, < r < Rp,

a = 0 for 0

Paper No

Ro Ro+x,

Subcooled Boling Flow

Paper No

n+2 .j r 1

n- yRb,

From Eqs. (2), (16) and (17), we can obtain the distribution

parameter for boiling flow in a round tube analytically as

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

to those employed in previous works, and using a Co,

values of 1.03, 1.04, and 1.05 corresponding to a Do/Po

parameters of 0.3, 0.5 and 0.7 respectively (see Fig. 2b)).

Consequently, the newly developed distribution parameter

correlation can be expressed, in a more condensed way, as,

8)

Sa

nd

n- 1-- (n+2)-2

S Rp21

n1- 1- x

R{ 1 j

Figure 4 shows the modification factor obtained from Eq. (1

and the numerical integration of Eqs. (2), (10) and (12) as

function of the distribution parameter for the round tube a

Do/Po values of 0.3, 0.5 and 0.7. In order to facilitate

use, the modification factor has been approximated to

polynomial function given by

4.23-10.n4,' ,,,, +14.28 ,,,,2-10.34COsh3+2.89COshl4 for Do/Po

A= 3.06-4.54COish+4.113 ,, 2-1.88COl3 +0.34COihl4 for Do/Po = 0.

2.41-0.90CO,ihn-4.32COshl2+6.55COish3-2.68CO,ishn4 for Do/Po =0.

0.0 0.2 0.4 0.6 0.8 1.0 1.2

Distribution Parameter, COIsh [-]

Figure 4. Modification factor A.

It is possible to observe that the influence of the

non-dimensional rod diameter value on the modification

factor is more important for low C0,, values lower than

0.5. The effect of the channel geometry in the distribution

parameter has a larger impact for low void fraction

conditions as reported by Hibiki et al. (2003). If higher

CO h, values are considered the differences are

insignificant. This fact can be easily explained by the

small dependency of the Co, parameter with Do/Pf (see

Fig. 2b)). If all the CO range is considered, the

modification factor equation for DO/Po = 0.5 can be used

for a non-dimensional rod diameter range from 0.3 to 0.7

assuming a 9% error (calculated with respect to a Do/Po

value of 0.5).

The results obtained by the use of the Eqs. (13) and (19)

have been fitted to a equation with a similar functional form

1-03 0_03 1

1.04 004 1

1-05 0-05 (I

V^

e 26 3(.)7" ) for DoPo

-e-21 2( 62 ) for Do/P

e-341(a>925) for DoPo

0.3

: 0.5 (20)

0.7

Vg, Modeling Considering Confined Channel Effect

Its In order to obtain the drift velocity in the sub-channel, the

a constitutive equation developed by Ishii (1977) for

distorted-particle regime will be considered. This

correlation has been chosen since it is a simple expression

= 0.3 in which all the parameters needed are usually known and

5 that has been successfully tested against different databases

(Ishii, 1977). More sophisticated expressions canbe found

in literature, even developed for rod bundle sub-channels

(19) (Tomiyama et al., 2003). Though, these expressions need

some input parameters such as some experimental data

usually unavailable in rod bundle experiments like bubble

diameter and aspect ratio. In this work, the expression

developed by Ishii (1977) has been modified by considering

the bubble size factor, and the final expression is given as

(1-())1.75

}Pf2

where o Ap, g and Bf are the surface tension,

density difference between the phases, gravitational

acceleration and the bubble size factor respectively. The

bubble size factor, Bf should be included in order to

consider the rod wall effect in the bubble rising velocity. In

this work, the reduction factor proposed by Wallis (1969)

will be used since it has been successfully used in rod

bundle geometries (Carlucci, et al., 2004).

1 Db for- <0.6

0. 9Lm Lma

B {f = 0 b 2 (22)

0.12Db_ for Db 0.6,

SLmax ) max

where Db is the bubble equivalent diameter and L_ is

defined as,

L = ( 0P -2R0) (23)

The Lm parameter has been chosen to replace the standard

pipe diameter since it corresponds to the maximum bubble

size. The bubble larger than L,_ will be deformed by the

sub-channel walls. Its value depends on the

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

non-dimensional rod diameter and it corresponds from 1.3 to

3.6 times the hydraulic diameter if a range of Do/Po

between 0.2 and 0.6 is considered.

Experiment

Figure 5 shows a schematic diagram of the SNU (Seoul

National University) steam-water boiling loop used in the

experiments. The test facility consists of a stainless steel

pump, a preheater, a by-pass line, a flow metering system, a

pressurizer system, a secondary heat removal system and a

storage tank. The heat exchange system allows working

with a subcooling temperature range, ATb, between 2 and

11 K. The distilled water was used as a coolant in a closed

loop. The water flow rate is measured by a conner-tap

orifice plate flowmeter with an accuracy of +4.5 %. The

mass flow rate range is between 250 and 522 kg/(m2s). The

test section is composed of a 49.8 mm x 49.8 mm

cross-section, 2000 mm high channel. Figure 6a) depicts the

cross-section view of the test channel and it is composed of a

squared lattice of 3x3 heating rods with 8.2 mm diameter.

The maximum heating power of each rod is 185 kW/m2. The

hydraulic diameter of the system, DH, is 18.6 mm or 34.6

mm depending if the complete rod bundle (including the 9

rods and the outer channel walls) or only the measured

sub-channel sections are considered, respectively. The local

probes are located 1.6 m downstream of the inlet of the test

channel (z/DH =46 for D =34.6 mm (sub-channel based

hydraulic diameter)) or z/DH =86 for D =18.6 mm (whole

bundle based hydraulic diameter)). A traversing system is

used to displace the probes inside the test section. Figure 6b)

shows the coordinates of the local two-phase flow

measurement positions for every flow condition. All the

experiments were performed in bubbly flow regime at an inlet

pressure of 0.12 MPa.

Figure 5. Schematic diagram of the flow loop.

4.1 8.3 16.6

------- ----- --

units in mm

A x=5.6 x=6.9 x=8.3 B

_ t ,=o0

S...... .- ....

C -- --.---.---

---------

x4.9 ----

units in mm

+-jy=1.4

-jy=2.8

4-y=3.3

,.jy=4.2I

* y=4.9

-Jy=5.6

jy=-6.9

-Jy=-8.3

Figure 6. a) Cross-sectional view of the sub-channel test

section, b) probe locations inside the sub-channel.

Flow Characteristics of Subcooled Boiling Flow in 3

X 3 Bundle

In order to obtain a general view of the local flow parameter

idA ^ : ^ ; +-, 11- ,_ ,,- 1 ib i i h bf 3D 1-% ^ ^ ,,, 1 ;

UItsIIULIUII1 111 nLIC S~U-tlllllllC Ie, a ;Ct ULs I JJ JIUtL 1s sOIIUWI 11

Figure 7. In this way, the effect of the heat flux on the

two-phase flow local parameters can be studied. The 3D

plots have been obtained from the local data points shown in

Figure 6b and the mesh has been created from linear

interpolation between the data points. Due to the finite size

,r Compressor of the conductivity probe and Pitot tube, measurements were

SCityWater impossible close to the rod. With the aim of facilitating the

visual information the mesh has been mirrored using the

diagonal of the sub-channel as mirroring axis and the rod

Valve surface has been included. In both graphs, the void fraction,

SPump a interfacial area concentration, a, bubble interfacial

Flow Meter velocity, v, liquid velocity, vf and bubble Sauter mean

Q Thermocouple

Pressure diameter, D,,, are shown for different flow conditions.

0 Transmitter

In Figure 7 the effect of the increment in the heat flux, q,

with a constant mass flow rate, G, is shown. The values of

inlet temperature, T and liquid subcooling, AT,,, are

kept constant. It is possible to observe a bubble (or void)

layer around the rod surface where the bubble related

parameter values are more significant. An increment in the

heat flux of the heater rod produces an increment in the void

fraction layer peak value and thickness around the rod surface.

The maximum void fraction peak value is around 0.3 and the

maximum bubble (or void) layer thickness occupies the 80 %

of the sub-channel cross-sectional area. An increment in the

Paper No

Paper No

heat flux of the heater rod also produces an increment in the

interfacial area concentration. The interfacial area

concentration shows a similar trend to the void fraction, since

the interfacial area concentration is approximately in

proportional to the void fraction in bubbly flow regime. In

addition, the interfacial velocity shows a 30 % increment for

the evaluated q range and it has a similar value in all the

sub-channel. The liquid velocity is not affected by the

increment in the q value. That means that the change in the

liquid velocity produced by the density change near the rod

- -- "

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

surface is not significant compared with the bulk liquid

velocity or that this effect is located very near the rod. The

Sauter mean diameter is greatly increased with the q values,

ranging from 2 mm to 6 mm. For low q values, the

maximum D s occurs near the rod, however when the heat

flux is incremented, the Ds, shows a maximum in the center

of the sub-channel.

x- i"

x- axib. irnrnj 8

E

8 C

8 3 C (

c

4

2 y-axis,[mm]

04-

02-

0 o

8 0

4y-

2 y-axis, [mm]

1000-

500-

01

0

S 8 0

a axl tnlnlJ 8 n y -axiS,[mm]

4 --

x- axis, [mm]

2 4

4

x-axis [mm] 6

E

o 05

8 0

4

2

8 0 y axis [mm]

4 4

6 2

x-axis, [mm] 0 y axis, [mm]

E

04

S02

0

8 0

4

2 y-axis, [mm]

1000-

500

01

0

4 -

x-axs, [mm] 6

8 0

500

S

x- axis, [mm] 8 0

4

2

y- axis, [mm]

2 y-axis,[mm]

E

E .

1 -

x- axis, [mm] 8

8 0

~ E

E

Q"

C

J

8 S

2 y- axis [mm]

E

E

C5

8

8 *

6

x- axis, [mm]

x2 4

2

8 0 y axis [mm]

40

4

,xis, [mm]x- axis, [mm]

Sy- axs, [mm]

8 0

Figure 7. Evolution of two-phase flow parameters with q. First column: G=339 kg/m2s, q=25.0 kW/m2, T,,=103.4 OC,

ATb=5.0 OC; second column: G=317 kg/ m2s, q=56.2 kW/m2, T,,=101.9 OC, AT,,b=5.1 oC; third column: G=354 kg/ mZs,

q=158 kW/m2, T,,=102.5 OC, AT,,b=4.5 oC.

2

4

x- axs, [mm] 6

S4

2 y-axis,[mm]

4

x-axis [mm] 6

x-ax [mm

x- axis, [mm] 6

J

L, [, , J

ja L

, &,,- ,,,j

4

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Evaluation of Co Model

The area-averaged data obtained in the 3x3 rod bundle have

been used to check the prediction capabilities of the drift-flux

models described in Section 2. In addition, these data have

been compared with those obtained by other drift-flux models

commonly used in rod bundle geometries and mentioned in

Section 2.

Figure 8 a) shows the dependence of distribution parameter

values on the area-averaged void fraction. The distribution

parameter values in the tested conditions are always lower

than 1, which corresponds to the typical wall peaked void

fraction profile present in subcooled boiling flow. The

distribution parameter is about 0.8 at (a) = 0.02 and

gradually increases with (a) Since the distribution

parameter is zero at (a) = 0 in subcooled boiling flow, very

rapid increase in the distribution parameter is expected at

(a) < 0.02. The extrapolation of the distribution parameter

at higher (a) implies the distribution parameter about 1.04,

confirming the results given in Fig. 2a), since a

non-dimensional rod diameter value of Do/Po -0.5 was

used in Yun et al. experiments (2008). In addition, in Fig. 8

a) the distribution parameters obtained (i) in this work for a

sub-channel, Eq. (20), (ii) by Ishii's equation (1977) for a

round pipe, Eq. (5), and (iii) by Hibiki's equation (2003) for

an internally heated annulus and modified by Basar et al.

(2008), are presented by solid, broken and dotted lines,

respectively. The change in the flow channel geometry has a

profound impact in the slope of the distribution parameter for

void fraction values (a) lower than 0.1. In this way, the

slope for the annular channel is 6 times higher than the one of

the round pipe. The rod bundle sub-channel slope is between

both the annular and the round pipe ones. This fact seems

feasible, since the sub-channel flow geometry can be

considered as an intermediate case between the annulus and

the round pipe. The agreement between the distribution

parameter correlation developed in this work and the

experimental data seems acceptable. The average predicting

error is 8.01%, a remarkable value since no experimental

data was used in the modeling. Finally, two additional facts

need to be considered: (i) the data presents some uncertainties,

as pointed out in Section 4, especially for low void fraction

values due to the lack of measurements near the rod wall and,

(ii) there is only one available database (Yun et al., 2008) that

provide local data in a rod bundle sub-channel and that,

therefore, can be used to check the proposed distribution

parameter constitutive equation. It is recommended that the

validity of the proposed distribution parameter in the

sub-channel is readdressed by additional experimental data,

especially for low void fraction values to be obtained in a

future study.

1.4

1.2

S1.0

0.8

S0.6

S0.4

0.2

0.0

0.(

).L

o

S0.1

00

0.0 '

0.00

0.05 0.10 0.15 C

Void Fraction, [-]

0.05 0.10 0.15

Void Fraction, [-]

b)

Figure 8. Comparison of drift-flux model results with

experimental data: a) Distribution parameter and b)

area-averaged drift velocity.

Evaluation of Vg, Model

In Fig. 8b), the drift velocity obtained by the modification of

Ishii's equation by the wall effects (Eq. (21)) and the Ishii's

equation (1977) (Eq. (21) with B, = 1) are indicated by solid

and broken lines, respectively. Here, a constant value of Db

of 1.3 mm has been chosen in Eq. (22), since it is the averaged

value of the data used in this study. This fact generates a

source of error in the figure, but it is lower than a 10 % for all

the flow conditions. Consequently, the information given in

the figure should be taken for comparative purposes.

As shown in Fig. 8b), the void-fraction weighted-averaged

drift velocity shows a slight decrease with the void fraction as

previously reported in bubbly flow conditions (Ishii, 1977;

Hibiki and Ishii, 2002). The drift velocity constitutive

equations also show this dependence. The results obtained

by Eq. (21) provide a prediction error of 13.1%. The main

source of the error in Eq. (21) is due to the experimental

scattering observed in the data that is usual in drift velocity

measurements (Hibiki and Ishii, 2002). If the rod wall effect

is not considered in Eq. (24), Bf = 1, the prediction error is

enlarged to 19.6%, which is still acceptable prediction

accuracy.

Paper No

--

0, O Experimental data (Yun et al, 2008) -

S-- C given by Eq (20)

/ - - C given by Eq (5)

I' ...... Co given by (Hibiki et al, 2003)

O Experimental data (Yun et al, 2008)

--
- v given by Eq (21) and B l

oo

- O o

S0

0(9 .01 o

0o

0

Paper No

Evaluation of Drift-Flux Model

In this section, the prediction accuracy of the area-averaged

void fraction is discussed. As shown in Fig. 9, all existing

correlations underestimate the void fraction values, except the

one developed by Chexal-Lellouche. The compared results

by newly developed drift-flux model are highlighted by solid

symbols. The lowest prediction error (14.4%) is obtained

by the use of the distribution parameter obtained by the

bubble layer thickness model, Eq. (20) and the drift velocity

given by Ishii (1977) modified by considering the wall effect,

Eq. (21). If no rod wall effect is considered in the drift

velocity correlation, Bf = 1, the prediction error is 20.4%,

still lower than the published correlations considered in this

work. In all the correlations, the prediction accuracy is

improved for increased area-averaged void fraction where the

distribution parameter effect is more pronounced than the

drift velocity effect. The results obtained by the Bestion

correlations are remarkable since it is a quite simple

correlation that is applicable to the whole range of void

fractions. However, the Bestion and Chexal-Lellouche

correlations present high scattering for low void fraction

conditions. The predictions of Inoue et al. and

Chexal-Lellouche correlations are similar providing

area-averaged void fraction prediction errors lower than

40%. The Maier and Coddington correlations do not

provide reasonable predictions since the error in the drift

velocity estimation is very high.

1i'1

,-1 I

0.00

0.00

0.05 0.10 0.15

Void Fraction, [-]

0.20

Figure 9. Comparison of area-averaged void fraction with

experimental data.

Conclusions

In this paper, new constitutive equations for the drift-flux

model developed for subcooled boiling bubbly flow in a rod

bundle sub-channel are presented and analyzed. In the case

of the distribution parameter, its asymptotic value, C0,,, has

been obtained analytically. In addition, its dependence on

the exponent of the j- and a distributions and the

non-dimensional rod diameter value, Do/Po has been

considered and discussed. A correlation for the constitutive

equation for subcooled boiling flow in a sub-channel is

obtained from the bubble layer thickness model. In this

7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

derivation an existing constitutive equation for subcooled

boiling flow in a round pipe (Ishii, 1977) is modified by

taking account of the difference in the flow channel geometry

between the sub-channel and round pipe. In the case of the

drift velocity the expression given by Ishii (1977) for round

pipes is modified in order to consider the rod wall effects.

The area-averaged data obtained by integrating over the

whole sub-channel have been used to validate the distribution

parameter and drift velocity constitutive equations. In

addition, the area-averaged void fraction results provided by

the developed constitutive equations have been checked with

the most used correlations found in literature.

- Distribution parameter: The averaged relative prediction

error by the newly developed correlation based on the bubble

layer thickness model presents a remarkable low prediction

error of 8.01%. However, more experimental data,

especially for low void fraction values, is needed to make a

further evaluation.

- Drift velocity: The best prediction results are provided by

the Ishii's correlation modified in order to take into account

the wall effect with an averaged prediction error of 13.1 %.

If this effect is not considered the prediction error given by

the mentioned equation is 19.6%.

- Void fraction: the predicting errors provided by the existing

correlations are lower than 40% (except for the Maier and

Coddington correlation) and the best results among them are

obtained using the Bestion correlation with a prediction error

of 23.8%, however, this correlation presents major

scattering for low void fraction conditions. Using the

distribution parameter distribution developed in this work and

the drift velocity constitutive equation given by Ishii it is

possible to reduce the prediction error to 20.4%. Finally, if

the rod wall effects in the drift velocity are taken into account

the prediction error can be reduced to 14.4%.

References

Bestion, D., The physical closure laws in the CATHARE code,

Nucl. Eng. Des. 124, 229-245, (1990).

Carlucci, L. N., Hammouda, N., Rowe, D. S., Two-phase

turbulent mixing and buoyancy drift in rod bundles, Nucl.

Eng. and Des. 227, 65-84, (21 1'4).

Chexal B., and Lellouche, G., A full range drift-flux

correlation for vertical flows (Revison 1), EPRI report

NP-3989-SR, USA, (1986).

Chexal, B., Lellouche, G., Horowitz J., and Healzer, J., Avoid

fraction correlation for generalized applications, Prog. Nucl.

Ener. 27, 255-295, (1992).

Hibiki T., and Ishii, M., Distribution parameter and drift

velocity of drift-flux model in bubbly flow, Int. J. Heat Mass

Transfer 45, 707-721, (2002).

Hibiki T., and Ishii, M. One-dimensional drift-flux model for

two-phase flow in a large diameter pipe, Int. J. Heat Mass

Transfer 46, 1773-1790, (2003).

Hibiki, T, and Ishii, M. Thermo-fluid dynamics of two-phase

flow, Springer, New York, USA, (2006).

Present workEqs. (20) and (21)

S 0 Bestion (1990)

SInoue et al. (1993)

A A Chexal-Lellouche (1992)

V Maier and Coddington (1997)

Co

0 o o

o- -20o

Paper No 7th International Conference on Multiphase Flow

ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Hibiki, T, Situ, R., Mi Y, and Ishii, M., Modeling of

bubble-layer thickness for formulation of one-dimensional

interfacial area transport equation in subcooled boiling

two-phase flow, Int. J. Heat Mass Transfer 46, 1409-1423,

(2003)

Inoue, A., Kurosu, T., Yagi, M., Morooka, S., Hoshide, A.,

Ishizuka T., and Yoshimura K., In-bundle void measurement

of a BWR fuel assembly by a X-ray CT scanner: assessment

of BWR design void correlation and development of new void

correlation, in: Proc. of the ASME/JSME Nuclear

Engineering Conference, (1993).

Ishii, M. One-dimensional drift-flux model and constitutive

equations for relative motion between phases in various

two-phase flow regimes, ANL-77-47, USA, (1977).

Maier D., and Coddington, P., Review of wide range void

correlations against an extensive data base of rod bundle void

measurements, in: Proc. of ICONE-5, paper 2434, (1997).

Morooka, S., Inoue, A., Oishi, M., Aoki, T., Nagaoka K., and

Yoshida, H., In-Bundle void measurement of BWR fuel

assembly by X-ray CT scanner, in: Proc. of ICONE-1, Paper

38, (1991).

Ozar, B., Jeong, J. J., Dixit, A., Julia, J. E., Hibiki, T. and Ishii,

M. Flow structure of gas-liquid two-phase flow in an

annulus, Chem. Eng. Sci. 63, 3998-4011, (2008).

Tomiyama, A., Nakahara, Y, Adachi Y, and Hosokawa, S.,

Shapes and rising velocities of single bubbles rising though a

inner subchannel, J. Nuclear Science and Technology 40,

136-142, (2003).

Venkateswararao, P., Semiat R., and Dukler, A. E., Flow

pattern transition for gas-liquid flow in a vertical rod bundle

Int. J. Multiphase Flow 8, 509-524, (1982).

Wallis, G. B., One-Dimensional Two-Phase Flow, McGraw

Hill, (1969).

Yun, B. J., Park, G. C., Julia J. E., and Hibiki, T, Flow

structure of sub-cooled boiling water flow in a sub-channel of

3x3 rod bundles, Journal of Nuclear Science and Technology

45, 402-422, (2008).

Zuber N., and Findlay, J. A. Average volumetric

concentration in two-phase flow systems, J. Heat Transfer 87,

453-468, (1965).