Review of audit calculation activities on the applicability of CFD software to nuclear safety problems

The Korea Institute of Nuclear Safety (KINS) is tasked with the technical supporting for approving the safety of specific components or design modifications in nuclear power plant. When a licensee submits analysis to the Nuclear Safety & Security Commission (NSSC) for acceptance, KINS staff reviews this analysis and sometimes conducts independent audit calculations. Though recently licensing applications supported by using Computational Fluid Dynamics (CFD) software are increasing for nuclear safety problems, there is no CFD software which obtains a licensing from the domestic regulatory body until now. Additionally, there is no domestic regulatory guideline for the comprehensive evaluation of CFD software. Therefore, from the nuclear regulatory perspective, it is necessary to perform the systematic assessment and prepare the domestic regulatory guideline for checking whether valid CFD software is used for nuclear safety problems. In this paper, regulatory audit calculation activities on the applicability of CFD software to nuclear safety problems was briefly explained.


Introduction
The KINS is tasked with the technical supporting for approving the safety of specific components or design modifications in nuclear power plant. When a licensee submits analysis to the NSSC for acceptance, KINS staff reviews this analysis and at times conducts independent audit calculations. In some cases, CFD methods are used for these confirmatory studies which support safety evaluations, build confidence and reduce uncertainty in the final regulatory decisions. Typical examples are (1) thermal-flow analysis for high density spent fuel storage racks (2) estimation of moderator subcooling margin for Pressurized Heavy Water Reactor (PHWR).
Though recently licensing applications supported by using CFD software are increasing for nuclear safety problems, there is no CFD software which obtains a licensing from the domestic regulatory body until now. Additionally, there is no domestic regulatory guideline for the comprehensive evaluation of CFD software. Therefore, from the nuclear regulatory perspective, it is necessary to perform the systematic assessment and prepare the domestic regulatory guideline for checking whether valid CFD software is used for nuclear safety problems.
In this study, regulatory audit calculation activities on the applicability of CFD software to nuclear safety problems was briefly explained for the following items.
• Effect of asymmetric inlet flow rate on reactor coreinlet flow-distribution • Single phase thermal stratification by emergency core cooling system injection • Flow distribution inside a fuel assembly with twist-split type mixing vane 2 Case I: Effect of asymmetric inlet flow rate on reactor core-inlet flowdistribution

Overview
In general, an impeller of a Reactor Coolant Pump (RCP) installed in a Pressurized Water Reactor (PWR) is manually manufactured and therefore there may be a possibility of performance difference (about maximum 4% of design flow rate) per a RCP. Therefore, Kim et al. [1] measured reactor core-inlet flow-distributions under asymmetric inlet flow rate condition (case1 in Table 2) and compared it with those from symmetric inlet flow rate condition.
In this study, CFD analysis was conducted to verify the validity of the reactor internal flow distribution measured under asymmetric flow rate conditions by using ANSYS CFX and compared with the calculated results under the existing symmetric flow rate condition. In addition, we attempted to confirm whether asymmetrically entering flows due to the difference in RCP performance are properly mixed while passing through low support structures of a reactor, and provide core inlet flow rate distribution equivalent to symmetric flow rate conditions. The internal structures of the reactor model (e.g., flow skirt and upper/lower core structures) had almost the same shapes as those in the original APR+, and satisfied geometrical similarity [1,2]. The core-inlet flow-rate distribution could be obtained by measuring the differential pressure and discharge coefficients at the venturi region of each core simulator. A total of 257 core simulators, which corresponded to the fuel assemblies, were installed in the reactor model. The upper head of the reactor, and some core-bypass flow-paths were neglected in the reactor model because these parts were expected to have little influence on the core-inlet flowrate distribution. The criteria of the allowable data scattering for each core simulator inlet flow-rate distribution was ±1.5% [2].

