Thermal Hydraulics Analysis of the Distribution Zone in Small Modular Dual Fluid Reactor

The Small Modular Dual Fluid Reactor (SMDFR) is a novel molten salt reactor based on the dual fluid reactor concept, which employs molten salt as fuel and liquid lead/lead-bismuth eutectic (LBE) as coolant. A unique design of this reactor is the distribution zone, which locates under the core and joins the core region with the inlet pipes of molten salt and coolant. Since the distribution zone has a major influence on the heat removal capacity in the core region, the thermal hydraulics characteristics of the distribution zone have to be investigated. This paper focuses on the thermal hydraulics analysis of the distribution zone, which is conducted by the numerical simulation using COMSOL Multiphysics with the CFD (Computational Fluid Dynamics) module and the Heat Transfer module. The energy loss and heat exchange in the distribution zone are also quantitatively analyzed. The velocity and temperature distributions of both molten salt and coolant at the outlet of the distribution zone, as inlet of the core region, are produced. It can be observed that the outlet velocity profiles are proportional in magnitude to the inlet velocity ones with a similar shape. In addition, the results show that the heat transfer in the center region is enhanced due to the velocity distribution, which could compensate the power peak and flatten the temperature distribution for a higher power density.


The Small Modular Dual Fluid Reactor
The molten salt reactor is one of six Generation IV reactor types, which adopts molten salt as fuel and has many unique characteristics compared to solid fuel reactors. The study on two-fluid molten salt reactor can be traced back to 1966, when a two-fluid molten salt breeder reactor (MSBR) was designed at Oak Ridge National Laboratory (ORNL) [1]. The load-following capability of this MSBR system was studied for various ramp rates of the power demand [2]. Later a fast breeder reactor of about 2 GW thermal output using molten chlorides both as fuel and coolant was proposed [3].
A new design of the two-fluid MSBR called the Dual Fluid Reactor (DFR), which employs liquid lead instead of molten salt for the secondary circuit as coolant, was proposed by the researchers at the Institute for Solid-State Nuclear Physics (IFK) [4]. Differing from Taube's design, the high volumetric heat capacity of the liquid lead enables a significant heat removal capability and, thus, allows operation with a much higher power density than that of conventional reactors. A series of studies covering the neutronic characteristics, sensitivity performance, coupled calculations, emergency drain tank estimation, depletion as well as hydraulic simulations can be found in references [5][6][7][8][9][10]. It has to be noticed that previous studies on the DFR focused mainly on a large type of reactor design with a thermal power of 3 GW. Yet, a smaller type of less than 0.3 GW would be more flexible and could be applied to a wider range of industrial processes, such as hydrogen production, water desalination, etc.
In this work a Small Modular Dual Fluid Reactor (SMDFR) is proposed with a nominal thermal power of 0.1 GW (100 MWt). The scheme is shown in Figure 1. The reactor core of the SMDFR, formed by hexagonally-arranged fuel tubes and surrounded by the reflector zone, is immersed at the center of a lead pool inside a cylindrical tank. As shown in Figure 2, the bottom and top of the core zone (from z = 0.2 m to z = 2.2 m) are joined by a distribution zone (from bottom, z = 0 m, to z = 0.2 m) and a collection zone (from z = 2.2 m to z = 2.4 m, symmetric with respect to the distribution zone), in order to distribute and collect the molten salt before and after flowing through the core region. Starting from the bottom part of the core region (Figure 2), the molten salt (depicted in golden color) goes into the distribution zone being distributed in the gaps between the hexagonally arranged coolant tubes (depicted in red) and then flows upwards inside the fuel tubes (depicted in golden color) through a plate with holes. After reaching the upper core part, it flows into the collection zone and then leaves the core region towards a chemical processing plant, which is employed for the online processing of the molten salt. During the operation, the salt coming from the collection zone is continuously passed to the chemical processing plant, in which the fission products are removed, and nuclear fuel components can be removed or added to maintain the optimal fuel composition. After being processed, it is pumped back into the bottom part for the next flow cycle through the core region. The liquid lead (Figure 2), acting as coolant in the secondary circuit, flows from the bottom center to the top center to cool down the core, which is heated by the nuclear fission reactions occurring in the molten salt. Upon leaving the core region, it is pumped to the heat exchanger located in the peripheral region of the tank and then transfers the heat to the tertiary circuit for utilization. Compared to the original design, the SMDFR is scaled down to 0.1 GW (100 MW), which makes it particularly suitable for providing energy for production facilities (e.g., water desalination [11], hydrogen production, or mines) in remote locations and more flexible, since it does not necessarily need to be hooked into a large power grid, and is able to be attached to other modules to provide increased power supplies if necessary. In addition, it also has changed from a pipe-type to a pool-type reactor in order to ensure the capability of establishing natural circulation during reactor transients or accidents. The hexagonal arrangement of fuel and coolant tubes realizes a dense lattice in the core and a high power density.  The technical data of the SMDFR are listed in Table 1. In this work a mixture of uranium tetrachloride and plutonium tetrachloride is chosen as fuel salt and the liquid lead is chosen as coolant. For the pipe walls silicon carbide is employed. The thermo-physical properties of these materials are listed in Table 2 according to [12][13][14][15]. The formulas in the reference are applied for each 100 K interval from 800 K to 1600 K (below the boiling point of molten salt: 1950 K [15]), which includes the operational temperature range. The corresponding interpolated values based on these temperature points are applied for the calculations.   [17] and Heat Transfer module [18] is used to perform the 3D pin-by-pin (for all tubes) thermal-hydraulics calculation. COMSOL Multiphysics is a multi-physics simulation code package which is used in many industry application, including fluid flow, and can provide a high accuracy numerical solution by using the Finite Element Method (FEM) to solve the corresponding set of partial differential equations. The software runs on a standalone PC with Intel R Core TM i7-7700 processor, 64 GB RAM and solid state disk.

