Thermal-hydraulics validation of CFD code for light water nuclear reactors against benchmark experimental results

The cooling of a nuclear reactor depends on a suitable fluid flow pattern among its fuel elements aiming the removal of heat produced in the fuel. In case of light water reactors, an excess of heat drives the fluid to change its phase from liquid to vapor, significantly reducing its capacity to remove heat and leading the reactor to a Loss of Coolant Accident. Numerical simulations using a CFD code is a suitable tool to address this kind of problem and explore the conditions that should be avoided during the reactor operation. The commercial CFD codes had proven to be reliable to simulate with a high accuracy and confidence the thermal-hydraulics of a sort of equipment and systems, avoiding spending efforts and financial resources in the development of new codes that, essentially, perform the same tasks. Despite of it, the CFD codes must be validated, such as against experimental results. To comply with this objective, a benchmark fuel element was purposed and experimentally essayed to provide experimental results for CFD codes calibration. The results of this essay are provided to the four types of subchannels for a 5x5 PWR fuel element, with results provided as density and void fraction. This work presentes the preliminary results obtained with CFD numerical simulations using the ANSYS-CFX code for the central subchannel with active rods for stead state operation. The results demonstrated that the ANSYSCFX is adequate to simulate with high accuracy the flow in this subchannel.


INTRODUCTION
In nuclear reactors, the cooling fluid develops a very important task, ensuring that the heat produced by the fissions on nuclear fuel is suitable removed. In Pressurized Water Reactors (PWRs), the cooling fluid should be kept in liquid phase to this, removing the heat efficiently from the fuel elements and its fuel rods. This liquid should flow through the spaces among the fuel rods; this space is named as subchannel [1]. By changing the flow conditions, that is, increasing its temperature, reducing its flow rate, increasing the power in fuel elements, etc, leads to form some vapor, increasing the quality of coolant. As the quality increases, a plug is developed, gradually blocking the flow of the coolant through the fuel element. This effect contributes to build up the heat that should be removed from fuel rods. However, the vapor phase has only a limited ability to remove heat transfer due to the low convection heat transfer coefficient and low thermal conductivity. As consequence, more heat is attained to the liquid phase, more vapor is formed and so on. Thus, the vapor formation should be avoided as much as possible. The PWR reactors are projected to avoid this vapor formation during its operation while Boiling Water Reactors (BWR) operates at this regime. Despite of PWR project avoids the vapor bubbles formation, it could easily be formed during transients, such as in case of accidents [1][2][3].
The PSBT is an experiment developed to investigate the limits of a reactor in several conditions, including overpressure, under-pressure, mass flow reduction, transient condition, etc. The experimental apparatus and parameters considered are available for four subchannel geometries, for a 5x5 fuel element in which the fuel rods are substituted by electrical resistances. The experimental results are available in form of void fraction, cooling fluid density, etc, which could be used as reference to calibrate numerical simulation codes [3][4][5].
The Computational Fluid Dynamics (CFD) numerical simulations allow to evaluate the fluid flow conditions, including multiphase flows, with detail level that often is not possible or is impracticable in experiments. As consequence, the reactor could be operated efficiently and within the safety margins [1][2][3]. The CFD simulations are performed by numerical codes, which solves the governing equations (mass, energy and momentum) using the finite volume method. Some of these codes have found in commercial licenses, such as the CFX  and the Fluent  , both of ANSYS  , while others could be specifically developed for research purposes [3][4][5][6][7].
The aim of this work is to present the initial part of the research under course in which CFD simulations are being conducted to evaluate the flow conditions in a 5x5 fuel element. The results obtained for the central-typical subchannel, with ¼ of active rod at each edge of the channel (or 1 active rod at the center of the channel) is presented and compared against the experimental results of PSBT experiment. A sensibility analysis regarding the influence of the bubbles mean diameter as well as the bubbles drag model is also provided, demonstrating that the bubbles diameter as well as the drag model affects significantly the results accuracy. The results demonstrate that the results obtained with the commercial code used and its models have a good accuracy, being a valuable tool to investigate the thermal-hydraulic phenomena without the need of construction of experimental apparatus. The rest of this work is organized as follows: section 2 describes the methodology and governing equations used in the present work. Section 3 details the boundary conditions. Section 4 presents the results and a discussion regarding them. Section 5 presents the conclusions and final remarks, followed by the acknowledges and references.

