Science and Technology of Nuclear Installations

Volume 2012 (2012), Article ID 746467, 20 pages

doi:10.1155/2012/746467

Research Article

Analysis of Void Fraction Distribution and Departure from Nucleate Boiling in Single Subchannel and Bundle Geometries Using Subchannel, System, and Computational Fluid Dynamics Codes

Taewan Kim,1,2 Victor Petrov,³ Annalisa Manera,³ and Simon Lo⁴

¹Nuclear Energy and Safety Research Department, Paul Scherrer Institute, 5232 Villigen, Switzerland

²Department of NPP Engineering, KEPCO International Nuclear Graduate School, Ulsan 689-882, Republic of Korea

³Department of Nuclear Engineering and Radiological Sciences, University of Michigan, Ann Arbor, MI 48109, USA

⁴CD-adapco, Trident House, Basil Hill Road, Didcot OX11 7HJ, UK

Received 30 April 2012; Accepted 20 July 2012

Academic Editor: David Novog

Copyright © 2012 Taewan Kim et al. 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.

Abstract.

In order to assess the accuracy and validity of subchannel, system, and computational fluid dynamics codes, the Paul Scherrer Institut has participated in the OECD/NRC PSBT benchmark with the thermal-hydraulic system code TRACE5.0 developed by US NRC, the subchannel code FLICA4 developed by CEA, and the computational fluid dynamic code STAR-CD developed by CD-adapco. The PSBT benchmark consists of a series of void distribution exercises and departure from nucleate boiling exercises. The results reveal that the prediction by the subchannel code FLICA4 agrees with the experimental data reasonably well in both steady-state and transient conditions. The analyses of single-subchannel experiments by means of the computational fluid dynamic code STAR-CD with the CD-adapco boiling model indicate that the prediction of the void fraction has no significant discrepancy from the experiments. The analyses with TRACE point out the necessity to perform additional assessment of the subcooled boiling model and bulk condensation model of TRACE.

1. Introduction

An international benchmark program, namely, the OECD/NRC NUPEC BWR Full-Size Bundle Test (BFBT) Benchmark [1], was organized by OECD/NEA and US NRC to encourage advancement in subchannel analysis methods. About 30 organizations from 15 countries had participated in the benchmark and an extensive validation of subchannel, computational fluid dynamic (CFD), and thermal hydraulic system codes had been carried out. Building on the success of the BFBT benchmark, OECD/NEA and US NRC have organized an additional international benchmark program for pressurized water reactor (PWR) conditions [2]. The benchmark, namely the OECD/NRC NUPEC PWR subchannel and bundle tests (PSBT) Benchmark, aims at assessing the capabilities of system codes, subchannel codes, and CFD codes to predict detailed void distributions and departure from nucleate boiling (DNB) in subchannels on the basis of experimental data measured at a full-scale prototypical PWR rod bundle in NUPEC test facility.

Meanwhile, the Paul Scherrer Institute (PSI) is engaged in strengthening its subchannel analysis capability for light water reactors in Switzerland. Following the experience from the European Union NURESIM and NURISP projects, PSI has opted for the subchannel analysis code FLICA4 developed at CEA in France [3]. Besides, the application of CFD for subchannel analysis is also considered as the new frontier of the state-of-the-art methodologies. In this context, PSI decided to participate in the PSBT benchmark with FLICA4 and a commercial CFD code, namely, STAR-CD [4]. In addition, analyses using the system code TRACE (TRAC/RELAP Advanced Computational Engine) [5] of the NUPEC experiments have been carried out in parallel to evaluate the applicability of TRACE for subchannel analysis.

This paper reviews the different analyses of the PSBT Benchmark that were conducted at PSI, and discusses and compares the results obtained with the different thermal hydraulic codes.

2. Benchmark Description

The OECD/NRC PSBT benchmark aims at encouraging advancement in subchannel analysis of fluid flow in rod bundles under the conditions typical for PWRs. The benchmark is aimed at assessing the capabilities of system-codes, subchannel codes, and CFD codes to predict void distributions, including DNB, in PWR rod bundle geometry on the basis of experimental data measured at the NUPEC test facility [2]. The NUPEC test facility depicted in Figure 1 consists of a high pressure and high temperature recirculation loop, a cooling loop, and instrumentation and data acquisition systems. The recirculation loop consists of a test section, circulation pump, preheater, steam drum (acting as a pressurizer), and a water mixer. Different test sections were constructed to represent a single subchannel and a complete rod bundle, respectively. The design pressure is 19.2 MPa and the design temperature is 362°C.

Figure 1: NUPEC test facility [2].

The benchmark consists of two phases: phase I for void distribution benchmark and phase II for DNB benchmark. Benchmark phase I includes four exercises: steady-state single subchannel exercise, steady-state and transient bundle exercises, and pressure drop exercise. Phase II consists of three exercises, which are steady-state fluid temperature exercise and steady-state and transient DNB exercises for bundle geometries.

3. Thermal Hydraulic Codes

Three different codes have been employed for the PSBT benchmark: the subchannel analysis code FLICA4 v1.10.8, the system code TRACE v5.0 Patch 2, and the CFD code STAR-CD v4.14.

3.1. FLICA4

