Development and Validation of Multiscale Coupled Thermal-Hydraulic Code Combining RELAP5 and Fluent Code

In nuclear reactor system research, the multiscale coupled thermal-hydraulic (T-H) system code and CFD code is one of the most prevalent research areas, and it could help improve simulation fidelity and optimize nuclear reactor design. Additionally, a new idea known as the function fitting method (FFM) for coupling parameter distribution has been newly proposed for exchanging data on the coupling interface, which uses math equations to present the velocity distribution characteristics at the coupling interface. This method could improve the simulation error and numerical instability. To verify and validate the abovementioned FFM, a comparison between the velocity function shape by FFM and real velocity distribution was completed. Besides, the Edwards pipe blowdown test results were used to verify the coupled code. The results showed good agreement with experiment results, and a better simulation accuracy compared to previous work. The current work will establish the ability to explore multiscale coupled thermal-hydraulic operation characteristics which permit precise local parameter distribution.


INTRODUCTION
Nuclear safety is a top priority for nuclear power application and expansion. Estimate tools for increasing safety analysis and evaluation requirements need to become better and more precise. In the past few decades, on a system scale, the best estimate codes such as RELAP5 (Allison et al., 1993), RETRAN (McFadden et al., 1981), CATHARE (Barre and Bernard, 1990), and MARS (Jeong et al., 1999) have dominated nuclear reactor operation, safety analysis, and severe accident analysis. However, these codes can only present one-dimensional transient system behaviors, which can not provide the local characteristics of the reactor. As computational resources have dramatically developed, component scale analysis codes like COBRA (Stewart et al., 1977), RELAP5-3D (RELAP5-3D Code, 2012), VIPRE (Stewart et al., 1989), and local scale codes such as Fluent (Rohde et al., 2007), CFX (Höhne et al., 2010), and Star-CCM+ (Cardoni, 2011) have emerged. These computational fluid dynamic (CFD) codes can provide three-dimensional features, which have been applied in pressurized thermal shock (PTS) (Egorov et al., 2004), boron dilution and distribution in reactor vessels (Muhlbauer, 2003;Scheuerer et al., 2005), and so on. The above CFD application has faced the challenge of the computational cost of transient safety analysis. To conquer these difficulties, a compromised method of developing coupled codes based on system and CFD codes has been promoted.
Ronnie Andersson (Andersson et al., 2004) utilized a system and CFD code to analyze turbulence intensity and dissipation on multiphase flow in the reactor. Daniele Martelli (Martelli et al., 2017) coupled RELAP5/MOD3.3 and Fluent to simulate a NACIE experiment loop and loss of flow accident, which was in good agreement with experiment data. J-J. Jeong (Jeong et al., 1999) developed an integral modular code MARS coupled RELAP5 and sub-channel COBRA code, which adopted a semi-implicit method and dynamic memory allocation method. W. L. Weaver Weaver, 2005) developed a series of studies coupling RELAP5-3D and CFX, D.L. Aumiller (Aumiller et al., 2001;Aumiller et al., 2002) further coupled RELAP5-3D and CFD-FLOW3D based on the parallel virtual machine (PVM) method, both studies used the Edwards blowdown test to conduct the verification and validation (V&V) process. Nolan Anderson (Anderson, 2006;Anderson et al., 2008) also developed RELAP-3D and Fluent coupled code using PVM, for Very High Temperature Reactor (VHTR) lower plenum analysis.
The existing coupled code development technologies mainly include PVM, dynamic link library (DLL) (Li et al., 2014), and boundary files modification methods. Challenges in the coupling process (Ivanov and Avramova, 2007) mainly concern the method of coupling, coupling approach, spatial mesh overlays, time step algorithms, and coupling numeric and convergence schemes. Spatial mesh mapping, especially the data exchange at the coupling interface plays a key role in simulation accuracy and numerical convergence. Therefore, in this paper, a new method called the function fitting method (FFM) for coupling parameter distribution was proposed for coupling RELAP5 and Fluent, aiming at providing precise data exchange. An Edwards pipe blowdown test was used to verify and validate the coupled code.