Research Objective
A novel design of the core (including the distribution zone, core zone and collection zone) for a small modular dual fluid reactor is proposed and the thermal hydraulics characteristics of the distribution zone are investigated. In order to obtain both systematic and local behaviors of the two fluids, a 3D computer model is built using COMSOL Multiphysics with the CFD and Heat Transfer modules and the corresponding numerical simulations are performed. In addition, the sensitivity analysis is conducted to investigate the system responses of various inlet velocities of both fluids.

Modeling and Simulation
The modeling and simulation process presented in this paper is divided into four parts: geometry extraction and building, setup of governing equations, meshing, selection of the numerical solver and performing the simulation.

Geometry
There are many more fuel tubes than those depicted in Figure 3, and a full core calculation requires very large computational resources. Because of the symmetrical hexagonal configuration of the core region, a sector of 30 • (Figure 4, right) of the lower part of the core with a height of 0.4 m (from z = 0 m to z = 0.4 m, Figure 4, left) around the central cylindrical axis, only one twelfth of the complete geometry, was selected to represent the full scale geometry by applying symmetric boundary conditions (mirror plane) [10]. The number of fuel tubes to be simulated was then decreased from 1027 to 100, among which 72 are complete tubes. A complete tube means that the whole tube was included in the selected 30 • sector and an incomplete tube means only one part of the tube was included. The tubes at the boundary were incomplete tubes, as shown in Figure 4 (right).
In Figure 4 (right) three colored domains can be observed: golden, red and grey, which represent respectively the fuel domain, the coolant domain, and the pipe wall domain. This geometric model contained the distribution zone with a height of 0.2 m and its upper part (from core zone) with a height of 0.2 m. In the fuel domain, the molten salt entered through the inlet pipe from the right side, flowed into the distribution zone, spread inside the gaps between the hexagonally-arranged coolant tubes, flowed upwards into the fuel tubes, and then left the core zone through the upper part. In the coolant domain, the liquid lead flowed into the core zone through the bottom coolant pipes, it then flowed upwards between the gaps along the fuel tubes, and, finally, it left the core region through the upper part.