FLICA4 is a three-dimensional (3D) two-phase flow analysis code developed for subchannel analysis by CEA in France. The two-phase flow model in FLICA4 is based on a 4-equation model, combined with a drift flux model to describe the relative velocity between phases. The drift-flux model developed by Chexal et al. [6] has been selected for the benchmark, based on a sensitivity analysis with cases randomly selected out of the benchmark database [7]. The Dittus-Boelter [8] and Jens-Lottes [9] correlations were employed to model the single-phase and two-phase heat transfer, respectively. The F3 correlation was employed to model the mass transfer between phases and proportional coefficient, KV0, in the correlation was decided by a sensitivity analysis [7] (note that the F3 model/correlation refers to the model/correlation used in FLICA3-M [10]). The friction coefficient in FLICA4 is calculated as a product of the single-phase friction factor, two-phase multiplier, and heated wall corrector. The single-phase friction factor was calculated by using the default F3 model, which employs the Blasius model [11], while the Friedel model [12] was used as the two-phase multiplier. The F3 model, which requires four user-defined coefficients, was employed for the heated wall corrector. The turbulence model based on a mixing length approach was applied in the calculation. This model includes four empirical coefficients in the correlations for turbulent viscosity and conductivity. The required coefficients were derived from a sensitivity analysis for several cases selected randomly from the benchmark database.

3.2. TRACE

TRACE is the latest best-estimate thermal-hydraulic system code developed by US NRC for analyzing steady-state and transient neutronic/thermal-hydraulic behavior of light water reactors. In TRACE, a two-fluid 6-equation model of the steam-water flow is employed. The Gnielinski correlation [13] is used as a single-phase heat transfer coefficient and two-phase multipliers by Aggour et al. [14] and Rezkallah and Sims [15] are employed for bubbly/slug regime. The subcooled boiling heat transfer is modeled by means of the Lahey’s mechanistic model [16]. As for wall friction, the Churchill formula [17] and a void-fraction-based two-phase multiplier [5] are employed for single-phase and two-phase conditions, respectively. The interfacial heat transfer in bubbly flow regime is modeled by means of the Ranz-Marshall correlation [18], with the interfacial area concentration by Ishii and Mishima [19]. For the cap bubble/slug flow regime, the Ranz-Marshall correlation with consideration of small bubbles is implemented. The interfacial heat transfer under subcooled boiling condition is modeled by using the Lahey and Moody model [20].

3.3. STAR-CD

STAR-CD is a CFD code developed by CD-adapco. The governing transport equations solved are the conservation laws of mass, momentum, and energy for each of the two phases. These equations are solved in 3 dimensions. The numerical algorithm used for this benchmark is IPSA (Interphase Slip Algorithm) [21–23]. It is fully implicit, using the pressure-correction-based method extended to multiphase flows. No flow regime map is used. Bubbly flow is assumed anywhere where the gaseous phase is present. Interfacial momentum transfer includes models for drag force, turbulent drag force, virtual mass, force and momentum transfer due to mass transfer. Interfacial mass transfer is due to evaporation or condensation as computed by the boiling model. Interfacial energy transfer is due to temperature difference between the two phases and the saturation temperature. Wall drag is calculated by standard wall function used in CFD codes. Wall boiling heat transfer is calculated by the “wall heat partitioning” model of Kurul and Podowski [24]. The lift force acting on the bubbles in radial direction was neglected. All the calculations were performed applying the liquid properties as a function of temperature and pressure.

Total wall heat flux is made up of three components, as follows:

where : convective heating, : quenching, and : evaporation.

For the nucleation site density, the Hibiki and Ishii model [25] with the following application range was used:

(i)pressure: 0.101–19.8 MPa;(ii)mass flux: 0–886 kg/m²sec;(iii)contact angle: 5°–90°;(iv)Number density: 1.0E+04 –1.51E+10 sites/m².

The bubble departure diameters were calculated using Kocamustafaogullari’s correlation [26], which is based on water experimental data at pressures from 0.0067 to 14.187 MPa. The computed bubble diameters in the flow were obtained from the Kurul and Podowski formulation using linear interpolation between two bubble diameters at a specified liquid subcooling [24].

The numerical algorithm and its parameters (i.e., relaxation factors, difference scheme, etc.) were kept the same for all geometries of the benchmark exercise and were not adjusted from run to run.

4. PSBT Benchmark Phase I: Void Distribution Benchmark

4.1. Single Subchannel Experiments

4.1.1. Description of Experiments

A series of steady-state experiments were performed to measure the void distribution in a single subchannel. The test section consists of four different geometries, which simulate various subchannels in a bundle as described in Table 1. The length of the heated section is 1.555 m and the void fraction was measured at 1.4 m from the bottom of the heated section. The cross-sectional view and the geometrical parameters of each subchannel test assembly are depicted in Figure 2 and listed in Table 2, respectively. The ranges of boundary conditions were as follows:

(i)pressure: 4.90–16.6 (MPa),(ii)mass Flux:  –  (kg/m²-sec),(iii)power: 12.5–90.0 (kW),(iv)inlet temperature: 164.1–345.0 (°C).

Table 1: Geometry and power shape for single-channel tests [2].

Item

Data

Assembly (subjected subchannel)

Subchannel type
Center (typical)
Center (thimble)
side
corner

Number of heaters
4 × 1/4
3 × 1/4
2 × 1/4
1 × 1/4

