Development of a Coupled Code for Steady-State Analysis of the Graphite-Moderated Channel Type Molten Salt Reactor

The molten salt reactor (MSR) is one of the six advanced reactor concepts selected by Generation IV International Forum (GIF) because of its inherent safety and the promising capabilities of TRU transmutation and Th-U breeding. In this study, a threedimensional thermal-hydraulic model (3DTH) is developed for evaluating the steady-state performance of the graphite-moderated channel typeMSR.The coupled code is developed by exchanging the power distribution, temperature, and fuel density distribution between SCALE and 3DTH. Firstly, the thermal-hydraulic model of the coupled code is validated by RELAP5 code. Then, the mass flow distribution, temperature field, keff , and power density distribution for a conceptual design of the 2MWt experimental molten salt reactor are calculated and analyzed by the coupled code under both normal operating situation and the central fuel assembly partly blocked situation.The simulated results are conductive to facilitate the understanding of the steady behavior of the graphite-moderated channel type MSR.


Introduction
The research on MSR concept can date back to the 1940s-1960s with the project of the aircraft reactor experiment (ARE) [1] and the molten salt reactor experiment (MSRE) [2,3] in Oak Ridge National Laboratory (ORNL).The successful design, construction, and operation of MSRE established the basic technologies, provided experimental data for MSR, and led to a conceptual design of the 1000MWe molten salt breeder reactor (MSBR) [4] at ORNL in the 1970s.After that, a series of studies on MSR for Th-U breeding or TRU transmutation were done by France [5], Japan [6], Russia [7], and China [8][9][10] since the 1980s.With the advantages of inherent safety, excellent neutron economy, and the capabilities of online reprocessing, MSR was selected as one of the six advanced reactors in GIF [11].In 2011, Chinese Academy of Sciences (CAS) launched the Thorium-based Molten Salt Reactor (TMSR) nuclear system project.It aims to build a 2MWt graphite-moderated liquid fuel molten salt reactor and to solve the key technical problems of TMSR [12].
Unlike the solid fuel reactors, the primary coolant of MSR is a liquid molten salt mixture of carrier salt and fissile material (such as U or Th), which serves as fuel at the same time.The graphite-moderated channel type MSRs generally adopt graphite as neutron moderator.In this type of MSR, most of the fission energy is released directly into the salt [13], and only a small part of energy is deposited in the graphite moderator due to gamma rays and fast neutron slowdown.The fuel salt flowing through the fuel channel transfers the fission energy out of the core and cools the graphite moderator and reflector.The thermal coupling of fuel salt in this type of MSR is achieved through heat conduction in the graphite moderator.Considering the above-mentioned characteristics, detailed understanding of the thermal-hydraulic behavior of MSR is important for its design.
Efforts on the development of suitable simulation tools for the graphite-moderated channel type MSR were made by different authors for different purposes.In ORNL, Engel et al. calculated the fuel salt and graphite temperature distributions of MSRE by an analytical method [14], in which the fuel assembly was equivalent to a hollow cylinder and the radial one-dimensional heat conduction equation was solved analytically.Zhang et al. performed the steady-state analysis for the MSRE by a multiple-channel analysis model, in which the distribution of the velocity and temperature was obtained [15].Guo et al. developed a multiple-channel analysis code and coupled it with MCNP4c to analyze the steady-state behavior of the MSRE [16].Křepel et al. developed the DYN1D-MSR [17] code and DYN3D-MSR [18] code for the dynamic analysis of graphite-moderated channel type MSR.The neutron diffusion model was adopted for neutronics calculation.And the one-dimensional singlephase flow model was adopted for the thermal-hydraulic calculation.Kópházi developed a three-dimensional analysis codes DT-MSR [19].The neutronic physics part of DT-MSR is based on the neutron diffusion model, which extended a term to consider the drift of delayed neutron precursors (DNPs).For the thermal-hydraulic part, the temperature of fuel salt was calculated by the one-dimensional heat convection equation.And the temperature field in the graphite moderator was calculated by the three-dimensional heat conduction equation on a structured mesh.In DT-MSR, the heat conduction in the graphite moderator was coupled to the temperature of fuel salt and the thermal coupling between fuel assemblies was initiatively realized.
Although many tools have been developed for the graphite-moderated channel type MSR, most of them only adopt the one-dimensional heat conduction model for the moderator and neglect the thermal coupling of the adjacent assemblies.Some tools like DT-MSR can simulate the threedimensional temperature field, but they are based on the structured mesh and need equivalent geometry modification.Therefore, one objective of the present work is to develop an accurate three-dimensional thermal-hydraulic model (named 3DTH hereinafter), which can calculate the temperature field in complex geometries on an unstructured mesh and consider the thermal coupling between different assemblies.Then, a coupled code is developed through exchanging the data between 3DTH and the reactor criticality safety analyses software SCALE [20] and is applied to the steady-state analysis of a conceptual design of the 2MWt experimental molten salt reactor (named 2MW-MSR hereinafter) [21].
In this paper, Section 2 will focus on the neutronic model, thermal-hydraulic model, and coupling code.In Section 3, the correctness of 3DTH will be validated.Section 4 will briefly introduce the 2MW-MSR.In Section 5, the distributions of mass flow, temperature, and power density of the 2MW-MSR under both normal operating situation and central fuel assembly partly blocked situation will be calculated and discussed.Finally, the conclusions are given in Section 6.