Analysis model
The test matrix consists of three flow conditions, i.e., the symmetric or asymmetric flow conditions for 4pumps operation, and the flow condition for 3-pumps operation. In this study, CFD simulation was conducted under the asymmetric flow conditions for 4-pumps operation (see Table 2). In Table 2, Asym. (case1~case3) correspond to asymmetric flow conditions for 4-pumps operation. In detail, case1 is real test condition and case2 & 3 are conditions for computation only. APR+ reactor internals are complex structures which support fuel assemblies, control rods and measuring instruments. The internal structures, especially those located in the upstream of the reactor core, may have a significant influence on the core-inlet flow-rate distribution; depending on both their shapes, and the relative distance between the internal structures and the core inlet [3]. Therefore, an exact representation of these internal structures is needed for CFD simulation of the core-inlet flow-rate distribution. However, such an approach requires a great deal of computing resources to analyze the real-flow phenomena inside a reactor.
In this study, as shown in Fig. 1, among the reactor internal structures located upstream of the reactor core, the real geometries of a flow skirt, Lower Support Structure Bottom Plate (LSSBP) and In-Core Instrumentation (ICI) nozzle support plate, were considered because these internal structures could significantly influence the flow-rate distribution at the core inlet. Meanwhile, to reduce total numbers of elements and thus minimize the required amount of computation, fuel assemblies and some internal structures (e.g., control-element guide tubes) were simply considered as each bulk volume (porous domain). Then, in order to reflect the velocity field and pressure drop occurring in the real-flow region; porosity and Isotropic Loss Models were applied to the porous domain.  Porosity is the ratio of the volume of fluid region to total volume; including both fluid and solid regions. It has an effect on flow acceleration in the porous domain. In this study, the porosity was determined by considering the real geometry of the reactor internal structures, and its magnitude was in the range of 0.5-0.75. A momentum source was used to model the momentum loss in the porous domain; which corresponds to a pressure drop in real reactor vessel. Loss coefficients were adjusted to match the magnitude of the pressure drop found in the porous domain, with those of the measurement.

Numerical modeling
The flow inside the scaled-down reactor model was assumed to be steady, incompressible, isothermal and turbulent. A high resolution scheme was used for both the convection-terms-of-momentum equations andturbulence equations. The solution was considered 'converged' when the residuals of the variables were below 310 -4 , and the variations of the target variables were small. Simulation was conducted with the commercial CFD software, ANSYS CFX R14.5.
The k- model, which is one of the most prominent Reynolds Averaged Navier-Stokes (RANS)-based turbulence models, was used to simulate the turbulent flow inside the scaled-down reactor. The reason is that this model has proven to be numerically stable and has offered a good compromise in terms of accuracy and robustness. In a previous study [4], turbulence models available in ANSYS CFX R13, for example k- model, Shear Stress Transport (SST) model, and SSG (Speziale, Sarkar and Gatski) Reynolds Stress model, were used to examine the turbulent flow inside the scaled-down reactor. Although the reactor internal-flow pattern differed locally; depending on the turbulence models used, the k- model showed the best agreement with the experimental data. Fig. 2 shows the grid system for the computational domain that had the same size as the test facility. A hybrid mesh, made up of tetrahedrons, pyramids and prisms, was generated to prevent the oversimplification of the geometry, and to have more efficient mesh distribution. Prism layers were used to get higher resolution in the near-wall region. Total numbers of nodes were 7.310 7 and maximum y + in the downcomer was about 305. By referring to the test conditions, as shown in Table  2, inlet flow-rates were imposed at each cold leg. Turbulence intensity at the inlet was assumed to be 5%.
Light water at 60 ℃ was used as the working fluid. The 'average pressure over the whole outlet' option; with a relative pressure of 0 Pa, was used at each hot leg as an outlet-boundary condition. A no-slip condition was applied at the solid wall. To model the flow in the nearwall region, scalable wall functions were applied.

