Abstract

This paper presents the assessment of CATHARE 3 against PWR subchannel and rod bundle tests of the PSBT benchmark. Noticeable measurements were the following: void fraction in single subchannel and rod bundle, multiple liquid temperatures at subchannel exit in rod bundle, and DNB power and location in rod bundle. All these results were obtained both in steady and transient conditions. Void fraction values are satisfactory predicted by CATHARE 3 in single subchannels with the pipe module. More dispersed predictions of void values are obtained in rod bundles with the CATHARE 3 3D module at subchannel scale. Single-phase liquid mixing tests and DNB tests in rod bundle are also analyzed. After calibrating the mixing in liquid single phase with specific tests, DNB tests using void mixing give mitigated results, perhaps linked to inappropriate use of CHF lookup tables in such rod bundles with many spacers.

1. Introduction

CATHARE 3 is a new two-phase thermalhydraulics system code developed at CEA Grenoble [1]. It has been designed to expand the capabilities of CATHARE 2 and to improve the simulation accuracy of light water reactor accidents. New features include additional field, like a droplet field or a bubble field, and coupled equations of turbulence transport for a continuous field or interfacial area transport for a dispersed field. Beside the unchanged choices for numerical schemes for time and space discretization, a numerical solver gathering the different modules of a circuit has been rewritten and improved compared to CATHARE 2 in order to allow new capabilities of coupling with external codes, for example, for neutronics, or detailed CFD. A preliminary version V1 needs a wide validation program. This paper deals with the 1D and 3D module validation of the code against various boiling experiments at subchannel scale.

Following the BWR Full-size Fine-Mesh Bundle tests (BFBT) benchmark, the PWR subchannel and bundle tests (PSBT) benchmark [2] is proposed by OECD/NRC. Both are based upon a NUPEC database obtained in full-scale subchannels and rod bundles and include detailed measurements of fluid temperature, void fraction, and critical power or DNB power in steady and transient conditions. These experiments are useful to check and validate the code closure laws in rod bundles, especially the turbulence dispersion coefficients for heat in single-phase flow and void in two-phase flow, the wall and interfacial friction coefficients, and the wall-to-fluid heat transfer models. The PSBT phase I exercises are devoted to the void fraction measurements, performed in single subchannels and in 5Γ—5 rod bundles in steady and transient conditions. In the first exercise, phase II features liquid temperature measurements in all subchannels of a heterogeneously heated rod bundle in steady conditions, and, in the following exercises, DNB measurements taken in various rod bundles in steady and transient conditions, that is, power value and location of first detected DNB.

Single subchannel experiments are simulated by the CATHARE 3 pipe 1D module while rod bundle cases are simulated with the CATHARE 3 3D module meshed at a subchannel scale, that is, one cell per subchannel in a horizontal cross cut. The 3D module for the rod bundle has been coupled with a 1D module in order to improve the inlet flow simulation along the downcomer.

Useful balance equations and closure laws are briefly presented in the following Section 2. Then, results of comparisons between simulations and measurements of void fraction for Phase I exercises are presented in Section 3. In Section 4 are presented the results of temperature and DNB simulations of Phase II exercises.

2. CATHARE 3 Balance Equations and Closure Laws

Both 1D and 3D modules of CATHARE 3 solve the same set of balance equations, except that the energy balances are written using enthalpy in the 3D module and internal energy in the pipe module. The closure laws remain identical as far as possible. A first-order donor cell scheme is used in both modules as far as space discretization is concerned. For time discretization, the pipe module calls a fully implicit scheme, while the 3D module uses a semi-implicit scheme.

Contrary to the preceding BFBT simulations [3], featuring high void two phase flow, most of the PSBT benchmark database remain in the low- or medium-void range and hence, simulations do not need an additional droplet field beside the standard 6-equation model because the boiling flow regime never changes towards an annular dispersed flow.

