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

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.


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.

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 systemcodes, 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.

Pressure
Heat flux

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.

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.

TRACE.
TRACE is the latest best-estimate thermalhydraulic 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 6equation 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 voidfraction-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].

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][22][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   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 q C : convective heating, q q : quenching, and q e : evaporation.
For the nucleation site density, the Hibiki and Ishii model [25] with the following application range was used: 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.

PSBT Benchmark Phase I:
Void Distribution Benchmark

Thermal Hydraulic Model
(1) 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, K t and M t , were set to 0.01 based on a sensitivity analysis.A value of 1.0E − 4 was assigned to the recondensation coefficient, KV0.As for the TRACE model, the single channel was modeled with a PIPE (onedimensional) component.Inlet mass flux and temperature and outlet pressure were imposed as boundary conditions.
(2) 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.
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.

Results
(1) 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 R 2 ) 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.
(2) STAR-CD.Figure 8 shows the radial void fraction distribution in the measuring crosssection 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.

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 nonheated 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.

Thermal Hydraulic Model. Analyses of the steadystate bundle experiments have been conducted by means of FLICA4 and TRACE only.
(1) 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 7.5E − 4 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, K t and M t , were set to 0.01 based on a sensitivity analysis.For assembly B7 that has a guide thimble at the center, K t and M t 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.
(2) TRACE.The bundles are described by using a 3D rectangular VESSEL component of TRACE as depicted in    7.
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 steadystate 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.

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 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.

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 Science and Technology of Nuclear Installations    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.

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 7.5E − 4 was imposed for the recondensation coefficient, KV0.The multipliers for turbulent diffusivity and viscosity, K t and M t , 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.

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.

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.
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.

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.

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.

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

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

Figure 9 :
Figure 9: Axial void profile for 4 runs with S1 type of geometry.
Top-down view of test section

Figure 12 .
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 Kfactors were prescribed at the relevant elevation in order to model the pressure drop introduced by spacers.No local Kfactor 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.

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

Figure 22 :
Figure 22: DNB Results with different CHF models.

Table 6 :
Mean absolute error of void fraction.

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 crosssectional view and the geometrical parameters of each subchannel test assembly are depicted in Figure2and listed in 4.1.Single Subchannel Experiments4.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
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 Figure13 indicate where smaller void fraction is expected.TRACE tends to overpredict the void fraction especially for values below 40%.In Figure15, 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