Analysis of the NUPEC PSBT Tests with FLICA-OVAP

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


Introduction
Based on NUPEC PWR subchannel and bundle tests (PSBT), an international benchmark has been promoted by OECD and NRC and has been coordinated by Penn State University [1].The aim of this benchmark is to encourage advancement and assessment of numerical models in subchannel analysis of fluid flow in rod bundles, which has very important relevance for the nuclear reactor safety margin evaluation.An important database of void fraction and critical heat flux measurements in steady-state and transient conditions has been carried out by NUPEC on a prototypical PWR rod bundle.Different types of subchannel or rod bundle geometries and a wide range of flow conditions at high pressure have been investigated (see Table 1) allowing the assessment of the key models and correlations in these conditions.
The Service de Thermohydraulique et de Mécanique des fluides (STMF) at the Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA), France, has been involved in the PSBT benchmark performing calculations with the FLICA-OVAP code [2].FLICA-OVAP is an advanced twophase flow thermal-hydraulics code based on a full 3D subchannel approach.It is designed to analyze flows in light water reactors cores such as PWRs, BWRs, and experimental reactors.To provide a relevant answer to different core concepts and multiple industrial applications, several models coexist in the FLICA-OVAP platform: the homogeneous equilibrium model, the four equations drift flux model, the two-fluid model, and finally, a general multifield model, with a variable number of fields for both vapor and liquid phases.For each model, an adapted set of closure laws is proposed concerning heat and mass transfer, interfacial and wall forces, and turbulence.

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

The Four-Equation Drift Flux
Model.The four balance equations describing two-phase flows in the drift flux model are, respectively, the mixture mass balance equation, the mixture momentum balance equation, the mixture energy balance equation, and the steam mass balance equation (porosities are omitted for the sake of simplicity): (i) mixture mass conservation where α k , ρ k , u k are the volume fraction, the density, and the velocity for the phase k.
(ii) mixture momentum balance where P is the pressure, g, the gravity, and F w , the friction forces.The tensor τ k represents the viscous and the Reynolds stress terms for the phase k.The mixture density ρ is defined as (iii) mixture energy balance where E k and H k are the total energy and the total enthalpy of the phase k, q k includes molecular and turbulent heat fluxes, and q w is the volumetric source term of thermal power.
(iv) steam mass The model is closed by a drift flux correlation and a general equation of state with the assumption that, in presence of liquid, the vapor is in saturation conditions.

The Drift Flux
Correlation.FLICA-OVAP includes several Zuber-Findlay type correlations in order to estimate the relative velocity u r between the vapor velocity u v and the liquid velocity u l .The general form of these correlations is: where C 0 is the distribution parameter, j = αu v +(1−α)u l is the area-averaged total volumetric flux, and V g, j is the voidweighted area-averaged drift velocity.The Chexal-Lellouche correlation [4] and a correlation derived from Ishii [5] are implemented in the code.The Chexal-Lellouche correlation covers a large range of pressure, diameters, flows whereas the Ishii correlation, firstly established for bubbly, slug, and churn turbulent flow in adiabatic conditions, allows taking into account nucleate boiling on heated walls.

Pressure Drop.
Friction forces F w are given by the sum of distributed F fric and singular pressure drops F sing Singular friction due to spacer or mixing grids or other pressure drops are given by where K sing is an antisymmetric tensor.Distributed friction at walls is instead accounted for by Science and Technology of Nuclear Installations 3  where the different friction terms f w are given by the product of the isothermal friction factor f iso , the heating wall correction f heat , and the two-phase flow multiplier f 2φ .The isothermal friction factor in turbulent flows is given by whereas the heating wall correction was estimated by an inhouse model already used in the FLICA-4 code [6].In this analysis the Chisholm correlation [7] was used for f 2φ , given by where φ 2 lo is the adiabatic two-phase frictional pressure drop multiplier, C φ , a parameter accounting for heat flux, D heat is the heated diameter, and q , the wall heat flux.In the runs of the PSBT benchmark C φ was set equal to 0 and thus f 2φ = φ 2 lo .

