Abstract
Heat transfer characteristics of fuel assemblies for a high flux research reactor with a neutron trap are numerically investigated in this study. Single-phase turbulence flow is calculated by a commercial code, FLUENT, where the computational objective covers standard and control fuel assemblies. The simulation is carried out with an inlet coolant velocity varying from 4.5 m/s to 7.5 m/s in hot assemblies. The results indicate that the cladding temperature is always lower than the saturation temperature in the calculated ranges. The temperature rise in the control fuel assembly is smaller than that of the standard fuel assembly. Additionally, the assembly with a hot spot is specially studied, and the safety of the research reactor is also approved.
1. Introduction
Research reactors are widely used to produce high neutron flux for research, training, education, and irradiation test. According to IAEA Research Reactor Data Base, there are 773 research reactors all over the world and 247 of them are being operated in different countries [1]. For certain purposes, larger central thermal neutron flux with a good quality factor is needed. For example, Chinese Mianyang Research Reactor (CMRR) in China Academy of Engineering Physics (CAEP) is such a research reactor with an inverse neutron trap. To obtain a high thermal neutron flux for generating fusion neutron through a neutron transformation target, the reactor must be reconstructed. Therefore, a new core configuration with a neutron trap should be put forward. At present, most existing research reactors are designed with an inverse neutron trap, for example, JRR-3M [2, 3] and CARR [4, 5]. JRR-3M is a 20 MW research reactor and the thermal neutron flux reaches the peak, 2.3 1014 n/cm−2·s, in the D2O reflector for the irradiation test. CARR is a 60 MW research reactor. In the D2O reflector, the peak thermal neutron flux reaches 8 × 1014 n/cm−2·s. Figure 1 illustrates the cross-section view of the CARR and JRR-3M, where the plate-type fuel assemblies are used in both reactors. They can be taken as references for designing new research reactors in this paper.
To realize a much higher central thermal neutron flux, a research reactor with a neutron trap is newly designed by our Reactor Engineering Analysis Lab (REAL) team at Tsinghua University. Nuclear physics design has been carried out by using a self-developed code, Reactor Monte Carlo (RMC) [6]. Figure 2 illustrates the core configuration of the new research reactor with a neutron trap. There are 18 standard plate-type fuel assemblies in yellow and 6 control fuel assemblies in grey or with a blue frame. Four irradiation tubes in blue are placed in different locations. A trap made of Be is located at the center, and a reflector made of Be and D2O surrounds the core. In the central location of Be trap, the peak thermal neutron flux reaches 4.04 × 1014 n/cm−2·s.
In this new reactor, the coolant, H2O, flows downstream the assemblies and cools the fuel plates with a forced circulation style. The core power is set as 20 MW. The fuel meat is U3Si2-Al and the fuel cladding is Al (6061). For safety concerns, the fuel meat temperature should keep lower than 400°C (673 K) and the cladding surface temperature must not exceed 220°C (493 K). No nucleate boiling is allowed throughout the assemblies, and meanwhile any kind of flow-induced vibration is avoided. This reactor is still in the design stage and heat-transfer calculation as well as safety analysis is required in advance.
With the development of computer technology, Computational Fluid Dynamics (CFD) method is a promising tool to solve nuclear thermal-hydraulic problems. Up to now, many researches have been designed under the aid of commercial CFD code, like STAR-CD, CFX, and FLUENT [7–9]. CFD work is very popular and effective especially when solving single-phase issues. Ikeda et al. [10] examined the applicability of the CFD model for predicting the critical heat flux position for a water DNB test. Single-phase simulation was proven to play a significant role in designing PWR spacer grids for improved CHF performance even though two-phase CFD model is not mature. Besides, many other studies have also proved the feasibility of CFD methods when studying the effects of the spacer and mixture grids [11–13]. The local information including flow, turbulence, and heat transfer can be of great help when designing and optimizing the rod bundles. In the past decades, Westinghouse has been using CFD code, STAR-CD, to model the precise flow with mixing vanes in PWR fuel assemblies. Meanwhile, CFX has already been used in CANDU reactor designs [14]. In some fields, CFD tends to be an irreplaceable and unique option. Salama and El-Morshedy [15] and Song et al. [16] have made detailed studies on flow blockage in MTR, which proves CFD a reliable tool for the problems of plate-type fuel assembly under both steady and transient conditions. With User Defined Functions (UDFs), power input and inlet conditions can be set freely to simulate different working conditions. Besides, CFD method has been used to analyze the flow distribution in System-Integrated Modular Advanced Reactor (SMART) [17].
All the investigations mentioned above demonstrate the feasibility and potential of CFD methods to solve thermal-hydraulic issues. For plate-type fuel assembly, it is very convenient and effective to construct the three-dimensional model and the boundary conditions. Until now, few researches concern heat transfer problems covering the whole plate-type assembly. On the contrary, most investigations only study the problems of the partial assembly. Meanwhile, fluid-structure interaction problems have not been studied sufficiently.
The aim of the present work is to provide guidance and instruction for designing the research reactor with a neutron trap. FLUENT is utilized for thermal-hydraulic calculations of standard and control fuel assemblies. A uniform distribution of inlet velocity and power input is first considered for hot assemblies to analyze the safety issue. Besides, a standard fuel assembly and a control fuel assembly with an average power are compared in detail. In addition, the assembly with a hot spot is analyzed based on the Monte Carlo calculation result.
2. Computational Setup
2.1. Computational Domain and Boundary Conditions
Figure 3 illustrates the schematic diagram of a standard and a control fuel assembly. In real situation, the flow rate in different channels of a fuel assembly has little difference. To simplify the calculation model, only the heat production parts colored by orange are considered in both assemblies. The height is 750 mm. Fuel meat, cladding, and coolant parts are all considered. There are 20 fuel plates in the standard fuel assembly and 16 fuel plates in the control fuel assembly. The fuel plate comprises the U3Si2-Al fuel meat and Al (6061) cladding. The width of coolant channel of the standard and control fuel assembly is 2.28 mm and 2.38 mm, respectively.
(a) Standard fuel assembly
(b) Control fuel assembly
Table 1 summarizes the basic thermal-hydraulic information of the research reactor. The working pressure is 0.152 MPa (abs) and the inlet coolant temperature is 35°C. The coolant flows downstream the assemblies and takes away the heat. The power distribution reads aswhere is the power peak factor and is the average power density. In this study, the power peak factor is 2.8, which is enough to provide adequate margin to ensure its safety.
Table 2 summarizes the computational boundary conditions and solving algorithm used in this study. The inlet condition of the computational domain is velocity inlet with different velocity magnitudes. The inlet turbulence intensity is set as 5% in all calculations. At the outlet, the pressure outlet is imposed and the outlet gauge pressure is set as 0. When pressure drop through the assembly is obtained, local pressure then can be calculated. The symmetry condition is applied on the side coolant planes while the surfaces of cladding and Al frame are set as no-slip wall condition. For solving the turbulence flow, standard model [8, 9] is chosen to solve the Reynolds-averaged Navier-Stokes equations.
In this study, the coolant parameters are not constant and they change against the temperature. The main parameters considered in this study are viscosity , thermal conductivity , and density , which read aswhere is the absolute temperature.
2.2. Governing Equations and Numerical Methods
In the present study, a three-dimensional, steady, incompressible, and turbulent flow issue is considered. Under the assumption of a single-phase Newtonian fluid, the Reynolds-averaged Navier-Stokes (RANS) equations are chosen to model the turbulent flow in the fuel assemblies of the research reactor with a neutron trap. In Reynolds averaging, the solution variables in the exact N-S equations are decomposed into the mean and fluctuation components. For the velocity components:where and are the mean and fluctuating velocity components ().
Substituting expressions of this form for the flow variables into the exact continuity and momentum equations and taking a time average yield the ensemble-averaged momentum equations. They can be written in Cartesian tensor form asThe Reynolds stresses must be modelled to close (5). Boussinesq hypothesis is employed to model the Reynolds stresses:where is the turbulence kinetic energy and is the turbulent viscosity. In standard model, two additional transport equations are solved, and is computed as a function of and :
In these equations, represents the generation of turbulence kinetic energy due to the mean velocity gradients. is the generation of turbulence kinetic energy due to buoyancy. represents the contribution of the fluctuating dilatation. , , and are constants. and are the turbulent Prandtl numbers for and , respectively.
For RNG model and realizable model, they have the similar equation forms, but they include some refinements compared with standard model in different ways. More information about the conservation equations, RANS equations, and standard values of the model constants employed here can be found in ANSYS handbook [8, 9].
In solid regions, the energy transport equation has the following form:where is the conductivity and is the enthalpy. is the volumetric heat source, which is defined by UDF in the paper.
The commercial CFD software, FLUENT, is used for solving the governing equations. The segregated solver is utilized and the Semi-Implicit Methods for Pressure Linked Equations (SIMPLE) algorithm is applied to treat pressure-velocity coupling.
2.3. Verification and Mesh Generation
Before carrying out the calculation, verification test needs to be done first, especially for CFD calculation. Jo et al. [18] have conducted a series of experimental investigations of convective heat transfer in a narrow rectangular channel. One case with an inlet velocity m/s and a thermal power kW is chosen to do the comparison. Figure 4 illustrates the calculation model where the width of the channel is 2.35 mm. Figure 4(a) shows the cross-section view of the experiment test section. Figure 4(b) is the CFD calculation model including coolant, insulation, and heater. The coolant in blue flows downward the narrow channel. The grey part is the insulation, whose thermal conductivity is set as 1 W/m·K. The purple part in the model is the heater, in which a uniform heat generation is applied.
(a)
(b)
Figure 5 shows the calculation differences resulting from different grid numbers. Structured hexahedral mesh is used to generate all the grids. As a result of increasing calculation grids, the outlet heater surface and coolant temperature as well as the pressure drop tend to be identical. Figure 6 shows the detailed comparison of temperature in the axial direction for the third case with 1080000 grid numbers. The calculation result agrees well with the experimental data [18]. For the cases with more girds, the difference between each other is rather small. Too many grids will take too much calculation resources, so meshing the model like the third case will be rational.
Based on the test above, similar meshing method is utilized for the fuel assemblies, which is illustrated in Figure 7. For both assemblies, the meshing methods are similar. The fuel meat in yellow and the cladding in grey are defined as solid, while the coolant in blue is defined as liquid.
is the dimensionless normal distance from the wall [8, 9]. In this paper, the range of is between 30 and 220, which is in the appropriate range of using standard wall function for treating the near-wall calculation. For additional verification and choosing a proper turbulence model, the calculation results with different k-ɛ models are compared in Figure 8. Some experiments have proven that the Dittus-Boelter (D-B) formula is available for both upward and downward flow in narrow-gap channels with a Reynolds number higher than 104 [2, 3, 19]. Jo et al. [18] have also proven it recently. Besides, Jo et al. [18] also proposed a new formula for heat transfer in narrow channels. Both empirical formulas are used to do the comparison and the results are shown in Figure 8. There is only a little difference between the curves when is larger than 2 × 104. As is larger than 2 × 104 in all the assembly calculations, all the models will be feasible. When carrying out thermal calculation of CARR, Yang [5] has proven the standard model is the best in solving narrow channels problems. Therefore, standard model is chosen as the turbulence model in the present study.
3. Results and Discussion
3.1. In-Core Temperature Distribution under Uniform Input Power and Velocity
In this study, a cosine-shape input power is utilized in height direction. However, for each fuel plate in the assembly, there is a uniform input power distribution. The average power density is 1.305 × 109 W/m3. The maximum peak factor is used in the hot assembly while only axial peak factor is considered in the average assembly.
Thermal-hydraulic calculations are carried out by forced-convection mode in the standard and the control fuel assembly. In this case, inlet velocity is the only variable to influence the distributions of temperature and pressure in the assemblies.
Figure 9 shows the temperature distribution in axial direction of a standard fuel assembly with an inlet velocity of 7.5 m/s. There are totally 20 fuel plates and 21 coolant channels, and they are numbered as shown in Figure 9(b). In the axial direction, the coolant temperature increases along the channels and takes away the heat produced by the fuel meat. The temperatures of fuel meat and cladding reach the peaks, 359 K and 352 K, at the location beneath the input power peak shown in Figure 9(d). Fortunately, the cladding temperature is much lower than the saturation temperature. Meanwhile, the fuel temperature is much smaller than the limit value even at the central part of the assembly. The temperature becomes lower from inside out and a temperature gradient exists on the cladding surface. The bulk coolant temperature is relatively low.
Figure 10 illustrates the temperature on the plane at the height of 326.26 mm where the peaks of fuel and cladding temperature occur. As all the fuel plates have the same power distribution, their temperature peaks only differ a little as depicted in Figure 10(c). Meanwhile, the temperature shown in Figure 10(b) is flat in each fuel plate since the heat source is the same in direction. At different heights, the profiles of temperature distribution are similar, where an obvious temperature gradient only exists near the wall.
Figure 11 depicts the local velocity and temperature near the cladding surface at a height of 326.3 mm. It should be noted that the coolant velocity in the center part of the channel is 8.1 m/s, which is a little higher than the inlet velocity. The black lines in Figure 11 represent the results calculated with much smaller grids. The comparison shows that refining the grids will show more detailed information in local parts. However, main parameters will change little, which proves the original meshing method is fine enough to get the key parameters. High velocity and temperature gradients can be found in the boundary layer.
As mentioned before, inlet velocity is the only variable determining the fuel, cladding, and coolant temperatures. Figure 12 shows the temperatures of cladding surface, coolant, and saturation point at the hot spot. Here the cladding surface temperature reaches the peak with different inlet velocities for the standard fuel assembly. Both temperatures of the coolant and the cladding surface decrease as the inlet coolant velocity becomes larger. Meanwhile, the saturation temperature becomes a little lower due to the lower local pressure. During the velocity range of 4.5 m/s to 7.5 m/s, the standard fuel assembly keeps safe since no local nucleate boiling occurs at the hottest spot.
To evaluate the cooling ability, number along the channels with different Reynolds numbers is calculated, as shown in Figure 13. With a higher Reynolds number, becomes a little larger and results in a stronger cooling ability. Furthermore, in different cases has similar profiles. is relatively large at the inlet due to the inlet effect and small temperature difference. Flowing downstream the channels, increases until the coolant arrives at the outlet.
Figure 14 illustrates the pressure and temperature distributions at the bottom of fuel plates against the inlet coolant velocity. Taking account of the pressure drop, the pressure at the bottom is the lowest. The corresponding saturation temperature reaches the lowest. From Figure 14(b), the exit coolant temperature is well below the corresponding saturation temperature. Therefore, the fuel assembly always keeps safe when the inlet velocity covers a wide range, 4.5 m/s through 7.5 m/s. For larger inlet coolant velocity will also be OK. But considering excessive inlet coolant velocity may lead to larger hydraulic load and then result in vibration. So it is not necessary to calculate more with larger velocity.
(a) Pressure
(b) Temperature
Figure 15 illustrates the calculation results of the standard and the control fuel assembly with an average power input. The two assemblies have the same fuel energy density. The cladding surface and coolant outlet temperatures are compared in this figure. The surface temperature is found to be much lower than that of the hot assembly. The temperature decreases as the inlet velocity becomes bigger. The temperature of the outlet coolant has a similar trend. The coolant temperature rise varies from 10 K to 17 K when the inlet velocity decreases from 7.5 m/s to 4.5 m/s.
3.2. Cases with Different Axial Power Densities in Different Fuel Plates
More precise calculations are carried out in this section. After neutronics calculations by the RMC code, the local power peak has been determined [20] as shown in Figure 16. The average power density in each fuel plate is plotted in Figure 16(b). The power data will be used as an input for the present calculations. The plate with the power peak is located on the side of the standard assembly. The power density decreases along the direction from the neutron trap to the other side. 20 UDFs are used to input the power distribution for each fuel plate. The average inlet velocity is set as 4.5 m/s.
(a)
(b)
Figure 17 illustrates the calculation results including the temperature distribution on the plane in Figure 17(a) where the fuel peak temperature exists, the peak temperature in different fuel plates in Figure 17(b), and the outlet coolant temperature in each channel in Figure 17(c). All the temperature profiles decrease as the plate or channel number becomes bigger. In this case, the maximum fuel peak temperature is 371 K, which is lower than that in the hot assembly with an inlet velocity of 4.5 m/s. It is reasonable to deduce that the assembly will work safely with larger inlet coolant velocities.
(a)
(b)
(c)
4. Conclusions
In the present study, heat transfer characteristics of a high flux research reactor with a neutron trap are numerically investigated. Standard model is adopted to simulate the turbulence. Verification has been done by comparing with experiment data. Both standard and control fuel assemblies are considered.
In the hot assembly, uniform distributions of power and inlet velocity are considered, and all the fuel and cladding peak temperatures are below the design limits. The fuel cladding surface temperature keeps lower than the saturation temperature within the inlet velocity range of 4.5 m/s to 7.5 m/s. With the same average power input, the temperature of control fuel assembly is lower than that of the corresponding standard fuel assembly. This is a result of a wider coolant channel for the latter. A precise case with a distributed power input in each fuel plate is simulated and also proves the assembly in a safe state.
Nomenclature
| Be: | Beryllium |
| : | Modeling coefficient |
| : | Generation of turbulence kinetic energy |
| : | Enthalpy |
| : | Diffusion flux of species |
| : | Turbulence kinetic energy |
| : | Nusselt number |
| : | Pressure |
| : | Prandtl number |
| : | Volumetric heat generation |
| : | Reynolds number |
| : | Volumetric heat source |
| : | Temperature |
| : | Velocity of fluid |
| : | Nondimensional wall distance |
| : | Contribution to the kinetic energy dissipation rate |
| : | Vertical distance from the bottom of the fuel assembly. |
| : | Rate of turbulence kinetic energy dissipation |
| : | Thermal conductivity of fluid |
| : | Dynamic viscosity of fluid |
| : | Density of fluid |
| : | Turbulent Prandtl number |
| : | Stress tensor. |
| avg: | Average value |
| : | Number. |
Highlights
(i)Heat transfer characteristics of plate-type fuel assemblies of high flux research reactor are numerically investigated with CFD method.(ii)Both standard fuel assembly and control fuel assembly are analyzed in the forced circulation mode.(iii)In hot assembly, the cladding surface temperature is always lower than the local saturation temperature and the fuel temperature is lower than the design limit.(iv)The temperature rise in the control fuel assembly is smaller than that of the standard fuel assembly.(v)Numerical results show that the research reactor with a neutron trap structure will keep safe under the inlet coolant velocity range between 4.5 m/s and 7.5 m/s.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work is partially supported by Project 2012AC01 of Key Laboratory of Neutron Physics, Chinese Academy of Engineering Physics.