Results and Discussion
Fig . 3 shows the circumferential distribution of velocity (y-direction) in the downcomer (-0.6 m from the center of the cold leg). In the downcomer, velocity showed a non-uniform distribution; with peak negative magnitude between cold legs. In case of asymmetric flow condition, as the flow rate in a cold leg CL2A decreased, asymmetry of y-directional velocity profile between CL2A and other cold legs increased.    4 shows the normalized, mass-flow rate at the core-inlet plane. The magnitude is the ratio of the massflow rate per each fuel assembly, to the average massflow rate at the core-inlet plane.
The core inlet flow rate range calculated from the asymmetric flow rate conditions (case1~case3) was 80% ~ 141% of the average flow rate of the fuel assembly, which was the same as the calculation result at the symmetric flow rate condition. For the measured values as well, the core inlet flow rate range in the asymmetric flow rate condition was 86% ~ 126% of the average flow rate of the fuel assembly, which was the same as the measured flow rate at the symmetric flow rate condition. Fig. 5 shows distribution of the normalized core-inlet mass-flow rate along the core centerline (A-A). For calculation, though there was a minor difference (maximum 0.07) for certain fuel assemblies (no.74, 110 etc.), core inlet flow rate showed overall similar distribution regardless of asymmetric magnitude in an inlet flow rate. For measurement, although there was a minor difference (below 0.01) in a fuel assembly no.47, core inlet flow rate showed overall same distribution regardless of asymmetry in an inlet flow rate. In summary, asymmetrically entering flows due to RCP performance differences were properly mixed while passing through the internal structures installed at the lower part of the reactor, confirming that the core inlet flow distribution is equivalent to that under the symmetric flow condition.
3 Case II: Single phase thermal stratification by emergency core cooling system injection

Overview
Pressurized Thermal Shock (PTS) phenomenon is characterized by multi-dimensional and non-equilibrium flow conditions. For example, in case Emergency Core Cooling System (ECCS) is operated during Loss Of Coolant Accident (LOCA) in a PWR, PTS phenomenon occurs as cooling water is injected into a cold leg, mixed with the hot primary coolant, and then entrained into a reactor vessel. Insufficient flow mixing may cause temperature stratification and steam condensation. In addition, flow vibration may cause thermal stresses in the surrounding structures, which reduce the life of the reactor vessel. Because of this importance, the effect of PTS has been actively studied either experimentally or numerically [5]. Among these, the ROSA project [6,7] undertaken by the OECD/NEA includes six types of the separate or the integral effect tests. In this study, the calculation with ANSYS CFX R17.2 was performed for the Test 1『Temperature stratification and coolant mixing in unsteady state during ECCS injection』and then the predicted results were compared with the measured data.

Analysis model
The objective of Test 1 is to clarify the phenomenon of temperature stratification on the cold legs and the upper region of downcomer of the reactor vessel, which are important for the evaluation of PTS during the coolant injection by the ECCS, and to obtain the measured data for the main parameters with the aim of validating the simulation code and the numerical modeling. Fig. 6 shows the schematic diagram of the analysis model.  As shown in Fig. 8, 9 and 10(a), to investigate the thermal stratification phenomenon in the cold leg, 21 thermocouples were installed in 3 columns and 7 rows at two cross sections (TE2, TE3) in the cold leg between the ECCS line A and the downcomer. The nominal accuracy of the thermocouple measurements is ± 2.75 K. In this study, the measured and the calculated results were compared for the cases where the cold coolant was injected through the ECCS line A during time intervals of 2,091 s ~ 2,187 s under the natural circulation condition after the main coolant pump stopped.

Numerical modeling
In this study, the turbulent flow field inside the ROSA/LSTF (Large Scale Test Facility) was calculated under unsteady, single-phase, incompressible and nonisothermal conditions by using the commercial computational fluid dynamics software ANSYS CFX R17.2. High resolution scheme with quasi-second-order accuracy was used for the convection-terms-ofmomentum and -turbulence equations. The second order backward Euler option, the default implicit time step method with the second order accuracy in ANSYS CFX R17.2, was used as the time discretization scheme. For the calculation of the unsteady flow, the calculation was performed during time intervals of 2,078 s ~ 2,218 s under the natural circulation condition with a time step of 0.1 s. It was determined that the solution was converged when the root mean square residuals of the individual equations were 10 -5 or less at each time step. Baseline Reynolds stress model was used to simulate the turbulent flow inside the ROSA/LSTF. To model the flow in the near-wall region, the automatic wall treatment was applied. Fig. 7 shows the grid system, generated by using ICEM-CFD software, for the computational domain that had the same size as the test facility. The grid shape was hexahedral, and the total number of elements used in the calculation was 1.99 × 10 6 . In the computational domain, the maximum y + was about 411, and its location was near the region where the inner wall of the pressure vessel and the hot leg were bordered. The measured unsteady flow rate and temperature were applied at inlet of the cold leg A/B and the ECCS injection line A. For reference, the temperature difference between the hot coolant at the inlet of the cold leg A/B and the cold cooling water at the inlet of the ECCS injection line A was 200 K or more. On the other hand, the turbulence intensity at each inlet was assumed to be 5.0%. The average static pressure condition was applied at the outlet. No-slip and adiabatic conditions were applied at all wall interfaces including the pressure vessel.  Cold cooling water injected at the inlet of the ECCS injection line A moved toward the lower part of the cold leg A and then flow mixing with the hot coolant occurred. The mixed flows proceeded along the cold leg A and then entered into the downcomer along the outer wall of the pressure vessel.