Diffusion Effects.
To account for viscous and turbulent diffusion effects, the tensor τ k is introduced in the momentum equation, given by where μ k M i j t,k is the turbulent viscosity, which is limited to the liquid phase.An anisotropic formulation is used for turbulent viscosity where Re = GD h /μ l is the Reynolds number, M  Table 3: Choice of X and K for the Katto and Ohno model.
Similarly, molecular and turbulent heat fluxes are given by k=v,l where h x = xh v + (1 − x)h l is the flow enthalpy based on the flow quality x.The turbulent conductivity is given by where K i j t0 , b K , Re t are parameters and f K ( f 2ϕ ) is a function of the two-phase flow multiplier.In this analysis, the value of M i j t0 and K i j t0 was varied as a function of the axial position.Two different values have been adopted for each parameter depending whether the considered axial position was downward a mixing or a spacer grid.
2.5.Wall Temperature.Wall temperatures are estimated on the basis of the bulk temperature and the heat transfer coefficient as The Nusselt number and the bulk temperature depend on the heat transfer regime.Four different regimes can be distinguished: single-phase convection heat transfer, subcooled nucleate boiling (SNB), saturated nucleate boiling (SANB) and post-critical-heat-flux (post-CHF) heat transfer.In single-phase heat transfer and SNB, the bulk temperature is equal to the liquid phase temperature, whereas in SANB it is equal to the saturation temperature.In forced convection conditions, the single-phase heat transfer coefficient is estimated by the Dittus-Boelter correlation.The onset of significant void (OSV), which is the transition between single-phase heat transfer and SNB can be predicted with the Forster and Grief correlation in low-pressure conditions or the Jens and Lottes correlation in high pressure conditions.In the present analysis, the Jens and Lottes correlation was used [8]    on the boiling characteristics, whether it is IAFB (inverted annular film boiling) or DFFB (dispersed flow film boiling).

The Mass
Transfer Term Γ v .The mass transfer term Γ v appearing in the steam mass balance equation is given by the sum of two contributions: the vapor generation on walls Γ wv and the mass transfer between the liquid and the vapor phase Γ vl .In subcooled nucleate boiling, only a portion χ v of the heat flux transferred from the wall to the mixture is used to vaporize the liquid phase, whereas the remaining part is used to heat the liquid phase up.
The vapor generation at walls is thus given by where χ v is a function of the saturation temperature, the liquid phase temperature, and the wall superheat demanded to have subcooled nucleate boiling, given by The mass transfer between the two phases Γ vl is instead given by where h v and h l are vapor and liquid enthalpy and q vl is the heat transferred between the two phases, given by where K v0 is a mass transfer parameter, G is the mass flux, Re = GD h /μ l is the Reynolds number, c * , the equilibrium quality based on the mixture enthalpy h = ch v + (1 − c)h l , f (P, ρ, μ l , u, u r ), a function depending on local conditions, D h , the hydraulic diameter, and μ l , the liquid viscosity.Re 0 is a parameter of the model.

Prediction of the Boiling Crisis
To predict boiling crisis conditions, several models are currently available within FLICA-OVAP.The W3 correlation [9] is appropriate to predict departure from nucleate boiling (DNB) phenomena of interest for pressurized water reactors.However, boiling crisis experienced in NUPEC tests can be close to dryout conditions.Two models have been therefore implemented in the code: the model of Shah [10] and the model of Katto and Ohno [11].Both models can predict critical heat flux values for both DNB and dryout conditions.An accurate description of them is given in the following sections.
3.1.The Shah Model.The Shah model [10] consists of two separate correlations to determine the boiling number Bo, defined as The first correlation covers conditions where the critical heat flux depends on the upstream conditions, named UCC (upstream conditions correlation).The second, named LCC (local conditions correlation), depends only on local quantities.
with In the previous equations, z eff is the effective tube length and x i,eff is the effective inlet quality, defined as where z crit is the distance from the inlet section and the location of the boiling crisis, where the critical heat flux is  calculated, and z sat is the boiling length, which for uniformly heated tubes is given by For water flows, when Y ≤ 10 4 , n = 0, otherwise it is given by for Y ≤ 10 6 , for Y > 10 6 . (27)

The LCC Correlation. The LCC correlation is expressed by Bo
The entrance factor F E is given by When the previous correlation gives  where P r is the reduced pressure given by P r = P/P c where P c is the critical pressure, that is 22.064 MPa for water.The value of F X depends on the quality at the location of the boiling crisis x crit .When x crit > 0, the following equation is used When P r > 0.6, c = 1, otherwise c = 0. F 3 is instead given by When x crit < 0, F X is given by As in the previous case, when P r > 0.6, b = 1, otherwise b = 0. F 1 is given by , Finally, F 2 is given by (35)