COUPLING METHOD
For multiscale coupled thermal-hydraulic code, the kinds of variables that are transferred through the coupling interface must be considered as priority. The RELAP5 code sets the boundary conditions through time-dependent control volume (TMDPVOL) and time-dependent junction (TMDPJUN); while Fluent has a predetermined velocity-inlet, pressure-outlet, and outflow boundary conditions, etc.
Assuming that the control volume Ri (RELAP5-interface) is connected to the coupling interface between Fluent and RELAP5, the mass and energy conservation equations are solved in the control volume Ri. The junction j is also connected to Fluent and RELAP5, then the momentum conservation equation is solved in junction j, and the vector variables are stored in j (as shown in Figure 1). In fact, the RELAP5 interface control volume Ri and the junction j do not have realistic components in the RELAP5 system; its main function is to be used as a coupling interface, which is also known as a "ghost cell." This coupling interface is an overlapped computational domain.
Under this scenario, the downstream Fluent computational mesh is regarded as the downstream control volume of the RELAP5 portion, and the abovementioned mass, energy, and momentum conservation equations can be derived in the following forms: where, P is the pressure and v is the fluid velocity. While, the subscript g and f stand for gas and fluid, respectively. The superscript n and n + 1 represent the current and next time step. Coefficients b, g 1 , g 2 , f 1 , and f 2 are column vectors, these coefficient vectors are known variables under the current time step n. b is the source term in the equation, and the coefficient matrices g and f represent the convection effect. A′, B′, C′, and D′ are coefficients that contain only the current time step variable.
The exchange variables between RELAP5 and Fluent codes are listed in the following Table 1.
The data transfer between the coupling interface is bidirectional. While, Fluent passes parameters to RELAP5, which can be calculated by surface summation or averaging. For the mass flow rate, the sum of the Fluent surface cells is equal to that of the RELAP5 interface. where, W is the mass flow rate and N c stands for the Fluent cell numbers on the coupling interface.
Other variables like temperature, pressure, and void fraction, are represented by Φ, and can be described as follow: where, A i is the area of cell surface i. For parameters passed from RELAP5 to Fluent, onedimensional parameters obtained by the lumped parameter method of RELAP5 will be converted into the twodimensional distribution of the Fluent interface. In past research, these data were generally averaged on the interface, but this method introduces errors into each calculation iteration. While, in some studies, the coupling boundary is set further afield so that it is far enough to minimize uneven effects on the interface. However, these two methods do not solve the different dimensional transformation problem of the coupling parameters. Therefore, in this paper, we worked on solving this issue by proposing a function fitting method for interface data.

FUNCTION FITTING METHOD
For nuclear power systems and equipment, most fluid areas are round tubes, such as pipes and fluid machinery (pumps, valves, etc.,). For the application of the RELAP5/CFD coupled code, the primary system and safety system are modeled by onedimensional system codes to obtain the transient characteristics. However, the pressure vessel, lower plenum, and the core are modeled by the three-dimensional CFD code Fluent. The connection between the system and the local equipment are mostly long round tubes, and the fluid flow therein can be considered fully developed. The flow and heat transfer of single-phase and two-phase fluids for round tubes have been widely studied. Therefore, in this paper, the function fitting method (FFM) for coupling parameter distribution was proposed. This method can accurately convert the onedimensional lumped variables into two-dimensional ones that satisfy the corresponding physical laws, thereby effectively improving the error of calculation instability and accelerating calculation convergence.

