Abstract
The thermal hydraulic and neutronics coupling analysis is an important part of the high-fidelity simulation for nuclear reactor core. In this paper, a thermal hydraulic and neutronics coupling method was proposed for the plate type fuel reactor core based on the Fluent and Monte Carlo code. The coupling interface module was developed using the User Defined Function (UDF) in Fluent. The three-dimensional thermal hydraulic model and reactor core physics model were established using Fluent and Monte Carlo code for a typical plate type fuel assembly, respectively. Then, the thermal hydraulic and neutronics coupling analysis was performed using the developed coupling code. The simulation results with coupling and noncoupling analysis methods were compared to demonstrate the feasibility of coupling code, and it shows that the accuracy of the proposed coupling method is higher than that of the traditional method. Finally, the fuel assembly blockage accident was studied based on the coupling code. Under the inlet 30% blocked conditions, the maximum coolant temperature would increase around 20°C, while the maximum fuel temperature rises about 30°C. The developed coupling method provides an effective way for the plate type fuel reactor core high-fidelity analysis.
1. Introduction
The high-fidelity simulation for reactor core requires more physics field coupling due to the strong feedback existing in the nuclear reactor, such as the reactivity interaction between core neutronics transport and thermal hydraulics. The neutron flux distribution determines the core power distribution and further determines the temperature distributions of fuel and coolant, which would affect the neutron fission cross sections, leading to the variation of core power distribution. Due to the strong feedback effects, the thermal hydraulics and neutronics research should be coupled together to achieve the accurate performances of coolant and fuel assembly in nuclear reactor core, especially in some unexpected accident conditions due to the large temperature variation. In the early years, neutron transport was often calculated with simplified point reactor neutron kinetics model due to the limited computing capability, and the results were usually conservative. With the development of computational capability, it has been feasible to realize the multiphysics coupling analysis for nuclear reactor high-fidelity simulations.
Actually, several thermal hydraulic and neutronics coupling analyses have been performed in the past several decades. Ji et al. coupled Monte Carlo code MCNP5 and system code RELAP5 to explore the feasibility of the pseudomaterial method for Doppler feedback in the VHTGR. In this study, the cross section libraries were interpolated to obtain the neutron cross sections at arbitrary temperature [1]. Hu and Uddin developed a new explicit coupling scheme using User Defined Function (UDF) to couple MCNP5 and Fluent, and some promising results were achieved [2]. Cardoni realized MCNP5 and STAR-CCM+ coupling to calculate 3D PWR fuel pin cell model and 3 × 3 model and obtained high-fidelity results of coolant temperature and density, fuel temperature, and power distribution, successfully demonstrating the coupling effects [3]. Yan et al. coupled STAR-CCM+ and deterministic code DeCART to calculate PWR 3 × 3 model with new mesh mapping methods [4]. Hoogenboom et al. built a flexible thermal hydraulic and neutronics coupling scheme with Python [5]. Sjenitzer et al. proposed a new method to perform transient coupling with Monte Carlo code [6]. Ward et al. used coupled system code RELAP5 and 3D spatial kinetics code PARCS to simulate safety performance of the I2S-LWR plant during accident conditions [7]. The thermal hydraulic and neutronics coupling algorithm in transient problems for the high-temperature gas-cooled reactor simulator TINTE was evaluated and developed [8]. Results indicated that the proposed coupling algorithms Picard and JFNK were better compared with the original semi-implicit coupling algorithm in TINTE.
In terms of plate type fuel reactor core analysis, Tian et al. developed a thermal hydraulic analysis code for China Advanced Research Reactor (CARR). The heat transfer and flow distribution characteristics in the reactor core were studied [9]. Gong et al. found that all fuel and cladding temperatures were below the design limits and remained safe with the inlet velocity ranging from 4.5 m/s to 7.5 m/s for the hottest assembly in 20 MW plate type fuel reactor [10]. Lu et al. analyzed the blockage accident of one assembly channel in 10 MW plate type fuel reactor, and the results showed that the obstructed channel causes temperature rise in adjacent channels. In addition, it was found that there was no boiling in obstructed channel except for the hot channel because of lateral heat conduction of adjacent channels [11]. Bousbia-Salah et al. calculated the neutron flux and power distribution of 10 MW MTR using MCNP5 code, and it showed that the MCNP5 was reliable for the plate type fuel core simulation and the calculation results were in good agreement with previous study [12]. Xoubi et al. investigated the impact of enrichment on neutron flux in the in-core facility of 10 MW MTR with OpenMC and found the importance of flux trap calculation while considering the conversion of reactor core from HEU to LEU [13]. It also demonstrated the feasibility of the Monte Carlo method in analyzing plate type fuel reactor.
It can be seen that most of the thermal hydraulic and neutronics coupling studies are performed for PWR rod bundle type fuel in the literatures, while the coupling analysis on plate type fuel reactor core is rare. Therefore, an innovative thermal hydraulic and neutronics coupling method was proposed based on the Fluent and Monte Carlo code through the UDF module for the plate type fuel reactor core in this paper. A typical plate type fuel assembly was analyzed with the coupling code, and the thermal hydraulic and neutronics features were studied and analyzed in detail.
2. Mathematic Model
The CFD method is widely utilized in the nuclear reactor thermal hydraulic analysis currently, especially for the local detailed three-dimensional flow and heat transfer process, such as the coolant mixing in the reactor core [14]. The basic of CFD method is the fundamental governing equations of fluid dynamics, including the mass, momentum, and energy conservation equations. For single-phase flow, the mass and momentum equations are given as follows:where the u is the velocity, ρ is the density, F is the body force per unit of fluid mass, and σ is the stress tensor.
For incompressible Newtonian fluid, the Navier–Stokes equations are given as follows:where the T is temperature, p is pressure, cp is the specific hear capacity at constant pressure, λ is thermal conductivity, is volumetric heat source, and is the dissipation function. For the steady calculation in the paper, the time derivative term is eliminated.
For solid domains like fuel and cladding, heat conduction equation can be simplified from equation (4) as follows:
For fuel, the volumetric heat source is fission heat or the heat generated by γ-rays in cladding. In this paper, the heat in cladding is ignored and is zero in cladding domain.
The direct numerical solution (DNS) method directly solves Navier–Stokes equations, and huge computing cost is required to simulate turbulence effect, which is unfavorable in reactor scale simulation. The computing cost could be decreased by applying turbulence models. For steady-state simulation, the Reynolds-averaged Navier–Stokes-based (RANS-based) models, which describe time-averaged motion of fluid flow, are widely used in the engineering fields. The typical RANS-based models include the Reynolds stress model (RSM) and the models that introduce eddy viscosity hypothesis such as the k-ε model and k-ω model. The standard k-ε model is the most common turbulence model used in CFD, and its equation is described below as an example of RANS-based model [15]:where μt is turbulent viscosity and it is modelled as , Gk is turbulent kinetic energy produced by laminar velocity gradient, Gb is turbulent kinetic energy produced by buoyancy, and C1ε, C2ε, Cμ, σk, and σε are model constants.
The Monte Carlo method is a stochastic algorithm based on probability statistics. The neutron transport in nuclear reactor is random with certain statistical properties, which fits the scope of application of the Monte Carlo method well. For the application of the Monte Carlo method in neutron transport, it simulates several histories of neutron from generating to disappearing continuously, instead of solving neutron transport equations directly.
With large numbers of neutron samples, many average parameters such as neutron flux distribution, keff, energy in the reactor, and the confidence interval are obtained according to central limit theory. For critical problem in Monte Carlo code, the fission source begins with a typical history and is converged after enough neutron generations.
3. Coupling Method
The coupling method between thermal hydraulic code and neutronics code is mainly divided into external coupling and internal coupling. The internal coupling integrates the two codes together and the data transfer process is achieved inside the integrated code. Although this method has relatively high accuracy and performs well in parallelism, it requires massive modification on the codes. The external coupling is achieved through user interface module to transfer data between two codes without modification on the codes. Therefore, the two codes are able to remain independent during external coupling process. In this paper, the UDF, which is a powerful user interface of Fluent for secondary development [16], is compiled to couple Fluent and neutronics codes and realize the transfer data between the two codes. The Fluent meshes are mapped with Monte Carlo code cells through UDF module. The Monte Carlo code provides heat source term to the Fluent meshes, while the Fluent results renew the temperature and density for Monte Carlo code, leading to the update of cross sections in the input file. The detailed coupling scheme is presented in Figure 1. The fuel domain in Fluent is provided with power distribution and the flow field is initialized, and then the CFD calculation begins. When the residual of continuity in Fluent is under 10−4, the calculated temperature and density are extracted. The Monte Carlo code input file is updated with the data and the cross section libraries are renewed. Then, the Monte Carlo code calculation begins. After Monte Carlo calculation, the fission energy deposition in the output file is extracted and transformed into power density, in which case the Fluent source term would be updated. The coupling code finally determines whether the iteration converges by monitoring keff.
In this paper, the geometry of plate type fuel assembly is simple and the one-to-one node mapping method optimized for plate type fuel model is adopted during the code coupling scheme. The mapping scheme is achieved through UDF traversing the coordinates of models in Fluent and Monte Carlo code and matching the cell ID, which makes it available for mapping considerable number of meshes, providing high accuracy at the cost of a little longer time. The mapping strategy also makes it convenient to modify the code to fit new plate type fuel model. The heat source term in Fluent is calculated by Monte Carlo code, and thus Monte Carlo code would not provide the power distribution directly. The tally module in Monte Carlo code is necessary for the fission energy deposition, and it is transformed into power density according to the following equation:where qv,i is the power density of mesh i, P is the total fuel power, Di is the average fission energy deposition in cell i, and Vi is the volume of cell i.
The initial heat source is provided by Monte Carlo code, and the new temperature and density are updated in the Monte Carlo code input file at each time step of Fluent calculation convergence. Then, the cross section and thermal S (α, β) libraries are updated according to new temperatures. There are mainly three methods to update the libraries: (a) update the libraries corresponding to the exact temperatures; (b) generate libraries with temperature interval of 2∼5 K with NJOY in advance and perform grouping approximation based on the calculated temperature; and (c) generate libraries with larger temperature interval and perform interpolation based on the calculated temperatures. The calculation cost of method (a) is too high. Method (b) costs less but still requires high memory. Compared to method (a), the computing cost of method (c) is less and the deviation of keff is just about 30 pcm, and finally method (c) is selected. The cross section libraries with 25°K temperature interval and thermal S (α, β) libraries with 50°K temperature interval are generated using NJOY. The cross section at any temperature for each cell is obtained by interpolating the libraries of adjacent temperature points. The details of the interpolation method are shown as follows:where TH is the highest temperature of selected interval, TL is the lowest one, Σ is the macrosection, and fL is the portion of lower temperature library.
4. Thermal Hydraulic and Neutronics Coupling Analysis
An active zone model of U3Si2-Al plate type fuel assembly in a typical reactor core is built in the paper. The model structure and cross sections are presented in Figure 2. The detailed parameters are presented in Table 1, and properties of fuel pellet and cladding are presented in Table 2 [17]. The variations of density, heat conduction coefficient, and coolant dynamic viscosity with temperatures are presented as follows:
(a)
(b)
The unit in equations (9) to (11) and Table 2 is K. All the temperature-dependent parameters are hooked by UDF.
Some reasonable assumptions are taken for simplifications: (a) symmetrical boundary conditions are set except the inlet and outlet; (b) the inlet and outlet planes are constant pressure surface; (c) no gap exists between fuel pellet and cladding; and (d) the impact of oxidation film on the fuel plate is ignored. The operation pressure is 0.683 MPa, and the inlet temperature is 308 K. Prior to the coupling analysis, grid and turbulence model independent sensitivities are performed and the core neutronics model is validated.
Totally five grid schemes are generated to perform the constant power steady-state calculation. The variations of pressure drop and temperature at outlet with grid numbers are presented in Figure 3. It can be seen that the scheme with 95256 meshes could be regarded as the grid-independent solution.
Three turbulence models, realizable k-ε, standard k-ω, and RSM are selected for sensitivity analysis to determine the most suitable turbulence model. The calculated coolant temperature and pressure are presented in Figure 4. It can be seen that the turbulence model selection has very limited influence on the calculation results, and finally the realizable k-ε model is selected in the paper.
(a)
(b)
The enrichment of 235U is 19.75% ± 0.20%, and uranium density in the fuel pellet is 4.3 gU/cm3. Nuclei component and design value of atomic density are presented in Table 3.
The meshing model in Monte Carlo code and Fluent is similar. The inlet and outlet planes are set as the vacuum boundary, and other boundaries are set as reflection. Totally 120 batches are run at one iteration, and each batch has 10000 neutron histories. Before counting, 20 neutron batches are skipped for convergence of fission source to reduce the final deviation. The libraries of neutron fission cross sections in fuel and neutron scattering cross sections in coolant are shown in Tables 4 and 5.
The reactivity coefficient calculation is performed to validate the feasibility of built neutronics model. The eight groups of calculations are carried out to determine the coolant temperature coefficient and fuel temperature coefficient. Group details are presented in Table 6, and the results with 95% confidence interval are presented in Table 7. The liner fitting is done after transforming keff to reactivity coefficient, as shown in Figures 5 and 6. The average fuel temperature coefficient is −2.086 × 10−5 (K), and average coolant temperature coefficient is −5.33 × 10−5 (K). The two coefficients documented in the literature are −2.2745 × 10−5 (K) and −8.0734 × 10−5 (K), and it can be seen that they are close to the calculated value in this paper. It demonstrates that the neutronics model and cross section settings are reliable.
4.1. Normal Operation Condition
The coupling and noncoupling simulations were performed, respectively, and the results are compared. For noncoupling simulation, cosine power distribution is applied as the heat source. In this case, the inlet velocity is set to 7.0 m/s, and the total assembly power is set to 7.85 MW. keff, outlet temperature, and pressure of each iteration are shown in Figure 7. It can be seen that after 5 iterations, the calculation is generally converged.
Figure 8 shows the coupled power density distribution in the whole assembly. The power distribution in a single plate along the thickness is treated as uniform distribution since its dimension is much smaller compared to that in the height and width direction. Figure 9 shows the comparison of power distribution in the middle plate between coupling analysis and noncoupling analysis. The calculated power with coupling analysis code is larger near the edge in the width direction compared to that with noncoupling analysis method due to the space self-shielding effect, which enlarges the fission rate near the edge. Figure 10 shows the variation of power distribution with the iterations, from the initial shape to the 5th iteration. The powers increase at the both ends, and power peak is lowered with the coupling analysis. It can be seen that the power is a little higher in the inlet side, while it is lower in the outlet side, and the power peak moves to the inlet side after coupling. It is because the coolant temperature is lower and the density is higher in the inlet side, where the moderation effect is stronger than that in the outlet side, making the thermal neutron flux higher and fission rate greater. The variation of power distributions before and after coupling shows that the coupling code introduces the feedback effect between thermal hydraulic and neutronics and the feasibility of coupling code is proved.
(a)
(b)
Figure 11 shows the comparison of coolant temperature distributions in the assembly before and after coupling analysis. Figure 12 provides the quantitative comparison of temperatures of coolant, cladding, and fuel pellet between the noncoupled and coupled analysis. The temperature peak with coupled analysis is lower than that with noncoupled analysis, indicating that cosine power distribution assumption is conservative. Figure 13 shows the variation of average heat transfer coefficient in the height direction. It increases at the beginning rapidly and descends along the flow direction gradually.
(a)
(b)
4.2. Blockage Condition Analysis
For the plate type fuel reactor, under certain accident conditions such as debris flowing into reactor and fuel blistering, the flow channel blockage will happen, causing temperature to increase sharply inside the fuel plate and leading to large temperature gradient along the plate, which may induce structure rupture of plates and cause severe consequences [18]. The fuel assembly operation features under the blockage conditions are analyzed with the coupling code in this section.
Totally three positions at the inlet of assembly are modelled as the solid partly to realize the fuel assembly blockage conditions. The cross sections of inlet blockage model are presented in Figure 14, and the 30% blockage of channel is set. The material is set as steel, and the corresponding Monte Carlo code coolant cells are changed into steel cells.
(a)
(b)
(c)
Figure 15 presents the coolant temperature in the whole assembly, and Figure 16 shows the temperature contours at z = 0.285, just behind the obstructed location. It can be seen that blockage causes temperature rise in both coolant and fuel plates. The temperatures of fuel plates along x = −0.025 are presented in Figure 17. The third blockage condition causes higher temperature rise due to the reflective boundary set in the Monte Carlo code. Figure 18 shows the x-direction temperature of fuel plates pointed with the black arrow in Figure 16. The large temperature gradient could be observed in all three conditions, which would threaten the fuel mechanical behaviors.
As thermal hydraulic and neutronics coupling effect is introduced in the calculation, the effect of blockage on local heat transfer can be more accurately obtained. The thermal hydraulic parameters in highlighted position (as shown in Figure 19) are studied. The local temperatures of fuel pellet, cladding, and coolant under three blockage conditions and normal conditions are presented in Figure 20. For the 30% blockage in one channel, the coolant velocity decreases, causing local heat transfer coefficient to decrease greatly. As shown in Figure 21, the heat transfer coefficient generally decreases by 5 × 103 W/(m2·K), and the maximum decrease is over 104 W/(m2·K), leading to coolant temperature increase around 20°C at the outlet and fuel maximum temperature increase about 30°C.
5. Conclusion
In this paper, a thermal hydraulic and neutronics coupling scheme was proposed for plate type fuel reactor core. The coupling platform is based on the Fluent and Monte Carlo code through the UDF module. Then, the thermal hydraulic and neutronics coupling analysis for the plate type fuel assembly was performed in detail. The mesh-independent analysis and turbulence model sensitivity analysis were carried out to determine the number of meshes and turbulence model firstly. Results show that the accuracy of power distribution is well improved and temperature distributions of fuel pellet, cladding, and coolant are refined due to the feedback effect introduction. The power density distribution in the width direction is not uniform, and the power near the edge is higher than that at inner position. Following the normal operating condition study, the blockage condition analysis with the coupling code was performed. Under the conditions of 30% blockage in one channel, the heat transfer coefficient decreases obviously. The maximum coolant temperature would increase around 20°C, and the maximum fuel temperature rises about 30°C. This work provides a promising analysis tool for the plate type nuclear reactor core high-fidelity simulation.
Data Availability
The access to data is restricted due to the commercial confidentiality.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (no. 11705139).