Axial heated length (mm)
1555
1555
1555
1555

Axial power shape
Uniform
Uniform
Uniform
Uniform

■: subjected subchannel, white circle: heated rod, gray circle: thimble rod.

Table 2: Geometrical parameters of each subchannel [2].

Subchannel type

Typical (S1)
Thimble (S2)
Side (S3)
Corner (S4)

Flow area, mm²
107.098
107.098
68.464
42.592

Heated perimeter, mm
29.845
22.384
14.923
7.461

Wetted perimeter, mm
54.645
54.645
44.923
33.161

(a)

(b)

(c)

(d)

Figure 2: Cross-sectional view of subchannel test assembly [2].

4.1.2. Thermal Hydraulic Model

FLICA4 and TRACE

Identical nodalizations were used for the single-channel analyses carried out with FLICA4 and TRACE. The single-channel test section was nodalized by means of a one-dimensional pipe divided into 32 axial nodes. Four different geometries were prepared, to take into account the difference in cross-section and hydraulic diameter of each case. The Chexal-Lellouche drift-flux model was employed for the FLICA4 calculations. Multipliers for the turbulent diffusivity and viscosity, and , were set to 0.01 based on a sensitivity analysis. A value of was assigned to the recondensation coefficient, KV0. As for the TRACE model, the single channel was modeled with a PIPE (one-dimensional) component. Inlet mass flux and temperature and outlet pressure were imposed as boundary conditions.

STAR-CD

In order to reduce the needs in computational resources, all subchannel geometries were modeled utilizing subchannel symmetry. 1/8 symmetry could be employed for the central subchannel, while the half symmetry could be exploited for the central (thimble) side and corner subchannels. Heated rods were not modeled explicitly, instead heat fluxes were applied as boundary conditions on the channel side walls. The height of the computational model is the same as the heated length of the experimental test section (1555 mm). A flat velocity profile was used for the inlet boundary condition, since no additional data about the inlet manifold geometry was made available to the benchmark participants. A summary of the boundary conditions applied to the CFD model is shown in Figure 3 (uniform velocity profile at the inlet; pressure outlet; constant heat flux at heated wall; adiabatic wall for unheated section; symmetry planes). Boundary values were specified according to the benchmark specifications [2].

Computational cells distribution and mesh size were kept the same for each run in the test series. The CFD meshes used are shown in Figure 4:

(i)53600 hexahedral cells for S1;(ii)214400 hexahedral cells for S2;(iii)114800 hexahedral cells for S3;(iv)62200 hexahedral cells for S4.

Figure 3: CFD model.

(a)

(b)

(c)

(d)

Figure 4: Hexahedral computational mesh.

The number of axial cell layers was set to 100 for each subchannel geometry. Mesh sensitivity studies performed on the 1/8 symmetry sector model of the type S1 geometry showed that there is no significant effect on the section averaged void fraction value after reaching 250 cells per axial cell layer.

4.1.3. Results

FLICA4 and TRACE

The results of the void fraction for S1 to S4 are depicted in Figure 5. Generally, both FLICA4 and TRACE could predict the experimental results well. However, both codes underpredicted the void fraction, taking into account a measurement accuracy of 4.0% (absolute void fraction), especially for the corner subchannel type (S4).

The bias of the calculation results from FLICA4 and TRACE has been assessed by using the linear regression method. Figure 6 shows the results from the linear regression for the void fraction calculated by FLICA4 and TRACE. Each result is fitted by a linear function with a high correlation coefficient (adjusted ) of 0.95. The linear regression analyses indicate that the void fraction results from FLICA4 and TRACE are biased from the experiments by 2.89% and 3.85% on average, respectively. The mean absolute errors of the FLICA4 and TRACE calculations are found to be 4.46% and 4.95%, respectively, which are slightly higher than the measurement error, assessed at 4.0% absolute void fraction. Both codes therefore reproduce the experimental data reasonably well.

(a)

(b)

(c)

(d)

Figure 5: Void fraction predictions for single-channel exercise.

(a)

(b)

Figure 6: Results from linear regression analyses.

STAR-CD

Figure 7 shows the comparison between experimental and calculated cross-section averaged void fraction for different subchannel geometries ((a) S1, (b) S2, (c) S3, and (d) S4).

The red lines indicate error bars of ±4%. Figure 8 shows the radial void fraction distribution in the measuring cross-section for selected cases; quantitative comparisons were not performed since only graphical data were available for the radial void fraction distributions. Figure 9 shows the computed axial void distribution profiles for selected cases. In general, for all geometries no significant discrepancy was found between computed and experimental cross-section averaged void fractions. However, a slight overprediction of the void fraction has been observed especially for cases with low void fractions. Further verification of the boiling model aimed at assessing the prediction of radial void fraction distributions is required; however, suitable experimental data was unfortunately not available within the PSBT benchmark.

(a) Case S1

(b) Case S2

(c) Case S3

(d) Case S4

Figure 7: Section average void fraction results.

(a)

(b)

(c)

(d)

Figure 8: Void fraction distribution in a measuring section (CFD versus experimental data).

Figure 9: Axial void profile for 4 runs with S1 type of geometry.

4.2. Steady-State Bundle Experiments

4.2.1. Description of Experiments

