Multiphysics Multiscale Coupling Modeling for Nuclear Reactor and Its Uncertainty Quantification
View this Special IssueResearch Article  Open Access
Linrong Ye, Mingjun Wang, Xin’an Wang, Jian Deng, Yan Xiang, Wenxi Tian, Suizheng Qiu, G. H. Su, "Thermal Hydraulic and Neutronics Coupling Analysis for Plate Type Fuel in Nuclear Reactor Core", Science and Technology of Nuclear Installations, vol. 2020, Article ID 2562747, 12 pages, 2020. https://doi.org/10.1155/2020/2562747
Thermal Hydraulic and Neutronics Coupling Analysis for Plate Type Fuel in Nuclear Reactor Core
Abstract
The thermal hydraulic and neutronics coupling analysis is an important part of the highfidelity 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 threedimensional 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 highfidelity analysis.
1. Introduction
The highfidelity 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 highfidelity 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 STARCCM+ coupling to calculate 3D PWR fuel pin cell model and 3 × 3 model and obtained highfidelity results of coolant temperature and density, fuel temperature, and power distribution, successfully demonstrating the coupling effects [3]. Yan et al. coupled STARCCM+ 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 I^{2}SLWR plant during accident conditions [7]. The thermal hydraulic and neutronics coupling algorithm in transient problems for the hightemperature gascooled reactor simulator TINTE was evaluated and developed [8]. Results indicated that the proposed coupling algorithms Picard and JFNK were better compared with the original semiimplicit 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]. BousbiaSalah 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 incore 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 threedimensional 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 singlephase 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, c_{p} 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 steadystate simulation, the Reynoldsaveraged Navier–Stokesbased (RANSbased) models, which describe timeaveraged motion of fluid flow, are widely used in the engineering fields. The typical RANSbased 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 RANSbased model [15]:where μ_{t} is turbulent viscosity and it is modelled as , G_{k} is turbulent kinetic energy produced by laminar velocity gradient, G_{b} is turbulent kinetic energy produced by buoyancy, and C_{1ε}, C_{2ε}, 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, k_{eff}, 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 k_{eff}.
In this paper, the geometry of plate type fuel assembly is simple and the onetoone 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 q_{v,i} is the power density of mesh i, P is the total fuel power, D_{i} is the average fission energy deposition in cell i, and V_{i} 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 k_{eff} 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 T_{H} is the highest temperature of selected interval, T_{L} is the lowest one, Σ is the macrosection, and f_{L} is the portion of lower temperature library.
4. Thermal Hydraulic and Neutronics Coupling Analysis
An active zone model of U_{3}Si_{2}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 temperaturedependent 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 steadystate 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 gridindependent 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 ^{235}U is 19.75% ± 0.20%, and uranium density in the fuel pellet is 4.3 gU/cm^{3}. 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 k_{eff} 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. k_{eff}, 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 selfshielding 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 5^{th} 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 xdirection 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 × 10^{3} W/(m^{2}·K), and the maximum decrease is over 10^{4} W/(m^{2}·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 meshindependent 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 highfidelity 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).
References
 J. Conlin, W. Ji, J. C. Lee, and W. R. Martin, “Pseudo material construct for coupled neutronicthermalhydraulic analysis of VHTGR,” Thermal Hydraulics of NextGeneration Nuclear Reactors, vol. 92, pp. 225–227, 2005. View at: Google Scholar
 J. Hu and R. Uddin, “Coupled neutronics and thermalhydraulics simulations using MCNP and FLUENT,” Transactions of the American Nuclear Society, vol. 98, pp. 606–608, 2008. View at: Google Scholar
 J. N. Cardoni, Nuclear Reactor MultiPhysics Simulations with Coupled MCNP5 and STARCCM+, University of Illinois, Champaign, IL, USA, 2011.
 J. Yan, B. Kochunas, M. Hursin, T. Downar, Z. Karoutas, and E. Baglietto, “Coupled computational fluid dynamics and MOC neutronic simulations of Westinghouse PWR fuel assemblies with grid spacers,” in Proceedings of the 14th International Topical Meeting on Nuclear Reactor Thermalhydraulics (NURETH14), Toronto, Canada, September 2011. View at: Google Scholar
 J. E. Hoogenboom, A. Ivanov, V. Sanchez, and C. Diop, “A flexible coupling scheme for Monte Carlo and thermalhydraulics codes,” in Proceedings of the International Conference on Mathematics and Computational Methods Applied to Nuclear Science And Engineering, Rio de Janeiro, Brazil, May 2011. View at: Google Scholar
 B. L. Sjenitzer, J. E. Hoogenboom, J. Jiménez Escalante, and V. Sanchez Espinoza, “Coupling of dynamic Monte Carlo with thermalhydraulic feedback,” Annals of Nuclear Energy, vol. 76, pp. 27–39, 2015. View at: Publisher Site  Google Scholar
 A. M. Ward, M. J. Wang, M. D. Neumann, M. Memmott, A. Manera, and T. J. Downar, “A simulation of I2SLWR selected transients,” Annals of Nuclear Energy, vol. 145, 2020. View at: Publisher Site  Google Scholar
 H. Zhang, J. Guo, J. Lu, F. Li, Y. Xu, and T. J. Downar, “An assessment of coupling algorithms in HTR simulator TINTE,” Nuclear Science and Engineering, vol. 190, no. 3, pp. 287–309, 2018. View at: Publisher Site  Google Scholar
 W. X. Tian, S. Z. Qiu, Y. Guo, G. H. Su, and D. N. Jia, “Development of a thermalhydraulic analysis code for CARR,” Annals of Nuclear Energy, vol. 32, no. 3, pp. 261–279, 2005. View at: Publisher Site  Google Scholar
 D. X. Gong, S. F. Huang, G. B. Wang, and K. Wang, “Heat transfer calculation on platetype fuel assembly of high flux research reactor,” Science and Technology of Nuclear Installations, vol. 2015, Article ID 198654, 13 pages, 2015. View at: Publisher Site  Google Scholar
 Q. Lu, S. Qiu, and G. H. Su, “Development of a thermalhydraulic analysis code for research reactors with plate fuels,” Annals of Nuclear Energy, vol. 36, no. 4, pp. 433–447, 2009. View at: Publisher Site  Google Scholar
 A. BousbiaSalah, H. Benkharfia, N. Kriangchaiporn, A. Petruzzi, F. D’Auria, and N. Ghazi, “MTR benchmark static calculations with MCNP5 code,” Annals of Nuclear Energy, vol. 35, no. 5, pp. 845–855, 2008. View at: Publisher Site  Google Scholar
 N. Xoubi, S. A. Darda, A. Y. Soliman, and T. Abulfaraj, “An investigative study of enrichment reduction impact on the neutron flux in the incore fluxtrap facility of MTR research reactors,” Nuclear Engineering and Technology, vol. 52, 2020. View at: Publisher Site  Google Scholar
 H. Ju, M. Wang, C. Chen et al., “Numerical study on the turbulent mixing in channel with Large Eddy Simulation (LES) using spectral element method,” Nuclear Engineering and Design, vol. 348, pp. 169–176, 2019. View at: Publisher Site  Google Scholar
 V. H. Kaarle and M. Weeratunge, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Pearson Education, London, UK, 2007.
 X. Zhao, M. Wang, C. Chen et al., “Threedimensional study on the hydraulic characteristics under the steam generator (SG) tube plugging operations for AP1000,” Progress in Nuclear Energy, vol. 112, pp. 63–74, 2019. View at: Publisher Site  Google Scholar
 A. M. Zhang and Y. L. Kang, “Design of U3Si2Al platetype fuel element for China advanced research reactor,” in Proceedings Of the 18th International Conference On Nuclear Engineering (ICONE18), Xi’an, China, May 2010. View at: Google Scholar
 L. Li, D. Fang, D. Zhang et al., “Flow and heat transfer characteristics in platetype fuel channels after formation of blisters on fuel elements,” Annals of Nuclear Energy, vol. 134, pp. 284–298, 2019. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Linrong Ye et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.