Hindawi Publishing Corporation

Science and Technology of Nuclear Installations

Volume 2012, Article ID 436142, 18 pages

doi:10.1155/2012/436142

Research Article

Analysis of the NUPEC PSBT Tests with FLICA-OVAP

Matteo Bucci and Philippe Fillion

CEA, DEN, DM2S/STMF, 91191 Gif-sur-Yvette, France

Correspondence should be addressed to Philippe Fillion, philippe.fillion@cea.fr Received 2 May 2012; Accepted 2 July 2012 Academic Editor: Alessandro Petruzzi

Copyright © 2012 M. Bucci and P. Fillion. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This paper discusses the results of a computational activity devoted to the prediction of two-phase flows in subchannels and in rod bundles. The capabilities of the FLICA-OVAP code have been tested against an extensive experimental database made available by the Japanese Nuclear Power Energy Corporation (NUPEC) in the frame of the PWR subchannel and bundle tests (PSBT) international benchmark promoted by OECD and NRC. The experimental tests herein addressed involve void fraction distributions and boiling crisis phenomena in rod bundles with uniform and nonuniform heat flux conditions. Both steady-state and transient scenarios have been addressed, including power increase, flow reduction, temperature increase, and depressurization, representative of PWR thermal-hydraulics conditions. After a brief description of the main features of FLICA-OVAP, the relevant physical models available within the code are detailed. Results obtained in the different tests included in the PSBT void distribution and DNB benchmarks are therefore reported. The relevant role of selected physical models is discussed.

1. Introduction

Based on NUPEC PWR subchannel and bundle tests (PSBT), an international benchmark has been promoted by OECD and NRC and has been coordinated by Penn State University [1]. The aim of this benchmark is to encourage advancement and assessment of numerical models in subchannel analysis of fluid flow in rod bundles, which has very important relevance for the nuclear reactor safety margin evaluation. An important database of void fraction and critical heat flux measurements in steady-state and transient conditions has been carried out by NUPEC on a prototypical PWR rod bundle. Different types of subchannel or rod bundle geometries and a wide range of flow conditions at high pressure have been investigated (see Table 1) allowing the assessment of the key models and correlations in these conditions.

The Service de Thermohydraulique et de Mécanique des fluides (STMF) at the Commissariat a l'Energie Atomique et aux Energies Alternatives (CEA), France, has been involved in the PSBT benchmark performing calculations with the FLICA-OVAP code [2]. FLICA-OVAP is an advanced two-phase flow thermal-hydraulics code based on a full 3D subchannel approach. It is designed to analyze flows in light

water reactors cores such as PWRs, BWRs, and experimental reactors. To provide a relevant answer to different core concepts and multiple industrial applications, several models coexist in the FLICA-OVAP platform: the homogeneous equilibrium model, the four equations drift flux model, the two-fluid model, and finally, a general multifield model, with a variable number of fields for both vapor and liquid phases. For each model, an adapted set of closure laws is proposed concerning heat and mass transfer, interfacial and wall forces, and turbulence.

2. The FLICA-OVAP Code

A detailed description of the FLICA-OVAP code can be found in [3]. In this section we give a summary description of the four-equation drift flux model and the relevant closure laws adopted in this analysis. Details of the models adopted to predict boiling crisis in the PSBT benchmark will be also presented extensively.

2.1. The Four-Equation Drift Flux Model. The four balance equations describing two-phase flows in the drift flux model are, respectively, the mixture mass balance equation, the

Table 1: Operating conditions of the NUPEC PWR test facility.

Quantity Range

Pressure 4.9-16.6 MPa

Mass flux 550-4150 kg/m2/s

Inlet coolant temperature 140-345°C

Wall heat flux 0.37-1.86 MW/m2

mixture momentum balance equation, the mixture energy balance equation, and the steam mass balance equation (porosities are omitted for the sake of simplicity):

(i) mixture mass conservation

d X akpk) + ^ akpk uJ = 0, (1)

at \k=v,l J \k=v,l J

where at, pk, U are the volume fraction, the density, and the velocity for the phase k.

(ii) mixture momentum balance

