CFD Analysis of a Void Distribution Benchmark of the NUPEC PSBT Tests : Model Calibration and Influence of TurbulenceModelling

The paper presents CFD calculations of the void distribution tests of the PSBT benchmark using ANSYS CFX-12.1. First, relevant aspects of the implemented wall boiling model are reviewed highlighting the uncertainties in several model parameters. It is then shown that the measured cross-sectionally averaged values can be reproduced well with a single set of calibrated model parameters for different test cases. For the reproduction of patterns of void distribution cross-sections, attention has to be focussed on the modelling of turbulence in the narrow channel. Only a turbulence model with the capability to resolve turbulent secondary flows is able to reproduce at least qualitatively the observed void distribution patterns.


Introduction
Based on NUPEC PWR Subchannel and Bundle Tests (PSBT), an international benchmark has been promoted by OECD and NRC and coordinated by Penn State University (PSU).The benchmark includes void distribution and departure from nucleate boiling exercises.In the first exercises the void fraction distribution was investigated in a steady state subchannel grade benchmark.In the present paper some of the single-channel steady state void fraction measurements are analysed to investigate the capabilities of present CFD modelling of wall boiling.
For engineering calculations, currently the most widely used CFD approach to model two-phase flows with significant volume fractions of both phases is the Eulerian two-fluid framework of interpenetrating continua (see, e.g., [1][2][3]).In this approach, balance equations for mass, momentum, and energy are written for each phase, that is, gas and liquid, separately and weighted by the so-called volume fraction which represents the ensemble averaged probability of occurrence for each phase at a certain point in time and space.Exchange terms between the phases appear as source/sink terms in the balance equations.These exchange terms consist of analytical or empirical correlations, expressing the interfacial forces, as well as the heat and mass fluxes, as functions of the average flow parameters.Since most of these correlations are highly problem-specific, their range of validity has to be carefully considered and the entire model has to be validated against experiments.
For the case of boiling flows, where heat is transferred into the fluid from a heated wall at such high rates that vapour is generated, additional source terms describing the physics of these processes at the heated wall have to be included.A CFD wall boiling model implemented in CFX following the lines of Kurul and Podowski [4,5] was calibrated and validated by several authors, for example, Krepper et al. [6] against experimental results of Bartolomej and Chanturiya [7].In these tests, subcooled flow boiling of water at high pressure flowing upwards in a vertical pipe heated from the outside was investigated and measurements of the axial development of void fraction, wall temperature and cross sectionally averaged liquid temperature were provided.
The aim of this work is to investigate the applicability of the CFX models to the PSBT tests.In Section 1 the PSBT tests for the void fraction distribution and the parameters of the selected tests are described briefly.In Section 2 the most important model details used for the simulations are described and critically reviewed.In Section 3 the sensitivity of the different parts of the models is investigated in view of the knowledge discussed in Section 2. The parameters are calibrated to only a single test.In the following Section 4 these calibrated models are applied to additional tests without any further change of model parameters.The cross sectionally averaged values of the selected tests are summarized for comparison to the calculations of other participants in the benchmark.In Section 5 the influence of the turbulence model on the void fraction distribution is investigated.

The OECD Benchmark Test
The PWR Subchannel and Bundle Tests (PSBT) were conducted by NUPEC (1987NUPEC ( -1993) ) within an extensive experimental campaign aimed at verifying the reliability of fuel assemblies used for commercial nuclear power plants [8,9].Void fraction measurements and departure from nucleate boiling (DNB) tests were performed under PWR thermalhydraulic conditions including steady states and transients such as power increase, flow reduction, depressurization, and temperature increase.The void fraction in each experiment was measured by gamma-ray transmission.These tests form the basis of an OECD benchmark for CFD and subchannel codes [10].
The subchannel test section, as shown in Figure 1, simulates a single subchannel of a PWR fuel assembly.The effective heated length is 1500 mm where the void measurement section is located near close to the top end at 1400 mm from the bottom of the heated section.
For the analysis different tests in different pressure regions were selected (see Table 1).