A series of steady-state experiments were performed to measure the void distribution in bundle geometry. The void fraction from the experiment is the one averaged over the four central subchannels. The three different 5 × 5 bundles as given in Table 3 were employed for the experiments. Bundles B5 and B6 have the same geometry but different axial power profiles: uniform and cosine shapes. Bundle B7 has a thimble rod at the center of the test section and cosine axial power shape. Two radial power distributions as depicted in Figure 10 were used for the experiments: type A and B, respectively. The two radial profiles are basically identical except for the non-heated thimble rod in type B. The cosine axial power profile applied in this experiment is given in Table 4. The bundles were equipped with three different grid spacers: spacer with mixing vane (MV), spacer without mixing vane (NMV), and simple spacer (SS), respectively. The number and location of each spacer are also listed in Table 3 and the bundle average pressure loss coefficients for the three types of spacers are provided in Table 5.

Table 3: Geometry and power shape for each bundle [2].

Item

Data

Assembly

Rods array
5 × 5
5 × 5
5 × 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 rods 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

Radial power shape
A
A
B

Axial power shape
Uniform
Cosine
Cosine

Number of MV spacers
7
7
7

Number of NMV spacers
2
2
2

Number of simple spacers
8
8
8

MV spacer location (mm)
471, 925, 1378, 1832, 2285, 2739, 3247

NMV spacer location (mm)

2.5, 3755

Simple spacer location (mm)
237, 698, 1151, 1605, 2059, 2512, 2993, 3501

White circle: heated rod, gray circle: thimble rod.

MV: mixing vane. NMV: no mixing vane.

Spacer location is distance from bottom of heated length to spacer bottom face.

Table 4: Cosine power profile [2].

Node
Relative Power

Cosine

(Bottom)

1
0.42

2
0.47

3
0.56

4
0.67

5
0.80

6
0.94

7
1.08

8
1.22

9
1.34

10
1.44

11
1.51

12
1.55

13
1.55

14
1.51

15
1.44

16
1.34

17
1.22

18
1.08

19
0.94

20
0.80

21
0.67

22
0.56

23
0.47

24
0.42

(Top)

Table 5: Bundle average spacer pressure loss coefficients [2].

Spacer type
Loss Coefficient

Simple spacer (SS)
0.4

Nonmixing vanes spacer (NMV)
0.7

Mixing vanes spacer (MV)
1.0

(a) Type A

(b) Type B

Figure 10: Radial power distribution [2].

4.2.2. Thermal Hydraulic Model

Analyses of the steady-state bundle experiments have been conducted by means of FLICA4 and TRACE only.

FLICA4

Since the geometry and the radial power distribution adopted for the bundle tests are fully symmetric, the symmetric boundary condition was employed so that each test assembly was described by means of a 1/8 symmetrical model as depicted in Figure 11. The symmetrical model consists of six subchannels subdivided in 100 axial nodes. A sensitivity study on the axial nodalization has shown negligible difference between 25 and 100 axial nodes. However, for better resolution, 100 axial nodes were employed in all the calculations.

The Chexal-Lellouche drift-flux model was employed and a value of was adopted for the recondensation coefficient, KV0, based on a sensitivity results with cases randomly selected out of the benchmark database. For all the bundle cases, except for assembly B7, coefficients for turbulent diffusivity and viscosity, and , were set to 0.01 based on a sensitivity analysis. For assembly B7 that has a guide thimble at the center, and were set to 0.05 also based on a sensitivity analysis. The pressure drop due to the spacers was taken into account by means of the singular pressure drop model in FLICA4.

Figure 11: 1/8 symmetrical model for FLICA4.

TRACE

The bundles are described by using a 3D rectangular VESSEL component of TRACE as depicted in Figure 12. This includes 100 heat structures, each representing a quarter of a single heater rod. The VESSEL component was nodalized with 25 axial levels and the K-factors were prescribed at the relevant elevation in order to model the pressure drop introduced by spacers. No local K-factor was used on the crossflow connections. Since the 3D VESSEL component cannot be connected directly to the FILL component (which provides the inlet boundary conditions) and the BREAK component (which provides the outlet boundary conditions), PIPE components were introduced at the inlet and outlet of the VESSEL components.

(a) Overview of TRACE Model

(b) Top-down view of test section

Figure 12: Nodalization for TRACE.

4.2.3. Results

The void fraction was measured at three different elevations: 2216 mm (lower), 2669 mm (middle), and 3177 mm (top) from the bottom of the test section. The calculated void fractions at the four central subchannels were averaged and compared with the experimental data. The results of both codes presented in Figure 13 indicate that, in general, both codes overpredicted the void fractions with respect to the experiments. In particular, TRACE overpredicted the void fraction in the lower and middle regions considerably. The linear regression method was applied to estimate the bias of each calculation quantitatively. As shown in Figure 14, the results from FLICA4 and TRACE are biased by +2.10% and +6.78% on average, respectively. The accuracies of the calculation represented by the mean absolute error are 4.09% and 8.71% for FLICA4 and TRACE, respectively. Considering the measurement error of 4.0% void fraction, the difference between the experiments and the TRACE results is relatively high. As confirmed by the mean absolute error of each region listed in Table 6, the discrepancy between TRACE and experimental data is mainly encountered for the lower and middle regions, where smaller void fraction is expected. TRACE tends to overpredict the void fraction especially for values below 40%. In Figure 15, the axial void fraction profile of a case from test series B5 is presented together with experimental data. The void fraction predicted by TRACE starts to rapidly increase at an elevation of about 1.8 m. This result points out the necessity to further investigate and validate the subcooled boiling model and bulk condensation model in TRACE. FLICA instead reproduces the experimental data reasonably well.