Theoretical Model
2.1.The Neutronic Model.SCALE [20] is widely used for calculating the criticality, depletion, and reaction rates distribution of reactor core.The TRITON computer code is a multipurpose SCALE control module for transport, depletion, and sensitivity and uncertainty analysis.In this work, the TRITON module is adopted for calculation of criticality and power distribution with a 238-group ENDF/B-VII crosssection data library.However, the TRITON module can only calculate the power of a material number.In order to obtain the power distribution of the fuel salt and graphite, the fuel assemblies are divided into several units in the axial direction, and each unit is labeled with a unique material number.Similar to the calculation of axial power distribution, the radial power distribution is also calculated with different material numbers.
It should be noted that the power distribution calculated by SCALE has not considered the drift of DNPs and decay heat.Since the fuel salt is circulated in the primary loop, the fission products (DNP and the products which generate decay heat) move towards the outlet of the reactor due to the fuel salt flow.For a reactor operating a long time, about 7% of the reactor power is generated in the form of decay heat.The drift of decay heat slightly changes the power distribution in the active core and causes a fraction of decay heat released in the primary loop [22], which will slightly influence the temperature rise and the temperature distribution in the core and primary loop.Because the fraction of delayed neutron (DN) is relatively small, the drift of delayed neutron precursors does not obviously affect the neutron fluxes and the power distribution under the steady-state condition [16,23].

