A Review of Molten Salt Reactor Multi-Physics Coupling Models and Development Prospects

Abstract

Molten salt reactors (MSRs) are one type of GEN-IV advanced reactors that adopt melt mixtures of heavy metal elements and molten salt as both fuel and coolant. The liquid fuel allows MSRs to perform online refueling, reprocessing, and helium bubbling. The fuel utilization, safety, and economics can be enhanced, while some new physical mechanisms and phenomena emerge simultaneously, which would significantly complicate the numerical simulation of MSRs. The dual roles of molten fuel salt in the core lead to a tighter coupling of physical mechanisms since the released fission energy will be absorbed immediately by the molten salt itself and then transferred to the primary heat exchanger. The modeling of multi-physics coupling is regarded as one important aspect of MSR study, attracting growing attention worldwide. Up to now, great efforts have been made in the development of MSR multi-physics coupling models over the past 60 years, especially after 2000, when MSR was selected for one of the GEN-IV advanced reactors. In this paper, the development status of the MSR multi-physics coupling model is extensively reviewed in the light of coupling models of N-TH (neutronics and thermal hydraulics), N-TH-BN (neutronics, thermal hydraulics, and burnup) and N-TH-BN-G (neutronics, thermal hydraulics, burnup, and graphite deformation). The problems, challenges, and development trends are outlined to provide a basis for the future development of MSR multi-physics coupling models.

1. Introduction

Molten salt reactors (MSRs) are liquid-fueled reactors characterized by their capabilities of online refueling, reprocessing, and helium bubbling to remove sparingly soluble fission products. The high utilization of thorium can be achieved in MSRs since ²³³Pa, which is an intermediate nuclide in the evolution chain of ²³²Th to ²³³U and has a half-life of 27 days, can be extracted from the core in a timely manner to reduce its neutron absorption [1]. Liquid-fueled MSRs are, hence, often associated with the ²³²Th-²³³U fuel cycle. Th-based fuel cycles are well known to have lower equilibrium radiotoxicity than U-based fuel cycles due to the much lower transuranic production from ²³²Th than from ²³⁸U. Therefore, MSRs also have the potential to greatly reduce the inventory of long-lived high-level radioactive wastes [2]. Other merits derived from the liquid fuel, including low excess reactivity over the core operation cycle and exemption from fuel cell manufacture, can also be provided by MSRs. Because of these unique features, MSRs have attracted growing attention worldwide, especially after they were selected as one of the GEN-IV reactor concepts in 2002 [3]. A series of MSR concepts, sorted by the neutron spectrum (thermal/epithermal/fast), molten salt (fluoride/chlorine), moderator (graphite/heavy water/ZrH), and power scale (micro/small-modular/large), have been proposed since the molten salt reactor experiment (MSRE), with a thermal power of 8 MW, was built and operated successfully in the 1960s by the Oak Ridge National Laboratory (ORNL) [4,5,6,7,8,9,10,11,12].

As in solid-fueled reactors, MSR engineering is a complex intersection of multiple physics, encompassing neutronics, thermal hydraulics, chemistry, and materials science. Besides this, it possesses other unique physical mechanisms, including fuel flow and continuous fission-product removal using physical (helium sparging) and pyrochemical processes, etc. [13,14,15,16]. A closer mutual interaction among the multiple physical mechanisms compared with solid-fueled reactors can also be observed in MSRs since fuel salt acts as both fuel and coolant [5]. A multi-physics coupling simulation is, hence, of great importance in MSR numerical calculation, directly influencing the core design and safety analysis. Since the 1960s, great efforts have been made regarding the development of an MSR multi-physics coupling model, especially in the aspects of neutronics and thermal hydraulics (N-TH) coupling. The point neutronic kinetics coupling with lumped parameter thermal hydraulics was the first N-TH coupling model developed by ORNL in the 1960s [17]. The delayed neutron precursor flowing out from the core was analytically modeled, based on the core residence time and external loop transit time. Afterward, the spatial effects of the MSR parameters were gradually taken into account, along with the approach employed in the neutronics and thermal hydraulics of nuclear reactors [18,19,20]. Furthermore, since 2002, various coupling solution strategies, methods, and solution methods were introduced to the MSR N-TH coupling model, causing the N-TH analysis model to approach a real MSR core operational condition with little approximation necessary [21,22,23,24,25]. Meanwhile, coupling other physical mechanisms, such as online refueling and reprocessing, helium bubbling, and dimensional changes of graphite, to the N-TH model is also currently being explored, thanks to a great improvement in computing technologies [26,27,28,29]. Nevertheless, the existing multi-physics coupling models of MSR should be further improved in the aspects of completeness and calculation cost.

In this paper, research activities regarding the development of an MSR multi-physics coupling model are extensively reviewed, aiming to identify the research trends and challenges. The paper is expected to provide a foundation/reference for those scholars who are interested in an MSR multi-physics coupling study. A comprehensive view of MSR research and development (R&D) worldwide and the main multi-physics mechanisms of MSR are provided in Section 2. The MSR multi-physics coupling models are then reviewed in Section 3. Section 4 gives the prospects for the development of an MSR multi-physics coupling model. Finally, our conclusions are drawn in Section 5.

2. R&D and Main Features of MSR

2.1. MSR R&D Worldwide

MSR is an old concept; its R&D activities have continued for around 70 years and have undergone ups and downs but have gained a great deal of attention worldwide since the beginning of the year 2000. An overview of the historical and current status of MSR R&D activities worldwide is shown in Figure 1. In the 1950s, ORNL built its first MSR, known as the aircraft reactor experiment (ARE), under the support of the US Aircraft Nuclear Propulsion program. This MSR adopted molten fluoride salt, NaF-ZrF₄-UF₄ (53-41-6 mol%), as the fuel, with beryllium oxide (BeO) as the moderator. It achieved criticality in 1954 and operated successfully for 100 h at a steady state [30]. Researchers then turned to MSR civilian applications for electricity production in the 1960s. A molten salt reactor experiment (MSRE) with a thermal power of 7.3 MW was built by ORNL in 1965 and was successfully operated for over 4.5 years by sequentially using ²³⁵U, ²³⁹Pu, and ²³³U fuels [6,31]. In parallel with the R&D of MSRE, the design activities of breeder MSRs were carried out by ORNL with the aim of breeding ²³³U from ²³²Th. From 1960 to 1968, a series of molten salt-breeder reactor (MSBR) concepts were proposed. A two-region two-fluid graphite-moderated and reflected thermal reactor concept with an electrical power of 1000 MW was the first proposed design of MSBR. In this design, two FLiBe salt fluids, dissolved with ²³³U and ²³²Th, respectively, were separated by the graphite structure, forming an inner core driver zone and an outer core blanket zone [32]. For both fluids, continuous fuel reprocessing and off-gas systems were planned. The Th-U breeding ratio (BR) could achieve around 1.06. Thereafter, several MSBR designs, including a modular MSBR (MMSBR), MSBR (Pa-Pb), a molten-salt epithermal breeder (MOSEL) (Pa-Pb), and a single-stream core-breeder (SSCB) were developed. The MMSBR is a small, modular, two-region two-fluid MSBR with an electrical power of 250 MW. Four modules of MMSBR form one plant. The plant availability factor can be improved since the stoppage of a fuel pump shuts down only one-quarter of the station capacity [33]. The MSBR (Pa-Pb) online extracts the Pa from the blanket salt and uses direct-contact cooling of the molten-salt fuel with molten lead within the reactor vessel, since lead is immiscible with molten salt. The fissile inventory external to the reactor core can be significantly reduced. MOSEL (Pa-Pb) is similar to MSBR (Pa-Pb), but its core is free from graphite and has an intermediate-to-fast neutron energy spectrum. Regarding SSCB, the fissile and fertile materials are contained in the fuel stream, and this fuel stream is surrounded by a blanket of thorium-containing salt [33].

In order to simplify the core structure, ORNL proposed a single-fluid MSBR with both fissile and fertile materials, incorporated in an FLiBe molten salt, in 1968. The core has an inner core zone, 13% of which is occupied by fuel salt, and an outer core zone with the fuel salt volume fraction at 37%, to enhance Th-U breeding in this region [34]. This MSBR program was closed down in 1972, when the USA consolidated the fuel breeder research into sodium-cooled fast reactors, to pursue higher U-Pu breeding. Nevertheless, modest MSR research activities on core design, materials development, fuel chemistry, and reprocessing were continued at ORNL until the early 1980s. During this period, another MSR design, the denatured molten salt reactor (DMSR), was proposed to address nuclear proliferation concerns [6]. DMSR is a direct outgrowth of the single-fluid MSBR and inherits many merits from MSBR. However, it presents many differences, including using enriched uranium (19.75 wt % ²³⁵U), adding ²³⁸U to denature the fuel, and does not perform continuous fuel reprocessing to comply with the rules of nuclear nonproliferation [35].

Apart from the USA, from the 1960s to the 1980s, the United Kingdom, Japan, and China have also made efforts in terms of MSR designs and experiments, especially the UK, in which the MSR program was initiated in 1964 and lasted for almost ten years. A 2.5 GWe lead-cooled molten salt fast reactor (MSFR) concept, with the core filled with chloride fuel molten salt, was proposed by the UK Atomic Energy Research Establishment (AERE) [36]. Japan started research into MSR in the 1980s and proposed a thorium molten salt nuclear energy synergetic system (THORIMS-NES) by the research group of Furukawa. In this system, the molten salt reactor (MSR-FUJI), which has a thermal power of 350 MW and contains three core regions to flatten power distribution, combined with the accelerator molten salt breeder (AMSB) to close the Th-U fuel cycle [6,8]. In China, the Shanghai Institute of Applied Physics and the Chinese Academy of Sciences (SINAP) worked towards building a 25 MWe MSR and established a critical facility in the 1970s, but this endeavor soon gave way to the Qinshan PWR project [37,38].