Table 6: Mean absolute error of void fraction.

Code
Lower
Middle
Upper
Overall

FLICA4

    Averaged error (%)
1.18
3.54
1.59
2.10

    Mean absolute error (%)
2.95
5.29
4.04
4.09

TRACE

    Averaged error (%)
9.54
7.45
3.35
6.78

    Mean absolute error (%)
10.1
9.49
6.49
8.71

(a)

(b)

(c)

(d)

Figure 13: Void fraction of bundle exercise.

(a)

(b)

Figure 14: Results from linear regression analysis.

Figure 15: Axial void fraction profile of case B5.2442.

4.3. Transient Bundle Experiments

4.3.1. Description of Experiments and Thermal Hydraulic Model

The transient experiments were conducted with assemblies B5, B6, and B7. Four different transient experiments were performed for each assembly: a power increase, a flow reduction, a depressurization, and a temperature increase. The initial conditions for each transient test are summarized in Table 7.

Table 7: Initial conditions for transient bundle tests [2].

Initial conditions

Test series
Assembly
Pressure (kg/cm² a)
Mass flux (10⁶ kg/m² h)
power (kW)
Inlet temperature (Celsius)
Transients

5T
B5
154.2
11.95
2282.0
300.4
Power increase

153.8
11.93
2244.0
301.2
Flow reduction

153.0
11.92
2236.0
300.4
Depressurization

152.5
11.94
2230.0
301.7
Temperature increase

6T
B6
158.2
11.55
2621.0
288.1
Power increase

158.4
12.03
2574.0
288.8
Flow reduction

154.6
12.02
2556.0
288.2
Depressurization

157.2
11.92
2603.0
288.8
Temperature increase

7T
B7
158.2
12.02
2500.0
291.9
Power increase

158.1
12.04
2405.0
292.0
Flow reduction

155.0
11.99
2577.0
291.8
Depressurization

158.8
11.99
2496.0
290.2
Temperature increase

The transient experiments were analyzed by means of FLICA4 only. The FLICA4 models for the transient experiments are the same as the ones used for the steady-state analyses. As for the steady-state cases, since the void fraction is averaged over the four central subchannels in the experiment, the void fraction of the subchannel 6 in the FLICA4 model can be considered as the averaged void fraction.

4.3.2. Results

The analysis has been carried out for all transient cases specified as part of the benchmark. As an example, the calculation results from transient exercise 5T are reported in Figure 16 together with the experimental data. In general, FLICA4 is found to reproduce the experimental data reasonably well even in fast transients such as power increase and flow reduction. However, a large discrepancy was observed for the temperature increase case. It is presumed that the discrepancy is originated from a delay in the experimental inlet temperature variation. In fact, since the inlet temperature measurement was located between the preheater and the inlet nozzle of the test vessel, the actual temperature increase at the inlet of the test section would happen with some delays with respect to the temperature variation at the measuring point. Because the distance between measuring point and the test section inlet is unknown, it is not possible to estimate the delay in inlet temperature variation appropriately. Nevertheless, when shifting the FLICA4 results for test case 5T by 6.0 sec with respect to the original simulation time (see Figure 17), very good agreement with the experimental data is obtained. The time shift by 6.2 sec due to the mass flow rate difference from case 5T produces consistent results for test cases 6T and 7T. Therefore, it can be reasonably concluded that the transient exercise was well predicted by FLICA4.

(a)

(b)

(c)

(d)

Figure 16: Axial void fraction of transient exercise 5T.

(a) 5T—temperature increase

(b) 6T—temperature increase

(c) 7T—temperature increase

Figure 17: Axial void fraction with original and modified time.

5. PSBT Benchmark Phase II: DNB Benchmark

5.1. Description of Experiments

PSBT benchmark phase II aims at developing and assessing mechanistic models for DNB prediction. Both steady-state and transient DNB exercises are included in phase II and all the tests were carried out with bundles. In the experiments, the occurrence of DNB was detected by a sudden increase of the surface temperature measured by thermocouples attached to the heater rods. The heating power was increased in fine steps to the vicinity of DNB power estimated by preliminary analysis and experience. In the experiment, a sudden surface temperature increase of more than 11°C confirmed the occurrence of DNB and the corresponding DNB power was defined as the power just before the sudden temperature increase.

As listed in Table 8, eight different assemblies are employed for phase II. All assemblies consist of a 5 × 5 subset of a typical 17 × 17 fuel assembly type, except for assembly A3 which is 6 × 6. As for the assemblies A8 and A12, a thimble rod is located in the center of each assembly. Three different types of spacer are included in the test assembly and the loss coefficient for each spacer is the same as the one used in phase I of the benchmark. Two axial power profiles are considered in this benchmark: uniform and cosine shapes, respectively. The cosine axial power profile used in phase II is the same as the one used in phase I. In total four radial power distributions are employed in this exercise, as depicted in Figure 18.