d (X akpk Uk ) + v- ( akpk Uk ® Uk

+ VP -V- (X akT± ) = pg + Fw,

\k=v,l J

where P is the pressure, g, the gravity, and Fw, the friction forces. The tensor Tk represents the viscous and the Reynolds stress terms for the phase k. The mixture density p is defined as

P = X akpk.

(iii) mixture energy balance

d hakpkEk )+v-(k^5„

akpkHk Uk

V-(X ak qk

\k=v,l )

qw + pg - u,

where Ek and Hk are the total energy and the total enthalpy of the phase k, qk includes molecular and turbulent heat fluxes, and qw is the volumetric source term of thermal power.

(iv) steam mass d

— (avpv) + V ■ (avpvuj - V ■ (Kc Vc) = Tv, (5)

The model is closed by a drift flux correlation and a general equation of state with the assumption that, in presence of liquid, the vapor is in saturation conditions.

<1 0.6 T3

i i i i i i i i i i i i i i i i i +20"exp / /<', ' /// ' 'Z-y :

/'/ 1 2 \ \ \

/ //, /À/ /¿7/ //M ¿B

-'M V /a □ □

gâ ■ , , , , ■ ■

0 0.2 0.4 0.6 0.8 1

Experimental void fraction

a Chexal-Lellouche □ Ishii

Figure 1: Comparison of void fractions predicted by the Chexal-Lellouche and the Ishii drift correlations against experimental data (subchannel S1 series).

2.2. The Drift Flux Correlation. FLICA-OVAP includes several Zuber-Findlay type correlations in order to estimate the relative velocity ur between the vapor velocity uv and the liquid velocity u;. The general form of these correlations is:

Uv = Co(j) + ( (vg,j)) = Co(j) + (yg,j),

where Co is the distribution parameter, (j) = auv+(1-a)u; is the area-averaged total volumetric flux, and Vg,j is the void-weighted area-averaged drift velocity. The Chexal-Lellouche correlation [4] and a correlation derived from Ishii [5] are implemented in the code. The Chexal-Lellouche correlation covers a large range of pressure, diameters, flows whereas the Ishii correlation, firstly established for bubbly, slug, and churn turbulent flow in adiabatic conditions, allows taking into account nucleate boiling on heated walls.

2.3. Pressure Drop. Friction forces Fw are given by the sum of distributed Ffric and singular pressure drops Fsing

Fw Fsing + Ffr

Singular friction due to spacer or mixing grids or other pressure drops are given by

1 pK. ||u||u,

2r =sing" 11 '

where K. is an antisymmetric tensor. Distributed friction

at walls is instead accounted for by

Table 2: Range of parameters for the Shah model covered by the NUPEC tests.

G (max/min)

x¡ (max/min)

xcrit (max/min)

Series 0 Series 2 Series 3 Series 4 Series 13 Series 8

0.009711 0.009711 0.009989 0.009711 0.009711 0.008867

376.67 376.67 366.21 376.67 376.67 412.55

4944.4/1408.3 4769.4/316.7 4702.8/1361.1 4725.0/566.7 3861.1/1361.1 4816.7/575.0

-0.062/-0.96 -0.053/-0.97 -0.045/-0.87 -0.058/-0.98 -0.057/-0.50 -0.058/-0.97

0.49/-0.20 1.00/-0.06 0.58/-0.05 1.00/0.06 0.50/0.20 1.00/0.11

£ 0.6

+2ffexp Á /</< ' / // ' //,'

/ / T// S / / z, uexp

/ /¿X, ' / x.í ' ////% /

y <M. •Jy 0 □ o v

yga-og &

0.2 0.4 0.6 0.8

Experimental void fraction

□ Series 1 a Series 2

o Series 3 Series 4

elat 400

200 400 600 800

Experimental density (kg/m3 )

Series 1 Series 2

Series 3 Series 4

Figure 2: Calculated void fractions (a) and densities (b) versus experimental data for subchannel series 1 to 4 with the Chexal-Lellouche drift correlation.

where the different friction terms fw are given by the product of the isothermal friction factor fjso, the heating wall correction /heat, and the two-phase flow multiplier f2,. The isothermal friction factor in turbulent flows is given by

fiso = 0.194Re

whereas the heating wall correction was estimated by an in-house model already used in the FLICA-4 code [6]. In this analysis the Chisholm correlation [7] was used for f2,, given by

/2, = 1 + (,2o - 1)( 1 + Df V) , (ID

where ,2o is the adiabatic two-phase frictional pressure drop multiplier, C,, a parameter accounting for heat flux, Dheat is the heated diameter, and q", the wall heat flux. In the runs of the PSBT benchmark C, was set equal to 0 and thus f2, = ,2o.

2.4. Diffusion Effects. To account for viscous and turbulent diffusion effects, the tensor Tk is introduced in the momentum equation, given by

„(1+«;i)(|k+dxj-3 s <12)

3, ^ dxi

l=x,y,z

where pkM'j is the turbulent viscosity, which is limited to the liquid phase. An anisotropic formulation is used for turbulent viscosity

Mj = M&Re -Ref , (13)

where Re = GDh/pi is the Reynolds number, M'/0, bM, Ret are parameters, and /M(/2lf) is a function of the two-phase flow multiplier.

~i-1-1-1-1-1-r

I 0.3 -

Subcooled 1 Saturated

J_■ ■■

: = Kv0 increase

t-1-1-r

0 0.5 1 1.5 2 2.5 3 Elevation (m)

-"-i Experimental data

KV0 = 1 X 10-4 --- kv0 = 1.5 X 10-4

...... kv0 = 2.5 X 10----kv0 = 5 X 10- 4

act 0.3

< i , y

i. / / ■; ■kts increase //

; ' //?'

! / i / ' .■/

-X-l-L'-

Subcooled i Saturated

i , / .7

! Hi! i i i.i

! ¡4 i'

' 'i i

j---1----l

/ /:' ¿/■A

0 0.5 1 1.5 2 2.5 3 Elevation (m)

|-"-H Experimental data

---kts = mts = 0.01

--- kts = mts = 0.025

■■■ kts = mts = 0.045

■ - kts = mts = 0.06

Figure 3: PSBT run 5.2442. Sensitivity analysis for Kv0 (a), kts and mts (b) on the void fraction profile in the central subchannel.

Table 3: Choice of X and K for the Katto and Ohno model.

R < 0.15

R > 0.15

< X2 X — X1

X1 <X5 X — X1

X1 >X2 X2 < X3

X1 >X2 X2 > X3

X — X2

X — X3

X1 > X5 X5 > X4

X1 > X5 X5 < X4

X — X5

X — X4

K1 > K2 K — K\

K\ < K2 K — K2

K1 > K2 K — K1

K1 < K2

K2 < K3

K1 < K2 K2 > K3

K — K2

K — K3

Similarly, molecular and turbulent heat fluxes are given

X akqk — (1 + K( l) Vhx,

k—v,l

where hx = xhv + (1 - x)hi is the flow enthalpy based on the flow quality x. The turbulent conductivity is given by

Kj — K0 (Re -Re t^frf),

where Kg,, bK, Ret are parameters and fK(f2f) is a function of the two-phase flow multiplier. In this analysis, the value

t0 and Ko was varied as a function of the axial position. Two different values have been adopted for each parameter depending whether the considered axial position was downward a mixing or a spacer grid.

2.5. Wall Temperature. Wall temperatures are estimated on the basis of the bulk temperature and the heat transfer coefficient as

Tw — Tb + \ — Tb + Nul/Dh.

The Nusselt number and the bulk temperature depend on the heat transfer regime. Four different regimes can be distinguished: single-phase convection heat transfer, subcooled nucleate boiling (SNB), saturated nucleate boiling (SANB) and post-critical-heat-flux (post-CHF) heat transfer. In single-phase heat transfer and SNB, the bulk temperature is equal to the liquid phase temperature, whereas in SANB it is equal to the saturation temperature. In forced convection conditions, the single-phase heat transfer coefficient is estimated by the Dittus-Boelter correlation. The onset of significant void (OSV), which is the transition between single-phase heat transfer and SNB can be predicted with the Forster and Grief correlation in low-pressure conditions or the Jens and Lottes correlation in high pressure conditions. In the present analysis, the Jens and Lottes correlation was used [8], which allows estimating the minimum wall superheat ATsat demanded to achieve net vapor generation

ATsat — 7.91 ( qr

Vapor generation starts when wall temperature estimated with Dittus-Boelter correlation exceeds this value. Finally, in post-CHF conditions, the choice of correlation depends

Experimental void fraction Experimental void fraction

□ Series 5 Series 6 o Series 7

□ Series 5 A Series 6 o Series 7

Experimental void fraction

□ Series 5 A Series 6 o Series 7

Figure 4: Steady-state rod bundle exercises. Comparison of calculated void fraction profiles in the central subchannel against experimental data: lower (a), medium (b), and upper (c) elevations.

Science and Technology of Nuclear Installations Table 4: Range of parameters for the Katto and Ohno model covered by the NUPEC tests.

Z [m] D [m] Z' R' (max/min) W' (max/min)

Series 0 3.658 0.009711 376.67 0.159185/0.070185 1.28 X 10-6/3.97 X 10-8

Series 2 3.658 0.009711 376.67 0.158337/0.030672 5.04 X 10-5/4.13 X 10-8

Series 3 3.658 0.009989 366.21 0.156657/0.070549 1.33 X 10-6/4.28 X 10-8

Series 4 3.658 0.009711 376.67 0.156156/0.031256 1.35 X 10-5/4.17 X 10-8

Series 13 3.658 0.009711 376.67 0.156156/0.096048 6.76 X 10-7/8.28 X 10-8

Series 8 3.658 0.008867 412.55 0.156323/0.031191 1.48 X 10-5/3.99 X 10-8

Table 5: Geometry and power shape of rod bundle test assembly (from [1]).

Assembly

(subjected subchannel)

Subchannel type Number ofheaters Axial-heated length (mm) Axial power shape Number ofruns

Center (typical) 4 X 1/4

1555 Uniform 43

Center (thimble) 3 X 1/4

1555 Uniform

Side 2 X 1/4

1555 Uniform 20

Corner 1 X 1/4

1555 Uniform 20

on the boiling characteristics, whether it is IAFB (inverted annular film boiling) or DFFB (dispersed flow film boiling).

2.6. The Mass Transfer Term rv. The mass transfer term rv appearing in the steam mass balance equation is given by the sum of two contributions: the vapor generation on walls rwv and the mass transfer between the liquid and the vapor phase rv;. In subcooled nucleate boiling, only a portion \v of the heat flux transferred from the wall to the mixture is used to vaporize the liquid phase, whereas the remaining part is used to heat the liquid phase up.

The vapor generation at walls is thus given by

Xvg" 4

hlv Dheat'

where \v is a function of the saturation temperature, the liquid phase temperature, and the wall superheat demanded to have subcooled nucleate boiling, given by

Tw,ic Tsat AT Si

Ti - AT s

If the mixture is in saturation conditions, it is \v = 1.

The mass transfer between the two phases rvi is instead given by

hv - hi'

where hv and hi are vapor and liquid enthalpy and qvi is the heat transferred between the two phases, given by

qvi = Kvo

log(1 +Re /Reo)

f(P, p, pi, u, Ur

pc(c* - c)

where Kv0 is a mass transfer parameter, G is the mass flux, Re = GDh/pi is the Reynolds number, c*, the equilibrium quality based on the mixture enthalpy h = chv + (1 - c)hi, f (P,p,pi, u, ur), a function depending on local conditions, Dh, the hydraulic diameter, and pi, the liquid viscosity. Re0 is a parameter of the model.

3. Prediction of the Boiling Crisis

To predict boiling crisis conditions, several models are currently available within FLICA-OVAP. The W3 correlation [9] is appropriate to predict departure from nucleate boiling (DNB) phenomena of interest for pressurized water reactors. However, boiling crisis experienced in NUPEC tests can be close to dryout conditions. Two models have been therefore implemented in the code: the model of Shah [10] and the model of Katto and Ohno [11]. Both models can predict critical heat flux values for both DNB and dryout conditions. An accurate description of them is given in the following sections.

3.1. The Shah Model. The Shah model [10] consists of two separate correlations to determine the boiling number Bo, defined as

The first correlation covers conditions where the critical heat flux depends on the upstream conditions, named UCC (upstream conditions correlation). The second, named LCC (local conditions correlation), depends only on local quantities.

Table 6: Geometry and power shape of rod bundle test assembly.

Item Data

B5 B6 B7

Rods array 5 X 5 5 X 5 5 X 5

Number of heated rods 25 25 24

Number of thimble rods 0 0 1

Heated rod outer diameter (mm) 9.50 9.50 9.50

Thimble rod outer diameter (mm) — — 12.24

Heated rod pitch (mm) 12.60 12.60 12.60

Axial-heated length (mm) 3658 3658 3658

Flow channel inner width (mm) 64.9 64.9 64.9

Pattern A Pattern B

Radial power shape

0.85 0.85 0.85 0.85 0.85

0.85 1.00 1.00 1.00 0.85

0.85 1.00 1.00 1.00 0.85

0.85 1.00 1.00 1.00 0.85

0.85 0.85 0.85 0.85 0.85

0.85 0.85 0.85 0.85 0.85

0.85 1.00 1.00 1.00 0.85

0.85 1.00 0.00 1.00 0.85

0.85 1.00 1.00 1.00 0.85

0.85 0.85 0.85 0.85 0.85

Axial power shape

Number of mixing vane spacers

Number of no mixing vane spacers

Number of simple spacers

Mixing vane spacer location (mm)

(Loss coefficient 1.0 [1])

No mixing vane spacer location (mm)

(Loss coefficient 0.7 [1])

Simple spacer location (mm)

(Loss coefficient 0.4 [1])

Uniform

Cosine

Cosine

471, 925, 1378, 1832, 2285, 2739, 3247 from the beginning of the test section 2.5, 3755 from the beginning of the test section 237, 698, 1151, 1605, 2059, 2512, 2993, 3501 from the beginning of the heated section

Table 7: Test series for void fraction measurements. Table 8: Average absolute errors in steady-state DNB tests (%).

Test series Assembly Test mode 1D Shah 1D Katto 3D Shah 3D Katto

Steady-state Transient and Ohno and Ohno

5 B5 74 runs 4 runs Series 0 15.5 10.9 13.7 5.6

5T Series 2 13.5 16.0 12.3 18.3

6 B6 74 runs 4 runs Series 3 8.5 10.9 10.0 15.0

6T Series 4 14.8 14.6 13.7 19.1

7 B7 74 runs 4 runs Series 8 11.5 15.6 14.3 12.9

7T Series 13 4.8 16.4 4.6 19.6

8 B5 31 runs

3.1.1. The UCC Correlation. The UCC correlation is

( D\ 0 89 /104 Bo = 0.124^ff i^M (1 - xff, (23)

" GDCpjr - ' G2 ' 0.4

[ ¿I \ _PfgD_ .Pg.

In the previous equations, zeg is the effective tube length and x;>eff is the effective inlet quality, defined as

Xi,eff = Xj, Zeff = Zcrit, for Xj < 0, Xj,eff = 0, Zeff = Zsat, for Xj > 0,

where zcrit is the distance from the inlet section and the location of the boiling crisis, where the critical heat flux is

Flow reduction

31 32 33 34 35 36 37 38 Time (s)

Time (s)

Exp. lower Exp. middle Exp. upper

- Calc. lower

Calc. middle - Calc. upper

Depressurization

Exp. lower - Calc. lower

Exp. middle Calc. middle Exp. upper - Calc. upper

Temperature increase

0.6 0.5 0.4 0.3 0.2 0.1 0

20 40 60

Exp. lower Exp. middle Exp. upper

80 100 120 140 160 178 Time (s)

- Calc. lower

Calc. middle - Calc. upper

90 100 Time (s)

110 120

Exp. lower Exp. middle Exp. upper

- Calc. lower

Calc. middle Calc. upper

Figure 5: Transient T5: comparison between calculated and measured void fractions.

calculated, and zsat is the boiling length, which for uniformly 3.1.2. The LCC Correlation. The LCC correlation is expressed

heated tubes is given by

zsat _ Zcrit, xi

D = D 4Bo ■

For water flows, when Y < 104, n = 0, otherwise it is given by

(1 - Xiff0

for Y < 106, ■ for Y > 106

Bo = FEFXBo0 The entrance factor FE is given by

Fe = 154 - 0^032 (D) ■

When the previous correlation gives FE < 1, it is used FE = 1.

100 120 Time (s)

Exp. lower Exp. middle Exp. upper

- Calc. lower

Calc. middle Calc. upper

100 120 Time (s)

Exp. lower Exp. middle Exp. upper

Calc. lower Calc. middle Calc. upper

Figure 6: Temperature increase transient T6: original (a) and delayed (b) inlet temperature boundary conditions.

Bo0 is defined as the boiling number at xcrit = 0, given by Finally, F2 is given by

the maximum value obtained by the following equations

Bo0 — 0.082 Y-a3(1 + 1.45Pr4.°3), Bo0 — 0.0024Y-a105(1 + 1.15Pr139),

where Pr is the reduced pressure given by Pr = P/Pc where Pc is the critical pressure, that is 22.064 MPa for water. The value of FX depends on the quality at the location of the boiling crisis xcrit. When xcrit > 0, the following equation is used

Fx — F3

F3-a29 - 1)(Pr - 0.6) 035

When Pr > 0.6, c = 1, otherwise c = 0. F3 is instead given by ' 1.25 X 105 \a833xclit

F3 —

When Xcrit < 0, Fx is given by

FX — F1

(1 - F2)(Pr - 0.6) 0.35

As in the previous case, when Pr > 0.6, b = 1, otherwise b = 0. Fi is given by

F1 = 1 + 0.0052 (-x0rit8) Ya41 for Y < 1.4 x 107,

F1 = 1 + 0.0052(-x0^8t8) (1.4 x 107)a41 for Y > 1.4 x 107.

F2 — F- for Fi < 4, F2 — 0.55 for F1 > 4.

3.1.3. Choice between UCC and LCC Correlation. For water, the UCC correlation is used when Y < 106. When Y > 106, the correlation giving the lower value of the boiling number Bo is used, with exception of cases where zeff > 160/Pr114, for which the UCC formulation is always adopted. In this analysis, the UCC correlation was used in all calculations.

3.1.4. Range Covered by the Shah Correlation. The Shah correlation was tested against 62 experimental databases with 23 different fluids, covering the following operating conditions:

(i) 0.315 x 10-3 <D< 37.5 x 10-3 m,

(ii) 1.3 < z/D < 940,

(iii) 4 <G< 29051 kg/m2/s,

(iv) -4 <xi< 0.85,

(v) -2.6 <Xcrit < 1.

Only 15 tests over the whole database are not included in the range covered by the Shah correlation. In Table 2 the values of experimental quantities relevant for the Shah model are reported.

Bo0 — 15Y-a612

1 . 1 1 1 1 1 1 1 ' f ' y / / & /C <§> /eP / o y / y / » <' //

+ 10% 0 <3/ t§>/ r. y ' /o°/ -10% / ' / /

® 0 OB ¿ 0 m J2 y Or y ?'// / / y y

/ / y '' <y // w / y 1 1D Shah model |

/// /// y/' y, s , , ,

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 0

1 1 1 1 1 1 - i " o O 0 0 5/o / / ' / S s / ' 0 / '' -/ / / /•

+1 0 0% é° / ' ' / ' <y / y y y ' -10% .

y GD ,"©/ / y y / ' y y

/ y^ // y ' / ' á' \y

\ 3D Shah model |

// y/y / - jT y , , ,

Steady-state DNB-series 0

1 1 1 1 1 1 , , ✓ / 1 y • y / °y /0 /<*> / ° . / oy <sy 0/

+ 10% Q.'' / A?/ & A y 0 0 <6 y / J ✓ y y -10%

© ®as 0 jéf /b<y y G -G y y y y y y / y y y

/ 0/7 0„ //, 0 O ' y ' F y y

\ 1D Katto model 1

/// y/y y ■ , , , , , , , ,

500 1000 1500

Experimental (kW/m2) (b)

Steady-state DNB-series 0

1 1 1 1 1 1 1 1 1 / p/o —1—/ /&> / / 0 / ■ / / / y G / G X ^ ■ / / y /

+10 0 / / / % / /Os ,'oy y 0/ y y y -10%

O/Jf py / M' y y

/ á/-y y y A//

3D Katto model

/yy y/y y , , ,

1000 1500 2000 2500 0 500 1000 1500

Experimental (kW/m2) Experimental (kW/m2)

(c) (d)

Figure 7: Computed versus experimental critical heat fluxes (Series 0).

3.2. The Katto and Ohno Model. The correlation of Katto and Ohno [11] provides the value of the critical heat flux as

^hf = XG(htv + K AhSub,i), (36)

where Ahsub,; is the subcooling inlet enthalpy. The terms X and K are functions of three-dimensionless terms:

R' = ^, Pf

aPf (G2z) ■

Five different values of X must be determined

CW'° 043

X1 = X2 =

X3 = X4 = X5 =

Z' 0 1R

'0 133 W'a 333

1 + 0 ■ 0031Z' 0 ■ 098R'a 133 W 'a 433Z'a 27 1 + 0 ■ 0031Z ' ,

0 ■ 0384Rft 6W 'ft 173

1 + 0 ■ 28 W'0233Z' ,

0 ■ 234R'0 513 W 'ft 433Z'°-27 1 + 0 ■ 0031Z '

Steady-state DNB-series 2

/ / / / o/ Z / 0/

+1 / 0% / / / / £ / / o' 3 ° -10%

AÁ / / o ✓ /TV /./ ® ''' oG O 0 o

<2 / / . / / / A OT O 00 o

3D Shah model

//y //fa / V

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 2

, , , , , , , , , / / / 1 / v / / / / / / / ✓ / /é

+1 / / 0% / / / / / / A / ® Vo g . -10%

/ / / nS / A ex / / / AA** // ° fe QD

G / / / / / / / /A 0* < ° o

3D Katto model

/GTV //O / , , ,

500 1000 1500

Experimental (kW/m2)

Figure 8: Computed versus experimental critical heat fluxes (Series 2).

3 1000

The value of C in the relationship of X1 is given by Three values of K must be determined

C = 0.25 for Z' < 50,

KZ' < 150, (39)

1.12[1.52W 'a233 + ( 1/Z')]

K = cW^, (40)

0.833[0.0124 + (1/Z ')] (41)

, , K2 R'0.133 W'0.333 , ( )

C = 0.25 + 0.0009(Z' - 50) for 50 < Z' < 150, (39) R W

C = 0.34 for Z' > 150. K3 =-R'0.6W'0.173-. (42)

1 1 1 1 1 1 ! c/°/ / / / :

+1 / o/ //a A -10% ■ / a/ ' À°A

/ / A Y A 0 O s

/y // y y/y o //'' 0 yy° 8 | 1D Shah model |

// /// / ■ JT y , , ,

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 3

1 ! 1 1 1 1 1 1 1 1 1 1 fO 1 > / / / / . / /4 , / n/o O 0O- / r 0

/ / +10% / / 0 / / Os' -10%

/ / <y A> / /> ' y 'VA / G ' O

/ / A / /Ü, / / ' ' / ' //'' * 0

3D Shah model

/y S ///

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 3

i . i i i 1 1 1 / / / V ' J / / / / . / / / / / y -> / ✓ /¿Io '' /ODD,/

+1 / / 0% '''/¿¡F '/rjf0 / CX3 o '' 0 /6 o - -10%

/ //' y/,® //& /A> o

/ /AA //y sê 0

1D Katto model

//s //s y /, s , , , , , , , , ,

500 1000 1500 2000

Experimental (kW/m2) (b)

Steady-state DNB-series 3

i i i I / / , 7 , , / / . -a y / O / ^OO

+1 / / / y 0% / / s / s ' y '„ ' / y o /<®aé o / o - D ° -10% 0

/ / A / / y A/ A AA/o0 y°yw / S io

// A' 0D®> o S>

3D Katto model

// /y/ 7 s , , ,

500 1000 1500

Experimental (kW/m2)

Figure 9: Computed versus experimental critical heat fluxes (Series 3).

5 1000

The appropriate values of X and K must be determined according to Table 3.

3.2.1. Range Covered by the Katto and Ohno Correlation. The Katto and Ohno correlation covers the following operating conditions:

(i) 0.01 < z < 8.8 m,

(ii) 0.001 <D< 0.038m,

(iii) 5 < Z' < 880,

(iv) 0.0003 <R' < 0.41,

(v) 3 X 10~9 <W' < 2 X 10~2.

All experimental tests are included in the range covered by the Katto and Ohno correlation (see Table 4).

4. Void Distribution Results

4.1. Subchannel Exercises. Four series of measurements of void fraction were performed in sections representative of subchannel types in a PWR assembly. Table 5 gives the

Steady-state DNB-series 4

500 1000 1500 2000

Experimental (kW/m2)

1 . .........)' y ' / ' / / / &■ /■ / s / / tí O' X <

+ 10%/' /Qc ' / / <po <8 so® ° " -10% ■

/ / / / ' / ' /A / O / QíC /./. <§> %

/ QÚ / / \ //Or 'AÏ' /

1 1D Katto model |

.......// /// / > , , ,

500 1000 1500 2000

Experimental (kW/m2) (b)

Steady-state DNB-series 4

d 1000

i i i 1 1 . 1 1 1 1 1 1 y o / y / 0 ¿ f 'o/ /Á\ /% y'' -

+1 / / 0% / /0 x / y''® O / y o o -10% o

/ / dy / 7 / / ' / ' // '/ y /o V° ¿CP

/ / / / ' / ' / / / y6 ■O o o

\ 3D Shah model |

//y /// / Jr / , y , , , , , , , , , , , ,

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 4

/ / / / ■ y / y / y / ; y / y y / y ' / y.' / /

+1 / / y y y' O O 0% y' / /o ° o -10% 10%

/ /'/y '/y' é> y/y' o ✓ ' o° ° 0 0 0 o ; . d

o £ 'o o o 3D Katt o model |

500 1000 1500

Experimental (kW/m2)

Figure 10: Computed versus experimental critical heat fluxes (Series 4).

5 1000

geometry (center, typical, and thimble subchannels, side and corner subchannels) and power shape for these series. The heated length is 1555 mm, and the void measurement section is located at 1400 mm from the bottom of the heated section.

According to the given test matrix (power value and boundary conditions), the flow is saturated at the location where the void fraction is measured for the most runs. This situation allows only partial assessment of the subcooled boiling models. In contrast, drift flux model plays a major role to obtain good agreement with experimental data.

Chexal-Lellouche and Ishii correlations have been compared for series 1 (Figure 1). It is shown that the Ishii correlation tends to underpredict the void fraction more than the Chexal-Lellouche one. Even if it is not obvious to extend this conclusion to other configurations, the Chexal-Lellouche correlation has been thus adopted for other subchannel and bundle calculations.

Calculated densities and void fractions are compared against the experimental data in Figure 2. Experimental densities measured by CT scan are declared to be affected

1 1 1 1 1 1 1 1 1 1 1 1 / 1 y 1 J % / /°r>é> • y ✓ y / JfS> / / / / /

+1 / £ rs 0%®°/ am> J> y6 '' b/m y y ' -10% -

A A /5 / „/ / / w y /// / / Y

/ / / / / / / ' / ' /

\ 1D Shah model |

// /// / , y , , ,

500 1000 1500 2000

Experimental (kW/m2)

Steady-state DNB-series 13

500 1000 1500

Experimental (kW/m2)

Steady-state DNB-series 13

/ / / / / / / / / / / / / / / /

+ / 10% / / / / //Y'*» / "i» G eP -10% 33 °o

/ / / / / y'/iy'b / ' cP

/ / / / / ' / ' ' / ' / // /

1D Katto model

/// / jp y

500 1000 1500 2000

Experimental (kW/m2) (b)

Steady-state DNB-series 13

500 1000 1500

Experimental (kW/m2)

Figure 11: Computed versus experimental critical heat fluxes (Series 13).

5 1000

by an uncertainty of 15kg/m3 (ffexp = 15kg/m3 in Figure 2(b)) [12]. Experimental void fraction resumed basing on measured densities are presumed to be affected by an uncertainty of 3% of void fraction (ffexp = 3% in Figure 2(a)) [12]. It happens that densities obtained by the code match well the experimental ones, whereas void fractions are slightly underestimated for the lowest value of this quantity. Indeed, physical properties and correlations adopted to convert density to void fraction could be the cause for this discrepancy.

4.2. Subchannel Exercises. A partial section of the full length 17 x 17 type PWR fuel assemblies was considered. The rod bundle test is a 5 x 5 square array. The heated section is 3658 mm high, and density measurements are set at 2216 mm (Lower), 2669 mm (Middle), and 3177 mm (Upper). In Table6 are shown the geometry and power shape of the rod bundle test assemblies considered. Three assembly configurations were considered: B5, B6, and B7. The main difference between these three configurations is in the power distribution. Assembly B5 has a uniform axial

Steady-state DNB-series 8

1 . 1 1 1 i i i o-/ y ? y /6000,'

+1 / / / /^O // s // s S re(G ' ß 9' / /<9 yó -10% o- 0 g g

/ / y / y / g( s /y ,'m> '/Â ° ' o g ° ° O

/ / / / , G

1 1D Shah model |

/// /// /// y o , , ,

500 1000 1500 2000

Experimental (kW/m2)

Steady-state DNB-series 8

1 1 1 1 1 1 1 1 1 ■ çefo / ■ o o g,«* q/ , ° é / / a/ "À /

+1 y y §y g / 0% /V 0 / ✓ ? S S 0 -10% g

/ °y' A //0, / y* ' ' y ' y /<§>/ / S g.® y

/% / / / ' y ' y

\ 3D Shah model |

y y y/y //y y/' ■ / / , y , , ,

500 1000 1500

Experimental (kW/m2)

500 1000 1500 2000

Experimental (kW/m2) (b)

Steady-state DNB-series 8

1 . / / / / 1 / v / / / >6 ■ y™ Qf / / / °

+1 / / 0% y gCa yyy /^0 0 ,0 0 c -10% 0 00 0 0 D

/ y / y / X> / / / / //vé y / & 0 G

/ / y / / s S S / / / // y 0

1 3D Katto model | '

/// //> /// / >

500 1000 1500

Experimental (kW/m2)

Figure 12: Computed versus experimental critical heat fluxes (Series 8).

5 1000

power distribution and pattern A radial power distribution. Assembly B6 has the same radial power distribution of assembly B5 but cosine axial power distribution. Assembly B7 has a cosine axial power distribution and pattern B radial power distribution. For this radial power distribution, the central rod is a thimble rod (i.e., no power is generated inside this rod) of larger diameter with respect to the fuel rods. The gamma-ray transmission method (chordal averaged) was adopted in the tests performed to measure the density and then converted to the void fraction of the gas-liquid two-phase flow. The declared uncertainties are 4% of void

fraction (ffexp = 4% in Figures 3 and 4) for steady-state bundle tests and 5% of void fraction (ffexp = 5% in Figures 5 and 6) for transient bundle tests. In Table 7 are shown the test series for void fraction measurements.

4.2.1. Steady-State Exercises. For the results at the lower and middle elevations, subcooled nucleate boiling models (condensation and wall heat transfer) play a major role for the prediction of void fraction. At the upper elevation, the fluid is generally saturated, and the drift model becomes important. Moreover, diffusion terms play also a major

Table 9: Comparison between calculated and experimental boiling crisis detection time.

Boiling crisis Series 11 Series 12

detection time (s) EXP CALC EXP CALC

Power increase 106.7 108.7 86.6 72.2

Flow reduction 52.9 48.8 55.0 42.9

Depressurization 88.8 90.8 143.8 140.9

Temperature decrease 140.6 143.3 128.8 121.4

role, and turbulent diffusion promoted by the presence of mixing grid spacers must be taken into account. Standard parameters Kt'0(z) = ktb and M'0(z) = mtb (the turbulent Prandtl number is assumed equal to 1) adopted for bare bundles (without any spacer) have to be increased in the region upstream mixing grid spacers in order to reproduce the enhancement of turbulent mixing. In (15) and (13), K'0 and M'/0 have been therefore adapted for the PSBT assemblies by means of piecewise functions: it was assumed that K(o(z) = kts and M'0(z) = mts downstream a mixing grid spacer, instead of taking the standard values ktb and mtb. Despite the values of ktb and mtb are reliably assessed to be around 0.010 for PWR type bundles [6], the appropriate values to be adopted for kts and mts are not known for these particular mixing grid spacers and thus, to better identify reasonable values for these coefficients, a sensitivity analysis has been performed. A sensitivity analysis for the coefficient Kv0 in the condensation model was also conducted. A sample result is shown in Figure 3, where combined effects of these parameters are reported for run 5.2442. In the subcooled and at the beginning of the saturated region, void fraction is strongly reduced as Kv0 increases, but this parameter does not affect results at the end of the bundle, where the flow is fully saturated. On the other hand, the effect of the turbulent diffusion coefficients kts and mts is significant when the void fraction increases, since it determines the flatness of the void distribution profile in the bundle. Detailed experimental density maps would have been useful to assess a reliable value for the turbulent diffusion coefficients, but these maps are not readily exploitable. As far as we could deduce basing on the void fraction profile in the central subchannel, a value for kts and mts was chosen equal to 0.045. The value of Kv0 was instead taken equal to 1.5 x 10~4, which is in the standard range of values adopted in thermal-hydraulic analysis performed with this model. Figure 4 shows the calculated void fraction against the experimental values for series 5, 6, and 7, with the following parameters: Kv0 = 1.5 x 10~4, kts = mts = 0.045, ktb = mtb = 0.01.

4.2.2. Transient Exercises. For the transient simulations, three different configurations were taken into account as shown in Table 7. Test series 5T considers the same assembly conditions used for steady-state series 5, with a uniform axial power distribution, and pattern A of radial power distribution.

For transient test series 6T the same assembly conditions used for steady-state series 6 were applied, which means a

cosine axial power distribution, and pattern A radial power distribution. Test series 7T considers the same assembly conditions used for steady-state test series 7, which are cosine axial power distribution, and pattern B radial power distribution.

For each transient analysis, four different scenarios were considered: power increase, flow reduction, depres-surization, and temperature increase. Boundary conditions were provided for each scenario, using table where the dependent variables (i.e., power, pressure, mass flux, and inlet temperature resp.,) were specified as function of the simulation time. Same parameters and models as for steady-state bundle computations were used.

Figure 5 shows complete results for transients 5T, also representative of the other transient series 6T and 7T. FLICA-OVAP presents an overall agreement against the experimental void fraction for the transients, although discrepancies with experimental data can be observed, in particular for the temperature increase transient. For the power increase transient, the evolution of the void fraction follows the data, but the void fraction is underestimated at the middle and upper locations. The results concerning the depressurization have more discrepancy as time increases, but are consistent with the experimental evolution.

For the temperature increase transient, a shift of calculated void fractions profiles with respect to experimental values can be noticed. This discrepancy has been noticed by several participants to the first PSBT/OECD workshop in Pisa. Fluid temperature probe is located in a pipe between the preheater and the inlet nozzle of the test section [12]. It is therefore reasonable to guess that a delay of several seconds occurs between the measurement point and the inlet of the heated section. This point has been investigated by JNES [12] but the distance from the measurement location to the inlet of the heated section (where the temperature boundary condition has to be applied in thermal-hydraulics codes) was not clarified. As far as we can say basing on the result obtained by FLICA-OVAP, if a delay of 6 s is assumed for the inlet temperature with respect to other boundary conditions (i.e., the initial temperature is maintained during 6 s and then the temperature variation is applied), the agreement between calculated and experimental void fraction profiles is improved. It can clearly be seen in Figure 6, where these profiles are compared for the calculations with and without the time delay of 6 seconds.

5. DNB Results

The two models adopted in this work, the Shah and the Katto and Ohno models, have been developed and tested mostly in single-channel uniform heat flux conditions.

In this work, the applicability of them to bundle geometries with nonuniform heat flux profile has been investigated. Each experimental series of the NUPEC database has been calculated by the two models, adopting a twofold approach.

In a first step, homogenized one-dimensional calculations have been performed, adopting a bundle-averaged description.

In a second step, three-dimensional calculations have been performed, adopting a subchannel description. In this case, the two boiling crisis models are applied at the subchannels scale. The fundamental difference with the first approach consists in having different hydraulic diameters depending on whether the selected subchannel is a central, a side or a corner subchannel. However, some special approach has to be taken for side and corner subchannels, since CHF models tend to underestimate their critical heat flux. Moreover, as observed by Tong and Tang [13, Section 5.4.5.2, Figure 5.54], for a given inlet enthalpy, critical heat fluxes in the presence of cold walls are even higher than for internal subchannels. For the sake of simplicity, the possibility to have boiling crisis in side and corner subchannels was thus excluded.

In the following paragraphs, for each series, results obtained via the bundle-averaged and the subchannel approach will be shown for the two models. The performances of the different models are also summarized in Table 8.

5.1. Steady-State Series

5.1.1. Series 0, Series 2, and Series 3. Series 0, Series 2 investigate boiling crisis in rod bundle (rods array: 5 x 5) and uniform heat flux. The two series differ only for the position of mixing and spacer grids. In Series 3, the bundle is also modified, adopting an array of 36 rods (6 x 6). The axial power shape is uniform for all series, but peripheral rods have a slightly lower power with respect to central rods: the radial power distribution factor is 0.85 for peripheral rods and 1.0 for central rods.

In Figures 7, 8, and 9, calculated and experimental critical heat fluxes are compared, respectively, for Series 1, Series 2, and Series 3. It is ascertained that the Shah model tends to predict higher critical heat flux values than the Katto and Ohno model, for both the bundle-averaged and the subchannel approach. Moreover, in this particular configuration with uniform heat flux and a moderate radial offset, 1D or 3D calculations give very similar results.

5.1.2. Series 4 and Series 13. Series 4 and Series 13 introduce the effect of nonuniform axial profile. A cosine shape was adopted on the same bundle geometry used for Series 2. To deal with nonuniform heat flux, a correction factor was adopted, as proposed by Tong and Tang [13]. In Figures 10 and 11 calculated and experimental critical heat fluxes are compared, respectively, for Series 4, Series 13. As for previous tests, the Shah model predicts higher critical heat flux values than the Katto and Ohno model, providing a better agreement with experimental values. Minor differences are experienced between the homogenized and the subchannel description.

5.1.3. Series 8. Series 8 includes the presence of a central thimble rod, representing guide tube for control rods. Con-trarily to the previous series, the results of the homogenized and the subchannel approaches are different (see Figure 12). As expected, in the presence of a radial heterogeneous power shape, a better description is achieved with the subchannel

approach. In general the Shah model gives a better estimate of critical heat fluxes.

5.2. Transient Series. Two experimental campaigns have been consecrated to detect boiling crisis in transient conditions. Different transient conditions have been investigated by NUPEC, including power increase, inlet temperature decrease, flow reduction, and depressurization. Two series of data have been proposed, Series 11 and Series 12, involving different bundle configurations. Series 11 is based on the bundle configuration adopted for Series 4 (cosine axial power shape, rods array 5 x 5), whereas Series 12 involves a central thimble rod. In the following tables, results obtained by the Shah model with the subchannel approach are listed. Boiling crisis is relatively well predicted for Series 11, whereas a larger discrepancy is ascertained for Series 12 (see Table 9).

6. Conclusions and Perspectives

This paper has presented the current capabilities of the FLICA-OVAP code in predicting void distribution and boiling crisis phenomena. The NUPEC database released in the frame of the OECD/PSBT international benchmark has been addressed.

Void fraction measurements in the subchannel configurations are of major interest for the validation of mass transfer model and the OSV criterion. In this analysis, the attention was focused on the mass transfer model, but further improvement of this work could include the analysis of other OSV criteria. Results obtained by FLICA-OVAP with a set of standard coefficients for the different models show a good agreement of the calculated densities and a slight underestimation of the void fraction at the measurement location, mainly located in the saturated regime for the considered runs.

For the steady-state bundle tests, the Kv0 coefficient and diffusion coefficients kt and mt are the key parameters to fit the void fraction. In particular, mixing grid spacers play a major role, since they enhance the turbulence in the downstream flow. It was found that a reasonably good agreement between calculated and experimental void fraction profiles is achieved when the turbulent diffusion associated to mixing grid is taken equal to kts = mts = 0.045. Adopting the same set of parameters, a reasonably good agreement between calculated and experimental void fraction was also ascertained for transient tests, even if a systematic underestimation of void fraction at the middle and upper measurement locations has been found.

Discrepancies noticed on temperature increase transients were corrected by applying a delay of 6 s of the original inlet temperature boundary conditions, in order to simulate the residence time of the fluid between the temperature probe where the boundary conditions are given and the inlet of the test section.

Results obtained from the subcooled region suggest possible directions to be pursued in order to improve the current modeling: further developments and validation will involve the OSV criterion and the modeling of heat flux and heat flux partitioning in subcooled nucleate boiling, but

also the modeling of inter-phase mass transfer and turbulent diffusion terms.

To predict boiling crisis, two different models have been tested: the Shah and the Katto and Ohno model. Both a bundle-averaged and subchannel approaches have been investigated. The Shah model, as implemented in the FLICA-OVAP code allows achieving a better prediction of critical heat fluxes in the different investigated configurations. Nevertheless, when the details of the radial power shape become important, better agreement was obtained with the subchannel approach permitting a local description of bundle thermal-hydraulics. Encouraging results were also obtained from the analysis of transient tests aimed at predicting the onset of boiling crisis. Future activities are also planned for future developments of the code, including the implementation of Groeneveld's look-up table [14] and other models based on local approaches.

Journal of Heat and Mass Transfer, vol. 27, no. 9, pp. 16411648, 1984.

[12] H. Utsuno, "OECD/NRC benchmark based on NUPEC PWR subchannel and bundle tests (PSBT). Assembly specification and Benchmark Database. Volume I Addendum," Tech. Rep JNES/SAE-TH08-0019, JNES, 2010, (Addendum).

[13] L. S. Tong and Y. S. Tang, Boiling Heat Transfer and Two-Phase Flow, Taylor & Francis, 2nd edition, 1997.

[14] D. C. Groeneveld, J. Q. Shan, A. Z. Vasic et al., "The 2006 CHF look-up table," Nuclear Engineering and Design, vol. 237, no. 15-17, pp. 1909-1922, 2007.

References

[1] A. Rubin, A. Schoedel, M. Avramova, H. Utsuno, H. Bajorek, and S. Velasquez-Lozada, "OECD/NRC benchmark based on NUPEC PWR subchannel and bundle tests (PSBT). Volume I: Experimental Database and Final Problem Specifications," Tech. Rep NEA/NSC/D0C(2010)1, NEA/NSC, November 2010.

[2] P. Fillion, A. Chanoine, S. Dellacherie, and A. Kumbaro, "FLICA-OVAP: a new platform for core thermal-hydraulic studies," Nuclear Engineering and Design, vol. 241, no. 11, pp. 4348-4358,2011.

[3] M. Bucci, A. Charmeau, and P. Fillion, "FLICA-OVAP: elements of validation for LWRs thermal-hydraulics studies," in Proceedings of the International Congress on Advances in Nuclear Power Plants (ICAPP '11), Nice, France, May 2011.

[4] B. Chexal, G. Lellouche, J. Horowitz, and J. Healzer, "A void fraction correlation for generalized applications," Progress in Nuclear Energy, vol. 27, no. 4, pp. 255-295, 1992.

[5] M. Ishii, "One dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow," Tech. Rep ANL-77-27, ANL, 1977.

[6] E. Royer, S. Aniel, A. Bergeron et al., "FLICA4: Status of numerical and physical models and overview of applications," in Proceedings of the 11th International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-11), Avignon, France, October 2005.

[7] D. Chisholm, "Friction during the flow of two-phase mixtures in smooth tubes and channels," Tech. Rep NEL-529, National Engineering Laboratory, Department of Trade and Industry, 1972.

[8] W. Jens and P. Lottes, "Analysis of heat transfer burnout, pressure drop and density data for high pressure water," Tech. Rep ANL-4627, Argonne National Laboratory (ANL), 1951.

[9] N. E. Todreas and M. S. Kazimi, Nuclear Systems I, Taylor & Francis, 1993.

[10] M. M. Shah, "Improved general correlation for critical heat flux during upflow in uniformly heated vertical tubes," International Journal of Heat and Fluid Flow, vol. 8, no. 4, pp. 326-335, 1987.

[11] Y. Katto andH. Ohno, "An improved version of the generalized correlation of critical heat flux for the forced convective boiling in uniformly heated vertical tubes," International

Copyright of Science & Technology of Nuclear Installations is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.