The R&D activities of MSR almost stagnated in the 1990s against the background of the nuclear industry depression induced by the Chernobyl nuclear accident. This trend, however, was reversed in the 2000s. Under the support of projects such as ISTC#1606 and ISTC#3749, Russia proposed a 2400 MWt molten salt actinide recycler and transmuter (MOSART), oriented to transuranic (TRU) transmutation. The core of this reactor is filled with the molten salt dissolved with TRU fuel and has a fast neutron spectrum [39]. In parallel, the European Union put forward a 3000 MWt two-flow fluid molten salt fast reactor (MSFR), under the EURATOM framework programs. The outer core region of MSFR is filled with the blanket salt fluid (LiF-ThF₄), functioning to reduce the neutron leakage and to breed ²³³U as well. The fuel salt, composed of LiF and (Th + ²³³U/TRU/LEU)F₄, circulates around the central core of the MSFR for power generation [12,40]. In the 2010s, the R&D activities of MSR were further enhanced worldwide. China launched a “Thorium Molten Salt Reactor Nuclear Energy System” project in 2011 to develop MSR technologies within 20–30 years, to realize effective thorium energy utilization [6,41]. The 168 MWe small modular liquid-fueled thorium molten salt reactor, TMSR-LF, was proposed, which adopts graphite as the moderator. Mature technologies adopted in this design are expected to facilitate its rapid construction and deployment in the near future [5]. Small modular MSRs have also drawn a great deal of attention from energy companies in the UK, Canada, and the United States. In 2015, the UK Moltex Energy Company proposed a 375 MW stable salt reactor (SSR), which is a hybrid reactor with a liquid fuel salt mixture contained in the fuel assemblies, cooled by NaF-KF-ZrF₄ salt [42]. The companies of Dual Fluid Energy Inc. and Terrestrial Energy in Canada, respectively, put forward a Dual Fluid Reactor (DFR), which adopts two separate liquid cycles, with one for fuel (molten salt fuel) and the other for the coolant (lead) [43], and an integral molten salt reactor (IMSR) with a core similar to DMSR [44]. An LFTR company from the United States developed a liquid fluoride thorium reactor (LFTR), with a core containing two-flow fluids, moderated by the graphite [45]. The Transatomic Power Corporation, from the USA, proposed a trans-atomic power (TAP) MSR, which uses zirconium hydride as the moderator, with the aim of achieving a higher burnup by using the same uranium fuel as that in the light water reactors [46]. The USA-based Thorcon company proposed a ThorCon MSR, based on the existing technologies of MSRE, to allow rapid deployment [9]. TerraPower LLC from the USA developed a molten chloride fast reactor (MCFR) with the core in the fast neutron spectrum, filled by molten chloride fuel salt [47]. The Kairos Power company from the USA proposed a fluoride salt-cooled high-temperature reactor (FHR), which is a novel advanced reactor technology that integrates the merits of the TRISO fuel in pebble form and the low pressure of fluoride salt coolant [48]. Copenhagen Atomics in Denmark and SINAP in China, respectively, proposed a 50 MW and a 2250 MW heavy water-moderated molten salt reactor (HWMSR), which uses the heavy moderator to address the management issue of the high-level radioactive-discharged graphite moderator [11,49].

Although many different types of MSRs were proposed worldwide, the graphite-moderated MSRs and MSFR are two main MSRs that are expected to be built in the near future and have gained much additional attention regarding core design, the nuclear fuel cycle, and safety analysis. Hence, this review focuses on studies of multi-physics coupling models for these two types of MSRs.

2.2. Main Operational and Physical Mechanisms of MSRs

Figure 2 demonstrates a typical graphite-moderated MSR system. This configuration can also represent the scenario of an MSFR, except for the free flow of fuel salt in the core, which is stronger in turbulence since the guidance offered by graphite tubes does not exist in the core. During core operation, the liquid fuel salt goes through the MSR core, where nuclear fission occurs, with the released heat absorbed by the molten salt. The heated fuel salt then flows to the heat exchangers, with the heat transferred to the secondary loop for electricity generation or other applications such as H₂ production and sea water desalination. To enhance neutronic economy, one bypass flow is guided to a helium bubbling system, where the fission gases and sparingly soluble noble metal elements entrained in the helium bubbles, including Kr, Xe, Se, Nb, Mo, Tc, Te, etc., are removed by porous carbon absorbers. Meanwhile, helium is continuously bubbled to the core [50]. Another bypass leads to the reprocessing system for the online recycling of heavy metals (HMs) and for removing soluble fission products. Fluorination is first performed to recover isotopes of U, Np, and Pu in fluoride form. Part of the residual HMs are thereafter separated from the molten salt by a reductive extraction process and are then stored in a stockpile for several months, to let the ²³³Pa decay to ²³³U. Fluorination is performed again to collect the decayed ²³³U, part of which, together with other HMs (such as Np, Am, and Cm), depending on choice, are then returned to the core to maintain the core critical operation. The fluoride salt and thorium are finally recycled via distillation and electrodeposition, respectively, leaving fission products (FPs) for storage and disposal [7,51]. In a severe accident, the fuel salt will be discharged quickly into the fuel salt drain tanks by thawing the freeze valve, to mitigate the consequences of the accident.

The unique operational processes of MSRs, as discussed above, lead to more tight and complicated couplings among physical mechanisms, compared with the solid-fueled reactors, as presented in Figure 3. The liquid fuel makes neutronics and thermal hydraulics coupling tighter in MSRs. Thus, the temperature rise in the core of MSRs would be more rapid in accidents concerning reactivity insertion. The flowing of fuel salt, on the one hand, causes the decay heat to spread through the primary loop, influencing the profile of fuel salt temperature rise and distribution in an accident. On the other hand, it results in the delayed neutron precursor flowing out from the core, reducing the delayed neutron fraction of the core, and subsequently affecting the safety control. Furthermore, the flowing of fuel salt leads to the nuclides being entrained in the molten salt that is continuously circulating through the core, where the nuclides experience neutron flux irradiation and nuclide decay, and the external core loop, including the pipes, pumps, and heat exchangers, in which only nuclide decay occurs. The nuclide evolution law in MSRs is, hence, significantly different from that in the solid-fueled reactors. Our previous studies demonstrated that the concentration of ¹³⁵Xe, one of the main core neutron poisons, is reduced by more than half when the flowing effect of fuel salt is considered [13]. Another neutron poison, ¹⁴⁹Sm, also presents a different behavior from solid reactors in nuclide concentrations during core transience, under the condition of fuel salt flow, such as core power rising and a scram procedure [14].

Graphite-related issues, including graphite deformation, graphite heat transfer, and graphite permeation by the fission gases and fuel salt are also important physical mechanisms that should be considered. The graphite in the core acts as the moderator and also the predominant structural material. It must endure high neutron flux irradiation under high-temperature conditions. Lattice displacement occurs and is aggravated by burnup, which subsequently influences the strength of the extension, thermal expansion, Young’s modulus, etc. The dimension of graphite would experience a change with the burnup in such a way that it first shrinks and then expands, complying with the function of the influence of fast neutrons and core operating temperature. Studies conducted by ORNL in the 1960s demonstrated that graphite shrinkage in the core can reach around 4.5–7%Δ v/v [52]. This significant deformation of the graphite in the core will change both the neutronic performance, by affecting the neutron moderation capability, and the thermal-hydraulic performance, by influencing the flow field. Energy is deposited in the graphite due to neutron moderation and gamma irradiation. Removing this heat by the fuel salt forms a heat transfer mechanism between the graphite and the fuel salt. This is of importance to the thermal-hydraulic calculation and directly determines the graphite temperature. Fission gas and fuel salt penetration of the graphite is another important mechanism and is also a predictable event, due to the porous properties of graphite, significantly influencing the neutronic performance of the core. This mechanism is accompanied by helium bubbling, making the simulation more complicated.

Online refueling and reprocessing can be regarded as the external mechanisms that change the intrinsic burnup course, influencing the composition contained in the fuel salt and core reactivity. As fuel elements are added to the core, the concentrations of corresponding nuclides will increase abruptly, which subsequently changes the concentrations of their daughter nuclides. Online reprocessing is an inverse process of online refueling. Nuclides, including fission products/actinide nuclides, that are extracted online from the fuel salt would decrease the concentrations of the related nuclides.

The above mechanisms, derived from graphite, fuel salt, both in the active core and external loop, and the online refueling and reprocessing couple with each other, introducing great challenges in the physical simulation of MSRs.

3. Development of MSR Multi-Physics Coupling Models

As performed in solid-fueled reactors, MSR multi-physics coupling refers to coupling among neutronics, thermal hydraulics, burnup, and mechanics/materials, but meanwhile, this necessitates taking into account the unique physical mechanisms of MSRs, as shown in Figure 4. These are coupled via power and temperature distribution, and we must consider the loss of the delayed neutron precursor due to the flowing of fuel salt. Neutronics, meanwhile, couples with the burnup by providing a neutron flux, while obtaining the nuclide composition from burnup. The mechanisms of fuel flowing, online refueling and reprocessing, and helium bubbling, which derail the nuclide evolution from its intrinsic course, should be taken into account in this coupling. Neutronics couple with the mechanism of the graphite cell through graphite deformation. Thermal hydraulics couples burnup and the graphite cell by offering flow-field and temperature distribution, while the graphite cell couples with burnup via nuclide permeation and evolution in the graphite. Until now, many studies have been conducted on N-TH coupling, N-TH-BN (burnup) coupling, and N-TH-G (graphite cell) coupling, but a complete multi-physics coupling, i.e., a N-TH-BN-G coupling, has not been reported. Hence, more efforts are needed for the development of a multi-physics coupling model (MPCM). The rest of this section gives an overview of the historical and current status of MPCM development, aiming to outline the problems and challenges of MPCMs.

3.1. N-TH Coupling

Since MSRs were proposed by ORNL in the 1960s, remarkable progress has been made in the study of MSR N-TH coupling. Models with dimensionality extending from 0D (0-Dimension) to 3D (3-Dimension) have been proposed sequentially, with the modeling accuracy significantly improved.

Table 1. Summary of the current N-TH coupling models.

Dim.	Codes/Developers/Platforms, Year	Neutronics Model	Thermal Hydraulics Model	Coupling Techniques	Application
0D	MURGATROYD, 1962 [17]; Sides, 1971 [53]; Lapenta et al., 2001 [20]; Dulla et al., 2004 [16]	Point kinetics	Lumped parameters	External	Thermal/fast MSRs
	Auwerda, 2007 [54]; Suzuki et al., 2008 [55]; Guo et al., 2013 [56]; Cai et al., 2014 [57]	Point kinetics	Single-channel	External	Thermal MSRs
	Improved RELAP5, 2016 [58]	Point kinetics	Multi-channel	External	Thermal MSRs
1D	Cinsf1D, 2003 [59]; DYN1D-MSR, 2005 [60]	1D neutron diffusion	Single-channel	External	Thermal MSRs
2D	Zhang et al., 2009 [61]	2D cylindrical neutron diffusion	Modeling heat transfer between the fuel salt and boundary graphite	External	Fast MSRs
	Improved SIMMER-III, 2006 [62]; Nicolino et al., 2008 [63]	2D cylindrical neutron diffusion	2D CFD	Internal	Fast MSRs
	Cammi et al., 2011 [64]	2D cylindrical neutron diffusion	2D CFD	External	Fast MSRs
	Yamamoto et al., 2006 [65]	2D cylindrical neutron diffusion	Porous medium	External	Thermal MSRs