Table 8: Test assemblies for phase II [2].

Assembly
Reference fuel type
Rods

array
Type of cell
Power distribution

Radial
Axial

A0
17 × 17 M
5 × 5
Typical cell
A
Uniform

A1
Typical cell
C
Uniform

A2
Typical cell
A
Uniform

A3
6 × 6
Typical cell
D
Uniform

A4
5 × 5
Typical cell
A
Cosine

A8
Thimble cell
B
Cosine

A11
Typical cell
A
Cosine

A12
Thimble cell
B
Cosine

(a) Type A

(b) Type B

(c) Type C

(d) Type D

Figure 18: Radial power distribution [2].

5.2. Thermal Hydraulic Model

FLICA4 is employed for the analyses of PSBT benchmark phase II. Thermal hydraulic models for the FLICA4 calculations are generated on the basis of information on the geometry and the power distribution. Unlike in phase I, it was not always possible to adopt 1/8 symmetry in the models for phase II due to the radial power distribution C, which allows implementing 1/2 symmetry only. Therefore, for consistency, all models have been generated by using a 1/2 symmetry, as depicted in Figure 19. The models were nodalized axially with 100 nodes.

(a) 5×5 assembly

(b) 6×6 assembly

Figure 19: 1/2 symmetry model.

The Chexal-Lellouche drift-flux model was employed and a value of was imposed for the recondensation coefficient, KV0. The multipliers for turbulent diffusivity and viscosity, and , were set to 0.01 based on the results from phase I. The pressure drop by the spacers was considered by means of the singular pressure drop model in FLICA4. In the FLICA4 analysis, the parameter used as indicator for the occurrence of DNB was the minimum DNB ratio (MDNBR), defined as the ratio of the critical heat flux (CHF) predicted by a given correlation to the heat flux at each axial node of each heater rod. DNB is confirmed when the MDNBR is less than unity. For the CHF prediction, the W3 correlation [27] and the Groeneveld look-up table [28] with a correction factor for the diameter as given in (2) were employed,

However, due to limitation in the application ranges, the W3 correlation could not be used for all the cases in phase II. Sensitivity studies carried out with both correlations, which will be discussed in more detail in Section 5.5, indicate that the DNB power predicted by the Groeneveld look-up table is slightly lower than the one from the W3 correlation. However, no significant discrepancy was observed. Considering robustness and conservatism in the analysis, it was decided to retain the Groeneveld look-up table for the phase II DNB analysis.

5.3. Results for Steady-State DNB Tests

Selected cases from the test series highlighted in Table 9 are analyzed for benchmark purpose. The results from all the calculations are plotted in Figure 20. In general, FLICA4 predicts lower DNB powers with respect to experimental data, that is, the results are conservative. The same trend was observed from most of the other participating organizations [29]. The mean absolute error of the DNB power prediction by FLICA4 is 10.1%. Consequently, a comparison with the results obtained by other participants indicates that FLICA4 predicts the DNB power slightly more conservatively but not less accurately than other state-of-the-art subchannel analysis codes.

Table 9: Test series for exercise 2 of PSBT benchmark phase II [2].

Test series
Test

section
Assembly
Test mode
Measurement

Steady

state
Transient
DNB
Fluid

temperature

0

A0
Y

Y

1
5 × 5
A1
Y

Y

2

A2
Y

Y

3
6 × 6
A3
Y

Y

4

A4
Y

Y

8

A8
Y

Y

11T
5 × 5
A11

Y
Y

12T

A12

Y
Y

13

A4
Y

Y

Figure 20: Results of steady-state DNB exercises.

5.4. Results for Transient DNB Tests

The transient DNB benchmark was conducted for test series 11T and 12T, which include four transient scenarios in each test as indicated in Table 10, a power increase, a flow reduction, a depressurization, and a temperature increase, respectively. The parameter of interest in this benchmark exercise is the time of DNB occurrence.

Table 10: Test series for transient DNB exercise [2].

Test

series
Assembly
Initial conditions
Transients

Power

(MW)
Mass flux (10⁶ kg/m²h)
Pressure

(kg/cm²a)
Inlet temperature (°C)

1 1T
A4
2.50
11.18
156.2
291.0
Power increase

2.50
11.19
156.1
293.1
Flow reduction

2.52
11.28
156.3
291.7
Depressurization

2.48
11.04
154.6
291.6
Temperature

increase

12T
A8
2.51
11.40
156.1
291.3
Power increase

2.51
11.71
156.3
292.5
Flow reduction

2.50
11.42
156.2
290.6
Depressurization

2.50
11.38
155.8
291.2
Temperature

increase

Results from the transient DNB benchmark exercise are depicted in Figure 21. The plot includes the calculated times of DNB which were shifted by 6.5 sec and 6.3 sec for the temperature increase transients of test series 11T and 12T, respectively, in order to take into account the location of the inlet temperature measurement and the mass flow rate of each test series, as done for the transient void fraction benchmark. The results indicate that FLICA4 predicts DNB earlier than in the experiment. This is consistent with the results from the steady-state DNB benchmark where lower DNB power is predicted. A comparison with the preliminary results by other participants [30] where the same earlier DNB occurrence was predicted indicates that FLICA4 produces acceptable results also for transient DNB calculations.

Figure 21: Results of transient DNB exercise.