Results and Discussion
To identify the thermal stratification region, wall surface temperature distribution was shown in Fig. 9. The strongest thermal stratification was found at the location where the cold cooling water injected at the inlet of the ECCS injection line A contacted with the lower part of the cold leg A. In addition, as previously explained in Fig. 8, because the mixed flow moved along the lower part of the cold leg A and then entered into the downcomer along the outer wall of the pressure vessel, the thermal stratification phenomena occurred in these regions.   Fig. 8 and 9. Therefore, it is expected that the relatively weaker flow mixing between the cold cooling water injected at the inlet of the ECCS injection line A and the hot coolant occurs in comparison with the section (TE3), and, as a result, the measured temperature at three positions of the section (TE2) showed the obvious temperature gradient at three locations in the cross-section TE2 during the ECC water injection period. Though there were a little difference between the measured and computational results, the predicted fluid temperature was relatively in good agreement with the measured data during the ECC water injection period. Due to the non-disclosure limitations of measurement data, quantitative values of fluid temperature corresponding to the vertical axis of the graph ( Fig. 10(b)~(d)) were not specified.

Case III: Flow distribution inside a fuel assembly with twist-split type mixing vane 4.1 Overview
In a PWR, the appropriate heat removal from the surface of fuel rod bundle is important for thermal margins and safety. A spacer grid that supports the fuel rods in a fuel assembly is equipped with mixing vanes that play a role in improving the heat transfer from the hot surfaces of the fuel rods to the coolant flow as the turbulenceenhancing devices. As shown in Fig. 11, flow patterns generated by a spacer grid with mixing vanes generally consist of both the swirl flow inside subchannel and the cross flow at the fuel rod gaps. The swirl flow improves the heat removal at the fuel rod surface by mixing the hot water near the fuel rod surface with the relatively cold water at the subchannel center, while the cross flow contributes to mitigating the hot peaking of subchannels by exchanging the enthalpy between subchannels [9]. Therefore, the geometrical shape and arrangement of the mixing vane are important factors that determine the performance of the mixing vane. Because a spacer grid may cause rigorous mixing as well as greatly increased local turbulence levels inside the sub-channel, prediction of sub-channel flows, even in isothermal condition, is very difficult. In general, subchannel analysis codes such as COBRA or VIPRE have been used to predict the flow and enthalpy distributions within fuel assemblies. However, these sub-channel codes rely on geometrically dependent mixing factors and empirical correlations to close the governing equations. The advantage of a CFD software for subchannel flow predictions is that it does not rely to the same extent on these empiricisms. Therefore, CFD results have the potential for wider applicability to capture the essential features of the turbulent structures downstream of the spacer grid.
In this study, in order to examine the flow distribution inside a fuel assembly with twist-split type mixing vanes as a part of benchmark simulation for a Coordinated Research Project (CRP) on the application of CFD codes to Nuclear Power Plant (NPP) design, simulations were conducted with the commercial CFD software, ANSYS CFX R18.1. The predicted results were compared with the measured data from the Omni Flow Experimental Loop (OFEL) test facility.

Analysis model
As shown in Fig. 12, OFEL test facility [10] was used to measure the flow distribution inside a fuel assembly with twist-split type mixing vane. Test rig consisted mainly of a water storage tank, a centrifugal pump, a flow meter, and a test section. Test section was made up of a square housing, fixing plate, rod bundle, rod support, and spacer grid with twist-split type mixing vane. Length and outer diameter (D) of a fuel rod, and rod-to-rod pitch (P) were 2,000 mm, 25.4 mm, and 34.29 mm, respectively. Therefore, P/D was 1.35 which corresponded to the regular rod-bundle configuration. As shown Fig. 13, velocity distribution inside the subchannel was measured at several cross-sectional planes by using either PIV (Particle Image Velocimetry) or LDV (Laser Doppler Velocimetry). Inlet temperature and pressure condition of the working fluid were measured at 35℃ and 0.1 MPa, respectively. The bundle-average axial velocity was estimated to be 1.5 m/s.