For a given generation of steam along a single heated channel, the local void fraction is governed by wall and interfacial friction. In a 3D flow inside a rod bundle, cross-flows between adjacent subchannels lead to void dispersion. Also turbulent dispersion or diffusion may affect the temperature map in the single phase region. The void dispersion phenomena can be modelled by a mixing term in the momentum balance equations. The temperature dispersion (caused by nonrandom flow from one subchannel to a neighbour) and diffusion (caused by random fluctuations of flow between adjacent subchannels) are modelled by a single term in the liquid energy balance equation. The velocity diffusion is presently neglected in the momentum equation; its implementation had no effect on results of several tests, either in single phase or two phase flow in bundles.

2.1. Momentum Balance

Consider the following equation: π›Όπ‘˜πœŒπ‘˜ξ‚ƒπœ•π‘‰πœ•π‘‘π‘˜+π‘‰π‘˜β‹…βˆ‡π‘‰π‘˜ξ‚„=βˆ’π›Όπ‘˜βˆ‡π‘ƒ+π‘π‘–βˆ‡π›Όπ‘˜+(βˆ’1)π‘˜πœπ‘–+πœπ‘π‘˜+π›Όπ‘˜πœŒπ‘˜π‘”.(1) Wall frictions πœπ‘ of both phases are calculated using the Blasius friction coefficient π‘“π‘˜ multiplied by a phase-dependent multiplier π‘π‘˜πœπ‘π‘˜=πœ’π‘π‘˜π‘“π‘˜πœŒπ‘˜||π‘‰π‘˜||π‘‰π‘˜2,with𝑐𝑔=𝛼1.25,𝑐𝑙=(1βˆ’π›Ό)πœŒπ‘™(1βˆ’π›Ό)πœŒπ‘™+π›ΌπœŒπ‘”.(2) Interfacial friction in bubbly, slug, and churn vertical flow is given by πœπ‘–=πΎπ‘™πœŒπ‘™+πΎπ‘”πœŒπ‘”πΏπ›Ό(1βˆ’π›Ό)3.6ξ€Ίπ‘‰π‘”βˆ’π‘‰π‘™ξ€»2.(3)𝐿 is the maximum bubble size, limited by the Laplace length β„“ and the hydraulic diameter 𝐾𝑔=29𝐾𝑙=ξ€·πΉπœ‡ξ€Έ0.25𝑓𝑙withπΉπœ‡=πœ‡π‘™βˆšπœŒπ‘™πœŽβ„“,𝑓𝑙𝐿=2.81+34π·β„Žξ‚Ά5ξ‚΅6βˆ’5πΏπ·β„Žξ‚Ά.(4)πœπ‘– and πœπ‘ are unchanged compared to CATHARE 2 6-equation model.

The mixing term 𝑝𝑖 is calculated from an assessment of the turbulent kinetic energy.

For a single-phase flow in a tube or a subchannel (far from a spacer grid), the turbulent kinetic energy π‘˜π‘™ can be assessed by π‘˜π‘™=0.0367𝑉2𝑙Reβˆ’1/6(5) (see [4]). The associated turbulent viscosity can be assessed as πœˆπ‘‘=0.5π·π»βˆšπ‘˜π‘™(6) and the dispersion term 𝑝𝑖: 𝑝𝑖=0.4πœ‡π‘‘π‘‰π‘™π·π».(7) The coefficient 0.4 comes from an order of magnitude for the velocity gradients between subchannels; the velocity difference is evaluated at 40% of the axial velocity, which is close to the velocity module.

At the end, it comes 𝑝𝑖=0.038πœŒπ‘™π‘‰π‘™2Reβˆ’1/12.(8) The coefficient 0.5 in the πœˆπ‘‘ formula has been adjusted so as to better match the void fraction measurements in the PSBT Phase I tests, and the temperature measurements in the Phase II, given a Pr𝑑 equal to 1 (see (9) below). This coefficient appears to be several orders of magnitude above the figure calculated using simple turbulence [4, 5]; the void dispersion due to cross flows (and not only diffusion) seems to be the main driving phenomena (see [6]).

2.2. Continuous Liquid Energy Balance