5.5. Assessment of CHF Models in FLICA4

FLICA4 includes three models to predict the CHF: the W3 correlation, the Groeneveld look-up table, and the SUDO correlation [31]. However, since the SUDO correlation was developed for rectangular channel, an assessment of this correlation is not conducted in this analysis.

The W3 correlation is one of the most widely used correlations for evaluation of DNB in PWRs and it is applicable to circular, rectangular, and rod bundle geometries. The correlation has been developed for axially uniform heat flux, with a correction factor for nonuniform flux distribution. In addition, local spacer effects can be taken into account by specific correction factors [32].

The Groeneveld look-up table was developed jointly by AECL (Canada) and IPPE (Russia) and has a very wide range of applications. Compared against the combined AECL-IPPE CHF database, it is known that the Groeneveld look-up table can predict the CHF data with an overall root-mean-square (RMS) error of 7.82%.

Selected cases from test series 0 and 8 of the benchmark are employed to assess the two CHF models. Each test series represents tests with or without a guide thimble at the center. The calculation results presented in Figure 22 indicate that, in general, the W3 correlation predicts slightly higher DNB power than the Groeneveld look-up table, that is the Groeneveld correlation produces more conservative results than the W3 correlation. The RMS errors are 5.2% and 6.1% for analyses with the Groeneveld look-up table and the W3 correlation, respectively. In addition, due to the limitations in the range of applicability, 3 out of 13 cases in test series 8 could not be analyzed with the W3 correlation. Therefore, considering all together the accuracy, the conservatism, and the range of applicability, the Groeneveld look-up table was deemed the most appropriate to estimate the CHF for this benchmark.

(a) Test series 0

(b) Test series 8

Figure 22: DNB Results with different CHF models.

6. Conclusion

In order to assess the range of validity and the accuracy of the subchannel analysis code FLICA4 and the commercial CFD code STAR-CD including a boiling model recently developed by CD-adapco, PSI has participated in an international benchmark based on the NUPEC PWR subchannel and bundle tests (PSBT) organized by OECD/NEA and US NRC. The tests have been analyzed by employing the US NRC code TRACE as well, in order to assess the applicability of TRACE to a subchannel analysis.

The results from the void distribution benchmark indicate that a reasonable agreement with the experimental data is obtained with FLICA4. The void fraction prediction by STAR-CD shows no significant discrepancy for the single Subchannel experiments. It is worthwhile mentioning that all the benchmark cases were calculated without any tuning of the boiling model and numerical algorithms. This fact demonstrates that the boiling model used in STAR-CD is able to predict void fractions over a wide range of void fraction values with acceptable accuracy. TRACE instead tends to overpredict the void fraction, especially for values lower than 40%. The analysis of the axial void fraction profile reveals that overprediction by TRACE is caused by an earlier increase of the void fraction along the axis of the channel, pointing out the necessity of additional assessments for the subcooled boiling model and bulk condensation model currently implemented in the TRACE code.

The DNB benchmark exercises have been analyzed with FLICA4 only. The steady-state benchmark exercise results indicate that FLICA4 slightly underpredict the DNB power. However, considering the accuracy of Groeneveld look-up table and the uncertainties in the experimental data, it can be concluded that the prediction of the DNB power by FLICA4 is acceptable and furthermore conservative. The transient DNB prediction reveals that FLICA4 predicts DNB earlier than in the experiments, which is consistent with the result from the steady-state DNB benchmark. In addition, an assessment of the CHF models of FLICA4 has been carried out by using the benchmark data. It is found that the Groeneveld look-up table predicts more conservative DNB powers with higher accuracy than the W3 correlation.

Acknowledgments