Modelling of Boiling at a Heated Wall: The General Model
Structure.In boiling, heat is transported from the hot wall to the fluid by several different mechanisms.On parts of the wall, where no bubbles reside, heat flows directly to the subcooled liquid in the same way as in single-phase flow.On parts of the wall where bubbles grow, heat is consumed by the generation of vapour which occurs at the so-called nucleation sites.Moreover, there is a liquid mixing mechanism due to the bubbles which leave the wall.As a consequence of the recirculation around the detaching bubbles, cold liquid from the bulk of the flow is brought into contact with the hot wall which leads to additional cooling.This mechanism is termed quenching.The total supplied heat flux is accordingly expressed as the sum of three contributions as where Q C , Q Q , and Q E denote the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively.The individual components in this heat flux partitioning are then modelled as functions of the wall temperature and other local flow parameters.Once this is accomplished, (1) can be solved iteratively for the local wall temperature T W , which satisfies the wall heat flux balance.Denoting the fraction of area influenced by the bubbles as A W , the heat flux components are expressed as discussed in the following.The turbulent convection heat flux is calculated in the CFX model version in much the same way as for a pure liquid flow without boiling, but multiplied by the fraction of area unaffected by the bubbles, that is, Here, h C is the heat transfer coefficient which is written using the temperature wall function T + (y + ) known from Kader [11] as where nondimensional variables (indicated by superscript "+") and the friction velocity u τ are defined as usual.Note that ( 2) and ( 3) may be evaluated at any location y, provided it is used consistently whenever a variable depends on position.Q Q is represented in terms of the quenching heat transfer coefficient h Q : A grid-independent solution for Q Q is obtained by evaluating the non-dimensional temperature profile at a fixed value of y + .The evaporation heat flux Q E is obtained via the evaporation mass flux at the wall: where the generated vapour mass ṁW is expressed in terms of the bubble diameter at detachment d W , bubble generation frequency f , and nucleation site density N as In terms of bubble detachment diameter and nucleation site density, the wall area fraction A W influenced by vapour bubbles is given by Here, a is the so-called bubble influence factor, for which a value of 2 is commonly used [4,5].Since A W = 1 corresponds to the case where the whole surface is under the influence of bubbles, A W as calculated by (7) has to be limited to values smaller than one.Moreover, it should be kept in mind that already as A W approaches 1 the assumptions of the model are not really satisfied anymore.
Correlations for the yet undetermined quantities used in the CFX wall boiling model are discussed in the following.

Bubble Size at Detachment.
The bubble size at detachment depends on the liquid subcooling.Also the liquid properties which depend on the system pressure, the flow rate, and the heat flux have an influence.An investigation of the bubble size at detachment was performed by Tolubinsky and Kostanchuk [12] for water at different pressures and subcoolings.The observed dependence on the liquid subcooling at atmospheric pressure can be fitted to a correlation:

Nucleation Site Density.
The situation concerning data on nucleation site density is much less clear.Most of the time, correlations are expressed in the form of power laws depending on the wall superheat as However, a recent compilation [13] shows that vastly different parameter values are required to match different data  sets.A likely reason for this fact is that nucleation site density is highly dependent on the microscale topography of the boiling surface, which in turn depends strongly on the processes that were used to finish the surface.These processes are very diverse and in most boiling experiments not specifically controlled or characterized.
For the tests of Bartolomej and Chanturiya [7] N ref = 0.8 × 10 6 m −2 and ΔT ref N = 10 K were found to yield the best results in the model framework.

Further Model Aspects.
The bubble detachment frequency f is given according to Cole [14] as a function of the detachment size d W .The quenching heat transfer coefficient is calculated as suggested by Mikic and Rohsenow [15].
In the bulk, vapour is assumed to be at saturation condition.Where the liquid is subcooled, that is, T L < T sat , vapour is condensing and calculated with the transfer rate: Here, h LG is the interfacial heat transfer coefficient, calculated according to Ranz and Marshall [16] and A I is the interfacial area.To close the phase transition model in the bulk bubbly flow with a mean bubble diameter d B , Kurul and Podowski [4] and also Anglart et al. [17] proposed to calculate the bubble diameter d B locally as a linear function of liquid subcooling T sub : Reference subcooling conditions for typical nuclear energy applications have been given as d B,1 = 0.1 mm at T sub,1 = −13.5K and d B,2 = 2 mm at T sub,2 = 5 K in [17].
In a two-phase turbulent flow gaseous bubbles have an influence on the liquid turbulence.This effect is described by increasing the turbulent viscosity according to Sato et al.  [18].Conversely the liquid turbulence flow structures influence the gas distribution, which is described by a turbulent dispersion force according to Burns et al. [19].
For momentum exchange between the phases, finally, the Ishii and Zuber [20] drag law was used.Furthermore, a lift force according to Tomiyama et al. [21], a turbulent dispersion force according to Burns et al. [19], and a wall force according to Antal et al. [22] were included.
3.5.Evaluated Quantities.Quantities requested for the benchmark have been evaluated as follows.The mixed density ρ MIX was calculated according to and the mixed enthalpy H MIX according to with the temperature-dependent static enthalpies for liquid H L and gas H G .The quality X was then calculated according to with the enthalpies at saturation for liquid H LSat and gas H GSat .unheated part of 0.5 m length was simulated to model flow development.For the liquid and gas properties the water steam tables according to IAPWS IF97 were applied.In Section 3 the turbulence of the liquid phase was modelled by a Shear Stress Transport (SST) model [23].In Section 4 two different Reynolds Stress Turbulence models were applied for comparison.The influence of the turbulence model is investigated in Section 5.