It is written using internal energy 𝑒𝑙 for the 3D module as follows: πœ•ξ€·π›Όπœ•π‘‘π‘™πœŒπ‘™π‘’π‘™ξ€Έξ€·π›Ό+βˆ‡β‹…π‘™πœŒπ‘™π‘’π‘™π‘‰π‘™ξ€Έ=π‘žπ‘™π‘–+πœ’π‘π‘žπ‘π‘™βˆ’Ξ“π»π‘™π‘ξ‚Έβˆ’π‘ƒπœ•π›Όπ‘™ξ€·π›Όπœ•π‘‘+βˆ‡.𝑙𝑉𝑙𝛼+βˆ‡β‹…π‘˜ξ‚΅πœ†π‘™π‘‡π‘™+πœŒπ‘™πœˆπ‘‘π‘™Prπ‘‘βˆ‡π‘’π‘™.ξ‚Άξ‚Ή(9) The molecular diffusion is neglected compared to the turbulent diffusion term.

To summarize, the turbulence is modelled here by two different algebraic terms: one described just above in the energy balance and the 𝑝𝑖 grad 𝛼 in the momentum equations. No additional transport equation of turbulence quantity was solved in this study.

The Departure from Nucleate Boiling appears on a hot wall when the heat flux towards the fluid exceeds the so-called β€œCritical Heat Flux,” which is assessed in six-equation model of CATHARE 2 and CATHARE 3 using home-made polynomials interpolating CHF lookup tables, given the local values of mass flux, pressure, and steam quality. The tables are based on the 1995 Groeneveld tables [7], using a simple rod bundle coefficient but no effect of the spacer grids.

3. Overview of the PSBT Benchmark

The PWR Subchannel and Bundle Tests (PSBT) benchmark is proposed by OECD/NRC. Pennsylvania State University (PSU) under the sponsorship of U.S. Nuclear Regulatory Commission (NRC) prepared the specification and organized the benchmark with the Japan Nuclear Energy Safety (JNES) Organization. The Nuclear Power Engineering Corporation (NUPEC) released a database including various single subchannel and full-scale rod bundle tests in boiling conditions, with detailed void distribution and DNB measurements. Both system codes and CFD codes can match their results against the averaged (macroscopic data at subchannel scale) or fine experimental results.

Two phases are proposed: The first one is devoted to the void distribution and includes four exercises. Exercise 1β€”steady-state single subchannel benchmark, Exercise 2β€”steady-state bundle benchmark, Exercise 3β€”transient bundle benchmark, Exercise 4β€”pressure drop benchmark. The second is devoted to DNB prediction and includes three exercises. Exercise 1β€”steady-state fluid temperature benchmark, Exercise 2β€”steady-state DNB benchmark, Exercise 3β€”transient DNB benchmark.

This benchmark, especially through its accurate measurements, is a very good opportunity to assess the capabilities of system codes such as CATHARE 3 to simulate boiling flows in PWR core geometry.

4. Void Fraction in PSBT Phase I Exercises

4.1. Exercise I.1 Steady State in Single Subchannels

Four different geometries of single subchannels, for central, side, and corner locations (side and corner are relative to a rod bundle in a square box) have been tested, resulting in void fraction measurements at 1400 mm level inside a 1555 mm long heated subchannel (Figure 1). They fit to standard PWR rod bundle subchannel geometry (except the heated length shorter than a real reactor core), with an additional test section corresponding to a central subchannel heated by 3 and not 4 contributing rods, one rod being replaced by a thimble. The axial power distribution is uniform. The CT scanner gives for every steady run a detailed void fraction array through the measuring section. Results can be compared with CFD simulations or averaged over the cross section for comparisons with 1D module simulation by system codes.

A set of 39 tests is proposed in the benchmark among a large database of 126 tests. The range of flow pressure is from 50 to 170 Bars and the range of mass flux is 500 to 4200 kg mβˆ’2 sβˆ’1. We calculated the whole database and compared the simulation results versus the measurement data.

Calculations were performed using a quasiuniform 31 cell meshing.

The void fraction has been measured by a Ξ³-ray attenuation CT scanner measurement in 1 mm wide beams, giving an array of values throughout the cross section. These values can be integrated in the whole cross section, giving a single void fraction figure associated to an experimental error bar of +/βˆ’3%.

The results are gathered in Figure 2. One can see a good coherence, but yet a slight bias, negative in any series, and a medium dispersion of the predicted values. Some rare points are outside the range +/βˆ’10%. The result statistics is presented in Table 1, given in % of void for the difference β€œpredicted minus measured fraction.”