3D	DYN3D-MSR, 2007 [66]; MOREL, 2015 [67]; TMSR-3D, 2022 [68]	3D neutron diffusion	Multi-channel	External	Thermal MSRs
	DALTON and THERM, 2009 [69]	3D neutron diffusion	Multi-channel and 3D heat transfer in graphite	External	Fast MSRs
	Nagy et al., 2014 [70]	3D neutron diffusion	Porous medial	External	Fast MSRs
	DALTON and HEAT, 2012 [23]; PHANTON-SN and DGFlows, 2020 [24]	3D neutron diffusion	3D-CFD	External	Fast MSRs
	OpenFOAM, 2014 [21]	3D neutron diffusion	3D-CFD	Internal	Fast MSRs
	COMSOL, 2020 [71]	3D neutron diffusion	3D two-phase flow	Internal	Fast MSRs
	Moltres under MOOSE framework, 2017 [25]	3D neutron diffusion	3D-CFD	Internal	Thermal MSRs
	PROTEUS-NODAL and SAM under MOOSE framework, 2022 [72]	3D neutron transport	Multi-channel for core and 1D for external components	External	Fast/Thermal MSRs

gives an overview of the N-TH coupling models.

3.1.1. Point Neutronic Kinetics (0-Dimension N) Coupling with Lumped TH Parameters (0-Dimension/1-Dimension/2-Dimension TH)

A point neutronic kinetics model is widely used in the field of reactor safety analysis due to its significant advantages in terms of computation cost. It describes the time behavior of a reactor system with point neutronic kinetics equations, in which the spatial shape of the neutron flux is assumed to be unchanged with time:

$$\frac{dn\left(t\right)}{dt}=\frac{\rho \left(t\right)-\beta }{\mathsf{\Lambda }}n\left(t\right)+{\displaystyle \sum }_{i=1}^N{\lambda }_i{C}_i\left(t\right)$$

(1)

$$\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_i}{\mathsf{\Lambda }}n\left(t\right)-{\lambda }_i{C}_i\left(t\right)$$

(2)

where $\rho$ is the core reactivity, which is determined by the temperatures of the fuel, coolant and moderator; $\beta$ is the effective fraction of delayed neutrons; $\mathsf{\Lambda }$ is the neutron generation lifetime; ${\lambda }_i$ is the DNP decay constant of family i; ${\beta }_i$ is the delayed neutron fraction of family i; N is the total number of DNP family, which is usually set to be six. In Equations (1) and (2), $n\left(t\right)$ is always written as power P(t), since the core power is linearly related to the magnitude of the neutron flux. In addition, a term correlating $C_i\left(t\right)$ with the fuel salt residence time in the core and external loop should be added in Equation (2) to account for the effect of DNP loss, caused by the flowing of fuel salt in an MSR.

Usually, a loose coupling scheme is adopted for the coupling of the point neutronic kinetics model with a TH model, as presented in Figure 5. The point neutronic kinetics model provides power to the thermal-hydraulic model during one time step, then the thermal-hydraulic model feeds back the temperature distribution to the point neutronic kinetics model at the next step, and so forth until the preset time steps come to an end. This coupling scheme significantly reduces the computational cost, compared with the tight coupling scheme, in which the iteration between neutronics and thermal hydraulics is continued over one step until it converges.

A point neutronic kinetics model was adopted by ORNL for MSR calculation in the 1960s. It was coupled with the lumped parameters of the TH model, forming the first MSR dynamic calculation tool, MURGATROYD [17], in which the same differential equations of point neutronic kinetics model as are used in solid fuel reactors are adopted, i.e., Equations (1) and (2), and the flow effect of fuel salt on delayed neutron fraction is simply considered by revising the fraction of delayed neutrons of the ith group, $v_i,$ in the core:

$$v_{i }=1-\frac{1-e^{-{\lambda }_i{\tau }_C}}{{\lambda }_i{\tau }_C}·\frac{1-e^{-{\lambda }_i{\tau }_L}}{1-e^{-{\lambda }_i({\tau }_C+{\tau }_L)}}$$

(3)

where ${\tau }_{\mathrm{C}}$ and ${\tau }_{\mathrm{L}}$ are the residence time of fuel salt in the core and external loop, respectively, which are calculated based on the constant neutron flux and slug flow under the steady state of the core. Therefore, this simple revision cannot accurately describe the time behavior of the delayed neutron in the core, especially in the case of transients and accidents. The TH model of MURGATROYD divides the core into a fuel salt block and a graphite block. The heat transfer between these two blocks is calculated. The coupling between the point neutronic kinetics model and thermal-hydraulic model is then realized by feeding back the temperature of fuel salt and graphite:

$$k_{\mathrm{eff}}=1+{\rho }_{in}+b\cdot t-\left|\frac{\partial k_{\mathrm{eff}}}{\partial T_{\mathrm{f}}}\right|\left(T_{\mathrm{f}}-T_{\mathrm{f}0}\right)-\left|\frac{\partial k_{\mathrm{eff}}}{\partial T_{\mathrm{g}}}\right| \left(T_{\mathrm{g}}-T_{\mathrm{g}0}\right)$$

(4)

where $k_{\mathrm{eff}}$ is the effective multiplication of the core and is related to the core reactivity $\rho$ in Equation (1), with a function of $\rho$ = ($k_{\mathrm{eff}}$ − 1)/$k_{\mathrm{eff}}$; ${\rho }_{in}$ is the initial step reactivity input, $b$ is the initial ramp reactivity input; $T_{\mathrm{f}}$and $T_{\mathrm{g}}$ are the temperatures of fuel salt and graphite, respectively; $T_{\mathrm{f}0}$ and $T_{\mathrm{g}0}$ are the temperatures of fuel salt and graphite in the steady state, respectively.

In 1971, in order to analyze the safety performance of MSBR, the dynamic calculation model was improved in terms of both the aspects of neutronics and thermal hydraulics, with the ability to model the transient behavior of the whole reactor system. The improved model accounts for the flowing effect by adding the terms of the delayed neutron precursor flowing out from the core (the third term on the right side of Equation (5)) and entering the core from the external loop (the fourth term on the right side of Equation (5)) to the delayed neutron precursor time-dependent differential equation:

$$\frac{d{C}_i\left(t\right)}{dt}=\frac{{\beta }_i}{\mathsf{\Lambda }}P\left(t\right)-{\lambda }_i{C}_i\left(t\right)-\frac{1}{{\tau }_C}{C}_i\left(t\right)+\frac{e^{-{\lambda }_i{\tau }_{\mathrm{L}}}}{{\tau }_C}{C}_i\left(t\right)\left(t-{\tau }_L\right).$$

(5)

However, the DNPs are divided into only two groups in the improved model. The TH lumped-parameter model was extended from the core to the whole reactor system by using the nodal method, as shown in Figure 6. The graphite and fuel salt in the core are modeled to be several nodal blocks, based on which heat transfer between graphite and fuel salt and heat delivery via fuel salt flowing was simulated. The temperatures of the fuel salt and graphite in the core can then be obtained by averaging the temperature of the nodal blocks. In the same way, the temperatures of the working fluid in the secondary loop and energy conversion system were calculated, in terms of time. The control system was modeled by a set of equations describing the adjustment of the flow rate and control rod behavior to control the steam temperature and reactor outlet temperature during transient events [53]. It being capable of simulating the transient behavior of the whole system is one main merit of this improved model, which has been widely referenced by the subsequent MSR dynamic codes.

The neutronics aspect of the MSR point dynamic model was improved by Lapenta et al. in 2001 [20]. Neutrons in the different core spatial positions for different energy groups will induce nuclear fission with different probabilities since nuclear fission cross-sections vary with neutron energy and a strong neutron leakage exists near the boundary of the core. Neutrons in the center of the core usually have higher importance than those in the peripheral core region, in terms of power generation. This neutron importance can be described by the adjoint neutron flux, which can be calculated by solving the adjoint neutron flux conservation differential equations. To account for the importance of fission neutrons, including both prompt and delayed neutrons, a normalization factor (F) was introduced to amend the parameters involved in the point neutronic kinetics model:

$$F={\displaystyle \sum }_{i=1}^R{\displaystyle \sum }_{n=1}^G⟨{\mathsf{\Phi }}_n^{†}|{\chi }_{i,n}{\lambda }_i{C}_{i,0}⟩+\left(1-\beta \right){\displaystyle \sum }_{n=1}^G{\displaystyle \sum }_{g=1}^G⟨{\mathsf{\Phi }}_n^{†}|{\chi }_n{(\nu {\mathsf{\Sigma }}_f)}_{g,0}{\mathsf{\Phi }}_{g,0}⟩$$

(6)

where ${\mathsf{\Phi }}_n^{†}$ is the adjoint neutron flux; $R$ and $G$ are the groups of the delayed neutron precursor and neutron flux, respectively; ${\chi }_n$ and ${\chi }_{i, n}$ indicate the probability of fission neutrons in the n-th energy group and neutrons released from the i-th delayed neutron precursor group in the n-th energy group, respectively; $\nu$ is the average neutrons generated in nuclear fission. Taking the DNP of family i as an example, the parameters of the point kinetics model are amended thus:

$${\tilde{\beta }}_i=\frac{1}{F}{\displaystyle \sum }_{n=1}^G⟨{\mathsf{\Phi }}_n^{†}|{\chi }_{i,n}{\lambda }_i{C}_{i,0}⟩.$$

(7)

In 2004 [16], a quasi-static method was employed to account for the spatial shape change of the neutron density and DNP concentration with time in the point neutronic kinetics model. The quasi-static method adopts the classic factorization projection scheme to solve the time-dependent neutron and adjoint neutron flux equations. Weighting vectors, in terms of neutron flux and delayed neutron precursor concentration, are obtained to correct the parameters during an interval step time, Δt. In this way, transients with a distortion of spatial neutron density shape and delayed precursor concentration can be simulated by using the point kinetics model, but this improvement is limited to those cases with a large change in neutron density shape.

In subsequent studies of MSR transient calculations associated with the point dynamic model, more attention was paid to the development of thermal-hydraulic models. A single-channel model (1D model) for thermal-hydraulic calculation was developed to couple point neutronic kinetics model by Auwerda in 2007 [54], Suzuki et al. in 2008 [55], Guo et al. in 2013 [56], and Cai et al. in 2014 [57]. It was then extended to a multi-channel model (2D model) by Shi et al. in 2016, via improving RELAP5, which is a point dynamic code designed for solid-fueled reactors and has the ability to perform system transient analysis by using the nodal method. The loss of a delayed neutron precursor due to the flowing of fuel salt was modeled by revising the DNP differential equations, as expressed in Equation (5) [58].