Governing Equations
Heat transfer between two fluids separated by a solid wall was the main thermal process to be simulated, with both fluids the molten salt and the liquid lead having high Reynolds numbers (molten salt: 1.17 × 10 4 , liquid lead: 1.28 × 10 4 ). Therefore, the physical COMSOL Multiphysics-interface "Conjugate Heat Transfer: Turbulent Flow, k-" is chosen for the calculation, which combines the "Heat Transfer in Solids and Fluids" interface and the "Turbulent Flow, k-" interface and, thus, can be used to simulate the coupling between heat transfer and fluid flow.
Three compressibility options are available in COMSOL Multiphysics: compressible flow (Ma < 0.3), weakly compressible flow, and incompressible flow. The option of compressible flow makes no assumptions for the system and takes into account any dependency that the fluid properties may have on the variables. The equations of weakly compressible flow look the same as these of compressible flow except that the density is evaluated at the reference absolute pressure. For the incompressible flow option, the density is regarded as constant and evaluated using the reference temperature and pressure. [19] Since the density of molten salt has a strong dependency on its temperature, the weakly compressible Reynolds-Averaged Navier-Stokes (RANS) equations were adopted as a compromise between computational cost and accuracy: where K is the viscosity term taking into account the interactions between the fluctuating parts of the velocity field. In order to close the above equations, two additional transport equations with two dependent variables, the turbulent kinetic energy k and the turbulent dissipation rate , are introduced: where µ T and P k are the turbulent viscosity and the production term and are given by: The model's constants used are: For the heat transfer in solids and fluids, the following equations are adopted: and the heat flux between the fluid with temperature T f and the wall with temperature T w is calculated by: Since there is no recommended turbulent Prandtl number model for liquid lead flow, the Kays-Crawford model is adopted: where the Prandtl number at infinity is Pr T ∞ = 0.85.

Meshing
The mesh was manually built by defining swept meshes and free tetrahedral meshes of various sizes for different regions. Since the flows of molten salt and liquid lead vary gradually along their channels, inside and outside the pipes, the finite elements can be quite stretched in the z-axial direction. For this reason, swept meshes were adopted in order to reduce the complexity of the mesh. The remaining parts of the model were relatively irregular and, thus, they were meshed by free tetrahedral mesh elements. The completely generated mesh ( Figure 5) has 599,333 elements with an average element quality of 0.6259, and its 3D-plot of mesh element quality is shown in Figure 6. Around 96.9% of the model volume was filled by the mesh elements with qualities greater than 0.1, which is considered accurate enough [20]. From the numerical point of view such a mesh was acceptable in terms of accuracy and stability of the solution.
Six boundary layers were added for both the molten salt and the liquid lead with a stretching factor of 1.1. To achieve a high accuracy the value of the wall resolution in viscous units was set to lie between 11.6 and 100, which corresponded to 4.12 × 10 −4 m for the molten salt and to 6.645 × 10 −5 m for the liquid lead.

Numerical Solver
Two fundamental classes of solvers, direct and iterative solvers, were provided by COMSOL Multiphysics. The direct solvers available within COMSOL Multiphysics are "PARDISO", "MUMPS", and "SPOOLES", as well as "Dense Matrix" solvers. The Dense Matrix Solver can be only used for Boundary Element Method models, and therefore not suitable for the model in this paper. Contrary to direct solvers, iterative solvers approach the solution gradually and consume less memory, rather than in one large computational step. Three direct solvers, "PARDISO", "MUMPS", "SPOOLES", and one iterative solver, "GMRES", were selected for solving the same model, and the convergence criterion used for all variables was a relative tolerance of 0.005. Their performance comparison is listed in Table 3. "PARDISO" was the fastest solver while consuming the highest physical and virtual memories except for "SPOOLES". The iterative solver "GMRES" consumed less memory and was slower then the direct solvers "PARDISO" and "MUMPS". All these solvers converged to the same results. Since the sizes of physical and virtual memories were not a bottleneck for the model to be solved, "PARDISO" was considered as the most efficient solver and was selected for solving all the simulations shown in this paper.