Some examples of void axial profiles simulated by CATHARE 3 in various flow conditions in the standard central subchannel are shown on Figure 3, as well as the measured value at 1.4 m elevation.

4.2. Exercise I.2 Steady State in Rod Bundles

Several types of rod bundles were tested, most of them including a 5Γ—5 matrix of rods settled with 17 spacer grids of 3 different types, with uniform or cosine power profile and with or without a central thimble instead of a heated rod. The heated length was 3658 mm; the rod diameter and pitch were 9.5 and 12.6 mm (Figure 4).

The heated part of the rod bundle was modelled by a 3D grid, with 6Γ—6 cells in the π‘₯-𝑦 directions at the subchannel scale, and 66 cells in the 𝑧 direction, taking one axial short cell for every 17 spacer grid, with 3 axial cells between two adjacent spacers. The 3 different kinds of spacer grids (Mixing Vane spacer, NonMixing Vane spacer, simple spacer) are described by their porosities, hydraulic diameter, and pressure loss coefficient in every subchannel.

An array of void fraction values in the different subchannels was measured at 3 different levels along the upper part of the heated length, reconstructed by 6 chordal averaged values in π‘₯ and 6 in 𝑦 directions. Here on Figure 5 is shown an example of the sensitivity of the void fraction distribution to the 𝑝𝑖 grad 𝛼 term.

Only the averaged value of the 4 central subchannel void fractions was available in the benchmark database, and this is the compared data versus the void calculated by CATHARE 3 hereafter on Figure 6 for all the tests of the benchmark.

The points are more dispersed than for the single subchannel tests. The statistics of the results (difference: computed minus measured void fraction given in absolute %) is presented in Table 2.

The series number correspond to 3 different bundles; the series 8 is tested with the same bundle as series 5 as repeated cases, which appear to be less satisfying.

4.3. Exercise I.3 Transient in Rod Bundles for Void Fraction Prediction

A set of 12 transient tests is proposed in the benchmark, including power increase (PI), flow reduction, temperature increase (TI), and depressurization in each of the 3 same tested bundles as in the steady tests of the exercise 2. The void fraction was also measured at the same 3 elevations during the transient.

The flow parameters (pressure, flow rate, and inlet temperature) were measured outside of the main vessel, near the inlet nozzle (referred as Coolant Inlet on Figure 7) and should be imposed as boundary conditions in the simulation. As the inlet conditions remain in quasi-incompressible liquid subcooled conditions, pressure and flow values remain more or less uniform all along the inlet pipes and devices, and the parameter variations are not delayed between the inlet nozzle and the bottom of the heated length in the bundle.

Hence, except for the β€œTemperature Increase” transients, as the inlet temperature remains quasiconstant during the transient, the same computational domain was considered as in the previous steady tests, that is, the heated length in the rod bundle only and the requested inlet conditions were applied at the bottom of the heated length.

An example of comparison is given below (Figure 8) for the test 7TPI, in the bundle B7 (with a central thimble and a cosine power profile) and given a linear increase of power, keeping constant the flow rate, pressure, and inlet temperature. The void fraction is slightly underpredicted at the upper location but is satisfactory in front of the 2 lower void measurement elevations. The other tests show less satisfactory results, the upper and medium level void often remaining underpredicted. This behaviour is consistent with the conclusions of the preceding exercise for steady tests, where the higher void fractions are underpredicted while the lower void fractions are more satisfactory.

For this 7TPI test, the location of the inlet temperature measurement is not sensitive because the temperature remains more or less constant during the transient. However, for Temperature Increase tests, this location must be at the boundary of the computational domain. Otherwise, a temperature delay would induce a bias in the simulation. Hence, for these TI transients, another domain has been set up, adding the downcomer as an axial module upstream the 3D module simulating the whole rod bundle.

An example of simulation is presented on Figure 9.

The comparison is less satisfying than for the other test 7TPI. All maximum void fractions are underpredicted and while the time of void take off is well predicted, the time of maximum void is delayed and the void curves seem to be widened. This seems to be the consequence of a too large axial diffusion of void and perhaps also of the temperature step in the inlet part of the domain.