The authors would like to thank the benchmark organizers and colleagues participating in the PSBT benchmark for the cooperation and meaningful discussion. This work has been performed within the STARS project (http://stars.web.psi.ch/) of Paul Scherrer Institute and was partially funded by the Swiss Federal Nuclear Safety Inspectorate (ENSI).

References

B. Neykov, A. Aydogan, L. Hochreiter et al., “NUPEC BWR full-size fine-mesh bundle test (BFBT) Benchmark, volume I: specifications,” NEA/NSC/DOC(2005)5, 2006

A. Rubin, A. Schoedel, M. Avramova, H. Utsuno, S. Bajorek, and A. Velazquez-Lozada, “OECD/NRC benchmark based on NUPEC PWR subchannel and bundle tests (PSBT), volume I: experimental database and final problem specifications,” NEA/NSC/DOC(2010)1, 2010

I. Toumi, A. Bergeron, D. Gallo, E. Royer, and D. Caruge, “FLICA-4: a three-dimensional two-phase flow computer code with advanced numerical methods for nuclear applications,” Nuclear Engineering and Design, vol. 200, no. 1, pp. 139–155, 2000.

CD-adapco, CCM User Guide: STAR-CD Version 4.06, 2008.

US NRC, TRACE V5.0 Theory Manual: Field Equations, Solution Methods, and Physical Models, U.S. Nuclear Regulatory Commission, Washington, DC, USA, 2009.

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.

T. Kim, “Collaboration with CEA saclay for OECD-NEA PSBT benchmark—single channel benchmark with FLICA4,” PSI Memorandum SB-RND-ACT-003-10.002, 2010.

F. W. Dittus and L. M. K. Boelter, “Heat transfer in automobile radiators of the tubular type,” University of California Publications in Engineering, vol. 2, no. 13, pp. 443–461, 1930, Reprinted in International Communications in Heat and Mass Transfer, vol. 12, pp. 3-22, 1985

CEA, FLICA4 User Guide – Reference Manual of Modules and Procedures, 2009.

P. Raymond, B. Spindler, and R. Lenain, “Pressurized water reactor thermal-hydraulic core analysis with the FLICA computer code,” Nuclear Engineering and Design, vol. 124, no. 3, pp. 299–313, 1990.

F. P. Incropera and D. P. Dewitt, Introduction to Heat Transfer, John Willey & Sons, New York, NY, USA, 3rd edition, 1996.

L. Friedel, “Improved friction pressure drop correlations for horizontal and vertical two phase pipe flow,” in Proceedings of the European Two Phase Flow Group Meeting, Ispra, Italy, June 1979.

V. Gnielinski, “New equations flow regime heat and mass transfer in turbulent pipe and channel flow,” International Chemical Engineering, vol. 16, pp. 359–368, 1976.

M. A. Aggour, M. M. Vijay, and G. E. Sims, “A correlation of mean heat transfer coefficients for two-phase two-component flow in a vertical tube,” in Proceedings of the 7th International Heat Transfer Conference, vol. 5, pp. 367–372, 1982.

K. S. Rezkallah and G. E. Sims, “An examination of correlations of mean heat transfer coefficients in two-phase two-component flow in vertical tubes,” AIChE Symposium Series, vol. 83, no. 257, pp. 109–114, 1987.

R. T. Lahey, “A mechanistic subcooled boiling model,” in Proceedings of the 6th International Heat Transfer Conference, vol. 1, pp. 293–297, Toronto, Canada, 1978.

S. W. Churchill, “Friction factor equations spans all fluid-flow regimes,” Chemical Engineering, vol. 84, no. 24, pp. 91–92, 1977.

E. Ranz and W. R. Marshall, “Evaporation from Droplets: part I and II,” Chemical Engineering Progress, vol. 48, pp. 141–180, 1952.

M. Ishii and K. Mishima, “Study of two-fluid model and interfacial area,” Argonne National Laboratory Report ANL-80-111 (NUREG/CR-1873), 1980.

R. T. Lahey and F. J. Moody, The Thermal-Hydraulics of a Boiling Water Nuclear Reactor, ANS Monograph, 1977.

D. B. Spalding, “The calculation of free-convection phenomena in gas—liquid mixtures,” in Heat Transfer and Turbulent Buoyant Convection Studies and Applications for Natural Environment, Buildings and Engineering Systems, D. B. Spalding and N. Afgan, Eds., vol. 2, pp. 569–586, Hemisphere, New York, NY, USA, 1976.

D. B. Spalding, “Numerical computation of multi-phase fluid flow and heat transfer,” in Recent Advances in Numerical Methods in Fluids, C. Taylor and K. Morgan, Eds., vol. 1, pp. 139–168, Pineridge Press, Swansea, UK, 1980.

D. B. Spalding, “Developments in the IPSA procedure for numerical computation of multi-phase flow phenomena with interphase slip, unequal temperatures, etc,” in Numerical Methodologies in Heat Transfer, T. M. Shih, Ed., Proceedings of the Second National Symposium, pp. 421–436, Hemisphere, New York, NY, USA, 1983.

N. Kurul and M. Z. Podowski, “Multi-dimensional effects in subcooled boiling,” in Proceedings of the 9th Heat Transfer Conference, Jerusalem, Israel, 1990.

T. Hibiki and M. Ishii, “One-demensional drift-flux model for two-phase flow in a large diameter pipe,” International Journal of Heat and Mass Transfer, vol. 46, no. 10, pp. 1773–1790, 2003.

G. Kocamustafaogullari, “Pressure dependence of bubble departure diameter for water,” International Communications in Heat and Mass Transfer, vol. 10, no. 6, pp. 501–509, 1983.

L. S. Tong, “Heat transfer in water-cooled nuclear reactors,” Nuclear Engineering and Design, vol. 6, no. 4, pp. 301–324, 1967.

D. C. Groeneveld, L. K. H. Leung, P. L. Kirillov et al., “The 1995 look-up table for critical heat flux in tubes,” Nuclear Engineering and Design, vol. 163, no. 1-2, pp. 1–23, 1996.

A. Rubin, “Comparative analysis of submitted participants’ preliminary results for exercise II-2,” in Proceedings of the 2nd PSBT Workshop, Stockholm, Sweden, April 2011.

A. Rubin, “Comparative analysis of submitted participants’ preliminary results for exercise II-3,” in Proceedings of the 2nd PSBT Workshop, Stockholm, Sweden, April 2011.

Y. Sudo, K. Miyata, H. Ikawa, M. Kaminaga, and M. Ohkawara, “Experimental study of differences in DNB heat flux between upflow and downflow in vertical rectangular channel,” Journal of Nuclear Science and Technology, vol. 22, no. 8, pp. 604–618, 1985.

N. E. Todreas and M. S. Kazimi, Nuclear Systems: Vol. 1, Thermal Hydraulic Fundamentals, Hemisphere, New York, NY, USA, 1990.