The Thermal-Hydraulic Model.
A three-dimensional thermal-hydraulic model 3DTH is developed for the graphite-moderated channel type MSR.The temperature distribution in the solid region which includes the graphite of fuel assembly, reflector, and reactor vessel is calculated by the three-dimensional heat conduction model.The steady heat conduction equation can be written as and T  are the thermal conductivity and temperature of the solid region, respectively; Q  is the power density deposited in the solid region.The outer wall of the solid region and the inner wall of fuel channel employ the adiabatic boundary condition and convective heat transfer boundary condition, respectively.The fuel salt in the channel is described by the onedimensional single-phase flow model.Under steady-state condition, the governing equations are momentum, mass, and energy conservation equations: v, , T  , C  are the velocity, density, temperature, and heat capacity of the fuel salt, respectively; A and P  are the flow area and perimeter of the fuel channel, respectively; P  is the frictional pressure drop; Q  represents the volumetric heat released in the fuel salt; Q ℎ denotes the convective heat transfer at the wall of the fuel channel.The Darcy formula is adopted to calculate frictional pressure drop, which can be written as f and D  are the friction coefficient and diameter of the fuel channel, respectively.For the laminar flow f can be written as The thermal coupling of fuel salt and solid region is built by the convective heat transfer boundary condition at the inner wall of the fuel channel.The heat flux Q ℎ (in (4)) can be described as and Nu are the thermal conductivity of fuel salt and Nusselt numbers, respectively.The empirical correlation provided by ORNL's experiment is adopted in the 3DTH to calculate the Nusselt numbers [24].For the laminar flow, the empirical correlation can be written as Re, Pr, and L are the Reynolds number, Prandtl number, and length of the channel, respectively.In (8),   and   are the dynamic viscosity of fuel salt calculated using fuel salt temperature and wall temperature, respectively.
In the 3DTH, the mass flow in each fuel channel is calculated based on the assumption of equal pressure drop over all channels.The line integral of momentum conservation equation (2) along the axial direction of fuel channel can be written as ΔP  is the acceleration pressure drop along the fuel channel; ΔP  and ΔP  are the friction and gravity pressure drop of the fuel channel, respectively; ΔP  and m  are the total pressure drop and the mass flow for the fuel channel numbered with i, respectively; ΔP   is the local pressure drop at the inlet of each fuel channel, which describes the pressure drop caused by the resistance of the flow distribution plate in the bottom plenum; ΔP   is the local pressure drop at the outlet of each fuel channel, which describes the pressure drop caused by the form resistance of top plenums.ΔP   and ΔP   have the form of is resistance coefficient, which is inputted for each channel.The average pressure drop of reactor is defined as In the process of mass flow calculation, the mass flow of each fuel channel will be adjusted according to the deviation between ΔP  and ΔP  before the condition |Δ  − Δ V |/Δ V <  is met [25], as shown in Figure 1.
In present work, the temperature field of the solid region is calculated with the commercial software COMSOL Multiphysics [26].Based on the module "LiveLink for Matlab" [27] provided by COMSOL Multiphysics, the geometric modeling and the thermal coupling between the three-dimensional heat conduction and the one-dimensional single-phase flow are achieved through Matlab software.The 3DTH developed in the current study enables us to perform calculations on both unstructured and structured mesh, which is capable of accurately modeling a complex geometry.

The Coupling Method. The coupling of SCALE and 3DTH
is based on the exchange of power density distribution, temperature, and density data.As each unit (a region with a material number) used in neutronics calculation contains many meshes for 3DTH, the average temperature and fuel density of each unit are passing to SCALE.And the meshes for 3DTH calculation located in the same unit get the same power density from neutronics calculation.The schematic overview of the coupled code is presented in Figure 1.At the initialization step, both neutronics and thermal-hydraulic model will be built.In this step, uniform fuel salt temperature and density and uniform graphite temperature are set in the input file of SCALE.And uniform mass flow distribution is adopted in 3DTH as the initial condition.In the process of the iterative calculation, the information set at the initialization step is used for the first time of neutronics calculation.After the execution of neutronics, the power distribution calculated by SCALE is passed to 3DTH.And the mass flow distributions and temperature field corresponding to the power distribution are calculated.Then the temperature and density data are updated for the next step of neutronics calculation.As shown in Figure 1, the iteration is continued till the k  is converged or the maximum number of iterations is reached.

Validation of the 3DTH
To verify the reliability of the 3DTH model for the channel type MSR, the result calculated by the RELAP5 code [28] and 3DTH is compared.The RELAP5 code adopts onedimensional heat conduction model for the solid region, neglecting the axial heat conduction.A round pipe with thin tube wall cooled by the coolant is chosen as the calculation model.For the temperature distribution in the tube wall, the axial heat conduction can be neglected, and only the one-dimensional heat conduction in radial direction is dominant.In this case, the three-dimensional heat conduction in the tube wall of the round pipe can be accurately approximated to one-dimensional heat conduction.
The main parameters of the round pipe and the physical property of the coolant are listed in Table 1.In the RELAP5 code, the round pipe is equally divided into 20 nodes in the axial direction, and tube wall is uniformly divided into 9 intervals in the radial direction.Just like in the MSR, the power is predominantly released in the coolant, as listed in Table 1.Since the coolant is single-phase turbulent flow in the round pipe, the Dittus-Boelter correlation is used to calculate the heat transfer coefficient in the RELAP5 code [29], and this correlation is also added to the 3DTH.
The radial temperature distributions in tube wall at three positions (node 19, node 10, and node 1) are compared in Figure 2(a).The three positions are chosen from the top, middle, and bottom of the tube, respectively.Due to the negligible effect of the axial heat conduction, the calculated results at different positions show a relatively good agreement between the 3DTH and the RELAP5 code.The coolant temperature profiles calculated by the two codes are compared in Figure 2(b).It is found that the coolant temperatures calculated by the 3DTH are consistent with the RELAP5 code.The temperature profile of the pipe calculated by the 3DTH is shown in Figure 2(c).The tube wall temperature and coolant temperature calculated by the 3DTH are consistent with those calculated by the RELAP5 code, which proves the correctness of 3DTH.