Numerical modeling
The flow inside the fuel assembly was assumed to be unsteady, incompressible, isothermal and turbulent. A high resolution scheme for the convection-terms-ofmomentum and -turbulence equations was used. 2nd Order Backward Euler scheme was used for the transient term. A time step of 0.001s was used with the maximum 10 iterations per time step. The solution was considered to be 'converged' when the residuals of variables were below 10 -5 and the variations of the target variables were small. Simulation was conducted with the commercial CFD software, ANSYS CFX R18.1. The ReNormalization Group (RNG) k- model was used to simulate the turbulent flow inside a fuel assembly. This model is based on renormalization group analysis of the Navier-Stokes equations. The transport equations for turbulence generation and dissipation are the same as those for the standard k- model, but the model constants differ. In general RNG k- model offers little improvement compared to the standard k- model. Fig. 14 shows the grid system for the computational domain that had the same size as the test facility. A hybrid mesh, made up of tetrahedrons, wedges, pyramids and hexahedrons, was generated to prevent the oversimplification of the geometry, and to have more efficient mesh distribution. Prism layers were used to get higher resolution in the near-wall region. Total numbers of nodes 4.010 7 . Because the flow area at the inlet of test section is different from that in the flow passage, uniform velocity, calculated by considering the above-mentioned condition, was used as an inlet-boundary condition. The 'average pressure over the whole outlet' option; with a relative pressure of 0 Pa, was used as an outlet-boundary condition. A no-slip condition was applied at the solid wall. To model the flow in the near-wall region, the scalable wall function method was applied. Fig. 15 shows the comparison of axial velocity profiles, measured and calculated along the z-axis (downstream of mixing vane) at both the subchannel center and the gap center of rod-to-rod. Both positions were indicated by the solid circles in Fig. 13. As shown in Fig. 15(a), defect (i.e. minimum) of the calculated axial velocity close to the mixing vane tip in the subchannel center agreed well with measurement. As the flow moved further downstream, the calculated axial velocity was slightly overestimated in the central sub-channel; compared with measurement. Additionally, at the gap center of rod-to-rod, though there was a little underestimation in the velocity magnitude, the calculated axial velocity showed qualitatively similar profile in comparison with measurement (see Fig. 15(b)).    to the mixing vane tip, a large elliptic vortex was generated in the subchannel center and two small secondary vortex appeared in the peripheral region near the upper and lower gaps between rods (see Fig. 17(a)). At the location Z/D = 3.0, major axis of a large elliptic vortex in the subchannel center rotated about 45 counter-clockwise and two small secondary vortex disappeared (see Fig. 17(b)). Finally, at the location Z/D = 6.0, the shape of central vortex became nearly circular (see Fig. 17(c)). The calculated velocity vector showed similar flow patterns in comparison with measurement.

Conclusions
In this study, regulatory audit calculation activities on the applicability of CFD software to nuclear safety problems was briefly explained. The major conclusions can be summarized as follows.
• Asymmetrically entering flows due to RCP performance differences were properly mixed while passing through the internal structures installed at the lower part of the reactor, confirming that the core inlet flow distribution is equivalent to that under the symmetric flow condition.
• The strongest thermal stratification was found at the location where the cold cooling water injected at the inlet of the ECCS injection line A contacted with the lower part of the cold leg A. Though there were a little difference between the measured and computational results, the predicted fluid temperature was relatively in good agreement with the measured data during the ECC water injection period. • As the flow passing through a mixing vane moved downstream, vortex shape in the subchannel center changed from elliptic to circular type. The calculated velocity vector showed similar flow patterns in comparison with measurement. • Although state of the art CFD software may solve reasonably these nuclear safety problems, it still has the limitation and a certain level of uncertainty in the calculation result. Main reasons may be grid quality, turbulence model and so on. • Therefore, lessons learned from these audit calculations are being used to prepare the domestic regulatory guideline on the applicability of CFD software.

Disclaimer
The opinions expressed in this paper are those of the author and not necessarily those of the Korea Institute of Nuclear Safety (KINS). Any information presented here should not be interpreted as official KINS policy or guidance.