Calibration of Boiling Model Parameters
4.2.1.Bubble Size at Detachment.Earlier investigations [24] have shown that the bubble size at detachment has a sensitive influence on the calculated gas fraction since it determines the part of heat transferred by evaporation.To demonstrate the influence of this parameter different values d ref (8) were investigated and the resulting cross sectionally averaged gas fraction for the test case 1.4324 is shown in Figure 2. Nevertheless, because of similar conditions as in the Bartolomej and Chanturiya [7] test cases, d ref = 0.6 mm was assumed for the further simulations.

Nucleation Site Density.
The nucleation site density has almost no influence on the liquid temperature, a small influence on the gas void fraction but a strong influence on the wall superheating T W − T sat .The common lack of information on the nucleation site density can be compensated by adjusting N ref to match the measured wall temperature.Unfortunately, the wall temperatures are not available in the PSBT data.For the test case 1.4324 different values of the reference nucleation site density N ref (9) were assumed to calculate the circumferentially averaged nucleation site density N (see Figure 3), bubble influence area A W (see Figure 4), wall temperature T W (see Figure 5), and the cross sectionally averaged gas volume fraction (see Figure 6).Figure 7 shows furthermore for N ref = 8 × 10 5 m −2 the distribution of the area fraction occupied by gas A W on the heated wall.For axial heights > 1 m (Figure 4) and for different areas (Figure 7) A W > 1 is found which does not seem reasonable.Therefore, N ref was set to a smaller value of 8 × 10 4 m −2 for all simulated cases.For all simulated cases compliance with the condition A W < 1 was checked.

Results
The calculations for the cross-sectional averages listed in Table 2 were performed applying a Shear Stress Transport (SST) model for the liquid phase turbulence.This model blends a k −ε model in the bulk and a k −ω model in the near wall region [23].Furthermore, the corresponding Reynolds Stress model, the Baseline Reynolds Stress model (BSL-RSM) were investigated.This model has the same model structure as the SST model.In some cases the convergence of the BSL-RSM runs was questionable.The k−ω Reynolds Stress model (O-RSM) proved to be numerically more stable.The details of these approaches to turbulence modeling can be found in the ANSYS solver manual.Calculating the cross sectionally averaged values, no significant differences between the different turbulence models could be discerned (examples for selected tests see Figure 8).In Section 5, however, a strong influence of turbulence modelling on the gas fraction distribution is shown.
The experiments do not provide values on bubble sizes or on temperatures of the heated wall.Exploiting this information gap for calibration of the correlations for the bubble size at detachment and the nucleation site density, at least cross-sectionally averaged values for the mixed density and gas fractions can be calculated with good agreement to experimental data.

Influence of the Turbulence Model on the Gas Fraction Distribution
Considering the cross-sectional distribution of the gas fraction, remarkable differences between the calculations using different turbulence models were found.Applying a turbulence model with a scalar turbulent viscosity, like the SSTmodel, the maximum of the gas was found near the heated wall in all cases.Figure 9(a) shows the calculated distribution for run 1.2422.In contrast, the measured gas distributions published in the benchmark specifications show significant amounts of gas also in the core of the channel.This is more similar to the results obtained by a Reynolds Stress model (see Figure 9(b)).The reason for these different patterns is probably that the SST model is not able to capture the secondary flows (see Figure 10).Figure 11 shows the comparison of different runs with increasing inlet temperature (compare Table 1).With increasing inlet temperature the maximum gas fraction is Science and Technology of Nuclear Installations shifted from the corner zones towards the core regions of the channel.The figure clearly shows the influence of the liquid flow field on the gas distribution similar to the measurements.

Conclusion
Current CFD implementations of a wall boiling model are able to calculate the gas fractions in good agreement with measurements.The model review in Section 3 has shown which information is necessary to perform meaningful calculations.With suitably calibrated correlations for the bubble size at detachment and the nucleation site density, at least cross sectionally averaged values for the mixed density and gas fractions can be calculated with good agreement to experimental data.In the corresponding experiments besides the measurement of gas fraction the determination of the temperature at the heated wall is urgent necessarily.Further research to reduce the number of model parameters that have to be determined by model calibration would be highly desirable.The distribution of gas within a cross section has been shown to depend strongly on the modelling of turbulence.Indication has been given that common two-equation models are too simplified for subchannel geometries while anisotropic Reynolds Stress models show a better potential.In order to make a stringent qualification of models, however, much more detailed CFD-grade data are urgently needed.While all of our calculations were performed for steady-state conditions, numerical convergence was a subject requiring special attention which could indicate that the real solution might contain transient elements.
) A least squares fit gives parameter values are d ref = 0.0013 m and ΔT ref d = 53 K under these conditions.To match the tests of Bartolomej and Chanturiya [7] which were conducted at much higher pressures relevant for typical nuclear energy applications the values of d ref and ΔT ref d had to be adjusted to d ref = 0.6 mm and T ref d = 45 K (e.g.,[6]).
Exp d ref = 1.2 mm

Figure 9 :Figure 10 :
Figure 9: Cross sectionally averaged gas fraction at an axial distance of 1.4 m for run 1.2422.