Fitting Function
The original intention of this method was to transfer the onedimensional parameters into the two-dimensional surface ones through appropriate function fitting, and accurately reflect the real distribution. Therefore, for Newtonian viscous fluids, there are many mature empirical formulas for one-dimensional flow, and they can be used as one of the independent variables of the fitting function.
For the fully developed single phase turbulent flow in a round tube, the velocity distribution is relatively flat in the middle.
While, in the viscous bottom layer near the wall surface, the velocity distribution changes sharply and the velocity gradient is relatively large. For the velocity distribution in the tube flow, the most influential factor is the Reynolds number and position. The flow velocity distribution function fitting method mainly considers the influence by friction coefficient and relative position.
The empirical formulas for the widely recognized turbulent flow frictional resistance coefficient are listed in Table 2.
Rr in Table 2 is the relative roughness. The friction coefficient calculations for the Nikuradse and Colebrook models in the table are implicit, that is, f appears on both sides of the equation; therefore, either the equation is iteratively solved, or the solution is interpolated in the Moody diagram (Moody, 1944). These solution methods are very inconvenient (Haaland, 1983). Therefore, in the function fitting method, it is preferred to select the simple and accurate explicit friction coefficient calculation as one of the independent variables. After a large number of experiments, the Filonenko model was chosen because it is applicable in a wider range of Reynolds numbers. Besides, in related research, the experiments of scholars Romeo (Romeo et al., 2002), Fang (Fang et al., 2011), and Yıldırım (Yıldırım, 2009) have verified that this model is more accurate among the explicit friction relationship.
The function fitting method is also related to the position of the one-dimensional grid to the two-dimensional grid spatial mapping. Therefore, for the round tube coupling interface, the radial position of the round tube is selected as another independent variable. According to the characteristics of the flatness of the central cross-section, and the steep edge, the fitted mathematical function should also have the above characteristics.
The simplest and most straightforward method of the fitting function is to use piecewise polynomial fitting. However, the function obtained by this method requires multiple constraints, which is not suitable for practical application. In elementary functions, logarithmic function has the characteristics of rapid transition in the normalized (0-1) interval. Therefore, the fitting function must contain a logarithmic function term, which also requires the consideration of relative distance.
The term related to the friction coefficient is in the form of the polynomial and power function. The determination of the index needs to be verified by a large number of calculations, so that the obtained function conforms to the flow velocity distribution under different flow Reynolds numbers.
Finally, the conversion formula for the function fitting method is: where, F is the fitting function, r is the distance from the central axial line, and R is the inner radius. The f in Eq. 7 is the friction coefficient. Here, the Filonenko model is selected, namely: The coupling parameter velocity conversion from the system code RELAP5 to the CFD code Fluent can be expressed as: where, u CFD is the fluid velocity at the Fluent coupling interface, while u RELAP5 is fluid velocity at the RELAP5 coupling junction. In addition, for the coupling process between different thermal-hydraulic codes, attention should also be paid to the conservation of the coupling parameters. Especially for the coupling parameter distribution function fitting method, it is necessary to assure that the integral flow on the Fluent interface is equal to the total flow at the RELAP5 coupling junction, as follows: where, W j is the mass flow rate at the RELAP5 coupling junction and A is the area of the Fluent coupling interface.

Verification and Validation of Fitting Function
In order to verify and validate the accuracy of the function fitting method proposed in the previous section, Fluent was used to analyze the flow in a round tube under different conditions. Besides, the results have been compared with the fitting function under the corresponding Reynolds number. Since the independent and dependent variables of the fitting function are dimensionless, in the CFD verification, the corresponding size and speed are also normalized accordingly. A horizontal tube with a length of 500 mm and an inner diameter of 10 mm was chosen; and the velocity distribution at the center of the tube was compared with the fitting function. The selected pressure, temperature, flow rate, and corresponding Reynolds number were included within the pressurized water reactor (PWR) operation and accident conditions. The specific verification conditions are shown in Table 3.
In this V&V work, in order to accurately simulate the flow features in the boundary layer, the y plus value was confirmed as 30; the wall surface was divided into 20 boundary layers, the grid thickness of the first layer of the boundary layer was calculated to be 7 × 10 −4 mm, and the growth ratio was 1.1. The mesh cross section is shown in Figure 2.
In Figure 3, the normalized velocity distribution is compared between the fitting function and Fluent results. Under different pressures, temperatures, velocities, and Reynolds numbers, the fitting function represented accurate results to prove its applicability. The x-axis of the fitting function was in a relative position to the axial center, and the y-axis was the normalized velocity, which determined that this function was not constrained by specific size.

VERIFICATION AND VALIDATION OF THE COUPLED CODE
In this paper, the Edwards pipe blowdown test was chosen to validate the multiscale coupled one dimensional (1D) and three dimensional (3D) code.

