Abstract
Pressurized thermal shock (PTS) modelling has been identified as one of the most important industrial needs related to nuclear reactor safety. A severe PTS scenario limiting the reactor pressure vessel (RPV) lifetime is the cold water emergency core cooling (ECC) injection into the cold leg during a loss of coolant accident (LOCA). Since it represents a big challenge for numerical simulations, this scenario was selected within the European Platform for Nuclear Reactor Simulations (NURESIM) Integrated Project as a reference two-phase problem for computational fluid dynamics (CFDs) code validation. This paper presents a CFD analysis of a stratified air-water flow experimental investigation performed at the Institut de Mécanique des Fluides de Toulouse in 1985, which shares some common physical features with the ECC injection in PWR cold leg. Numerical simulations have been carried out with two commercial codes (Fluent and Ansys CFX), and a research code (NEPTUNE CFD). The aim of this work, carried out at the University of Pisa within the NURESIM IP, is to validate the free surface flow model implemented in the codes against experimental data, and to perform code-to-code benchmarking. Obtained results suggest the relevance of three-dimensional effects and stress the importance of a suitable interface drag modelling.
1. Introduction
The European Platform for Nuclear Reactor Simulations (NURESIM) Integrated Project aims at developing a common European Standard Software Platform for modelling, recording, and recovering computer simulation data for current and future nuclear reactor systems [1]. NEPTUNE [2] is the thermal-hydraulics tool of NURESIM, and is designed to simulate two-phase flow in all situations encountered in nuclear reactor power plants. The present work is aimed at the validation and benchmarking of NEPTUNE_CFD, the two-phase CFD tool of the NEPTUNE platform [3, 4].
Since PTS has been identified as one of the most important aspects related to nuclear reactor safety, some relevant PTS scenarios were chosen as reference test cases for CFD code validation within the NURESIM IP [5, 6]. A severe PTS scenario limiting the reactor pressure vessel (RPV) lifetime is the cold water emergency core cooling (ECC) injection into the cold leg during a loss of coolant accident (LOCA). Complex phenomena take place during this scenario, such as turbulent mixing in the jet region and downstream of the impingement zone, stratified two-phase flows, phase change at the steam water interface. This paper deals with the study of a stratified air-water flow experiment performed at the Institut de Mécanique des Fluides de Toulouse in 1985 [7, 8]; this flow configuration is likely to share common physical features with the chosen PTS scenario.
To validate the two-phase models implemented in NEPTUNE_CFD V1.0.6, numerical simulation results were compared with both experimental data and predictions from two commercial CFD codes, namely, ANSYS CFX 10.0 and FLUENT 6.1.
2. Experimental Facility and Tests
The experimental facility (see Figure 1) consists of a quasihorizontal channel with rectangular cross-section (0.1 m high, 0.2 m wide, and 13.0 m long), with an inclination of 0.1° with respect to the horizontal plane. The rectangular channel is connected to the water and air inlet and outlet volumes. Desired mass flow rates are granted by a recirculation pump and all the facility walls are adiabatic.
The facility is equipped with sensors located at 7.05 m, 9.10 m, and 11.10 m from the inlet section, which provide the measurements of mass flow rates, local instantaneous water height (by capacitance wire sensors), local mean and fluctuating values of horizontal and vertical velocity components, as well as Reynolds stress tensor (by a laser Doppler anemometer in water and by hot wire anemometer in air). Documentation is available [7] for three test cases conducted at ambient pressure and temperature, characterized by constant water mass flow rate and different air mass flow rates. This work deals only with one of these tests, namely, the “T250” experiment, in which water and air bulk velocities are 0.395 m/s and 3.66 m/s, respectively. Under this condition, mean water depth is measured to be 38 mm.
3. Experimental Test Simulation with Ansys Cfx and Fluent
3.1. Preliminary Results of Two-Phase Calculations
Since the width of the duct is large compared to the height, a two-dimensional velocity field has been assumed in order to perform preliminary calculations. These analyses, carried out with the ANSYS CFX 10.0 code [9], adopted a 1-mm wide one-cell thick grid since it does not allow assuming a real two-dimensional domain. A spatial discretization of the channel has been created using ANSYS ICEM CFD 10.0 code, counting ~20 k hexahedral elements (~40 k computational nodes). Elements refinement has been provided near the walls and at the expected height of interface between fluids; anyway it is worth noting that in more realistic applications, the interface is moving in the domain, so that such grid refinement could be obtained only with dynamic meshing adaptation.
The “inhomogeneous” two-phase flow model was selected, since some interface instability has been observed in all preliminary simulations. This model solves one velocity field for each phase (resulting in two different coexisting fields in the domain); while the “homogeneous” setting has been adopted for turbulence together with the shear-stress transport (SST) model, providing only one field shared by the phases. Moreover, the code default interface treatment model has been used, the so-called “standard free surface” [9].
A uniform velocity for both air and water has been assumed at the inlet. Air constant pressure and water hydrostatic profile have been imposed at the outlet section according to measured water height. Finally, adiabatic upper and lower walls with no slip and symmetric lateral faces of the domain have been selected.
Figure 2 shows the calculated velocity in the test section compared with experimental data. Air-water stratification has been correctly predicted, but relevant differences can be observed: water level is calculated about 24% lower, while the maximum air velocity is reduced by 10% and it is no longer placed between the top wall and the interface (66 mm), but closer to the wall (81 mm, ~20% higher). These relevant mismatches suggest that the modelling of the frictional drag between the phases is overestimated, and a deeper understanding of the experimental data is needed. In order to investigate these problems, single-effect analyses have been carried out as described in the following Sections, together with some sensitivity analyses on the most relevant parameter.
(a)
(b)
3.2. Experimental Data Understanding
All performed experimental tests assume the same value of water flow rate. Different water heights are measured into the channel depending on the different equilibrium conditions between the forces acting on the fluids, that is, the drag force between air and water, the longitudinal component of the gravitational force due to the inclination of the duct, and the friction forces acting on walls. Except for the gravitational force, the others depend on fluid velocity, thus changing their values flowing into the channel: the drag force increases while the interface friction force decreases up to reaching the equilibrium condition. An incorrect prediction of one of them can justify the mismatch between calculated profiles and experimental data.
Furthermore, it is worth noting that water average velocity resulting from measured data at probe section (calculated by means of the trapezes rule, thus underestimating the real value) is >4% greater than the value provided in [7] for the water bulk velocity. The same occurs for the air. Since the flow rate is constant along the channel, this mismatch could be due to the development of a 3D profile. In fact, in a real 3D geometry, the velocity profile measured at the symmetry plane (as shown in all the plots) is the maximum profile of a developed 3D flow (see Figure 3). It is possible to conclude that considering a 2D computational domain implies a loss of information related to the flow development.
3.3. Single-Phase Analysis
In this analysis, the computational domain was splitted into two separate parts along the experimental interface level. Spatial discretizations have been created for each phase channel in both two- and three-dimensional configurations. Different node distributions have been employed to evaluate the grid requirements for CFD simulation to correctly reproduce near-wall effect and flow developing. The most relevant grid details are reported in Table 1 for both FLUENT 6.1 [10] and CFX 10.0 [9] codes. These characteristics have been established following the main findings of preliminary mesh sensitivity investigations.
In single-phase calculations, the interface has been modelled as a wall moving at the expected velocity of the free surface. Since this value is not available, it has been imposed in the range 0.50–0.65 m/s, according to the top measured water velocity, 0.502 m/s, which is the available data closer to the interface. The boundary conditions have been imposed according to the preliminary calculation documented in Section 3.1.
Direct numerical simulation (DNS) and large eddy simulation (LES) of stratified flows [11] were used in the framework of the NURESIM Integrated Project to derive some closure laws for interfacial momentum, turbulence, and heat transfer. Future work is still necessary to implement these laws in CFD codes as well as to compare predictions with DNS-LES studies on the same flow conditions [12] and to validate them against experimental data. Anyway, this subject is beyond the aim of the present article.
3.3.1. Single-Phase Water Flow
In Figure 4, the obtained results are shown in terms of longitudinal water velocity profile at the test section located at 9.1 mm from the inlet. The experimental profile is correctly predicted from a qualitative point of view by both CFX and FLUENT codes. No relevant improvements are obtained varying the interface velocity except for the region closer to the moving wall in both two- and three-dimensional calculations. However, two-dimensional calculations underestimate the velocity values by about 10% than three-dimensional ones. Thus, domain three dimensionality has a great relevance on water velocity profile and cannot be neglected in the simulations. Moreover, grid sensitivity analysis has shown that limited improvement is obtained by increasing the number of cells.
(a)
(b)
(c)
Figure 5 shows the comparison between calculated turbulent kinetic energy and experimental data both for two- and three-dimensional simulations. The third dimension has great relevance on near-wall values, increasing them by about 20%. Calculated values in the region near the interface have a different behavior from experimental data owing to the presence of a solid wall. Turbulence produced by the contact between fluids at the interface is not included in the model. Although these differences in shape and local values, predicted turbulent kinetic energy is in good agreement with measurements, especially for three-dimensional calculations.
(a)
(b)
3.3.2. Single-Phase Air Flow
Figure 6(a) shows the transversal air velocity profiles at the test section in three-dimensional simulations calculated with CFX, which is predicted with relevant differences on both shape and values. Negligible effects on the results are obtained varying the interface velocity. Two-dimensional calculations show the same behavior as well as results obtained using FLUENT. Moreover, as shown in Figure 6(b), the experimental data provide a greater maximum velocity with respect to both two-dimensional (15.8% higher), and three-dimensional simulations (11.4% higher). This is a further confirmation of the three-dimensional flow structure supposed in Section 3.2. However, the systematic underestimation of the velocity needs further investigations.
(a)
(b)
Finally, from sensitivity analysis performed with FLUENT, no relevant advantages are obtained by increasing the cell number or changing the turbulence near-wall treatment.
4. T250 Experimental Test Simulation with Neptune_cfd
In the hypothesis of planar symmetry, the computational domain is constituted by a 2D section of the channel; it has been modelled with three successively refined 2D grids built up with the ANSYS ICEM code: 7530, 10040 and 15060 cells, respectively. All imposed boundary conditions were as for CFX calculations (Section 3.1); initially a stratified air/water flow was established and parabolic velocity profiles were imposed for both water and air flows at inlet. Calculations were run with NEPTUNE_CFD V1.0.6 by means of an input deck kindly provided by Mr. Pierre COSTE (CEA/Grenoble). The model was considered for both air and water turbulence together with interface turbulence production, “Pierre Coste Large Interface Model” [13] was selected for the drag coefficient.
As Figure 7 shows, water velocity profile is quite well predicted in all three cases (maximum error ~7%), while air velocity profile is appreciably underestimated, especially in the bulk region between wall and free surface (maximum error ~12% for the coarser grid and ~10% for the finer one). It seems that mesh refinement does not produce important improvements, except for the air velocity profile in the region near the interface.
(a)
(b)
Figures 8(a) and 8(b) show the turbulent kinetic energy profile at probe location for water and air flow, respectively. As in the previous case, the profile is qualitatively well predicted for the water flow, while for the air flow, it is significantly underestimated in the bottom region and overestimated in the upper one. It can be observed that calculated values with refined grids better match experimental data in the upper region of the water stream, with some underestimation near the wall (maximum error ~45%). It is worth noting that the code is able to catch the increase of water turbulence near the free surface due to the air stream, but not due to the presence of solid walls. On the contrary, air turbulent kinetic energy profile does not get significantly better with mesh refinement; results are underestimated near the interface (maximum error ~66%) and slightly overestimated elsewhere.
(a)
(b)
Calculations were also run considering the “separated phases model” for the interphase drag coefficient. The resulting velocity profiles seem to be very similar to that predicted by ANSYS CFX and shown in Figure 2. In both cases, the interface level is underestimated and the maximum air velocity is reached in the region near the upper wall instead of the air stream core. This could be due to an incorrect modelling of the drag coefficient. Taking into account results presented in Section 3.3, a 3D simulation was also set up considering a 1004020 grid and the “Pierre Coste Large Interface Model” for the drag coefficient. Unfortunately, only 29 seconds were calculated after one-week run on two processors, but preliminary results were encouraging.
5. Conclusions
A Computational fluid dynamic analysis of a stratified air-water flow experimental investigation performed at the Institut de Mécanique des Fluides de Toulouse in 1985 [7] was performed. The aim was the comprehension of the experimental data and of the role played by some fundamental parameters. The simulation has been performed by means of three different CFD codes: NEPTUNE_CFD V1.0.6, FLUENT 6.1, and CFX 10.0. The spatial discretizations have been modelled with GAMBIT 1.0 and ANSYS ICEM 10.0 softwares.
Preliminary results of two-phase CFD calculations with a two-dimensional domain suggested that three-dimensional effects are not negligible, so that 2D simulations are not suitable to correctly predict this stratified fluid flow. To better understand the physics of the problem, single-phase analyses were conduced with FLUENT and ANSYS CFX comparing 2D and 3D simulations for both air and water single-phase domain. As a result, relevant improvements of both water and air velocity profile were achieved with 3D simulations. It is worth noting that in such single-phase analyses, the water level was not calculated but fixed according to experimental data.
Two-phase simulations by means of NEPTUNE_CFD code, despite taking into consideration a 2D domain, showed better agreement with measured data when considering the new “Pierre Coste Large Interface Model” for the drag coefficient: water level was correctly predicted and error in velocity profiles decreased, even if some underestimation of the air velocity is still present. Moreover, CFX and NEPTUNE_CFD standard models gave similar results, putting in evidence the fundamental role played by the drag coefficient modelling. Nevertheless, a systematic underestimation of the air medium velocity suggests that further information on the experiment and boundary conditions is needed.