4.4. Exercise I.4 Pressure Drop in Rod Bundle

No data is available for code-to-data comparison for this exercise, except a single value given at the beginning of one transient test (7TPI) in nonboiling steady-state condition. The benchmark specifications recommended the pressure loss coefficients to be used for every type of spacer grid. Using these values, CATHARE 3 predicted the overall bundle pressure drop at 1.85 kg cmβˆ’2, while the measured pressure drop is 1.6 kg cmβˆ’2. This can be considered as a satisfying bias, considering the constant pressure loss coefficient with no dependency on the flow Reynolds number.

5. Departure from Nucleate Boiling in PSBT Phase II Exercises

5.1. Exercise II.1 Steady-State Fluid Temperature in Rod Bundles

This exercise is particularly useful to assess the code capabilities for turbulent dispersion and diffusion in single-phase flow.

In a 5Γ—5 rod bundle featuring a heterogeneous power distribution (Figure 10), a set of 36 thermocouples measure fluid temperature in every subchannel 50 cm above the top of the heated length.

Nine tests at high pressure (from 50 to 170 bars) are proposed for simulations in a wide range of mass fluxes (between 500 and 4700 kg/m2s). In Figure 11 are presented the compared W/E profiles of temperatures, averaged in the N/S direction. The π‘₯-axis numbers correspond to the subchannel columns 1 to 6. The shown temperature gradient is due to the power distribution and is governed by the diffusion and dispersion across the subchannels. The actual profile is more complex due to mixing vanes, which tend to swirl the flow but this effect is not modeled by CATHARE 3.

A first step of analysis allowed us to calibrate the turbulent viscosity used in the liquid energy balance (also in the void mixing term which is not useful in this exercise). Figure 11 shows the sensitivity of the turbulent viscosity value.

The profiles of half of the proposed tests are correctly predicted. The other tests show unbalanced measured temperatures at the outlet compared to inlet flow parameters and bundle power and hence, comparisons are not significant.

In the Table 3, one can see that the parameters of the two correct tests presented on Figure 12 are very close except the inlet temperature (which obviously may have a slight effect on the Reynolds number). However, the temperature profiles show unexplained different behaviours, which are not predicted by CATHARE 3.

Using the available correct tests, a satisfying value of the temperature dispersion parameter has been selected, and implemented in CATHARE 3 for the other void and DNB simulations.

5.2. Exercise II.2 Steady State DNB in Rod Bundles

As for the void fraction tests proposed in phase I, different bundles were tested. The DNB is detected both in experiments and simulation by a significant rise of the wall temperature (more than 11Β°C) when the bundle power is slowly increased.

We calculated 6 test series in 5 different bundles corresponding to several geometries and power profiles, as shown in Table 4. The DNB location is given only in the A4 and A8 bundles. All calculated tests of bundles 4 and 8 were run at the same pressure of 150 bars, while the range in series 0,2,3 spreads from 50 to 170 bars.

Some statistics of the simulation results (relative power: computed over data in % of experimental data) are given in Table 4.

One can see that the result for bundle A0 is satisfying results while the others show a significant bias.

The A3 bundle, featuring a 6Γ—6 rod array, is not better than the A2 bundle, while the side effects are weaker.

The results in the bundle A2 are weakened by two tests at very low flow rate (330 kg mβˆ’2 sβˆ’1), which are overpredicted contrary to the 9 other tests; this enlightens the large value of the standard deviation for this bundle. A similar behaviour exists in series 4 and series 8 where the 2 tests at very low flow rate show significant differences (larger DNB power) compared to the other tests. Generally speaking, the location of the first detected DNB matches better in the series 4 than in the series 8.

The general underprediction of the DNB power in rod bundles may be linked to the use of lookup tables in a 3D analysis; such tables can predict CHF or DNB given 3 parameters: mass flux, pressure, and steam quality. These tables were built using 1D analysis of numerous tests. But in a 3D analysis, the steam quality and the mass flux must obey a local definition with local void fraction and velocities and may display wrong values. As a consequence, the code computation of the local CHF may deviate from the recommended value. Better results can be expected when this point is improved.