Choice between UCC and LCC Correlation.
For water, the UCC correlation is used when Y ≤ 10 6 .When Y > 10 6 , the correlation giving the lower value of the boiling number Bo is used, with exception of cases where z eff > 160/P 1.14 r , for which the UCC formulation is always adopted.In this analysis, the UCC correlation was used in all calculations.

The Katto and Ohno Model. The correlation of Katto and
Ohno [11] provides the value of the critical heat flux as where Δh sub,i is the subcooling inlet enthalpy.The terms X and K are functions of three-dimensionless terms: Five different values of X must be determined  The value of C in the relationship of X 1 is given by C = 0.25 for Z < 50, C = 0.25 + 0.0009(Z − 50) for 50 < Z < 150, Three values of K must be determined The appropriate values of X and K must be determined according to Table 3.All experimental tests are included in the range covered by the Katto and Ohno correlation (see Table 4).

Subchannel Exercises.
Four series of measurements of void fraction were performed in sections representative of subchannel types in a PWR assembly.Table 5 gives the geometry (center, typical, and thimble subchannels, side and corner subchannels) and power shape for these series.The heated length is 1555 mm, and the void measurement section is located at 1400 mm from the bottom of the heated section.
According to the given test matrix (power value and boundary conditions), the flow is saturated at the location where the void fraction is measured for the most runs.This situation allows only partial assessment of the subcooled boiling models.In contrast, drift flux model plays a major role to obtain good agreement with experimental data.

Chexal-Lellouche and Ishii correlations have been compared
for series 1 (Figure 1).It is shown that the Ishii correlation tends to underpredict the void fraction more than the Chexal-Lellouche one.Even if it is not obvious to extend this conclusion to other configurations, the Chexal-Lellouche correlation has been thus adopted for other subchannel and bundle calculations.
Calculated densities and void fractions are compared against the experimental data in Figure 2. Experimental densities measured by CT scan are declared to be affected by an uncertainty of 15 kg/m 3 (σ exp = 15kg/m 3 in Figure 2(b)) [12].Experimental void fraction resumed basing on measured densities are presumed to be affected by an uncertainty of 3% of void fraction (σ exp = 3% in Figure 2(a)) [12].It happens that densities obtained by the code match well the experimental ones, whereas void fractions are slightly underestimated for the lowest value of this quantity.Indeed, physical properties and correlations adopted to convert density to void fraction could be the cause for this discrepancy.

Subchannel Exercises.
A partial section of the full length 17 × 17 type PWR fuel assemblies was considered.The rod bundle test is a 5 × 5 square array.The heated section is 3658 mm high, and density measurements are set at 2216 mm (Lower), 2669 mm (Middle), and 3177 mm (Upper).In Table 6 are shown the geometry and power shape of the rod bundle test assemblies considered.Three assembly configurations were considered: B5, B6, and B7.The main difference between these three configurations is in the power distribution.Assembly B5 has a uniform axial power distribution and pattern A radial power distribution.Assembly B6 has the same radial power distribution of assembly B5 but cosine axial power distribution.Assembly B7 has a cosine axial power distribution and pattern B radial power distribution.For this radial power distribution, the central rod is a thimble rod (i.e., no power is generated inside this rod) of larger diameter with respect to the fuel rods.The gamma-ray transmission method (chordal averaged) was adopted in the tests performed to measure the density and then converted to the void fraction of the gas-liquid two-phase flow.The declared uncertainties are 4% of void fraction (σ exp = 4% in Figures 3 and 4) for steady-state bundle tests and 5% of void fraction (σ exp = 5% in Figures 5 and 6) for transient bundle tests.In Table 7 are shown the test series for void fraction measurements.