MATERIALS AND METHODS
The simulations were conducted with the CFD code ANSYS-CFX  Release 19.2 [7]. The geometry refers to the central-typical subchannel of a 5x5 fuel element, with four quarters of active fuel rods located in each corner of the subchannel (the S1 geometry of the PSBT experiment). The main dimensions of the geometry as well as the mesh and boundary conditions are detailed in section 3 [4]. The simulations were run as follows: • Symmetry relations were used to save the volumes used in the mesh to discretize the geometry. In this manner, only a quarter of the geometry was modeled using the Design Modeler of ANSYS Workbench ® platform. No extended lengths were included. Despite of it, the simulations converged normally without any inter-occurrence; • Only the fluid domain was modeled since the objective of the present work is evaluate the thermal-hydraulics at the coolant only; • The mesh used is an extruded type mesh, applied to the face of subchannel cross section area. This kind of mesh is a 2D mesh, which is copied in layers to form a 3D mesh. The quantity of layers were adjusted to 500 layers, which is refined enough to capture suitable the phase change across the channel length. This type of mesh contributes to save the quantity of volumes required do discretize the geometry and thus, reduces the computational effort without compromise the accuracy of results.
• The boundary conditions that describes the coolant fluid were inputted in the CFX ® Preprocessor. Two fluids were used to make possible to describe the multiphase flow. These fluids are water and vapor (water vapor). The boundary conditions considered are available at [3,4,6,8,9] and detailed in section 3; • The RMS residuals and imbalances of governing equations was considered as convergence criterion. In this first part of the research a value of 10 -4 was considered, which is the minimum value recommended in the ANSYS-CFX  manual for a simulation be considered converged [10]. The simulations converge when the residuals are equal or below this target value, when the imbalances are as close as possible of 0% and when the parameters values of interest (void fraction, fluid density and outlet temperature) does not varies anymore.
• The simulations were run in the ANSYS ® Solver, which solves the governing equations using the boundary conditions and models settled in the pre-processor. The governing equations are given by equations (1) to (5) for mass (equation 1), energy (equation 2) and momentum conservation (equations 3 to 5), as given by [10]; • The results are extracted and evaluated in the ANSYS ® Post-processor. The results evaluated are the void fraction and density of the coolant fluid at the cross-sectional plane located 1400mm above the subchannel inlet (as in the experimental apparatus) [4].
• The results of each simulation were inserted in a Microsoft Excel  spreadsheet, allowing to analyze and compare them with the experimental data. They are available in section 4; • A mesh dependency study was conducted to evaluate the influence of mesh refinement over the results. Four meshes were used in this study, which results are available in section 4; (1) in which ρ is the coolant fluid specific mass, t is the time, x, y and z are the cartesian coordinates, u, v e w are the velocity components to x, y and z directions. (2) in which Cp is the coolant fluid constant pressure specific heat, T is the temperature, k is the fluid thermal conductivity and ST is the source term. (5) in which P is the pressure, μ is the fluid absolute viscosity, SMx, SMy e SMz are the source terms in x, y e z directions, respectively.

BOUNDARY CONDITIONS
The experimental apparatus taken as benchmark for the simulations performed is show schematically in Figures 1 and 2.  Source: [4] The void fraction and coolant density were evaluated in the cross-section plane located at 1400mm above the channel inlet (experimental and simulation). The quantity of volumes and nodes of the meshes simulated are available in Table 1 Table 2, and is based in the considerations and data of [4]. For the mesh dependency study was conducted considering the experiment #1.1223, chosen randomly among the other experiments. The simulations were run considering the steady state flow and turbulent regime. Transient simulations are subject of future works. The turbulence intensity was settled as high intensity (10%, using the values available in the code), after some preliminary simulations using with 1% (low intensity) and 5% (medium intensity). Only the values obtained using the high intensity turbulence is presented in the present work because its better agreement with experimental values. The turbulence intensity is defined according equation (6) [10].
The inlet boundary conditions were settled according [4], to each of the experiments considered.
These and other boundary conditions are summarized in Table 2. The channel walls were divided in three distinct regions, as shown in Figure 4. These regions are, namely: heater, symmetry, insulated walls.