Moreover, the better results for the A0 bundle seem point out that the CHF is better predicted with 13 spacers than with 17. The number of 13 is closer to the usual number of spacers in industrial bundles and CHF experiments used to build up the lookup tables. So, the DNB predictions would depend strongly on the spacer number, either through the CHF calculated with the lookup tables, or through the mixing effects simulated in the 3D computation. From this point of view, CATHARE 3 is not able now to predict the DNB power at subchannel scale with accuracy better than 20% in a new bundle.

Given the underprediction of the high-void fraction noticed in the preceding sections, the predicted DNB power should have been overpredicted because the local CHF increases when the local void fraction decreases. This also shows that these lookup tables are not convenient for the analysis of the PSBT tests. It has to be noted that the original purpose of the 3D module of CATHARE2 and 3 was an improvement of the code behaviour in very large volumes within the reactor vessel such as the downcomer and the lower plenum. Applying this module in core subassemblies at subchannel scale is beyond the usual scope of the module, and the BFBT and PSBT benchmarks were just opportunities to check the module capabilities.

5.3. Exercise II.3 Transient DNB in Rod Bundles

Several tests were proposed in the benchmark specification report, in the rod bundles A4 and A8 (see Table 4 for details). As an example, the Temperature Increase transients prescribed a linear increase of the inlet temperature while the three other flow parameters (flow rate, pressure, and power) remained unchanged. The simulation results show the same behaviour as in the steady tests: the DNB power occurred at 86% of the experimental value. This behaviour is also seen in other transients: Flow Reduction, Power Increase, and Depressurization.

6. Conclusion

The 1D and 3D modules of the CATHARE 3 system code were used for the simulations of PSBT benchmark tests. Results of void fraction in phase I and temperature measurements and DNB power measurements in phase II have been compared to calculation results. The void comparisons show that our models of wall and interfacial friction, coupled with void dispersion, lead to satisfactory results, with a slight bias towards void underprediction, for both single subchannels and full rod bundles.

The exercise 1 of phase II, devoted to single-phase mixing and cross flows in liquid phase, shows good results as far as the experimental heat balance of the tests remains satisfactory.

In the exercise 2 of benchmark phase II, steady DNB simulations in 5 different rod bundles show significant underprediction of the critical power in the whole bundle (20% bias). The main reason should be a poor local CHF assessment, more than a rough mixing model. The analysis of transient tests for exercise 3 confirmed this behaviour. The CHF lookup tables in CATHARE 3 underpredict the CHF values in the PSBT rod bundles, which have 17 spacer grids. Improving these features would require improving the modelling of the interaction between spacer grids, turbulence, and local CHF or using CHF correlations designed for specific bundles.

Nomenclature

Roman Letters
π‘βˆΆWall friction phase multiplier
𝑒:Internal energy
𝑓:Wall friction coefficient
𝑔:Gravity
𝐻:Enthalpy
π‘˜:Turbulent kinetic energy
𝐾:Interfacial friction coefficient
β„“:Laplace length
𝑃:Pressure
𝑝𝑖:Void dispersion coefficient
Pr:Prandtl number
π‘ž:Heat flux
𝑑:Time
𝑇:Temperature
𝑉:Velocity
Greek Letters
𝛼:Phase volume fraction
Ξ“:Boiling/condensation rate
πœ’:Friction area (or heated area) over control volume
πœ†:Molecular heat conductivity
πœ‡:Dynamic viscosity
𝜈:Kinematic viscosity
𝜌:Phase density
𝜎:Surface tension
πœπ‘–:Interfacial stress
πœπ‘:Wall stress
Indexes
𝑖:Interface
π‘˜:Any phase
𝑙:Liquid phase
𝑝:Wall
𝑑:Turbulent.

Acknowledgments

This work was made in the frame of the development of the NEPTUNE project, which is jointly developed by the Commissariat Γ  l’Energie Atomique et aux Γ©nergies alternatives (CEA, France) and ElectricitΓ© de France (EDF) and also supported by the Institut de Radioprotection et de SΓ»retΓ© NuclΓ©aire (IRSN, France) and AREVA-NP.