Results and Discussion
The results and discussion section is divided into four parts: verification, hydraulic characteristics, heat transfer analysis, and sensitivity analysis of inlet velocities.

Verification
Since the SMDFR design is a new concept and no experimental data can be used as a reference, the solution should be verified through the wall resolution in viscous units for the first layer, and also the consistence of results from different mesh scenarios for the mesh convergence.
The value of wall resolution in viscous units tells how far into the boundary layer the computational domain starts and should be lower than 500 to ensure that the logarithmic layer meets the viscous sublayer. As shown in Figure 7, for more than 90% of the fluid domains the value of the wall resolution in viscous units lay between 11.6 and 100, which means the first boundary layers approximately corresponded to the beginning of the logarithmic layer and the required accuracy of the results through the first boundary layers (next to the wall) was achieved. There was no need to further refine the boundary layer mesh. In order to investigate the mesh convergence, five mesh scenarios (M0 to M4) were selected, M0 was the original mesh with six boundary layers, M1 and M2 had three and eight boundary layers respectively, and M3 and M4 have coarse and refined meshes separately. The meshing parameters are shown in Table 4. In order to quantify the simulation results for comparison and analysis, the outlet surface was divided into 19 channels for molten salt, and 20 channels for liquid lead, along the radial direction from center to circumference, as shown in Figure 8. The mass flow rate averaged velocity and temperature of each channel were selected as the quantities to be compared and analyzed.
As shown in Figures 9 and 10, all the mesh scenarios delivered consistent results for velocity and temperature profiles of both molten salt and liquid lead, especially in the peripheral regions. However, in the center regions (next to the center-line/axis of the core), some divergences were observed for M1 (black) and M3 (green). In order to ensure a high enough accuracy in the whole domain, M0 (red), which had consistent results with M2 (magenta) and M4 (blue) but consumed less computational resources, was verified against the finer meshes and adopted for the final simulations.

Hydraulic Characteristics
From Figure 9 it can be seen that the molten salt velocity was relatively low in the central regions, which was caused by the high flowing resistance of the cross-flow through the coolant tubes. In order to get a better understanding of this phenomenon, a plane of the fuel inlet pipe was obtained by cutting the model by the blue plane (z = −0.1 m), as shown in Figure 11 (top). A majority of the molten salt tended to go upwards along the fuel tubes and only a small portion of the molten salt could pass the peripheral regions and reach the central core regions (Figure 11, bottom). Taking a look into the the velocity vector field (Figure 12), one cold and two hot vortices could be observed located at the lower left corner (cold), top right corner (hot) and lower right corner (hot). Special consideration must be given to these locations, since the risk of structure material failure was relatively high and vibration might occur in the regions with vortices. Figure 11. Top, plane of the fuel inlet pipe cut by the blue plane; bottom, Velocity vector field on the plane of the fuel inlet pipe of the molten salt. For liquid lead the situation was the contrary, starting from the first two channels it had a high velocity and after that the velocity became relatively uniform until the last two channels (Figure 9). Since in the region next to the circumference there was an additional channel surrounding the fuel channels, a certain portion of liquid lead was distributed to this channel and thus the velocity in peripheral region was reduced. However, the liquid lead in the central region (first two channels) accumulated due to the symmetry structure and thus had a relatively high velocity compared to others.
In order to quantify the energy loss during the flowing, Bernoulli's equation applied for a micro element of fluid is introduced: multiplying both sides by a micro flow element dm (dm = ρ 1 gu 1 dA 1 = ρ 2 gu 2 dA 2 ), and then integrating both sides by the corresponding cross section areas A 1 and A 2 : where h w is the energy loss in W, which determines the demand of pumping power. Selecting the inlet and outlet cross section areas of the fluids as A 1 and A 2 , the energy loss in the distribution zone was then obtained, which equals h w . The h w of molten salt in the distribution zone (for this one twelfth core sector) is 64 W, which was much lower if compared to that of the liquid lead. Since the heat removal function was mainly accomplished by the liquid lead, the flowing of molten salt is only needed for chemical processing and its mass flow rate does not have to be large. Due to its high mass flow rate the h w of liquid lead (for this one twelfth sector) was quite significant: 14,515 W. Considering the whole structure and including the collection zone with the assumption of complete symmetry, the h w of liquid lead was around 348 kW only for the distribution and collection zones. It means that the pump had to at least overcome this energy loss and it should have provided larger power when taking into account the energy losses in other parts of the secondary circuit, especially in the heat exchanger.