Figure 3: Cross section of the central subchannel with active rods -S1 geometry.
Source: adapted from [4] Different boundary conditions were considered to each of these regions, as following detailed.

Symmetry
• In ANSYS-CFX  there is only an adjust to this condition: symmetry. In this manner, the code consider that the phenomena that occur in one side of the symmetry plane, occur equally at the other side. Notwithstanding, this condition allow a virtual mass exchange (inlet and outlet of mass in the domain simulated through the symmetry plane);

Fluids
• Liquid is a continuum medium and vapor is the dispersed phase; • Eulerian approach; • Non-homogeneous total energy model; • Bubble mean diameter: 0.5mm and 1.0mm; • Buoyancy and gravitational forces; • Non-homogeneous multiphase flow; • Saturation temperature: according the pressure of each experiment (see Table 2); • Mass transfer due to phase change in the coolant; • Turbulence: non-homogeneous, zero equation; • Turbulent heat transfer enhancement; • Surface tension: constant, varying according the inlet temperature of the fluid; • Wall lubrication force: o Frank model. The coefficients and constants values used are those recommended by the ANSYS-CFX  code manual [10]; • Boiling models [10]:

RESULTS AND DISCUSSION
In this section the results obtained for the void fraction and density at the subchannel simulated, changing the mean bubble diameter, turbulence models and its turbulence intensity are presented in section 4.1. The section 4.2 presents the mesh dependency study.

Results void fraction and density
The selection of the bubble mean diameter was made after some initial simulations in which the sensibility of the results to this parameter was evaluated, within the range from 0.1mm to 2.5mm.
Based on these initial results, it was concluded that the diameter that best fits with the experimental results are diameter of 0.5mm and 1.0mm, which were used in the remaining simulations. The results obtained for the coolant void fraction and density, for bubbles with mean diameter of 0.5mm and 1.0mm, are available in Tables 3 to 6.    As could be observed, the models of ANSYS-CFX ® used in the present work, the parameters values adopted and boundary conditions, represents with high accuracy the void fraction and coolant density in the central-typical subchannel of a PWR fuel element in different operational conditions. It could be observed also that the results obtained for bubbles with 0.5mm has a higher accuracy than those obtained considering bubbles with 1.0mm of mean diameter, considering the absolute and relative errors, which are calculated according equations (7) and (8), respectively.
When the relative error is compared with the accuracy of the experimental measurements (accuracy: 3%; void fraction tolerance: 0.030; density tolerance: 15kg/m³; [4]), it is possible to see that the majority of simulations has errors within the experimental error margin. Thus, the simulation results are considered representative [3,9].
In which experimental value is the value obtained during the experiments and simulated values is the values obtained with the CFD simulations.
Besides the values of Tables 3-6  Source of image b): [6] As could be observed, the results obtained using the code ANSYS-CFX ® are physically consistent and very similar to the experimental void distribution. It is important to highlight that these results were obtained with a very coarse mesh and that they could be improved with a refined mesh. This evaluation is presented in the next section.

Mesh dependency study
For the geometry simulated (see Figure 4), a mesh dependency study was conducted to evaluate the influence of mesh refinement over the results and the simulation convergence. The results obtained for the coolant void fraction and density are available in Figures 6 a) and b) considering as case of study the experiment #1.1223.

CONCLUSION
In the present work the code CFD ANSYS-CFX  Release 19.2 was used to simulate a quarter of a central-typical subchannel (geometry S1) of the benchmark 5x5 fuel element. The results were compared with experimental data and demonstrated a good agreement since most part of results have absolute and relative errors within the experimental error. The best diameter to bubbles is 0.5mm, despite of the bubbles with 1.0mm had presented a good agreement too. Despite of it, additional efforts should be spent in investigation of the reasons that make the simulations results in some conditions fall out of the experimental error margins. Regarding the mesh refinement, it was demonstrated that no significant dependency exists since there is a small spready in results obtained for the void fraction and density in the coolant fluid, not justifying the additional computational effort associated to the refined meshes.

ACKNOWLEDGMENT
Authors wish to acknowledge UFABC, ANSYS and FAPESP to support this work.