Steady-State Exercises.
For the results at the lower and middle elevations, subcooled nucleate boiling models (condensation and wall heat transfer) play a major role for the prediction of void fraction.At the upper elevation, the fluid is generally saturated, and the drift model becomes important.Moreover, diffusion terms play also a major (z) = m tb (the turbulent Prandtl number is assumed equal to 1) adopted for bare bundles (without any spacer) have to be increased in the region upstream mixing grid spacers in order to reproduce the enhancement of turbulent mixing.In (15) and ( 13), K i j t0 and M i j t0 have been therefore adapted for the PSBT assemblies by means of piecewise functions: it was assumed that K i j t0 (z) = k ts and M i j t0 (z) = m ts downstream a mixing grid spacer, instead of taking the standard values k tb and m tb .Despite the values of k tb and m tb are reliably assessed to be around 0.010 for PWR type bundles [6], the appropriate values to be adopted for k ts and m ts are not known for these particular mixing grid spacers and thus, to better identify reasonable values for these coefficients, a sensitivity analysis has been performed.A sensitivity analysis for the coefficient K v0 in the condensation model was also conducted.A sample result is shown in Figure 3, where combined effects of these are reported for run 5.2442.In the subcooled and at the beginning of the saturated region, void fraction is strongly reduced as K v0 increases, but this parameter does not affect results at the end of the bundle, where the flow is fully saturated.On the other hand, the effect of the turbulent diffusion coefficients k ts and m ts is significant when the void fraction increases, since it determines the flatness of the void distribution profile in the bundle.Detailed experimental density maps would have been useful to assess a reliable value for the turbulent diffusion coefficients, but these maps are not readily exploitable.As far as we could deduce basing on the void fraction profile in the central subchannel, a value for k ts and m ts was chosen equal to 0.045.The value of K v0 was instead taken equal to 1.5 × 10 −4 , which is in the standard range of values adopted in thermal-hydraulic analysis performed with this model.Figure 4 shows the calculated void fraction against the experimental values for series 5, 6, and 7, with the following parameters: K v0 = 1.5 × 10 −4 , k ts = m ts = 0.045, k tb = m tb = 0.01.

Transient Exercises.
For the transient simulations, three different configurations were taken into account as shown in Table 7. Test series 5T considers the same assembly conditions used for steady-state series 5, with a uniform axial power distribution, and pattern A of radial power distribution.
For transient test series 6T the same assembly conditions used for steady-state series 6 were applied, which means a cosine axial power distribution, and pattern A radial power distribution.Test series 7T considers the same assembly conditions used for steady-state test series 7, which are cosine axial power distribution, and pattern B radial power distribution.
For each transient analysis, four different scenarios were considered: power increase, flow reduction, depressurization, and temperature increase.Boundary conditions were provided for each scenario, using table where the dependent variables (i.e., power, pressure, mass flux, and inlet temperature resp.,) were specified as function of the simulation time.Same parameters and models as for steadystate bundle computations were used.
Figure 5 shows complete results for transients 5T, also representative of the other transient series 6T and 7T.FLICA-OVAP presents an overall agreement against the experimental void fraction for the transients, although discrepancies with experimental data can be observed, in particular for the temperature increase transient.For the power increase transient, the evolution of the void fraction follows the data, but the void fraction is underestimated at the middle and upper locations.The results concerning the depressurization have more discrepancy as time increases, but are consistent with the experimental evolution.
For the temperature increase transient, a shift of calculated void fractions profiles with respect to experimental values can be noticed.This discrepancy has been noticed by several participants to the first PSBT/OECD workshop in Pisa.Fluid temperature probe is located in a pipe between the preheater and the inlet nozzle of the test section [12].It is therefore reasonable to guess that a delay of several seconds occurs between the measurement point and the inlet of the heated section.This point has been investigated by JNES [12] but the distance from the measurement location to the inlet of the heated section (where the temperature boundary condition has to be applied in thermal-hydraulics codes) was not clarified.As far as we can say basing on the result obtained by FLICA-OVAP, if a delay of 6 s is assumed for the inlet temperature with respect to other boundary conditions (i.e., the initial temperature is maintained during 6 s and then the temperature variation is applied), the agreement between calculated and experimental void fraction profiles is improved.It can clearly be seen in Figure 6, where these profiles are compared for the calculations with and without the time delay of 6 seconds.

DNB Results
The two models adopted in this work, the Shah and the Katto and Ohno models, have been developed and tested mostly in single-channel uniform heat flux conditions.
In this work, the applicability of them to bundle geometries with nonuniform heat flux profile has been investigated.Each experimental series of the NUPEC database has been calculated by the two models, adopting a twofold approach.
In a first step, homogenized one-dimensional calculations have been performed, adopting a bundle-averaged description.
In a second step, three-dimensional calculations have been performed, adopting a subchannel description.In this case, the two boiling crisis models are applied at the subchannels scale.The fundamental difference with the first approach consists in having different hydraulic diameters depending on whether the selected subchannel is a central, a side or a corner subchannel.However, some special approach has to be taken for side and corner subchannels, since CHF models tend to underestimate their critical heat flux.Moreover, as observed by Tong and Tang [13, Section 5.4.5.2, Figure 5.54], for a given inlet enthalpy, critical heat fluxes in the presence of cold walls are even higher than for internal subchannels.For the sake of simplicity, the possibility to have boiling crisis in side and corner subchannels was thus excluded.
In the following paragraphs, for each series, results obtained via the bundle-averaged and the subchannel approach will be shown for the two models.The performances of the different models are also summarized in Table 8.