Edwards Pipe Blowdown Test
The original intention of the Edwards pipeline blowdown test was to simulate the phase change process in the safety analysis of the water-cooled reactor, which is very similar to the coolant loss process of the PWR loss of coolant accident (LOCA). The experimental facility contained a heating pipe filled with water, and the pressure was maintained above the saturation point. Before the glass plate at the end of the pipe was broken, the pressure was adjusted to the required level. During the blowdown test, measurement gauge stations were set up along the pipeline, with which pressure, temperature, and density changes were measured.
Before the blowdown process, the horizontal pipe was filled with supercooled water (602.15 K, 7.0 MPa). The ambient temperature was 293.15 K, and the ambient pressure was atmospheric pressure. The pipe was 4.096 m long and had an inner diameter of 73.15 mm. The design pressure of the pipeline was 17.2 MPa, and the design temperature was 616.5 K. The initial experimental pressure ranged from 3.55 to 17.34 MPa, and the temperature ranged from 514.8 to 616.5 K. The experiment layout is shown in Figure 4. The end of the pipe was the blowdown section, sealed by a toughened glass disc, the diameter of the glass disc was 88.9 mm, and the thickness was 12.7 mm. At the beginning of the experiment, the glass disc was ruptured by the lead pellet from a compressed air gun, and the water in the pipe was discharged into the environment. The reduction in the spray area accounted for 13% of the pipe cross-sectional area. A total of seven measuring gauge stations (GS1-GS7) were set up along the axial direction of the pipeline, and each was equipped with a facility for fitting fast response pressure and temperature transducers. The pipe was heated electrically with sheathe band heaters, formed to the curvature of the pipe. They covered approximately 70% of the pipe circumference and each had a capacity of 700 watts. Heat losses were reduced by using asbestos insulation. In order to verify the RELAP5/Fluent coupled code, the Edwards pipe blowdown test was used as the benchmark problem in this paper. Meanwhile, the standalone RELAP5 code and other scholar's coupled code results (Li et al., 2014) were compared with our work.

Simulation Model
The standalone RELAP5 model nodalization is shown in Figure 5, in which the pipe component 003 PIPE stands for the experiment pipe, and the single junction component 004 SNGJUN connects the experiment pipe and rupture boundary. The time-dependent volume component 005 stands for the ruptured end of the pipe and also provides boundary conditions for the simulation. In RELAP5, the bubble radius change rate in the flash evaporation process is based on the Plesset-Zwick model (Plesset and Zwick, 1954), and the corresponding Nusselt solution is based on the Lee-Rypley model (Lee and Ryley, 1968): where, the subscript b represents bubbles.
In the coupled code analysis, the Edwards pipeline simulation model was divided into two parts along the axial direction. The upstream was modeled by Fluent, and the downstream to the end of the pipeline was simulated by RELAP5 (as shown in Figure 6). The downstream RELAP5 part contained 10 mass and energy control volumes, and nine momentum junctions among each control volume. Since the Edwards blowdown test is a strong transient process, the maximum time step of RELAP5 was set to 0.0001 s. In each time step, the inlet boundary of RELAP5 provided the pressure, temperature, and void fraction to the outlet boundary of Fluent; the Fluent returned the void fraction, pressure, temperature, and mass flow rate of the gas and liquid phase inlet boundary condition of RELAP5.
For the Fluent calculation part, the two-phase Euler model was adopted for the multiphase flow simulation, and a bubble diameter d bubble 1 mm. The thermal phase change model was adopted for the rapid and intense evaporation-condensation process (ANSYS, 2015). The water properties were based on the National Institute of Standards and Technology (NIST) compressible gas model, which was important to capture the spread of pressure waves.