2MW-MSR Description
In this work, we select a conceptual design of the 2MWt experimental molten salt reactor designed by TMSR center of CAS [21] as the research object for the coupled code.The vertical and horizontal view of the 2MW-MSR are presented in Figure 3(A) and Figure 3(B), respectively.The active core of the 2MW-MSR contains 85 fuel assemblies.Both height and diameter of the active core are 1.1 m. Figure 3(C) shows the hexagonal fuel assembly, the half-pitch and fuel channel radius of the assembly is 5.5 cm and 2.4 cm, respectively.In the 2MW-MSR, graphite is chosen as the moderator and reflector material, and Hastelloy-N alloy is used as the reactor vessel material and other structural materials.
The parameters of the 2MW-MSR are listed in Table 2.The fuel salt composition is LiF-BeF2-ThF4-UF4 (70.2 mol%-26.8mol%-0.073mol%-2.927mol%) with 99.95% enrichment of Li-7 and 16.2% enrichment of U-235.The inlet and outlet temperature of the core are 600 and 620 ∘ C, respectively.The rated mass flow is designed as 59.25kg/s.The fission heat   generated in the core is mainly removed by the part of fuel salt which flows through the active core.In the 2MW-MSR, a small part of fuel salt (the outer salt layer in Figure 3) which flows upward between the reflector and the reactor vessel is used to cool the structure material.
Because the fuel assembly is arranged in a hexagonal manner in the 2MW-MSR, the power difference between different assemblies in the same ring is small.Therefore, considering efficiency and accuracy of neutronics calculation, the fuel assemblies are divided into 11 units in the axial direction and 6 rings in the radial direction, where the height of each unit is 10 cm.In this work, 700 generations with 20000 neutrons per generation are adopted in the neutronic simulation.Due to the geometric symmetry of the 2MW-MSR (as shown in Figure 3(B)), only one-twelfth of the core is chosen as the calculation model for 3DTH.And only the active core, radial reflector, radial outer salt layer, and reactor vessel are considered in 3DTH, while the plenums and the reflector at the top and bottom are ignored.The crosssectional view of the calculation model is shown in Figure 4, and each fuel assembly in the model is numbered based on their radial positions.For the thermal-hydraulic calculation of the 2MW-MSR, the mass flow in each fuel assembly is calculated based on the assumption of equal pressure drop  over all the channels, and the mass flow in the outer salt layer is set as 5% of the rated flow.In the present study, the local pressure drops at the inlet and outlet of each fuel channel are ignored for the simulation.

Results and Discussions
The steady-state coupled calculations for the 2MW-MSR are performed under two different conditions.First, the normal operating state of 2MW-MSR is calculated.Then the central fuel assembly under partly blocked state is analyzed.