From the above, it can be concluded that point dynamic models have played an important role in MSR safety analysis and have attracted attention worldwide, even in recent years. Thermal-hydraulic models with 0–2 dimensions were also developed, and some efforts have been made to address the drawbacks when describing the shape change of the neutron flux in the core. Nevertheless, great errors and uncertainties would be inevitably introduced for the analysis of accidents with a large and abrupt shape change in terms of neutron flux and local reactivity insertion.

3.1.2. The 1-Dimension Neutronics Coupling with 1-Dimension Thermal Hydraulics

In parallel with the development of the MSR point dynamic model in the years after 2000, a 1D dynamic model was proposed with the aim of improving the accuracy of transient analysis. Two energy groups, i.e., the thermal and fast neutron energy groups, are generally adopted in the 1D time-dependent neutronic diffusion equations:

$$\begin{array}{c}\frac{1}{{\nu }_1}\frac{\partial {\mathsf{\Phi }}_1\left(z, t\right)}{\partial t}=\frac{\partial }{\partial z}\left(D_1\frac{\partial {\mathsf{\Phi }}_1\left(z, t\right)}{\partial z}\right)+\left(1-\beta \right)\left[{\nu }_1{\Sigma }_{f1}\left(z\right){\mathsf{\Phi }}_1\left(z, t\right)+{\nu }_2{\mathsf{\Sigma }}_{f2}\left(z\right){\mathsf{\Phi }}_2\left(z, t\right)\right]\\ -\left[{\mathsf{\Sigma }}_{a1}\left(z\right)+{\mathsf{\Sigma }}_s\left(z\right)+D_1{B}^2\right]{\mathsf{\Phi }}_1\left(z, t\right)+{\displaystyle {\displaystyle \sum }_{i=1}^M}{\lambda }_i{C}_i\left(z, t\right)\end{array}$$

(8)

$$\frac{1}{{\nu }_2}\frac{\partial {\mathsf{\Phi }}_2\left(z, t\right)}{\partial t}=\frac{\partial }{\partial z}\left(D_2\frac{\partial {\mathsf{\Phi }}_2\left(z, t\right)}{\partial z}\right)-\left[{\mathsf{\Sigma }}_{a2}\left(z\right)+D_2{B}^2\right]{\mathsf{\Phi }}_2\left(z, t\right)+{\mathsf{\Sigma }}_s\left(z\right){\mathsf{\Phi }}_1\left(z, t\right)$$

(9)

$$\frac{\partial C_i\left(z, t\right)}{\partial t}={\beta }_i\left[{\nu }_1{\mathsf{\Sigma }}_{f1}\left(z\right){\mathsf{\Phi }}_1\left(z, t\right)+{\nu }_2{\mathsf{\Sigma }}_{f2}\left(z\right){\mathsf{\Phi }}_2\left(z, t\right)\right]-{\lambda }_i{C}_i\left(z, t\right)-\frac{V\left(z, t\right)\partial C_i}{\partial z},i=1,\dots ,M$$

(10)

where the subscripts “1” and “2” indicate “fast” and “thermal”, respectively; $\mathsf{\Phi }$ is the neutron flux; $\nu$ is the speed of neutrons; $D$ is the neutron diffusion coefficient; ${\mathsf{\Sigma }}_f$ is the macroscopic fission cross-section; ${\mathsf{\Sigma }}_a$ is the macroscopic absorption fission cross-section; ${\mathsf{\Sigma }}_s$ is the macroscopic scattering cross-section from the fast group to the thermal group; $V$ is the flow speed of the molten salt in the core; M is the groups of delayed neutrons, which is generally set to be six. The terms in the right side of Equation (8) represent the time change of the fast neutron flux caused by neutron diffusion, prompt neutron (thermal and fast) fission, neutron absorption, neutron scattering, neutron leakage, and delayed neutron fission, respectively. Equation (9) describes the time change of the thermal neutron flux, in terms of neutron diffusion, neutron absorption, neutron leakage, and neutron scattering from the fast neutrons. The time change of the delayed neutron precursor is specified by the neutron groups, while the axial position is as expressed in Equation (10). With this model, this z-axial shape change of neutron density under transients and accidents can be simulated.

The 1D neutron kinetic model was generally coupled with a single-channel model of thermal hydraulics. Cinsf1D (Cinétique pour réacteur à sels fondus, 1D) proposed by Lecarpentier et al. in 2003 is one typical code of this model [59]. It adopts a time-dependent two-energy-group diffusion theory, as shown in Equations (8)–(10), for neutronic calculation. The effective neutronic cross-sections associated with the temperatures of the molten salt and graphite were prepared in advance by performing the neutronic transport calculation of fuel cells. The core was modeled as a single channel in which the graphite prism was treated as a cylinder, with the fuel salt located in the centerline for thermal-hydraulic calculation. The time-dependent equations, in terms of axial fuel salt temperature, graphite temperature, and heat transfer between them were solved to feed back the core temperature distribution to the neutronics calculation [59]:

$$\frac{\partial T_s\left(z, t\right)}{\partial t}=\frac{P_f\left(z, t\right)}{\pi R^2{\rho }_s{c}_s}-V\left(t\right) \frac{\partial T_s\left(z, t\right)}{\partial z}-\frac{2h}{R{\rho }_s{c}_s}[T_s\left(z, t\right)-T_g\left(r, z, t\right)]$$

(11)

$$\frac{\partial T_g\left(r,z, t\right)}{\partial t}=\frac{2h}{R{\rho }_s{c}_s}\left[T_s\left(z, t\right)-T_g\left(r, z, t\right)\right]$$

(12)

$$\frac{{\rho }_g{c}_g}{{\mathsf{\lambda }}_{\mathrm{g}}}\frac{\partial T_g\left(r, z,t\right)}{\partial t}=\frac{1}{r}\frac{\partial T_g\left(r,z,t\right)}{\partial r}+\frac{{\partial }^2{T}_g\left(r,z,t\right)}{\partial r^2}$$

(13)

where the subscripts “s” and “g” indicate fuel salt and graphite, respectively; $P_f$ is the power contained in the fuel salt; $R$ is the radius of the graphite cylinder; $\rho$ is the fuel salt density; $c$ is the specific heat capacity; $h$ is the heat convection coefficient between fuel salt and graphite; $\mathsf{\lambda }$ is the thermal conductivity of graphite. Another MSR 1D N-TH coupling dynamic calculation code reported in the open literature is DYN1D-MSR, developed by Křepel et al. in 2005 [60]. The main improvements of this code compared with Cinsf1D are in terms of thermal-hydraulic calculation. In DYN1D-MSR, the graphite moderator in the core is divided into several segments to more accurately model the temperature distribution of the graphite, while the mass and energy balance of the fuel salt are taken into account. Compared with the point dynamic model, the 1D dynamic model is able to capture the axial shape change of neutron density. However, significant simplification in both neutronics and thermal-hydraulic calculation would still introduce many errors and uncertainties into the simulation, leading to very limited efforts made on the development of this 1D N-TH coupling dynamic model.

3.1.3. 2-Dimension Neutronics Coupling with 2-Dimension Thermal Hydraulics

Coupling the 2D neutronic kinetics model with the 2D thermal-hydraulics model can accurately simulate the nuclear reactor transient events that occur in those types of MSRs that have a symmetric geometry in the X-Y plane, e.g., MSFR and MOSART. The differential equations of neutronic diffusion involved in the model have similar forms, with the scenario of 1D expressed in Equations (8)–(10), except for the terms associated with the neutron diffusion coefficient ($D$) and the flow velocity of fuel salt ($V$), which should be extended to the R-Z dimensions. Relatively large differences exist in the thermal-hydraulic model, compared with the scenario of the 1D model. In the 2D TH models, the mass, momentum, and energy balance of the fuel salt should be taken into account, as described in the following equations:

$$\frac{\partial {\rho }_s\left(r, z, t\right)}{\partial t}+\nabla ·\left[{\rho }_s\left(r, z, t\right)V\left(r,z,t\right)\right]=0$$

(14)

$$\frac{\partial [{\rho }_s\left(r, z, t\right)V\left(r,z,t\right)]}{\partial t}+\nabla ·\left[{\rho }_s\left(r, z, t\right)V^2\left(r,z,t\right)\right]=-\nabla P\left(r,z,t\right)+\nabla ·\tau \left(r,z,t\right)$$

(15)

$$\frac{\partial [{\rho }_s\left(r, z, t\right)h\left(r,z,t\right)]}{\partial t}+\nabla ·[{\rho }_s\left(r, z, t\right)h\left(r,z,t\right)V\left(r,z,t\right)]=Q$$

(16)

where $\tau$ is the stress tensor caused by friction and gravity, $h$ is the energy enthalpy of fuel salt, P is the pressure drop of fuel salt, and $Q$ is the energy produced in the fuel salt. In addition, the vorticity of fuel salt should be considered for the TH calculation in fast MSRs, while the heat transfer between fuel salt and the moderator should be taken into account for thermal MSRs.

To realize the coupling between neutronics and thermal hydraulics in the R-Z dimensions, different methods/schemes were proposed. In the 2D dynamic calculation model developed by Yamamoto et al. in 2006, the time-spatial dependent equations involved in the neutronics and thermal-hydraulics calculations were firstly discretized, based on the control volume and implicit methods. These discretized equations were then solved using Patankar’s SIMPLE algorithm and relaxation algorithm [65]. In the same year, Wang et al. made an improvement on a traditional 2D dynamic code, SIMMER-III, which was originally developed for the safety assessment of sodium-cooled fast reactors, to account for the effects of fuel flow. This code couples N and TH externally, for which the transient neutron flux distribution is calculated with the improved quasi-static method, and temperature distribution is evaluated using the Eulerian fluid-dynamics model [62]. In 2008, Nicolino et al. realized 2D N-TH coupling via solving the equations of neutronic diffusion and full-core computational fluid dynamics together, using a rigorous and fully implicit approach Jacobian Free Newton Krylov (JFNK) algorithm [63]. In 2009, Zhang et al. developed a 2D dynamic simulation code by using the alternating direction implicit TDMA algorithm to solve the equations of N and TH in an internal coupling approach, but the TH calculation only took into account the heat transfer of fuel salt [61]. In 2011, Cammi et al. applied a segregated solver of the finite element method, using quadratic Lagrangian elements to solve the overall equation systems of the reactor, based on the multi-physics coupling platform, COMSOL. This multi-physics method can take into account the turbulence and buoyancy effects, provide detailed information on the fields of temperature and velocity and is suitable for the simulation of MSRs with complex geometry [64]. Overall, 2D N-TH coupling offers an improvement in terms of accuracy in safety analysis but being unable to accurately simulate the local transients/accidents, or transients/accidents occurring in MSRs with an asymmetric core geometry, limits its development.

