Science and Technology of Nuclear Installations

Science and Technology of Nuclear Installations / 2012 / Article
Special Issue

Selected Papers from OECD-NEA PSBT Benchmark

View this Special Issue

Research Article | Open Access

Volume 2012 |Article ID 746467 |

Taewan Kim, Victor Petrov, Annalisa Manera, Simon Lo, "Analysis of Void Fraction Distribution and Departure from Nucleate Boiling in Single Subchannel and Bundle Geometries Using Subchannel, System, and Computational Fluid Dynamics Codes", Science and Technology of Nuclear Installations, vol. 2012, Article ID 746467, 20 pages, 2012.

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

Academic Editor: David Novog
Received30 Apr 2012
Accepted20 Jul 2012
Published11 Oct 2012


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.

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) [2123]. 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/m2sec;(iii)contact angle: 5°–90°;(iv)Number density: 1.0E+04 –1.51E+10 sites/m2.

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/m2-sec),(iii)power: 12.5–90.0 (kW),(iv)inlet temperature: 164.1–345.0 (°C).

Item Data

Assembly (subjected subchannel)

Subchannel typeCenter (typical)Center (thimble)sidecorner
Number of heaters4 × 1/4 3 × 1/4 2 × 1/41 × 1/4
Axial heated length (mm)1555155515551555
Axial power shapeUniformUniformUniformUniform

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

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

Flow area, mm2107.098107.09868.46442.592
Heated perimeter, mm29.84522.38414.9237.461
Wetted perimeter, mm54.64554.64544.92333.161

4.1.2. Thermal Hydraulic Model

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.

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.

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

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.

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.

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.

Item Data

Rods array 5 × 55 × 55 × 5
Number of heated rods252524
Number of thimble rods001
Heated rod outer diameter (mm)9.509.509.50
Thimble rod outer diameter (mm)12.24
Heated rods pitch (mm)12.6012.6012.60
Axial heated length (mm)365836583658
Flow channel inner width (mm)64.964.964.9
Radial power shapeAAB
Axial power shapeUniform Cosine Cosine
Number of MV spacers 777
Number of NMV spacers222
Number of simple spacers888
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.

NodeRelative Power


Spacer typeLoss Coefficient

Simple spacer (SS)0.4
Nonmixing vanes spacer (NMV)0.7
Mixing vanes spacer (MV)1.0

4.2.2. Thermal Hydraulic Model

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

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.

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.

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.


 Averaged error (%)1.183.541.592.10
 Mean absolute error (%)2.955.294.044.09
 Averaged error (%)9.547.453.356.78
 Mean absolute error (%)10.19.496.498.71

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.

Initial conditions
Test seriesAssemblyPressure (kg/cm2 a)Mass flux (106 kg/m2 h)power (kW)Inlet temperature (Celsius)Transients

5T B5154.211.952282.0300.4Power increase
153.811.932244.0301.2Flow reduction
152.511.942230.0301.7Temperature increase

6T B6158.211.552621.0288.1Power increase
158.412.032574.0288.8Flow reduction
157.211.922603.0288.8Temperature increase

7T B7158.212.022500.0291.9Power increase
158.112.042405.0292.0Flow reduction
158.811.992496.0290.2Temperature 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.

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.

AssemblyReference fuel typeRods
Type of cellPower distribution

A017 × 17 M5 × 5Typical cellAUniform
A1Typical cellCUniform
A2Typical cellAUniform
A36 × 6Typical cellDUniform
A45 × 5Typical cellACosine
A8Thimble cellBCosine
A11Typical cellACosine
A12Thimble cellBCosine

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.

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.

Test seriesTest
AssemblyTest modeMeasurement

15 × 5A1YY

36 × 6A3YY

11T5 × 5A11YY

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.

AssemblyInitial conditionsTransients
Mass flux (106 kg/m2h)Pressure
Inlet temperature (°C)

1 1TA42.5011.18156.2291.0Power increase
2.5011.19156.1293.1Flow reduction
12TA82.5111.40156.1291.3Power increase
2.5111.71156.3292.5Flow reduction

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.

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.

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.


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 ( of Paul Scherrer Institute and was partially funded by the Swiss Federal Nuclear Safety Inspectorate (ENSI).


  1. 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. View at: Google Scholar
  2. 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. View at: Google Scholar
  3. 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. View at: Publisher Site | Google Scholar
  4. CD-adapco, CCM User Guide: STAR-CD Version 4.06, 2008.
  5. US NRC, TRACE V5.0 Theory Manual: Field Equations, Solution Methods, and Physical Models, U.S. Nuclear Regulatory Commission, Washington, DC, USA, 2009.
  6. 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. View at: Google Scholar
  7. 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. View at: Google Scholar
  8. 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. View at: Google Scholar
  9. CEA, FLICA4 User Guide – Reference Manual of Modules and Procedures, 2009.
  10. 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. View at: Google Scholar
  11. F. P. Incropera and D. P. Dewitt, Introduction to Heat Transfer, John Willey & Sons, New York, NY, USA, 3rd edition, 1996.
  12. 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. View at: Google Scholar
  13. 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. View at: Google Scholar
  14. 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. View at: Google Scholar
  15. 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. View at: Google Scholar
  16. 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. View at: Google Scholar
  17. S. W. Churchill, “Friction factor equations spans all fluid-flow regimes,” Chemical Engineering, vol. 84, no. 24, pp. 91–92, 1977. View at: Google Scholar
  18. E. Ranz and W. R. Marshall, “Evaporation from Droplets: part I and II,” Chemical Engineering Progress, vol. 48, pp. 141–180, 1952. View at: Google Scholar
  19. M. Ishii and K. Mishima, “Study of two-fluid model and interfacial area,” Argonne National Laboratory Report ANL-80-111 (NUREG/CR-1873), 1980. View at: Google Scholar
  20. R. T. Lahey and F. J. Moody, The Thermal-Hydraulics of a Boiling Water Nuclear Reactor, ANS Monograph, 1977.
  21. 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. View at: Google Scholar
  22. 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. View at: Google Scholar
  23. 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. View at: Google Scholar
  24. N. Kurul and M. Z. Podowski, “Multi-dimensional effects in subcooled boiling,” in Proceedings of the 9th Heat Transfer Conference, Jerusalem, Israel, 1990. View at: Google Scholar
  25. 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. View at: Publisher Site | Google Scholar
  26. 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. View at: Google Scholar
  27. L. S. Tong, “Heat transfer in water-cooled nuclear reactors,” Nuclear Engineering and Design, vol. 6, no. 4, pp. 301–324, 1967. View at: Google Scholar
  28. 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. View at: Google Scholar
  29. A. Rubin, “Comparative analysis of submitted participants’ preliminary results for exercise II-2,” in Proceedings of the 2nd PSBT Workshop, Stockholm, Sweden, April 2011. View at: Google Scholar
  30. A. Rubin, “Comparative analysis of submitted participants’ preliminary results for exercise II-3,” in Proceedings of the 2nd PSBT Workshop, Stockholm, Sweden, April 2011. View at: Google Scholar
  31. 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. View at: Google Scholar
  32. N. E. Todreas and M. S. Kazimi, Nuclear Systems: Vol. 1, Thermal Hydraulic Fundamentals, Hemisphere, New York, NY, USA, 1990.

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.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.