Steady-State Series
5.1.1.Series 0, Series 2, and Series 3. Series 0, Series 2 investigate boiling crisis in rod bundle (rods array: 5 × 5) and uniform heat flux.The two series differ only for the position of mixing and spacer grids.In Series 3, the bundle is also modified, adopting an array of 36 rods (6 × 6).The axial power shape is uniform for all series, but peripheral rods have a slightly lower power with respect to central rods: the radial power distribution factor is 0.85 for peripheral rods and 1.0 for central rods.
In Figures 7, 8, and 9, calculated and experimental critical heat fluxes are compared, respectively, for Series 1, Series 2, and Series 3. It is ascertained that the Shah model tends to predict higher critical heat flux values than the Katto and Ohno model, for both the bundle-averaged and the subchannel approach.Moreover, in this particular configuration with uniform heat flux and a moderate radial offset, 1D or 3D calculations give very similar results.13.Series 4 and Series 13 introduce the effect of nonuniform axial profile.A cosine shape was adopted on the same bundle geometry used for Series 2. To deal with nonuniform heat flux, a correction factor was adopted, as proposed by Tong and Tang [13].In Figures 10 and 11 calculated and experimental critical heat fluxes are compared, respectively, for Series 4, Series 13.As for previous tests, the Shah model predicts higher critical heat flux values than the Katto and Ohno model, providing a better agreement with experimental values.Minor differences are experienced between the homogenized and the subchannel description.

Series 4 and Series
5.1.3.Series 8. Series 8 includes the presence of a central thimble rod, representing guide tube for control rods.Contrarily to the previous series, the results of the homogenized and the subchannel approaches are different (see Figure 12).As expected, in the presence of a radial heterogeneous power shape, a better description is achieved with the subchannel approach.In general the Shah model gives a better estimate of critical heat fluxes.

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

Conclusions and Perspectives
This paper has presented the current capabilities of the FLICA-OVAP code in predicting void distribution and boiling crisis phenomena.The NUPEC database released in the frame of the OECD/PSBT international benchmark has been addressed.
Void fraction measurements in the subchannel configurations are of major interest for the validation of mass transfer model and the OSV criterion.In this analysis, the attention was focused on the mass transfer model, but further improvement of this work could include the analysis of other OSV criteria.Results obtained by FLICA-OVAP with a set of standard coefficients for the different models show a good agreement of the calculated densities and a slight underestimation of the void fraction at the measurement location, mainly located in the saturated regime for the considered runs.
For the steady-state bundle tests, the K v0 coefficient and diffusion coefficients k t and m t are the key parameters to fit the void fraction.In particular, mixing grid spacers play a major role, since they enhance the turbulence in the downstream flow.It was found that a reasonably good agreement between calculated and experimental void fraction profiles is achieved when the turbulent diffusion associated to mixing grid is taken equal to k ts = m ts = 0.045.Adopting the same set of parameters, a reasonably good agreement between calculated and experimental void fraction was also ascertained for transient tests, even if a systematic underestimation of void fraction at the middle and upper measurement locations has been found.
Discrepancies noticed on temperature increase transients were corrected by applying a delay of 6 s of the original inlet temperature boundary conditions, in order to simulate the residence time of the fluid between the temperature probe where the boundary conditions are given and the inlet of the test section.
Results obtained from the subcooled region suggest possible directions to be pursued in order to improve the current modeling: further developments and validation will involve the OSV criterion and the modeling of heat flux and heat flux partitioning in subcooled nucleate boiling, but

Figure 2 :
Figure 2: Calculated void fractions (a) and densities (b) versus experimental data for subchannel series 1 to 4 with the Chexal-Lellouche drift correlation.
i j t0 , b M , Re t are parameters, and f M ( f 2ϕ ) is a function of the two-phase flow multiplier.

Figure 3 :
Figure 3: PSBT run 5.2442.Sensitivity analysis for K v0 (a), k ts and m ts (b) on the void fraction profile in the central subchannel.

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

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

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

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

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

Table 7 :
Test series for void fraction measurements.

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