3.1.4. 3-Dimension Neutronics Coupling with 2-Dimension/3-Dimension Thermal Hydraulics

As computing technologies constantly improve, the coupling of 3D neutronics with thermal hydraulics has attracted growing attention from MSR dynamics analysis. In 2007, Krepel et al. developed a 3D MSR dynamic code, termed DYN3D-MSR, by improving the light water reactor dynamic code, DYN3D. The 3D neutron diffusion equations were split into three 1D equations and solved using the nodal expansion method. A multi-channel thermal-hydraulic model was adopted to externally couple with the 3D neutronic model. With this developed code, channel-wise transients were successfully simulated. However, the heat insulation treatment on the edge of fuel channels in the multi-channel model made DYN3D fail to simulate the accidents that involve heat transfer between the adjacent channels, e.g., an accident due to a blockage in the fuel channel [66]. In 2009, Kophazi et al. developed a similar 3D dynamic code by coupling the 3D neutronic code, DALTON, and the in-house thermal-hydraulic code, THERM. This code is capable of modeling heat conduction in the moderators of graphite over the entire core, using the conjugate gradient method to solve the related equations. Every fuel channel in the core can, hence, be connected thermally through the graphite [69]. The multi-channel model was also adopted by Zhang et al. in 2015 [67] and by Cui et al. in 2022 to couple their respective in-house 3D neutronic diffusion solvers externally, forming the 3D dynamic analysis codes of MOREL and TMSR-3D, respectively. In these two codes, however, heat transfer between two adjacent graphite cells was not considered, as in DYN3D-MSR [68].

In the work of Nagy et al. in 2014, an in-house thermal-hydraulic model using the porous medium method was applied. In this TH model, the equations of fuel salt by Navier-Stokes, fuel salt temperature, and graphite temperature are solved to evaluate the change in temperature distribution and flow field. This couples with DALTON to perform dynamic analysis of graphite-moderated MSRs [70]. Based on the same neutronic code, DALTON, in 2012, Linden made a couple with an in-house CFD program, HEAT, in which the RANS (Reynolds-averaged Navier-Stokes) equation was solved with k − ε model, to perform dynamic analysis of MSFR [23]. In 2020, Tiberga developed a high-fidelity multi-physics simulation tool by coupling the in-house CFD parallel solver DGFlows and the neutron transport code PHANTON-S_N to accurately evaluate the behavior of the physical mechanisms of neutron transport, fluid flow, and heat transfer in a fast-spectrum MSR under transients/accidents [24].

In addition, efforts have been also made to realize the N-TH coupling on multi-physics platforms, to take advantage of advanced computing techniques. Based on the multi-physics coupling platform OpenFOAM, in the work of Cammi et al. in 2014, the differential equations of 3D neutron diffusion and the RANS model for mass and momentum conservation were solved together, via Picard iterations [21]. MOOSE (multiphysics object-oriented simulation environment) is a parallel computational framework oriented for solving the nonlinear equations of coupled systems by using the JFNK method and allows researchers to easily implement advanced finite element methods and to use niche basis functions. In 2017, Ridley developed a high-fidelity dynamic analysis model for MSRE by using one application of MOOSE, Moltres, to solve the governing equations of neutron diffusion and thermal hydraulics [25]. Later, in 2022, Yang et al. made a coupling of the neutron transport code PROTEUS-NODAL with the system thermal-hydraulic calculation code SAM, based on MOOSE [72]. COMSOL is a finite element solver and multiphysics simulation software that is oriented for the solution of partial differential equations of coupled systems. With COMSOL, in 2020, Bajpai et al. carried out a coupling of neutronics with thermal hydraulics by using a two-phase Euler-Euler model to evaluate the influence of helium bubbling [71].

3.1.5. N-TH Coupling for Steady-State Analysis

The steady state of the core is one important state necessitating deep analysis for core design. It is generally regarded as a condition that the performances of the reactor system are unchanged with time. The terms associated with differential time in the partial differential equations of neutron diffusion/transport and thermal hydraulics can be set as zero. The solution processes for the coupling calculation are, hence, significantly simplified. Most importantly, the coupling of N-TH in the dimension of time is removed, which substantially reduces the calculation cost and allows for applying more time-cost codes of neutronics and thermal hydraulics, such as Monte Carlo codes and CFD, to perform accurate core modeling and calculation. As early as 1962, Engel et al. made an N-TH coupling in 2D for MSRE steady-state analysis by using an analytical solution [73]. In 2008, Zhang et al. conducted a coupling of the deterministic neutron code, Dragon, with a multi-channel thermal-hydraulic model to perform 3D N-TH calculation for thermal MSRs [74]. Based on the Monte Carlo code MCNP, in 2013, Guo et al. [75] and Laureau [76] carried out a coupling with a multi-channel thermal-hydraulic model and CFD in 3D for a steady-state analysis of thermal MSRs and MSFR, respectively. In 2020, Deng et al. conducted a steady state analysis of MSFR by coupling the open-source codes of OpenMC (3D neutronics calculation) and OpenFOAM (3D thermal-hydraulic calculation) to accurately analyze the distribution of temperature, power, and flow field in the core [77].

3.2. N-TH-BN Coupling Development

Fuel burnup is a measure of how much energy is extracted from fuel, reflecting the fuel depletion and the consequent changes in fuel composition. One main task of fuel burnup analysis is to accurately evaluate the time evolution of all nuclides of fuel in the core, which is also called burnup calculation. Up until now, the calculation of burnup coupling with neutronics has been extensively studied by considering the above MSR unique mechanisms, while thermal hydraulics is rarely coupled, which necessitates more efforts in terms of future N-TH-BN study. The main features of current N-TH-BN coupling models are summarized in Table 2.

3.2.1. The N-TH-BN Coupling Model, Considering Online Refueling and Reprocessing

In a reactor core, the time evolution of the isotopes of fuel under a neutron flux can be generally described using the Bateman equations:

$$\frac{d{N}_i}{dt}={\displaystyle \sum }_j{N}_j\mathsf{\Phi }\mathsf{\Delta }{\sigma }_{j\to i}-{\displaystyle \sum }_j{N}_i\mathsf{\Phi }\mathsf{\Delta }{\sigma }_{i\to j}+{\displaystyle \sum }_j{N}_j{\mathsf{\lambda }}_j\mathsf{\Delta }{\sigma }_{j\to i}-N_i{\mathsf{\lambda }}_i$$

(17)

where N is the nuclide concentration, σ is the neutron reaction cross-section associated with temperature, and λ is the nuclide decay constant. In the equations, nuclide transmutation under neutron flux $\mathsf{\Phi }$ (the first and second term on the right side) and nuclide decay (the third and fourth term on the right side) are considered. To solve Equation (17), neutron flux $\mathsf{\Phi }$ and the temperature of the fuel should be provided in advance; in turn, the change of nuclide concentration would influence the neutron flux $\mathsf{\Phi }$ and temperature distribution, and so forth. Therefore, a coupling calculation among them is required [78,79,80]. For MSRs, the nuclides in the core that undergo online refueling and reprocessing would deviate somewhat from the law described in Equation (17). As a result, the traditional burnup calculation tools for reactors with solid fuel are no longer suitable for MSRs.

In the 1960s, to evaluate the fuel cycle performance of the circulating fuel reactors, ORNL developed a code, ROD, which couples the neutronics and burnup calculation while taking no account of thermal hydraulics. The time change of nuclide concentration in MSRs is described by two equations:

$$V_j\frac{\mathrm{d}{N}_{ij}}{\mathrm{d}t}=Q_{ij}^f+R_{ij}+F_{ij}+T_{ij}+D_{ij}-N_{ij}\left(t_{ij}+d_{ij}+q_{ij}-r_{ij}\right)=0$$

(18)

$$\left[N_{ij}^s-N_{ij}^{s-1}\right]{\nu }_i{\displaystyle \sum }_k{C}_{ijk}^f={\displaystyle \sum }_{ijk}{N}_{ij}^{s-1}{C}_{ijk}^a-{\displaystyle \sum }_{ijk}{N}_{ij}^{s-1}{C}_{ijk}^f{\nu }_i$$

(19)

where $V_j$ is the volume of stream j; $Q_{ij}^f$ is the feed rate of nuclide $i$ into stream j; $R_{ij}$ is the production rate of nuclide i in stream j, due to recycling from other streams; $F_{ij}$ is the production rate of fission fragment i in stream j; $T_{ij}$ and $D_{ij}$ indicate the production rate of nuclide i in stream j due to neutron absorption and the radioactive decay of other nuclides, respectively; $t_{ij}, d_{ij}$and $q_{ij}$ denote the rate coefficient for the loss of nuclide i in stream j because of neutron capture, radioactive decay, and processing removal, respectively; $r_{ij}$ is the production rate coefficient of nuclide i in stream j, resulting from recycling from stream j; $N_{ij}^s$ is the concentration of nuclide i in stream j at iteration s; ${\nu }_i$ represents the neutrons produced per fission in nuclide i; $C_{ijk}^f$ and $C_{ijk}^a$ indicate the fission and absorption reaction rate coefficient of nuclide i in stream j in region k, respectively. Equation (22) is an improved Bateman equation, considering online feeding and reprocessing, while Equation (23) is the conservation equation to overcome the deficiency/excess of neutron production at iteration s-1 by adding/removing fissile materials in iteration s. By coupling the 2D neutron diffusion calculation performed by the MODRIC module from ROD, the above two equations were iteratively solved by module ERC of ROD to realize the equilibrium nuclide concentration calculation [81].

In 2005, Nuttin et al. introduced a fictive time-decay constant, ${\lambda }_e,$ to model the fuel salt reprocessing in nuclide evolution calculation, by assuming that the amount of nuclide i extracted from inventory N_i during dt is proportional to the quantity N_i $\times \frac{\mathrm{dt}}{{\mathrm{T}}_{\mathrm{r}}}$, with extraction efficiency ${\epsilon }_e$:

$$d{N}_i=N_i\frac{dt}{T_r}{\epsilon }_e=N_i{\lambda }_{e}dt, {\lambda }_e=\frac{{\epsilon }_e}{T_r}$$

(20)

where $T_r$ is the time during which the fuel salt in the primary loop is reprocessed once. The nuclide evolution in MSRs can, then, be written as:

$$\frac{\mathrm{d}{N}_i}{\mathrm{d}t}={\displaystyle \sum }_j{N}_j\mathsf{\Phi }∆{\sigma }_{j\to i}-{\displaystyle \sum }_j{N}_i\mathsf{\Phi }∆{\sigma }_{i\to j}+{\displaystyle \sum }_j{N}_j{\lambda }_j∆{\sigma }_{j\to i}-N_i{\lambda }_i-N_i{\lambda }_e+F_i^{fed}$$

