Abstract
This paper deals with the validation of the two-phase flow models of the CFD code NEPTUNEC-CFD using experimental data provided by the OECD BWR BFBT and PSBT Benchmark. Since the two-phase models of CFD codes are extensively being improved, the validation is a key step for the acceptability of such codes. The validation work is performed in the frame of the European NURISP Project and it was focused on the steady state and transient void fraction tests. The influence of different NEPTUNE-CFD model parameters on the void fraction prediction is investigated and discussed in detail. Due to the coupling of heat conduction solver SYRTHES with NEPTUNE-CFD, the description of the coupled fluid dynamics and heat transfer between the fuel rod and the fluid is improved significantly. The averaged void fraction predicted by NEPTUNE-CFD for selected PSBT and BFBT tests is in good agreement with the experimental data. Finally, areas for future improvements of the NEPTUNE-CFD code were identified, too.
1. Introduction
The validation of the two-phase flow modelling capability of CFD codes, for example, NEPTUNE-CFD, is mandatory for its application in the design and safety evaluation of energy systems. The goal thereby is to demonstrate that the two-phase flow models of CFD codes are able to predict the most relevant flow regimes under pre- and postcritical heat flux (CHF) conditions. The focus is on the accurate prediction of the pressure drop, void fraction, critical power, departure from nucleated boiling (DNB), and so forth. The ways how the CFD codes are modelling the heat transfer between a solid and the coolant (wall heat transfer) and the liquid-vapour interphase heat transfer (bulk heat transfer) differ from code to code. The validation of the two-phase flow heat transfer models of NEPTUNE-CFD requires the understanding of the implemented mathematical-physical models in the code as well as their interaction with the heat conduction models of the solids and the fluid dynamics. In the medium term a combined application of CFD and system or subchannel codes will lead to a more realistic prediction of safety-relevant phenomena in nuclear reactors. The simulations performed during this work are focused on the void fraction prediction in rod bundles of light water reactors (LWR) using the NEPTUNE-CFD 1.0.8. In this paper, the investigations performed to validate the two-phase flow models of NEPTUNE-CFD under steady-state and transient conditions using the database provided by the PSBT and BFBT tests are presented and discussed. This research code is being developed and tested within the European NURISP and NURESAFE projects. First of all, the main features of NEPTUNE-CFD are presented in detail. Then, the validation work based on both PSBT and BFBT is discussed in detail. Finally, a summary and the main conclusions are given.
2. The NEPTUNE-CFD Code
The NEPTUNE-CFD solver is based on a pressure correction approach to simulate multicomponent multiphase flows by solving a set of three balance equations for each field (fluid and/or gas phase) in a Reynolds Averaged Navier Stokes (RANS) approach. These fields can represent many kinds of multiphase flows; among them also is bubbly flow. The solver is based on a finite volume discretization together with a collocated arrangement for all variables. The two-fluid models of the NEPTUNE-CFD are designed specifically for the simulation of two-phase transients in nuclear reactors. The models and closure laws used by NEPTUNE-CFD [1] presented in this paper are focused only on the description of the boiling phenomena and heat flux partitioning at the wall under boiling conditions.
When an averaging operation is performed, the major part of the local information at the interfaces and the physics governing the different types of exchanges at a microscale are lost. As a consequence a number of closure relations (also called constitutive relations) must be supplied to close the balance equations so that they can be mathematically solved. Three different types of closure relations can be distinguished: interfacial mass and heat transfer terms (i.e., the molecular and turbulent transfer terms) and the wall heat transfer terms. In the next paragraphs, a short description of the interfacial and wall transfer terms important for the description for the boiling phenomena will be presented. Details about the NEPTUNE-CFD models can be found in the theory and user manual [2, 3].
2.1. Heat and Mass Transfer
If the mechanical terms are neglected in comparison to the thermal terms in the averaged form of the energy jump condition, this condition reduces to This important relation (together with the mass jump condition: ) allows computing the mass transfer terms as functions of the interfacial heat transfer terms and the interfacial-averaged enthalpies : Since no information about the dependence of the interfacial-averaged enthalpies is known, two basic assumptions can be made: either the interfacial-averaged enthalpies are identified to the phase-averaged ones or the interfacial-averaged enthalpies are given by the saturation enthalpies. In NEPTUNE-CDF the first assumption is made. Each interfacial heat transfer term is the product of the interfacial area concentration by the interfacial heat flux density (): The interfacial heat flux density can be defined as (3), where denotes a heat transfer coefficient, the average temperature of phase , and the saturation temperature. The interfacial area concentration is expressed as , where is the void fraction, is the Sauter mean bubble diameter (SMD), and is the thermal conductivity of the liquid. Depending on the Jakob number given by (4), there are two possible scenarios: condensation () or evaporation (): The thermal capacity of the liquid is , and is the latent heat of vaporization. In case of condensation, the Nusselt number is where is the bubble Reynolds number, is the liquid Prandtl number, is the liquid kinematic viscosity, and is the thermal liquid diffusivity. In case of evaporation the Nusselt number is defined as follows: where Pe is the Péclet number and it is defined as follows: The heat transfer term between the vapour and the interface for the case of bubbles is written as where is the gas heat capacity at constant pressure and is a characteristic time given by the users. This relation simply ensures that the vapour temperature remains very close to the saturation temperature , which is the expected result for bubbly flows with sufficiently small bubbles (flow in a PWR core in conditions close to nominal ones).
2.2. Interfacial Momentum Transfer
The interfacial transfer of momentum appears in the balance momentum equation as a source term. It is assumed to be the sum of five forces: The five terms are the averaged drag, added mass, lift force, turbulent dispersion force, and wall lubrication force per unit volume. The drag force definition used in the case of bubbly flow is the one given by the Ishii and Zuber correlation [4], where the calculation of the drag coefficient is based on the local flow regime. The added-mass force takes into account the effect of the bubbles concentration according to Zuber [5] and Ishii [6]. The lift force describes the particular case of a weakly rotational flow around a spherical bubble in the limit of infinite Reynolds number according to Auton [7]. It has been empirically modelled by Tomiyama et al. [8]. The turbulent dispersion force tends to move the bubbles to locations with less void density; it is correlated by Lance and Lopez de Bertodano [9].
2.3. Interfacial Area Modelling
The interfacial area concentration (IAC) transport equation is given by (10) according to Yao and Morel [10]. This formulation is based on the model of Wu et al. [11]: The advantage of (10) is that the nucleation phenomenon clearly appears as a single term given by the product of the bubble number density source term with the surface of a nucleated bubble. If for this particular term the SMD () is replaced by the bubble detachment diameter (for wall nucleation), the IAC is changed accordingly due to the fact that the newly nucleated bubbles often have a smaller diameter than the SMD. The bubble detachment diameter is calculated using Unal's correlation ((17) and (15)). Closure relations must be proposed for the bubble number density source terms , coalescence , and breakup . An example of such models is proposed by Yao and Morel [10].
2.4. Wall Boiling Model
The nucleate boiling term () appears in (11) as a wall-to-fluid heat transfer term per unit volume and time. It is assumed that all applied heat is transferred to the liquid phase. Hence the contribution to the vapour () is zero. To model the wall-to-fluid heat transfer at nucleate boiling, a two-step approach is used. The two steps include the calculation of the condition for boiling incipience in terms of critical wall superheat according to the Hsu criterion [12] and calculation of heat flux partitioning. Following the analysis of Kurul and Podowski [13], the wall heat flux is split into three terms:(i)a single phase flow convective heat flux at the fraction of the wall area unaffected by the presence of bubbles;(ii)a quenching heat flux where bubbles departure bring cold water in contact with the wall periodically;(iii)evaporation heat flux needed to generate the vapour phase.Each of these three phenomena is expressed by a heat flux density (per unit surface of the heated wall) which is related to the volumetric heat flux by the following relation: where is the heated wall surface in contact to the cell having volume ; therefore, is expressed in W/m3 and as well as , , and are expressed in W/m2. The quantities , , and denote the heat flux densities due to liquid convective heat transfer, quenching, and evaporation. The liquid convective heat transfer per unit surface of the heated wall is written as where is the wall temperature and is a heat exchange coefficient which is given by in which is the wall friction velocity and is the nondimensional liquid temperature. The velocity is calculated from the logarithmic law of liquid velocity in the wall boundary layer. The nondimensional temperature follows a similar logarithmic profile.
The heat flux density due to quenching is written as where is the wall fraction occupied by bubble nucleation, is the bubble detachment frequency, is the quenching time, and al is the liquid thermal diffusivity. The two fractions and are given by Here is the active nucleation sites density (per unit surface of the heated wall) and is the bubble detachment diameter. The active nucleation sites density is modelled according to Kurul and Podowski [13], as a function of the wall superheating. The bubble detachment diameter is given by the correlation from nal [14]. This correlation is valid for subcooled liquid but has been extended to saturated liquid. The bubble detachment diameter is given by where is the pressure and , , and are given by the following relations: Here and denote the wall conductivity and thermal diffusivity, respectively, specifies the vapour density, and is the latent heat of vaporization. In the modified correlation is given by where is the norm of the liquid velocity and St is the Stanton number which is defined by The quantity appearing in (17) is given by The quenching time and the bubble detachment frequency are modelled as The third heat flux density used for evaporation is given by To ensure a grid independent solution, the liquid temperature in the wall boiling equations is calculated from the logarithmic temperature profile in a given nondimensional distance from the wall at . This solution is proposed by Egorov and Menter [15]. The reason is that at the centre of the wall-adjacent cell high temperature and void fraction gradients are expected which are strongly dependent on the nodalisation. Taking into account the self-similarity, the non-dimensional temperature profile in the wall boundary layer at y+ = 250 reads where the subscript WC denotes the wall-adjacent cell. This approach is valid only if the wall-adjacent cells remain in log region of the wall boundary layer ().
2.5. Boiling Model Extension for DNB Modelling
The basic wall heat flux partitioning model assumes that the amount of water on the wall is sufficient to remove heat from the wall to be used for evaporation. Superheating of the vapour that occurs at high void fractions is not modelled. Under these assumptions, the basic heat flux partitioning model cannot be used for CHF conditions. In order to take into account the phenomenon of temperature excursion at DNB conditions, the heat flux partitioning model can be generalized as follows: The fourth part of the wall heat flux, , is the diffusive heat flux used to preheat the gas phase: where is the wall heat transfer coefficient calculated from the temperature wall function for the vapour phase and is the vapour temperature at the centre of the wall-adjacent cell. is the phenomenological function, which depends on the liquid volume fraction and takes care for the numerically smooth transition between nucleate boiling regime and CHF regime. The generalized model assumes function in the following form: The extension of the wall-heat-flux-partitioning model was used to take into account the CHF condition. The local void fraction equal to 0.8 can be used as a criterion for the CHF occurrence. This value is close to the Weisman and Pei DNB criterion with the void fraction equal to 0.82 according to [16].
3. Contribution to the NEPTUNE-CFD Validation
The NEPTUNE-CFD code with the two-phase heat transfer models explained in Section 2 is used for the posttest analysis of the PWR and BWR-relevant bundle tests: PSBT and BFBT. The PSBT [17] and BFBT [18] provide void fraction measurements for conditions that are representative for LWR.
3.1. OECD/NRC NUPEC PWR PSBT Benchmark
3.1.1. Scope and Description of the Benchmark
The test bench of the PSBT experiment is shown in Figure 1. The geometry described is a single centred isolated subchannel with a pitch of 8 mm between the walls of the rods. The subchannel has 4 electrically heated walls with a heated axial length of 1.55 m and a uniform axial power distribution. Averaged void fraction data is provided over a cross-sectional area located at 1.4 m distance from the bottom of the heated section. The experimental data are collected by X-ray densitometer. For the comparison, six different experiments have been selected from the PSBT database for the simulation with NEPTUNE-CFD. These experiments belong to Exercise 1 (steady-state single subchannel benchmark) from Phase 1 (void distribution Benchmark). The main parameters of these tests such as power, inlet temperature, pressure, mass flow, and subcooled liquid temperature are specified in Table 1. The NEPTUNE-CFD code lacks a steady-state algorithm; thus a null transient is performed for these simulations.
3.1.2. Applied NEPTUNE-CFD Models
The turbulent transfer term applied is the - model for the liquid phase and the local equilibrium turbulence model [2] for the dispersed phase. A standard wall function adapted for two-phase flow is applied. The steam/water properties are based on the International Association for the Properties of Water and Steam (IAPWS) data [19]. The steam is set at saturation temperature, and it has a slip condition at the wall. The drag and nondrag forces explained in the Section 2.2 are applied to the simulation. The wall lubrication force is not estimated. In case of boiling at a wall, the model can consider an extra artificial quenching flux. Superficial heat flux has been imposed as a boundary condition at the heated walls. As initial conditions, a water mass flow rate at the inlet and a constant pressure at the outlet are imposed.
3.1.3. Space Discretization of the Studied Domain
The cross section of the subchannel geometry is illustrated in Figure 2. The design of the space discretization was performed with the tool ICEM, which is included as a part of the ANSYS-CFX package. The cells in the near wall region are thinner in order to describe the velocity and temperature gradients accurately. At those locations where the velocity is not expected to have steep gradients, the grid is coarser, for example, at the centre of the subchannel.
To catch the physical phenomena near the wall, for example, the void fraction by the numerical codes, a more refined discretization is necessary. Hence, four different spatial resolutions of the subchannel were tested, all of them consist of structured meshes in order to avoid diffusivity problems and to reduce the number of cells; see Figure 2. The nodalisations have the same number of axial levels (176). Geometrical aspects of the spatial discretization are summarized in Table 2.
3.1.4. Mesh Sensitivity Analysis
The case PSBT 1.6222 has been simulated using different meshes (M1 and M4). In Figure 3 the axial void fraction profile in the near wall region predicted by NEPTUNE-CFD for the two meshes is compared to each other. In Figure 4, the axial coolant temperature profiles predicted near the wall (T.near wall) and at the centre of the subchannel (T.center) for M1 and M4 are compared to each other. The largest difference can be observed for the void generation near the wall, where refined meshes generate more steam than coarse meshes. NEPTUNE-CFD predicts a fast axial void fraction increase using M1 (fine meshing) than using M4 (coarse meshing). For M1, the maximal void fraction is calculated at 0.55 m, while for M4 this value shifts to 1 m elevation.
The heat flux partitioning at the wall is commonly divided into 3 fluxes according to Kurul and Podowski [13]. This model releases the heat flux from the heated wall into the water phase, and then it is subdivided into evaporation, quenching, and convective flux. Hence, the water phase receives all the heat fluxes. In the refined mesh (M1) this is problematic for the water, which cannot dissipate the heat fast enough; hence some small peaks appear in the water temperature once the cell has been filled with vapour (Figure 4). This phenomenon is reproduced also in the coarse mesh but is not that severe. To avoid it, the wall heat flux is solved according to the four-flux model described in Section 2. The heat flux is redirected to the steam following a CHF criterion based on the void concentration on the near wall region. This measure can relax the water temperature, but if the heat flux is high enough, the problem is reproduced in the steam temperature. The run number 1.2211 has been selected to test the influence of the four different meshes on the NEPTUNE-CFD predictions. A constant bubble diameter of 0.1 mm has been assumed for all mesh analysis. The location selected to visualize the local values of temperatures and void fraction (VF) is illustrated in Figure 1(b); this is a 8.32 mm distance between heated walls at m. The temperatures near the wall reach saturation for all cases as illustrated in Figure 5. A coarse nodalisation generates flatter profiles than finer meshes, especially regarding the void distribution calculation as shown in Figure 6. There, the steam production is higher for the smaller cells.
For the comparison of four different meshes, NEPTUNE-CFD predicts large gradients of void fraction in the near wall region for the refined meshes (M1) and (M2). These void fraction differences affect the calculated average void fraction over the cross-sectional area at the measurement position. The water temperature evolution calculated for the different mesh configurations at two axial locations can be observed in Figure 7. The first location is at centre of the subchannel, the second is in near wall heated region (first cell near the wall). At the centre of the subchannel, there are no big differences in the temperatures; only close to the outlet there are some deviations. In the near heated wall region locally the difference can reach 6 degrees K between refined mesh (M1) and coarse mesh (M4).
The liquid temperature profile is the combination of several phenomena. The liquid phase near the wall is heated by the wall heat flux, which is divided into convective, evaporation, and quenching heat flux. Once the bubbles are generated, they migrate and condense within subcooled liquid in the core of the flow and hence heat the liquid. The molecular and turbulent heat fluxes inside the liquid phase also modify the temperature profile.
The results computed by NEPTUNE-CFD are summarized in Table 3. The experimental value for the case 1.2211 is 3.8%. The first three nodalisations overpredict the averaged void fraction and the last one slightly underpredicts it. Nevertheless, this experimental value has an error of 3% in the VF measured and is not clear which discretization provides better results. The computational time for the coarse M4 mesh is around two hours; by applying more refined nodalisation, this time increases, reaching around seven hours of calculation for the mesh M1. Some of the cases selected to be simulated have rather high VF concentration at the measurement location; for example, case 1.4326 has a 53.1% VF. Case 1.6222 has a 30.6% VF concentration, and as was illustrated in Figure 3, the first cell near the heated wall is quickly filled with steam. This amount of gas leads to overheating problems in the water and steam phase. The coarse mesh M4 can shift this negative effect in the water and steam temperatures to higher locations and in some cases avoid them.
The choice of the coarse mesh to perform the rest of the simulations has been made to preserve the numerical stability while solving the heat transfer problem. In addition, the coarse mesh (M4) is producing maximum values around 300, which makes it still valid for the application of the selected - turbulence model.
3.1.5. Bubble Size Sensitivity Analysis
In the previous simulations, a constant bubble diameter (0.1 mm) has been applied. Incidence of other bubble diameters or the selection of an interfacial area equation (IAE) for the simulation is discussed in this subchapter. For the mesh M4 and the case 1.2211 previously studied, three different configurations for the IAC are considered. The first is by applying the previous 0.1 mm constant diameter. In the second configuration the diameter is increased to 0.2 mm. The third applies one IEA, described by (10). The bubble diameter influences the simulation results. Therefore, the VF at the measurement location varies depending on this parameter and different results can be obtained. When the code especially is dealing with simulations where the main phenomena are subcooled boiling and low VF are expected.
For the three configurations proposed, the water temperature in the near wall region is illustrated in Figure 8. Here, rather good agreement between cases can be observed. The local bubble size calculated by the IAE is shown by Figure 9. The SMD calculated in this case is much lower, around 0.07 and 0.03 mm, compared with the other cases (0.1 and 0.2 mm). By calculating this small bubble size with an IAE, the IAC is higher compared to the other cases. As a consequence, the heat and mass transfer is higher and there are two important effects. First, the condensation into subcooled liquid is stronger, decreasing the VF within subcooled regions. Second, the boiling in the superheated near wall region is higher. These effects are illustrated in Figure 10. Here, the VF profile illustrates that steam bubbles are nucleated at the heated wall surface and condense in the subcooled liquid in the core of the flow. Larger bubbles (0.2 mm) do not condense so easily leading to a higher VF concentration in the bulk flow region. Concerning the impact of the bubble size in the steam velocity (Figure 11), slightly higher velocities are registered for the steam in the case of 0.2 mm bubble diameter in the bulk region. This can be explained as a consequence of a higher VF concentration in this region.
In Figure 12 the axial VF evolution in the centre of the subchannel is shown. In the centre of the subchannel the VF generated by the 0.2 mm bubbles is clearly higher due to the reasons previously explained. The void calculated by the IAE with bubbles of 0.03 mm of diameter is lower due to the strong condensation into subcooled liquid.
The pressure drop in the channel is calculated also for the different bubble sizes. The simulated results are presented in Figure 13. Only the case predicted with the IAE exhibits higher values. In this case the concentration of bubbles in the near wall region is bigger leading to an increment of the pressure drop.
By using an IAE there are more parameters to control. One of the most important measures is to clip the value of the minimum bubble diameter. Very small bubble size can lead to numerical instabilities regarding the nondrag forces applied like the added mass force or the turbulent dispersion force. For the simulations the smallest bubble diameter of 0.01 mm is allowed in the computational domain.
3.1.6. Selected Results
The VF calculated by NEPTUNE-CFD for the six cases selected from the PSBT database is illustrated in Table 4 together with the experimental measurements. These simulations have been performed with a single IAE and the mesh M4.
Three experiments are overpredicted (1.2223, 1.2237, and 1.2211) and the other three are underpredicted (1.4326, 1.4325, and 1.6222) by NEPTUNE-CFD. In Table 4 the relative error between the experimental data and the computed results is shown. The simulations performed show a variation of the relative error from −34% to 36% compared with the measured PSBT data. It must be noted that the two-phase flow models of CFD codes are still under development, and hence the deviations of the predictions from the measured data are still considerable. In Table 4, the predicted void fraction by ANSYS CFX 12.1 using the shear stress turbulence model (SST) and the respective deviations from the experimental data are also given. These values were taken from [20] and they confirm the trends predicted by NEPTUNE-CFD that the deviations may be huge depending on the modelling being used to describe the turbulence effects among others.
The axial VF profile for each case is illustrated in Figure 14; the experimental measure and its measurement error (3%) are also included.
3.2. The NUPEC BWR Full Size Fine Mesh Bundle Test (BFBT) Benchmark
3.2.1. Scope and Description of the Benchmark
The BFBT void distribution benchmark [18] was made available by the Nuclear Power Engineering Corporation (NUPEC). It is one of the most valuable databases identified for thermal hydraulic modelling. The NUPEC database includes subchannel void fraction and critical power measurements in a representative (full-scale) BWR fuel assembly (FA). The high resolution and high quality of subchannel void fraction data encourage advancement in understanding and modelling of complex two-phase flow in real bundles and make BFBT experiments valuable for the NEPTUNE-CFD multiphase models validation. The benchmark consists of two parts: void distribution benchmark (Phase I) and critical power benchmark (Phase II). Each part has different exercises including simulations of steady-state and transient tests as well as uncertainty analysis. An exercise from phase I has been selected for the validation of NEPTUNE-CFD. The transient tests performed in the frame of this benchmark represent the thermal hydraulic conditions that may be encountered during a postulated BWR turbine trip transient without bypass and a recirculation pump trip. From this postulated transient scenarios important thermal hydraulic parameters are derived for the test, such as the evolution of the pressure, total bundle power, mass flow, and radial and axial power profile, which serves as initial and boundary conditions for the CFD simulations.
The test section of the experiment is a full sized BWR FA (Figure 15), with sixty electrically heated rods (12.3 mm diameter and 16.2 mm rod pitch) and one water channel (34 mm diameter). The electrically heated section is 3708 mm high; the heaters are surrounded by an insulator (boron nitride) and by the cladding (Inconel 600). An X-ray densitometer is used to measure the averaged void fraction at three different axial levels from the bottom from the heated section ( m, m, and m). What the BFBT experimental data provides is the evolution in time of the averaged VF at the three axial levels mentioned.
To reproduce in the test bench the turbine trip without bypass and the recirculation pump trip, the BFBT database provides the evolution of the water mass flow rate, the system pressure, and the power. This data is used as boundary condition for NEPTUNE-CFD. In Figures 16, 17, and 18 the evolution of the outlet pressure, the mass flow rate, and the power during the 60 seconds transient is given for both scenarios. The experiment was performed with a uniform axial power shape. The radial power shape is described by Figure 19. The water inlet temperature remains constant at 552°K during the experiment.
3.2.2. Applied NEPTUNE-CFD Models
The numerical simulation is performed with the models explained in Section 2; - turbulence model for the liquid phase is used with the two-phase modified wall function. Local equilibrium turbulence model [2] is considered for the dispersed phase. The libraries for the liquid properties are provided by the IAPWS data [19]; furthermore steam close to saturation condition is assumed. For heat transfer through the gas/liquid interface, a thermal phase change model is applied. The heaters are modelled in terms of a heat flux boundary condition. At the inlet, a mass flow rate is set and the pressure is imposed at the outlet. The temporal evolution of mass flow rate, power, and pressure are taken from the experimental conditions (Figures 16, 17, and 18).
The time step is adaptive depending on the Courant number. The time step width ranges from 1 to 3 ms. The drag and nondrag forces, lift, added mass, and turbulent dispersion force, are computed for the simulation.
3.2.3. Modelling of the Rod Bundle
The nodalisation of the PSBT test section was done following the best practice guidelines for the use of CFD codes. The nearest wall heated region cell has a constant width of 0.3 mm. The mesh is composed of 135 axial levels and 12 cross cells in each subchannel. Globally the NEPTUNE-CFD nodalisation has 211928 cells. The maximum values are located close to the outlet. They oscillate between 300 and 400 depending on the transient conditions. These high values correspond to the high velocities locally achieved by the steam.
Taking into account the radial symmetry of the FA only 1/8 of the fluid domain is modelled with two symmetry planes (Figure 15). This decision is taken although being aware that the radial power distribution has no axial symmetry according to Figure 19, but larger models penalize the computational time. For the initialization of the simulations, the power is increased gradually from 0 to nominal values which are reached after 6 seconds. The simulations are conditioned by a large amount of void present in the domain. Numerical simulations find convergence problems in very refined meshes with boiling flow, especially with respect to the water and steam temperature close to the heated wall region. To set the heat flux defined by the experimental measurements (Figure 18) two different configurations are investigated. The first option is to locate the flux at the rod wall of the fluid domain using NEPTUNE-CFD as standalone, Figure 20(a). The second option is to model the clad and insulator of the test bench, where the heat flux is placed at the inner diameter of the insulator’s surface, Figure 20(b). In both cases the power is set as superficial heat flux (W/m2). NEPTUNE-CFD is not able to solve the thermal problem for the solid domain by itself. It has to be coupled with the heat conduction tool SYRTHES [21].
The averaged void fraction for the cross-sectional area of the FA is calculated each 0.02 seconds at the three axial levels ( m, m, and m) from the bottom of the heated section. Figure 15 shows the location of the measured sections. The original experimental measurements were considered not to be accurate enough and a correction factor has been proposed to be applied; see [22]. is the original experimental data in percentage of void, and is the resulting data after the application of the correction factor. The new corrected data is obtained according to the following:
3.2.4. Selected Results for the Turbine Trip without Bypass Experiment
Computed void fraction results corresponding to two configurations are illustrated. The first is by modelling the problem with NEPTUNE-CFD standing alone and the second configuration takes into account the wall thermal inertia effect by adding the SYRTHES code to NEPTUNE-CFD.
Figure 21 illustrates with different figures the local VF distribution calculated by NEPTUNE-CFD and SYRTHES in the domain for different time steps during the transient: at second 7, at second 11.5 (during the maximum void concentration), at second 20, at second 30, and at second 50, where the combination of power and mass flow decreases the void to the minimum values of the simulation. Furthermore the location of the measured axial levels is shown.
(a)
(b)
(c)
(d)
(e)
(f)
The comparison between the evolution of the calculated void averaged calculatedby NEPTUNE-CFD and SYRTHES and the experimental data is shown in Figure 22 for the three different axial levels; the experimental measurement error of the mean averaged VF (±2%) is also included.
At this point, it is important to remark that the models implemented regarding interfacial interactions like mass, heat transfer, and momentum are designed for bubbly flow only. The water is set as continuous phase and the vapour is set as the disperse phase; this condition is valid for bubbly flow. When the amount of void in the fluid exceeds 60%, the regime is no longer bubbly flow. Therefore, the results of the second and third elevation cannot be considered reliable and more developments in the modelling are necessary to properly describe the flow at those elevations. The description of different flow regimes in the same domain requires the definition of a change of continuous and dispersed phase at each location. Nevertheless the code is able to trace the tendency of the void generation in line with the experimental data even in those flow regimes.
Figure 23 illustrates the relative error of the computed void fraction against the measurements for the first axial level ( m) together with the normalized values of the pressure mass flow rate and power of the transient. In Figure 23, it can be observed how the code struggles especially when the mass flow rate decreases. The code has a constant relative error around −10% during the rest of the transient. At the end of the transient, the predicted void fraction agrees very well with the data at all three axial elevations.
The coupling of NEPTUNE-CFD with SYRTHES leads to a significant improvement due to the accurate steam temperature calculation. In addition SYRTHES contributes to relax the temperature calculation and the solid-liquid interface increasing the time step and accelerating the simulation.
Figure 24 illustrates the evolution during the turbine trip simulation of the local temperatures of the water and steam together with the saturation and clad temperature.
The temperatures from Figure 24 are taken from one point near the wall region of the rod 4 at the end of the heated section, where the temperatures are expected to be the highest in the domain. By applying NEPTUNE-CFD and SYRTHES during the power peak, the steam temperature remains as maximum at 3 degrees oversaturation, while in other simulations without solving the thermal wall problem (NEPTUNE-CFD standalone), this temperature can reach locally hundreds of degrees oversaturation which is physically not correct. On the other hand if the heat flux is defined at the insulator, the water and steam temperatures are calculated in a better way since SYRTHES is providing the wall temperature. Thus, the heat flux at the solid-fluid interface is calculated according to the following considering the wall and fluid temperature (, ) and a heat transfer coefficient: where is the wall friction velocity and is the nondimensional liquid temperature. Hence the heat flux from the wall to the steam phase is regulated avoiding excessive overheating. The wall temperature increases at those locations with high void fraction. In Figure 24, the local temperature jumps between the wall and the liquid oscillate from 20 to 30 degrees K during the transient, and when the power peak occurs, this difference rises up to 90 degrees K.
3.2.5. Selected Results for the Recirculation Pump Trip Experiment
Using the experience acquired from the turbine trip simulation, another exercise of the BFBT database is simulated. For this case, only the simulation combining NEPTUNE-CFD and SYRTHES is performed.
The comparison between the experimental data provided by the BFBT database and the computed void fraction of NEPTUNE-CFD/SYRTHES is illustrated in Figure 25.
For this case the prediction fits well with the experimental data during the void increase between seconds 11 and 13. At this time there is no power peak and the void rises due to the decrease of mass flow. At this point the simulation is underestimating the experimental data only at the first axial level. For the lapse of time between seconds 15 and 42, an overprediction occurs for all axial levels. At the last third of the transient when the mass flow increases up to nominal conditions, the void predicted is slightly underestimated for axial levels 1 ( m) and 2 ( m), while for axial level 3 ( m) a good agreement is achieved.
Figure 26 shows the local liquid temperature distribution for the 3 axial levels of the measurements at second 30 of the recirculation pump trip simulation. For the first axial level, the temperature is 3 degrees overheated in the near wall region and subcooled in the centre of the subchannels. In upper locations the steam concentration in the near wall region is above 80% (Figure 27), and all the heat flux is transferred into the steam relaxing the liquid temperature.
(a)
(b)
(c)
(a)
(b)
(c)
In Figure 28, the local temperature evolution is shown for two locations: the near wall region and the bulk of the fluid in the centre of one subchannel. Both locations are m high from the beginning of the heated section. The largest deviation from saturation takes place for the steam temperature in the near wall region. It can reach locally 17 K oversaturation. This steam temperature has different values by decreasing the time scale returning to saturation (8). This time scale should be lower than the time step to work properly. In this figure a time scale of 0.05 seconds is compared with a time scale of 1 ms. The temperature difference between those two settings is about 10 degrees. The lower time scale helps to maintain the steam temperature close to saturation in the near wall region. The temperature of the liquid in the near wall region is slightly oversaturated and it is not affected by the different time scales of the steam heat transfer coefficient. This parameter is not affecting significantly the steam generation or the temperature of the steam and the liquid at the bulk region which remains at saturation.
4. Summary and Conclusions
The validation basis of NEPTUNE-CFD has been extended by using LWR-relevant transient test data obtained in both a PWR-specific (NUPEC PSBT) and a BWR-specific (NUPEC BFBT) test facility. It is the first time that the two-phase flow models of NEPTUNE-CFD have been validated using transient data obtained for tests representing the LWR plant conditions of postulated transients like a turbine trip or a recirculation pump trip. This validation process has clearly shown the status of the two-phase flow modelling of NEPTUNE-CFD. It is possible to summarize the global conclusions about the NEPTUNE-CFD capabilities, identifying its strengths and weaknesses.
According to the interfacial exchange terms, the main flow regime currently implemented in NEPTUNE-CFD is the bubbly flow, which is good enough to describe the flow with a certain void concentration (<60%). However, in BWR the amount of void generated is beyond the limits of the bubbly flow, and other flow regimes (slug flow or annular flow) must be taken into account for an accurate description of safety-relevant phenomena. How to model the transition between flow regimes describing a whole flow map is an open issue for CFD codes.
A study of the sensitivity to the bubble size has been performed applying two constant bubble sizes (0.1 and 0.2 mm) and a variable size with an IAE for the PSBT experiments. The heat and mass exchange is strongly influenced by the IAC. Bubbles may rise close to the heated wall and eventually depart from it and migrate into subcooled liquid where they condense. Smaller diameters produce more IAC and more condensation. Hence, for smaller selected bubble sizes, reduced void fraction will appear in the domain. The smaller bubble size generated by the IAE produces more interfaces, and thus the condensation is stronger and the void predicted is less than the other two bigger bubbles selected (0.1 and 0.2 mm).
By using an IAE it is assumed that there is a single bubble size per cell. If all bubbles have locally the same size, they will condense at the same speed and their diameter will decrease. If a multisize model is applied, the small bubbles will decrease and collapse rapidly leaving the bigger ones, which increase the mean bubble size. Therefore, the assumption of a single bubble size can lead to an underestimation of the bubble mean size in this case, affecting the void generation.
NEPTUNE-CFD alone cannot solve heat conduction in solid domains. For this reason it is coupled with SYRTHES which is in charge of calculating the temperature in solid domains. The capabilities of the NEPTUNE-CFD coupled with SYRTHES have been demonstrated by the prediction of thermal inertia at the walls. By applying SYRTHES, the simulations of the turbine trip are in better agreement with experimental data compared with the application of NEPTUNE-CFD standalone. In addition, this heat conduction solver helps to reasonably control steam and water temperatures in the near wall region by computing a proper heat transfer coefficient at the solid/fluid interface. The use of SYRTHES has a positive contribution to the NEPTUNE-CFD prediction capabilities. However, it is not parallelized and hence it penalizes the computational time.
Due to the different mesh distribution applied for each code (tetra volumes for SYRTHES and hexa volumes for NEPTUNE-CFD), a careful mesh design is required to match the nodes at the solid/liquid interface and have a good energy balance at the interface.
The classical 3 fluxes formulation for the wall/fluid heat transfer assumes three fluxes: quenching, convection, and evaporation. If boiling occurs in this region and the amount of liquid decreases, the water enthalpy rises and the temperature can be several degrees above saturation. To avoid this problem, the four-flux model decomposition (see (25)) transfers all the heat flux into the steam when the void is above the 80% in the near wall region (condition frequently present near the wall). With this measure, the liquid temperature remains close to saturation, but now the steam can be locally overheated in the near wall region if the heat flux increases. The steam temperature problem can be fixed by increasing the near wall cell size. By this way, the y+ values increase; it is important to maintain those values in a range according to the turbulence model selected. The turbulence model selected in this case is the two-equation based -. The wall-adjacent cells must belong to the log region of the wall boundary layer () for this model to work properly. Other more sophisticated turbulent models have been tested, for instance, a 7-equation model, but no convergence has been obtained, mainly because they require a refined near wall region discretization to reach y+ values below 10.
The near wall cells size results as an agreement between the thermal and momentum modeling. This agreement consists in solving the heat transfer problem properly without penalizing excessively the values.
By setting a very refined axial nodalisation (more than 200 axial levels for the presented domains), the code struggles calculating the pressure field. It is recommended to enlarge the cells in the axial direction in the locations where large amounts of steam are generated.
According to the results obtained during the validation process, the subcooled and saturated boiling description can be considered sufficient enough in its range of applicability. Even if further developments are required, NEPTUNE-CFD has demonstrated to be a valid TH tool for the two-phase flow modeling in LWR applications.
List of Symbols
| Interfacial area concentration (m2/m3) | |
| : | Wall thermal diffusivity (m2/s) |
| : | Liquid thermal diffusivity (m2/s) |
| : | Thermal capacity of the liquid (J/(mol·K)) |
| : | Gas heat capacity at constant pressure (J/(mol·K)) |
| : | Sauter mean bubble diameter (m) |
| : | Bubble detachment diameter (m) |
| : | Bubble detachment frequency (Hz) |
| : | Interfacial-averaged enthalpies (kJ) |
| : | Heat transfer coefficient |
| : | Jakob number |
| : | Thermal conductivity (W/(m·K)) |
| : | Latent heat of vaporization (kJ/kg) |
| : | Interfacial transfer of momentum |
| : | Active nucleation sites density |
| : | Nusselt number |
| Pr: | Liquid Prandtl number |
| : | Péclet number |
| : | Convective heat flux (W/m2) |
| : | Evaporation heat flux (W/m2) |
| : | Quenching heat flux (W/m2) |
| : | Volumetric heat flux (W/m3) |
| : | Interfacial heat flux density |
| : | Bubble Reynolds number |
| : | Stanton number |
| : | Temperature of phase (K) |
| : | Nondimensional liquid temperature |
| : | Quenching time (s) |
| : | Wall friction velocity (m/s) |
| : | Phase velocity (m/s) |
| : | Wall-adjacent cell |
| : | plus value |
| : | Void fraction |
| : | Mass transfer condition (Kg) |
| : | Kinematic viscosity (m2/s) |
| : | Time scale returning to saturation (s) |
| : | Bubble nucleation source terms |
| : | Bubble coalescence source term |
| : | Bubble break-up source term |
| : | Wall conductivity (W/(m·K)) |
| : | Phase density (kg/m3). |
Abbreviations
| BFBT: | Boiling water reactor full bundle test |
| BWR: | Boiling water reactor |
| CEA: | Commissariat de l’Energie Atomique |
| CFD: | Computational fluid dynamics |
| CHF: | Critical heat flux |
| DNB: | Departure from nucleated boiling |
| EDF: | Electricité de France |
| IAE: | Interfacial area equation |
| IAC: | Interfacial area concentration |
| IRSN: | Institut de Radioprotection et de Sûreté Nucléaire |
| IAPWS: | The International Association for the Properties of Water and Steam |
| KIT: | Karlsruhe Institute of Technology |
| NUPEC: | The Nuclear Power Engineering Corporation |
| RANS: | Reynolds Averaged Navier Stokes |
| SMD: | Sauter mean bubble diameter |
| SST: | Shear stress turbulence model |
| OECD: | Organisation for Economic Cooperation and Development |
| PSBT: | Pressurized water reactor subchannel and bundle tests |
| PWR: | Pressurized water reactor |
| VF: | Void fraction. |
Conflict of Interests
The authors hereby declare that no conflict of interests is present between them and the commercial entities mentioned in the context of the paper.
Acknowledgments
The authors thank the Program “Nuclear Safety Research” of KIT for the financial support of the Research Topic “Multiphysics Methodologies for Reactor Dynamics and Safety” and the EU Project NURISP.