The Normal Operating State.
Figure 5 shows k eff as a function of iteration number.The black dot line represents the initial condition of 600 ∘ C which is set for the first neutronics calculation, and the red dot line indicates the initial condition of the melting point temperature (465 ∘ C) of fuel salt.The values of k eff are converged after a few iteration numbers.As shown in Figure 5, the converged values of k eff are not affected by initial temperature.
The power density distribution is one of the important factors influencing the energy output and the safety of the reactor.The distribution of power density in fuel salt is presented in Figure 6, which shows the peak of the power density located in the central region of the reactor core.Because some neutrons are reflected and moderated by the radial reflector, a small power peak exists in the channel nearby the reflector.Besides, the fractions of fission energy deposited in the graphite moderator, radial reflector, bottom plenum, and top plenum are 3.38%, 1.14%, 11.45%, and 12.87%, respectively.
Based on the energy fractions deposited in the plenums, the temperature rises in the bottom plenum and the top plenum are both about 2 ∘ C.During the thermal-hydraulic calculation, the inlet temperature of each fuel assembly has considered the temperature rise in the bottom plenum.
The temperature fields of both graphite and fuel salt for the 2MW-MSR are shown in Figure 7.The temperature in the solid region reaches the maximum (about 642 ∘ C) at the radial reflector, which is much lower than the boiling temperature of fuel salt.Because of the cooling effect of the outer salt layer, the reflector temperature at the outer edge is obviously lower than that of the central region.Figure 7(right) shows the one-dimensional temperature profiles of the fuel salt in each fuel channel, which illustrates that the fuel salt heats up as it flows along the axial direction.At normal operating steady state, since the graphite is cooled by fuel salt, the graphite temperature is higher than the salt temperature.As shown in Figure 7, the temperature difference between fuel salt and graphite is not huge, which reveals that the energy deposited in fuel and graphite is easily removed by the flowing salt.
To analyze the temperature field in detail, the temperature field and isotherms at the cut plane near the outlet of the fuel assembly are presented in Figure 8.The reflector temperature at the outer edge is lower than the reflector temperature at the inner edge, which indicates that the outer salt layer plays a more important role in cooling the reflector compared with the fuel channel.Because the energy deposited in the  graphite of fuel assembly is removed by the fuel salt, the isotherms densely circle around the wall of fuel channel, as shown in Figure 8(right).Besides, the shapes of isothermal line in different fuel assemblies are obviously different.The fuel assemblies (C1-C8) are not adjacent to the radial reflector.The isothermal lines in these fuel assemblies take the typical concentric circles distribution patterns, which indicates the heat flux at the interfaces of these assemblies is not obvious.The fuel assemblies (C9-C11) play an important role in removing part of energy deposited in the radial reflector.Due to the influence of the radial reflector, the shapes of isotherm in the C9-C11 depart from the circle shape.
The mass flow distribution is an important parameter for the hydraulic design of MSR. Figure 9 shows the mass flow in each fuel assembly.The mass flow distribution in the 2MW-MSR is relatively uniform.The maximum mass flow in the central fuel assembly C1 is only 1.12 times of the minimum mass flow in the C7 fuel assembly.The total power deposited in each fuel assembly is also presented in Figure 9.Because the local pressure drops for each channel at top plenum and bottom plenum are ignored and the density of fuel salt decreases with rise of temperature, more mass flow is distributed to the pipe deposited with more energy and the flow distribution nicely follows the power distribution.

The Central Fuel
Assembly Partly Blocked State.The fuel assembly blockage can be caused by the breakage of the pump or heat exchanger [30].During the fuel assembly blocked situation, the increase of temperature in the fuel salt and moderator might give rise to security risks to the structure material.In the 2MW-MSR, the partial blockage of the central fuel assembly is simulated by reducing its mass flow to 20% 1.1x10 7   9.0x10 6   7.5x10 6   6.0x10  of the average mass flow in one assembly and maintaining the rated mass flow in the core.k  as a function of iteration number is presented in Figure 5 (the blue dot line).It can be found that the converged k  is very close to the value under the normal operating state and is not obviously affected by the blockage.The reason is that the blocked fuel assembly has not seriously changed the average temperature and the temperature distribution in the reactor.
In the state of blockage, the power density distribution of fuel salt in three assemblies from the inner to the edge of the active core is shown in Figure 10.It can be found that the peak of the power density in this case is still located at the middle position of the reactor.Compared with that at the normal operating state, the power density distribution is not obviously influenced by the blockage of the central fuel assembly.
The temperature distribution in the solid region and the fuel salt in the state of blockage is displayed in Figure 11.The maximum temperature in the fuel salt and graphite is located at the outlet region of the blocked assembly.It can be found that the fuel salt temperature is substantially higher than the graphite temperature in the outlet of the central assembly, which means the direction of heat flux at the wall of fuel channel is reversed and a part of power deposited in the fuel salt is removed by the graphite.Based on the temperature difference of fuel salt and total power deposited in the blocked assembly, about 5.8% of total power in the central assembly is removed by the adjacent fuel assemblies through heat conduction.
Figure 12 gives a cross-sectional view of temperature field and isotherms at the cut plane near the outlet.The results show that the isotherms densely circle around the wall of fuel channel, especially for the blocked assembly.The temperature distribution and the shapes of isotherms in the fuel assemblies C3-C11 and in the radial reflector are similar to those presented in the normal operating state (Figure 8), which illustrates that the partial blockage of C1 only affects the temperature distribution of the adjacent fuel assembly C2.
The mass flow distribution and total power deposited in each fuel assembly are shown in Figure 13.The maximum mass flow in the fuel assembly C2 is 1.11 times of the minimum mass flow in the C7.Except for the blocked channel C1, the distribution of mass flow in this case is similar to that in the normal operating state and approximate to the power profile in the radial direction.