(21)

where $F_i^{fed}$ refers to nuclide refueling. This equation is solved based on the neutron flux provided by MCNP to perform the nuclide evolution calculation [7].

The above MSR burnup calculation method was widely adopted in the subsequent studies of MSR fuel cycle modeling. In 2008, to evaluate the fuel element compositions contained in the molten salt in a state of equilibrium, Becker et al. developed a code, MOCUP, by coupling MCNP5 and ORIGEN2, with the former used for neutronic calculation while the latter was used for nuclide evolution calculation [82]. Based on the same neutronic code, MCNP5, in 2011, Merle-Lucotte et al. made a coupling with an in-house nuclide evolution code REM for the burnup calculation of MSFR [83]. This coupled model was improved by Doligez et al. in 2014, to assess nuclide evolution in the reprocessing unit and fission gas storage beside it in the core [84]. SCALE and SERPENT are two other Monte Carlo neutronic calculation codes that are usually taken for N-BN coupling calculation. In 2012, Sheu et al. developed a special sequence based on SCALE6 for the N-BN calculation, which, however, can only account for batch reprocessing and the refueling of MSRs [85]. In the special MSR reprocessing sequence (MSR-RS) developed by Zou et al. in 2015, based on SCALE6, the function was extended to cover the modeling of online reprocessing [86]. In 2013, Aufiero et al. made an improvement to the SERPENT code to perform N-BN coupling calculation, considering reprocessing and refueling in MSRs [87]. Apart from the Monte Carlo method, deterministic neutronic methods were also reported as being used. In 2013, Fiorina et al. realized the N-BN coupling calculation for MSRs by improving a deterministic code, ERANOS [88].

In the above studies, the nuclide evolution calculation took no account of the thermal hydraulics coupling and might be significantly mis-estimated. To address these issues, recently, some efforts have been made regarding the N-BN-TH coupling calculation. In 2013, to accurately evaluate the nuclide evolution in an MSFR, Frima added a new term related to flow velocity, $\nabla ·\left[u_i\left(r,t\right)N_i\left(r,t\right)\right]$, to the space-time-dependent nuclide evolution equations, which is an extension of Equation (21), to account for flow effect. The in-house code, LOWFAT, was taken to perform the multi-group neutron diffusion calculation, while the thermal hydraulic calculation was carried out with the in-house code HEAT, via solving the Reynolds averaged Navier-Stokes (RANS) equations to model the turbulent flow [89]. In 2019, Graham et al. developed a multi-physics model in which the neutronics, thermal hydraulics, and burnup are coupled via the flow velocity and temperature, based on the virtual environment for reactor applications (VERA) [90]. The flowing effect and FPs deposited on the surface of structure materials were considered in the nuclide evolution equations (highlighted in bold in the following equation):

$$\begin{array}{c}\frac{\mathrm{d}{N}_i\left(r,t\right)}{\mathrm{d}t}={\displaystyle {\displaystyle \sum }_j}{N}_j\left(r,t\right)\Phi \left(r,E\right)∆{\sigma }_{j\to i}\left(T,E\right)-{\displaystyle {\displaystyle \sum }_j}{N}_j\left(r,t\right)\mathsf{\Phi }\left(r,E\right)∆{\sigma }_{i\to j}\left(T,E\right)\\ +{\displaystyle {\displaystyle \sum }_j}{N}_j\left(r,t\right){\lambda }_j∆{\sigma }_{j\to i}\left(T,E\right)-N_i\left(r,t\right){\lambda }_i-N_i{\lambda }_e+F_i^{fed}-\mathbf{\nabla }\\ ·\left[{\mathit{u}}_{\mathit{i}}\left(\mathit{r},\mathit{t}\right){\mathit{N}}_{\mathit{i}}\left(\mathit{r},\mathit{t}\right)\right]-\mathbf{\nabla }·{\mathit{J}}_{\mathit{i}}\left(\mathit{r},\mathit{T}\right)\end{array}$$

(22)

The neutronics calculation was performed by a module, MPACT, using a method of characteristics (MOC), while the thermal-hydraulic calculation was conducted via a subchannel module CTF. In 2022, Ronco et al. conducted an in-depth study on the modeling of FP deposit on the surface of the structure material by coupling neutronics and thermal hydraulics, considering FPs turbulent diffusion $D_{eff}$ [27,91]:

$$\frac{d{J}_i}{dt}+\nabla ·\left(u_i{J}_i\right)=\nabla ·\left(D_{eff}\nabla J_i\right)-{\lambda }_i{J}_i+y_i{\displaystyle \sum }_g{\mathsf{\Sigma }}_{f,g}{\mathsf{\Phi }}_g$$

(23)

where $y_i$ is the fission yield, ${\mathsf{\Sigma }}_{f,g}$ is the fission cross-section, and g is the energy group.

In general, the time step for the above nuclide evolution calculation is in a magnitude of days, while the time scale of transient analysis is in a magnitude of a second. To realize transient analysis covering the whole fuel life cycle, in 2022 Pathirana et al. proposed an MSR dynamics model which is dependent on depletion. The simulation time was factorized into transient time and depletion time. The depletion-dependent data were pre-calculated before the transient evaluation. The transient calculation was then executed at a certain time, using the constant depletion-dependent parameters [28].

3.2.2. N-TH-BN Coupling Model with Considering Helium Bubbling

As introduced in Section 2.2, helium bubbling functions as a way to remove the sparingly soluble FPs, including fission gases and noble metal elements, and subsequently influenced the evolution behavior of these nuclides. ¹³⁵Xe is one type of fission gas nuclides having a large neutron absorption cross-section and is regarded as a major neutron poison in the core. Modeling of its evolution behavior has been extensively studied by gradually coupling the neutronics, thermal hydraulics, and helium bubbling. With a given temperature and neutron flux, a simple model considering several mechanisms, including the ¹³⁵Xe interaction between fuel salt and graphite, ¹³⁵Xe decay in the external loop, and ¹³⁵Xe stripping in the pump bowl, was developed by ORNL in 1962 [92,93]. In 1967, this model was improved by including all of the important and potentially important behavior mechanisms, encompassing the ¹³⁵Xe generation rate, the decay rate in the salt, the burnup rate in the salt, the stripping rate, the migration rate to the graphite, and the migration rate to the circulating bubbles. The core was divided into 72 annular rings to accurately assess the ¹³⁵Xe behavior in the graphite. Many of the gas-transport constants involved in the governing equations of ¹³⁵Xe behavior were inferred from a pre-operational experiment run in the MSRE, with ⁸⁵Kr as the tracer [94]. This calculation model can reproduce reasonably well the observed core poisoning of ¹³⁵Xe at a steady state, but it could not adequately describe its transient behavior. In front of this background, a more elaborate mathematical model that specifically described the behavior of ¹³⁵Xe in the primary loop and considered the cover-gas solubility was developed by ORNL in 1971. As demonstrated in Figure 7, the behavior of ¹³⁵Xe is realistically simulated by dividing the primary loop into four calculation modules, along the flow of fuel salt. In each calculation module, the mechanisms of ¹³⁵Xe generation, disappearance, and interaction among gas bubbles, liquid fuel salt, and graphite (only for the core calculation module) were taken into account. A total of 31 time-related differential equations were applied to accurately model the transient behavior of ¹³⁵Xe in four calculation modules, as well as the connections among them [95]. However, in this detailed mathematical model, the values required for the parameters of xenon mass transport were significantly different from the predicted values; it failed to describe the transient observations in MSRE adequately. Hence, additional investigations would be needed to further elucidate the behavior of ¹³⁵Xe in MSRs. In 2020, to accurately simulate the internal behavior of ¹³⁵Xe in the graphite, Price et al. made an improvement to the model by adding a cylindrical diffusion equation for ¹³⁵Xe [96].

In the above ¹³⁵Xe evolution models, the neutron flux and temperature distribution are given in advance without coupling with N-TH, which might misestimate the positioning of ¹³⁵Xe in the core. In 2017, Wu et al. analytically derived the time-dependent equations of the nuclide evolution of ¹³⁵Xe contained in the flowing fuel salt by considering the stripping rate of helium bubbling. The neutron flux, calculated by a deterministic code (CITATION), was coupled to solve the nuclide evolution equations and to evaluate the behavior of ¹³⁵Xe regarding steady-state and transients [13,97]. In 2022, Caruggi et al. extended their previous OpenFOAM-based N-TH coupling model to evaluate the behavior of gaseous fission products in an MSFR [26]. The Euler-Euler solver for this two-phase calculation was employed to solve the governing equations that describe the mass exchange between bubbles and fuel salt. Besides this, the mechanisms of fission gas production resulting from fission, consumption due to neutron absorption, and removal by the off-gas system were modeled by coupling with N-TH:

$$S_{Xe}-{\alpha }_m\left[{\lambda }_{dec}+{\displaystyle \sum }_n\left({\sigma }_{c,n}{\mathsf{\Phi }}_n\right)\right]C_{Xe, m}+\frac{d{M}_{Xe, m}}{dt}=\frac{\partial {\alpha }_m{C}_{Xe, m}}{\partial t}+\nabla ·\left({\alpha }_m{u}_m{C}_{Xe, m}\right)-\nabla ·\left[\frac{\partial {\alpha }_m{u}_m}{S{c}_m}\nabla \left(\frac{C_{Xe, m}}{{\rho }_m}\right)\right]$$

(24)

$$-{\alpha }_g{\lambda }_{dec}{C}_{Xe, g}+\frac{d{M}_{Xe, g}}{dt}=\frac{\partial {\alpha }_g{C}_{Xe, g}}{\partial t}+\nabla ·\left({\alpha }_g{u}_g{C}_{Xe, g}\right)-\nabla ·\left[\frac{\partial {\alpha }_g{u}_g}{S{c}_g}\nabla \left(\frac{C_{Xe, g}}{{\rho }_g}\right)\right]$$

(25)

$$\frac{d{M}_{Xe, m}}{dt}=-\frac{d{M}_{Xe, g}}{dt}=\frac{6{\alpha }_g\mathrm{Sh}}{d_b^2}{D}_m\left(H{C}_{Xe, g}-C_{Xe, m}\right)$$

(26)

$$S_{Xe}=y_{Xe}\frac{m_{mol}}{N_{Av}}{\displaystyle \sum }_n({\mathsf{\Sigma }}_{f,n}{\mathsf{\Phi }}_n)$$

(27)