RESULTS AND ANALYSIS
The Edwards pipe blowdown test process can be divided into two periods: one is the rapid pressure discharge period caused by the single-phase water loss, and the other is the slower pressure discharge period by the loss of the two-phase mixture. The first  period happens almost instantaneously, at approximately 2 ms. The second period continues until the pressure is discharged close to the environment pressure. Figure 7 shows the pressure variation of the GS-5. After the blowdown started, the pressure dropped sharply due to the rapid loss of relatively high-density supercooled water until saturation conditions were reached. The pressure drop during the discharge phase of the two-phase mixture was compensated by vapor generation. The pressure results simulated by the coupled Fluent/RELAP5 code was in good agreement with the Edwards pipeline test results (Edwards and O'Brien, 1970). The maximum deviation was 0.26 MPa, which corresponds to an error of 18.4% from the experiment results, which is acceptable for two-phase transient flow. The pressure results simulated by the coupled code was higher than the experiment and RELAP5 results. It is worth mentioning that the prediction of coupled analysis is not restricted by independent RELAP5. Compared with the coupled code with a simple boundary condition treatment from other papers, our results were closer to the experiment results. Figure 8 shows the pressure variation at GS-7. The coupled code results were within a reasonable range compared to the experiment results. In 0-250 ms, the simulation results of the coupled code were close to the experimental results, after which the simulation results were higher than the experiment ones. The maximum pressure deviation between the coupled code and the test results was 0.41 MPa, and the corresponding error was 28%. The pressure drop process simulated by the Fluent/RELAP5 coupled code and the standalone RELAP5 code began earlier than in the experiment. This phenomenon was partly due to the thermal phase change flashing model in the Fluent code. Its flashing rate was determined by the temperature difference between each phase and their corresponding saturation temperature, which led to a faster pressure drop. Figure 9 shows void fraction variation at GS-5. In the first millisecond, the vapor generation rate was relatively slow, and the coupled code results were slightly lower than the experiment ones. As the blowdown process reached saturation conditions, the largest increase rate of the void fraction occurred in the middle of the whole process. In the following period, due to the equilibrium state with the environment, the vapor generation rate decreased. The error between the coupled code results and the experiment ones can be attributed to the assumption that the average bubble size was simplified to a fixed value (1 mm) in the Fluent model. The size of the pre-existing nuclei in the sub-cooled liquid was substantially smaller than the prescribed 1 mm. This indicated a significant under-prediction of the interfacial area density available for the initiation of flashing. While the steam bubbles grew rapidly during the flashing process, which resulted in a mean size much larger than 1 mm in a short period (Liao and Lucas, 2017). Nevertheless, the uncertainty was weakened with the increase of the void fraction. The effects of prescribed bubble sizes were studied in Liao's paper (Liao et al., 2013). Figure 10 shows the mass flow rate variation of the control volume at the coupling interface. Since there is no available experiment data in the Edwards pipe blowdown test, the results of the coupled code were compared with those of the standalone RELAP5 code. The results of the coupled Fluent/ RELAP5 code showed similar trends with those of RELAP5. The mixture mass flow rate by the standalone RELAP5 code was  larger than that of the coupled code, and the falling tendency of the mass flow rate variation simulated by the coupled code began earlier. Figure 11 shows the mass flow rate variation at the break of the pipe. After the rupture of the end of the pipe, the discharge flow rate quickly reached the maximum value. The mass flow rate of the break simulated by the Fluent/RELAP5 coupled code was in good agreement with that of the standalone RELAP5, and the      falling trend of mass flow also showed a similar relationship with that of the pressure. Although there is a lack of experimental results for the break mass flow rate, the results in this paper were close to those of the standalone RELAP5. It also indicated that the function fitting method encountered smaller errors in the simulation, which weakened the error induced by the coupling parameter.

CONCLUSION
In this paper, a multiscale coupled thermal-hydraulic method was studied, and a coupled RELAP5/Fluent code was developed. To solve the problem of exchanging different dimensional thermalhydraulic parameters through the coupling interface, the function fitting method was proposed. The physical distribution of different parameters was described by mathematical functions. The Edwards pipe blowdown test was used to verify the multiscale coupled thermal-hydraulic method. The results show that the FFM can help simulate the strong transient process accurately, and it can also improve the calculation accuracy when compared to the previous study that used uniform parameters distribution at the interface.

DATA AVAILABILITY STATEMENT
Due to the intellectual property limitations, the simulation data could not be published to the readers.