Heat Transfer Analysis
From Figure 10 it can be seen that the molten salt temperature was continuously increasing from the center to the external circumference and the liquid lead had a relative uniform temperature distribution ( Figure 13). This can be easily explained by the characteristics of the cross (lower part) and co-current flows (upper part). In the center region, there were two factors enhancing the heat transfer: the long flow path of molten salt crossing the fuel tubes ( Figure 14) and the large velocity difference between molten salt and liquid lead. The heat transfer was weak due to the low velocity difference in the peripheral region, which resulted in a high molten salt temperature. Since the radial distribution of power in the core follows a Bessel function of order zero (first period), the heat transfer enhancement in the center region could compensate the power peak and flatten the temperature distribution to achieve a higher operational temperature for both molten salt and liquid lead, which means the potential for a higher power density.
The exchanged heat between molten salt and liquid lead in the distribution zone resulted in their enthalpies changing and yielding a total thermal power transfer of 1.65 MW for this one twelfth sector. Considering the whole structure and including the collection zone, the heat power transferred increased to 39.6 MW, which was around 40% of the total thermal power. It means that the distribution zone as well as the collection zone not only distributed and collected the fluids, but also had a significant impact on the total heat removal process of the core, while the remainder thermal energy was transferred inside the core zone.

Sensitivity Analysis of Inlet Velocities
In order to investigate the influence of the inlet velocities on the thermal-hydraulic characteristics of the distribution zone, five inlet scenarios (Table 5: I0 for the original one and I1 to I4 for the changed  inlet velocities) were simulated and compared.   As shown in Figure 15, the magnitude of outlet velocity of molten salt was proportional to its inlet velocity and its profile's shape was similar. The same conclusion can be made for the liquid lead. The temperature profile of liquid lead (Figure 16) was not influenced too much by the inlet velocities thanks to its heat capacity and mass flow rate. However, the thermal behaviour of the molten salt was influenced by its inlet velocity: a higher molten salt inlet velocity resulted in a higher magnitude but a similar trend in the temperature profiles and vice-versa. The impact of the liquid lead's inlet velocity on the molten salt's temperature profile was negligible, which means that the heat transfer between the fluids in the distribution zone occurred mainly in the lower part (cross flow). It has to be pointed out that the heat transfer in the core zone (co-current flow) was definitely influenced by the velocity variation of both fluids due to the change of relative velocity between them, and a strong coupling effect should be considered during transient conditions with rapid velocitiy variation expected. The thermal hydraulic coupling effects between the molten salt and liquid lead, as well as the thermal hydraulics-neutronics coupling, will have to be considered for further simulations of the core, especially in the case of transients and accident situations.

Conclusions
The computer model built by COMSOL Multiphysics is verified by examining the value of the wall resolution in viscous units and by performing mesh sensitivity studies. The energy loss and heat exchange in the distribution zone are quantified, and it is found that the pumping power delivered to the liquid lead has to be much larger than that required by the molten salt in order to overcome the large frictional energy loss of the flowing liquid lead. An amount of thermal power corresponding to about 40% of the total reactor thermal power is transferred in the distribution and collection zones, which not only affects the hydraulic performance of both fluids but also has a significant influence on the heat removal process. Since the radial distribution of power in the core follows a Bessel function of order zero (first period), the heat transfer enhancement in the center region could compensate the power peak and flatten the temperature distribution to achieve a higher operational temperature for both molten salt and liquid lead, which means the potential for a higher power density. The outlet velocity profiles observed are proportional in magnitude to the inlet velocity ones with a similar shape.