$$\frac{\partial C_{Xe}}{\partial t}=y_{Xe}{\displaystyle \sum }_n({\mathsf{\Sigma }}_{f,n}{\mathsf{\Phi }}_n)-C_{Xe}{\displaystyle \sum }_n({\mathsf{\Sigma }}_{a,n}{\mathsf{\Phi }}_n)-\left({\lambda }_{dec}+{\lambda }_{bub}\right)C_{Xe}$$

(28)

where the subscript m and g indicate the molten salt and gas, respectively; $\alpha$ is the interfacial area of exchange per unit volume; $C_{Xe}$ is the concentration of ¹³⁵Xe; $\frac{d{M}_{Xe}}{dt}$ denotes the mass transfer of ¹³⁵Xe to or from the molten salt/gas; $Sc$ is the dimensionless Schmidt number; $Sh$ is the Sherwood number; $D$ is the diffusivity coefficient; $S_{Xe}$ is the production of ¹³⁵Xe, resulting from fission; $m_{mol}$ is the molar mass of ¹³⁵Xe; $N_{Av}$ is the Avogadro number; ${\lambda }_{bub}$ is the fictive time decay constant resulting from helium bubbling. Equations (24) and (25) describe the behavior of ¹³⁵Xe in the molten salt fuel and helium bubbles, respectively. Equation (26) models the interaction of ¹³⁵Xe between two phases. Equation (28) evaluates the time evolution of ¹³⁵Xe, considering helium bubbling. In the same year, Taylor et al. carried out a similar study based on the multi-physic platform VERA by using the liquid-gas model. Neutron flux, temperature distribution, and fluid-flowing field were provided by the N-TH calculation for simulating fission gas behavior. Both the above two models, however, neglect the migration of ¹³⁵Xe in the graphite, which necessitates an in-depth study in the future [98].

The evolution behaviors of noble metal fission products are significantly affected by helium bubbling. In 2021, Walker et al. developed a coupled NM (noble metal)-helium bubble model to simulate the behaviors of NM transport, deposition, and extraction. According to the behaviors of NM elements in the primary loop, the insoluble NM species were split into four categories: species dissolved in the bulk liquid, species deposited on a stationary structural surface of the loop, species transported either on the surface of or in helium bubbles and species accumulated in the pump bowl. For each sub-species, a governing differential equation was set up with respect to time, i.e., advective mass transport, due to a velocity field, source, or sink terms resulting from nuclide decay and deposit, and mass transfer terms among sub-species. These equations for all the species were solved by coupling the thermal-hydraulic calculation (CTF), while using a given neutron flux, which can be improved in the future [99].

Table 2. Development status of the N-TH-BN coupling models.

Codes/Developers/Platforms	N Model	TH Model	Burnup Model
			Nuclide Evolution	Reprocessing	Refueling	Helium Bubbling
ROD, 1962 [81]	2D diffusion	Given temperature	Equilibrium	Removal rate	Feed rate	-
ORNL, 1960s–1970s [92,93,94]	Given neutron flux	Given temperature	135Xe time-step evolution	-	-	135Xe migration among fuel salt, graphite and bubble
Nuttin et al., 2005 [7]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
MOCUP, 2008 [82]	Monte Carlo	Given temperature	Equilibrium	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Sheu et al., 2012 [85]	Monte Carlo	Given temperature	Time-step	100% removal efficiency	Feed rate	100% removal efficiency for fission gas and NM removal
Improved ERANOS, 2013 [88]	3D diffusion	Given temperature	Equilibrium	Fictive time decay constant	Feed rate	-
Frima, 2013 [89]	3D diffusion	CFD	Time-step	Fictive time decay constant	Feed rate	-
Improved SERPENT-2, 2013 [87]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Doligez et al., 2014 [84]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant, accounting for nuclide evolution in reprocessing unit and storage	Feed rate	Fictive time decay constant for fission gas and NM removal
MSR-RS, 2015 [86]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Wu et al., 2017 [13,97]	3D diffusion	Given temperature	135Xe time-step evolution	-	-	135Xe removal rate
VERA, 2019 [90]	3D neutron transport	Sub-channel	Time-step, accounting for interaction of each nuclide with system components	Removal rate	Feed rate	Removal rate for fission gases and NMs
Price et al., 2020 [95,96]	Given neutron flux	Given temperature	135Xe time-step evolution	-	-	135Xe migration among fuel salt, graphite and bubble, accounting for 135Xe diffusion in the graphite
Walker, et al., 2021 [99]	Given neutron flux	Sub-channel	Noble metal time-step evolution	-	-	Two-phase model, accounting for noble metal transport, deposition and extraction
Caruggi et al., 2022 [26]	3D neutron transport	CFD	Fission gas time-step evolution	-	-	Two-phase model, accounting for fission gases interaction between a liquid and gas
Ronco et al., 2022 [27,91]	2D neutron transport	CFD	FPs time-step evolution, accounting for FPs deposition on the wall	-	-	-
VERA, 2022 [98]	3D neutron transport	Sub-channel	Fission gas time-step evolution	-	-	Two-phase model, accounting for fission gases interaction between a liquid and a gas
Pathirana, et al., 2022 [28]	Point kinetics	Lumped parameters	Given nuclide evolution data over the burnup cycle	-	-	-

Table 2. Development status of the N-TH-BN coupling models.

Codes/Developers/Platforms	N Model	TH Model	Burnup Model
			Nuclide Evolution	Reprocessing	Refueling	Helium Bubbling
ROD, 1962 [81]	2D diffusion	Given temperature	Equilibrium	Removal rate	Feed rate	-
ORNL, 1960s–1970s [92,93,94]	Given neutron flux	Given temperature	135Xe time-step evolution	-	-	135Xe migration among fuel salt, graphite and bubble
Nuttin et al., 2005 [7]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
MOCUP, 2008 [82]	Monte Carlo	Given temperature	Equilibrium	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Sheu et al., 2012 [85]	Monte Carlo	Given temperature	Time-step	100% removal efficiency	Feed rate	100% removal efficiency for fission gas and NM removal
Improved ERANOS, 2013 [88]	3D diffusion	Given temperature	Equilibrium	Fictive time decay constant	Feed rate	-
Frima, 2013 [89]	3D diffusion	CFD	Time-step	Fictive time decay constant	Feed rate	-
Improved SERPENT-2, 2013 [87]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Doligez et al., 2014 [84]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant, accounting for nuclide evolution in reprocessing unit and storage	Feed rate	Fictive time decay constant for fission gas and NM removal
MSR-RS, 2015 [86]	Monte Carlo	Given temperature	Time-step	Fictive time decay constant	Feed rate	Fictive time decay constant for fission gas and NM removal
Wu et al., 2017 [13,97]	3D diffusion	Given temperature	135Xe time-step evolution	-	-	135Xe removal rate
VERA, 2019 [90]	3D neutron transport	Sub-channel	Time-step, accounting for interaction of each nuclide with system components	Removal rate	Feed rate	Removal rate for fission gases and NMs
Price et al., 2020 [95,96]	Given neutron flux	Given temperature	135Xe time-step evolution	-	-	135Xe migration among fuel salt, graphite and bubble, accounting for 135Xe diffusion in the graphite
Walker, et al., 2021 [99]	Given neutron flux	Sub-channel	Noble metal time-step evolution	-	-	Two-phase model, accounting for noble metal transport, deposition and extraction
Caruggi et al., 2022 [26]	3D neutron transport	CFD	Fission gas time-step evolution	-	-	Two-phase model, accounting for fission gases interaction between a liquid and gas
Ronco et al., 2022 [27,91]	2D neutron transport	CFD	FPs time-step evolution, accounting for FPs deposition on the wall	-	-	-
VERA, 2022 [98]	3D neutron transport	Sub-channel	Fission gas time-step evolution	-	-	Two-phase model, accounting for fission gases interaction between a liquid and a gas
Pathirana, et al., 2022 [28]	Point kinetics	Lumped parameters	Given nuclide evolution data over the burnup cycle	-	-	-

3.3. N-TH-BN-G Coupling Development

Graphite deformation is another important mechanism that closely couples with neutronics, thermal hydraulics, and burnup in a thermal MSR. Table 3 gives an overview of the development status of the N-TH-BN-G coupling models. In 1969, the volume deformation $\nu$ of graphite was evaluated using a parabolic curve, in terms of core temperature and neutron fluence by Scott and Eatherly [100]:

$$\nu =A\mathsf{\Phi }t+B{\left(\mathsf{\Phi }t\right)}^2$$

(29)

where A and B are the core temperature-dependent constants. Under neutron flux (especially fast neutron with energy > 50 keV) irradiation, the graphite deformation ν firstly decreases to a minimal value ν_m and then increases with time. To avoid introducing a large amount of reactivity during core operation, it is necessary to exclude the fuel salt and xenon from the graphite. As for the first requirements, the pores of the graphite should be no larger than ~1 μ in diameter. Regarding the xenon permeating to the graphite, an allowable permeability is in the order of 10⁻⁸ cm²/s. Because we lack the definitive data, it is usually assumed that the pore size and permeability requirements can be maintained during neutron irradiation until the original graphite volume is reattained. At this point, the time τ is defined as the graphite lifespan, and can be calculated as:

$$0=A\mathsf{\Phi }\tau +B{\left(\mathsf{\Phi }\tau \right)}^2.$$

(30)

Then, the graphite distortion volume ν can be rewritten as:

$$v=4v_m\frac{t}{\tau }\left(1-\frac{t}{\tau }\right).$$

(31)

The parameters$\mathsf{\Phi }$τ and ν_m were estimated to be $\left(9.36-8.39\times {10}^{-3}T\right)\times {10}^{22}$ n/cm² and $\left(-12.0+8.92\times {10}^{-3}T\right)\%$, respectively. The above equations were solved by coupling the core neutron flux and temperature distribution to evaluate the lifetime of graphite in the MSBR core. In this model, the mechanism of burnup is neglected, and the graphite temperature is simply calculated based on the heat transfer from the graphite to fuel salt, using an axial sine distribution of power. Meanwhile, a constant neutron flux (~5 $\times {10}^{14}$ n$·$cm⁻² sec⁻¹, with a neutron energy $\ge$ 50 keV) is taken by assuming the peak power density to be 63 W/cm³ [100].

In the study by Kasten et al. in 1969 [101], the mechanisms of graphite deformation were explored from the viewpoint of graphite stress under neutron flux irradiation at a given core temperature distribution. The corresponding stress-strain equations of graphite were given:

$$\begin{gathered} {ϵ}_i=\frac{1}{E}\left[{\sigma }_i-\mu \left({\sigma }_j+{\sigma }_k\right)\right]+\frac{1}{E}\left[{\sigma }_i-\frac{1}{2}\left({\sigma }_j+{\sigma }_k\right)\right]+\underset{0}{\overset{t}{{\displaystyle \int }}}k\Phi \left[{\sigma }_i-\frac{1}{2}\left({\sigma }_j+{\sigma }_k\right)\right]dt+\underset{0}{\overset{t}{{\displaystyle \int }}}gdt \\ +\alpha \left(T-T_0\right) \end{gathered}$$

(32)

where ${ϵ}_i$ is the total strain in the ith direction (I, j, k = r, θ, z), ${\sigma }_i$ indicates the stress in the ith direction, E is Young’s modulus, $\mu$ is Poisson’s ratio, $k$ is the constant of secondary creep induced by irradiation, $g$ is the time rate of radiation-induced dimensional changes, $\alpha$ is the differential dimensional change due to thermal expansion, and $T_0$ is the reference temperature. The right side of Equation (32) sums up the elastic strain, the saturated primary creep strain, the secondary creep strain, the imposed radiation-induced distortion, and the thermal strain. It was solved using an iterative process, where a zero-order approximation is first generated and then iterated by the subsequent first-order approximation, via updating the third term in the right side of Equation (32). The neutron flux and the temperature of the center channel in a thermal MSR were taken for the N-TH-G coupling calculation.

In 2020, Zhu et al. made an effort to couple the mechanisms of neutronics, thermal hydraulics, burnup, and graphite deformation to accurately evaluate the graphite lifespan in a small modular MSR. The neutronics calculation was carried out based on MCNP to provide a neutron flux. The nuclide evolution calculation, considering online refueling and reprocessing, was calculated using ORIGEN, and the thermal-hydraulic calculation was performed using a multi-channel model. The graphite dimensional change was then evaluated using a parabolic curve, as presented in Equation (29), based on the calculated neutron flux and temperature distribution. In turn, the dimensional change of graphite was fed back to the input of the geometry for the MCNP calculation but was not taken into account by the thermal hydraulics calculation. On the other hand, the temperature distribution of the core calculated by thermal hydraulics was also not fed back to the MCNP [102]. Later, in 2021, an attempt was made by Wang et al. to improve this model by employing Fluent (a CFD software) to provide more accurate temperature distribution. However, the improved model cannot perform burnup calculations and feed back the dimensional changes of graphite to the MCNP and CFD, mainly because of the technical challenges in coupling the calculation and a significantly high calculation cost [103].

In 2020, Stewart proposed an N-TH-G coupling model, termed GeN-foam-G, based on the multi-physics platform GeN-foam. The dimensional change of graphite under neutron flux irradiation was evaluated from the viewpoint of graphite mechanics. The model incorporates two graphite irradiated strain components. One of them accounts for the assessment of graphite dimensional change strain, while the other one is responsible for the creep strain analysis. These two irradiated strain components are performed, based on the neutron flux and temperature distribution provided by GeN-foam, in which the thermal-hydraulic calculation is implemented using the κ-ε turbulence porous media approach, while neutron calculation is conducted by adopting a neutronics subsolver with multigroup neutron diffusion. Feedback on the graphite dimensional changes to N-TH calculation, however, was not available, which is expected to be addressed in the future by the authors [29].

Table 3. Development status of the N-TH-BN-G coupling models.

Codes/Developers/Platforms	N Model	TH Model	BN Model	G Model
Scott and Eatherly, 1969 [100]	Given neutron flux	Lumped parameters	-	Parabolic curve
Kasten et al., 1969 [101]	Given neutron flux	Given temperature	-	Stress-strain
Zhu et al., 2020 [102]	Monte Carlo	Multi-channel	Time-step	Parabolic curve
Wang et al., 2021 [103]	Monte Carlo	CFD	-	Parabolic curve
Stewart, 2020 [28]	3D neutron diffusion	Porous media	-	Stress-strain

Table 3. Development status of the N-TH-BN-G coupling models.

Codes/Developers/Platforms	N Model	TH Model	BN Model	G Model
Scott and Eatherly, 1969 [100]	Given neutron flux	Lumped parameters	-	Parabolic curve
Kasten et al., 1969 [101]	Given neutron flux	Given temperature	-	Stress-strain
Zhu et al., 2020 [102]	Monte Carlo	Multi-channel	Time-step	Parabolic curve
Wang et al., 2021 [103]	Monte Carlo	CFD	-	Parabolic curve
Stewart, 2020 [28]	3D neutron diffusion	Porous media	-	Stress-strain

4. Improvements in the Future Work of the Multi-Physics Coupling Model Study

Calculation cost has become an important issue that greatly challenges the development of MSR multi-physics coupling models, as more mechanisms are coupled. In addition, orchestrating the solution processes of multiple, domain-specific codes running at different time scales is another issue that should be addressed in the future. Dedicated multi-physics platforms, such as MOOSE, can provide advanced computing techniques on parallel computing platforms and a coupled partial differential equations solution. Developing an MSR multi-physics coupling model on multi-physics platforms is an effective way to address the above two issues, attracting growing research attention. Another potential method is the machine-learning-based data-driven method, which can derive relationships or set of rules from data and solve more nonlinear complex problems [104]. It can act as a fast estimator or surrogate model for simulating some MSR mechanisms that have enough data available for training, such as neutronics and thermal hydraulics, to decrease the calculation cost. Up until now, many studies have been conducted on the application of machine learning in nuclear reactor neutronics, thermal hydraulics, and burnups calculation, which has formed a strong basis for the future improvement of MSR multi-physics coupling models [105].

The completeness of coupling is another issue that should be addressed for MSR multi-physics coupling model development. As introduced above, most of the MSR multi-physics coupling models only take into account part of the physical mechanisms or set a fixed value/distribution for certain physical mechanisms without feedback. Currently, there is a trend to integrate the mechanism of salt chemistry into the multi-physics coupling model, in view of the extraordinarily close connection between reactor physics and fuel salt chemistry in MSRs. In a technical report on ORNL in 2018, a framework coupling with the codes of Thermochimica (thermochemistry calculation), SCALE6 (nuclide evolution and power density calculation), and COBRA-TF (thermal-hydraulic calculation) was proposed [106]. Thermochimica is an open-source implementation of the CALPHAD method. It takes molar elemental combinations as the input and can calculate the lowest energy state, from which the equilibrium chemical combinations can be predicted. Some thermodynamic parameters, including heat capacity, enthalpy, the chemical potential of the system components, and the driving forces of metastable phases can be provided for thermal-hydraulic calculation. In turn, the thermal-hydraulic (COBRA-TF) calculation feeds back the temperature and pressure distribution to the chemistry evaluation (Thermochimica). SCALE6 is used to prepare the data of isotropic evolution and power density for Thermochimical and COBRA-TR, respectively. In 2021, Poschmann et al. realized the coupling between Thermochimical and ORIGEN (one module of SCALE6 was used to perform the isotropic evolution calculation), based on which, the time behavior of radio-toxic gases’ quantities in a molten salt fuel loop was explored [107]. In the same year, under the Nuclear Energy Advanced Modeling and Simulation (NEAMS) program, several NEAMS codes that are under development adopted the Gibbs energy minimization (GEM) solver to accurately evaluate the phase and chemical equilibria of FPs and actinides in molten salt systems, and to simulate the leaching and deposition of FPs by coupling with thermal fluids and electron-transport mechanisms [108]. From the view of completeness and the accuracy of simulation, coupling the mechanisms of fuel salt chemistry with the MSR reactor multi-physics is one aspect that should be seriously considered in the future.

5. Conclusions

Molten salt reactors are one type of liquid-fueled reactor that is able to perform online refueling, reprocessing, and helium bubbling to remove sparingly soluble FPs. The fuel utilization, safety, and economics can be significantly enhanced in an MSR, but the unique features of MSR would introduce a series of new physical mechanisms and phenomena, which causes great challenges in terms of numerical simulation. Compared with solid-fueled reactors, coupling between neutronics and thermal hydraulics is tighter since molten salt acts as both the fuel and the coolant. A delayed neutron precursor would flow out from the core, together with the flowing fuel salt, which decreases the fraction of delayed neutrons in the core and, in turn, influences the core safety. Online refueling, reprocessing, and helium bubbling can be regarded as the external physical mechanisms that change the evolution of nuclides in real-time, leading to a failure in the application of the conventional burnup calculation tools to MSR. The dimensional change of the graphite in the core under neutron irradiation is another important physical mechanism that influences the behavior of neutronics and thermal hydraulics. The above physical mechanisms and phenomena derived from the flowing of liquid fuel, graphite, and online refueling, reprocessing and helium bubbling couple with each other, necessitating the development of an MSR multi-physics coupling model.

In this study, the development status of the MSR N-TH, N-TH-BN, and N-TH-BN-G coupling models are extensively reviewed, to outline the problems and challenges that exist. The MSR N-TH coupling model is considered to be the key model of multi-physics coupling. Its development can be dated back to the 1960s, in which the first molten salt reactor, MSRE, was proposed and built by ORNL. Over around 60 years of development, the N-TH coupling model has been improved from the 0D model that couples point neutron kinetics with the lumped thermal hydraulics at the beginning of the 3D model. There is a trend that realizes N-TH coupling on multi-physics platforms, to take advantage of advanced computing techniques. The MSR N-TH-BN coupling model adds a dimension of burnup that must take into account the many unique mechanisms of MSR, including online refueling, reprocessing, and helium bubbling. Up until now, the calculation of burnup coupling with neutronics has been extensively studied, while thermal hydraulics is rarely coupled, which will necessitate more effort in any future N-TH-BN coupling study. The N-TH-BN-G coupling model additionally considers the mechanism of graphite dimensional change under neutron irradiation. Some efforts have been made for this coupling model, but it is still at the initial development stage, with some mechanisms being unilaterally coupled without feedback.

In the future, new methods and mechanisms can be considered to address the issues of calculation cost, the orchestration of the codes with the solution processes in different time scales, and the completeness of the multi-physics coupling model. Dedicated multi-physics platforms, such as MOOSE and VERA, can provide advanced computing techniques in parallel computing and the solution of coupled partial differential equations. They have recently been attracting growing attention and can be viewed as an effective way to decrease the calculation cost and orchestrate the codes with solution processes on a different time scale. Machine learning is another method that is expected to be used to accelerate the calculation of the MSR multi-physics model as it can solve nonlinear complex problems with a low calculation cost via establishing a fast estimator or surrogate model. Regarding the completeness of the MSR multi-physics coupling model, besides all the potential reactor’s physical mechanisms, molten salt chemistry is another mechanism that should be coupled since MSR is actually a chemical nuclear reactor, in which the reactor, physics and fuel salt chemistry have an extraordinarily close connection.