Conclusion
In this work, a special steady-state analysis code for the channel type MSR is developed by coupling the SCALE and the three-dimensional thermal-hydraulic model (3DTH).The 3DTH adopts a three-dimensional heat conduction model for the solid region and one-dimensional single-phase flow model for the fuel salt.Based on the coupled code, steady-state behaviors of the 2MW-MSR under the normal operating state and central fuel assembly partly blocked state are analyzed in terms of k eff , power distribution, mass flow rates, and temperature field.The main conclusions drawn from this study are as follows:   (iii) Under the normal operating steady state, the maximum temperature region is located at the radial reflector.Except for the fuel assemblies close to the radial reflector, the heat exchange through heat conduction between the fuel assemblies C1-C8 is not obvious.Besides, the mass flow distribution in the 2MW-MSR is relatively uniform, and it is similar to the power profile in the radial direction.
(iv) Under the central fuel assembly partly blocked state, the maximum temperatures in fuel salt and graphite are located at the outlet of the center assembly.In this case only the adjacent fuel assemblies are significantly influenced by the blocked assembly, and a part of the energy deposited in the blocked assembly is removed by the adjacent fuel assemblies.Except for the blocked assembly, the mass flow distribution in each fuel assembly is similar to that at the normal operating state and is not significantly influenced by the blocked assembly.

Figure 1 :
Figure 1: The flowchart of the coupled code.

Figure 2 :
Figure 2: The temperature distribution in the pipe.(a) The radial temperature distribution in the tube wall (line represents the 3DTH and square represents the RELAP5).(b) Coolant temperature distribution along the axial direction (line represents the 3DTH and square represents the RELAP5).(c) The temperature profile of the pipe (unit: K).

Figure 4 :
Figure 4: The cross-sectional view of the 1/12 core of the 2MW-MSR.

Figure 5 :
Figure 5: k eff as a function of iteration number (black dot line: initial temperature is set at 600 ∘ C during the normal operating state; red dot line: initial temperature is set at 465 ∘ C during the normal operating state; blue dot line: initial temperature is set at 465 ∘ C during the central fuel assembly partly blocked state).

Figure 7 :
Figure 7: Temperature distribution in the solid region (left) and fuel salt (right) (unit: ∘ C).

Figure 8 :
Figure 8: The cross-sectional view of temperature field (left) and isotherms (right) in the solid region near the outlet (z=1.0m)(unit: ∘ C).

Figure 9 :
Figure 9: The mass flow and power in each fuel assembly.

Figure 10 :
Figure 10: Power density distribution of fuel salt for the C1, C5, and C10 assemblies (line represents the situation of central fuel assembly partly blocked and square represents normal operating state).
(i) The correctness and reliability of the 3DTH code are verified by the RELAP5 code by a single pipe model.The tube wall and coolant temperature calculated by 3DTH are consistent with those calculated by RELAP5.(ii) k eff of the 2MW-MSR is converged after a few times of iteration.The converged values of k eff are not obviously affected by initial temperature and the partial blockage of the central fuel assembly.Besides, the power density distribution in the 2MW-MSR is not obviously influenced by the blockage.

Figure 12 :Figure 13 :
Figure 12: The cross-sectional view of temperature field (left) and isotherms (right) in the solid region near the outlet (z=1.0m)(unit: ∘ C).

Table 1 :
Parameters used for the calculation model.