Core and blanket thermal-hydraulic analysis of a molten salt fast reactor based on coupling of OpenMC and OpenFOAM

Core and blanket thermal-hydraulic analysis of a molten salt fast reactor based on coupling of OpenMC and OpenFOAM 1 2 3 first-author https://orcid.org/0000-0002-2325-5976 1 2 3 https://orcid.org/0000-0002-0748-8980 1 2 3 corresp Corresponding author. E-mail address: chenjg@sinap.ac.cn; caixz@sinap.ac.cn Corresponding author. E-mail address: chenjg@sinap.ac.cn; caixz@sinap.ac.cn chenjg@sinap.ac.cn https://orcid.org/0000-0002-7801-7895 1 2 https://orcid.org/0000-0002-5960-1682 1 2 https://orcid.org/0000-0002-0087-7991 1 2 https://orcid.org/0000-0002-5163-3903 1 2 3 1 2 3 corresp Corresponding author. E-mail address: chenjg@sinap.ac.cn; caixz@sinap.ac.cn Corresponding author. E-mail address: chenjg@sinap.ac.cn; caixz@sinap.ac.cn chenjg@sinap.ac.cn https://orcid.org/0000-0002-7108-2344 Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China CAS Innovative Academy in TMSR Energy System, Chinese Academy of Science, Shanghai 201800, China CAS Innovative Academy in TMSR Energy System, Chinese Academy of Science, Shanghai 201800, China University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China In the core of a molten salt fast reactor (MSFR), heavy metal fuel and fission products can be dissolved in a molten fluoride salt to form a eutectic mixture that acts as both fuel and coolant. Fission energy is released from the fuel salt and transferred to the second loop by fuel salt circulation. Therefore, the MSFR is characterized by strong interaction between the neutronics and the thermal hydraulics. Moreover, recirculation flow occurs, and nuclear heat is accumulated near the fertile blanket, which significantly affects both the flow and the temperature field in the core. In this work, to further optimize the conceptual geometric design of the MSFR, three geometries of the core and fertile blanket are proposed, and the thermal-hydraulic characteristics, including the three-dimensional flow and temperature fields of the fuel and fertile salts, are simulated and analyzed using a coupling scheme between the open source codes OpenMC and OpenFOAM. The numerical results indicate that a flatter core temperature distribution can be obtained and the hot spot and flow stagnation zones that appear in the upper and lower parts of the core center near the reflector can be eliminated by curving both the top and bottom walls of the core. Moreover, eight cooling loops with a total flow rate of 0.0555 m 3 ∙s −1 ensure an acceptable temperature distribution in the fertile blanket. Key words: Molten salt fast reactor Core and blanket thermal-hydraulic analysis Neutronics and thermal hydraulics coupling 1 Introduction 1 The molten salt reactor (MSR) was selected as one of six advanced reactor concepts by the Generation IV International Forum in 2002 because of its considerable advantages of inherent safety, sustainability, excellent neutron economy, and nuclear nonproliferation [ 1 1 ]. Most MSR concepts were designed as thermal systems in the early stage of research and development. In 1954, the first MSR in the world, the Aircraft Reactor Experiment [ 2 2 ], was built at Oak Ridge National Laboratory; it had a power of 2.5 MW. In the 1960s, the Molten Salt Reactor Experiment [ 3 3 ], with a power of 8 MW, achieved criticality, further demonstrating the feasibility of operating a liquid-fueled reactor. After decades of MSR development, renewed interest in MSRs focuses on advanced goals by developing reactor and fuel cycle technologies [ 4 4 – 7 7 4–7 ]. The sixth and seventh EURATOM Framework Programs (the ALISIA and EVOL projects) proposed the concepts of the molten salt actinide recycler and transmuter [ 8 8 ] and the molten salt fast reactor [ 9 9 ]. The former was developed to burn transuranium elements (TRU) and minor actinides, whereas the latter is used mainly for breeding thorium fuel. Compared to traditional graphite-moderated MSRs, the MSFR not only has better 232 Th– 233 U breeding and TRU burning capability, but is also inherently safer, because it can avoid the positive reactivity feedback caused by a graphite moderator [ 10 10 ]. The innovative MSFR was initially designed as a single-cylinder core with a radial fertile blanket and absorber (the benchmark model of the MSFR) [ 9 9 ]. Thermal-hydraulic (TH) analyses of the MSFR benchmark (see Figs.S1 and S2) show that a large recirculation flow of fuel salt is formed near the fertile blanket, causing a hot spot to appear in this region and thus degrading the heat transfer performance of the reactor core [ 11 11 , 12 12 ]. The corrosion rate of structural materials in the core is known to increase rapidly with increasing operating temperature. Thus, the core model should be further optimized to eliminate the hot spot and improve the core lifetime of the MSFR. For the channel-type MSR, the flow rate distribution is critical for optimizing the core temperature field because the location and peak value of local hot spots are determined by the fuel flow in the graphite channels [ 13 13 , 14 14 ]. Furthermore, TH studies of the MSFR indicate that it is essential to obtain a relatively flat flow field. Rubiolo et al. analyzed the core thermal hydraulics of the MSFR under nominal conditions [ 15 15 ], and Laureau et al. performed a steady-state neutronics and TH coupling simulation of an MSFR [ 16 16 ]. Both groups adopted an optimized MSFR model in which curves were added to the radial walls of the core cavity in the benchmark model. Although the recirculation flow reported in these studies is slightly improved, a hot spot still exists on the wall near the fertile blanket. That is, the core geometry could be further optimized. Rouch et al. proposed a preliminary core TH design for an MSFR with the main objective of optimizing the core cavity walls, which can increase the overall flow mixing and reduce the temperature peak in the fuel salt and on the core walls [ 17 17 ]. The heat source used in the computational fluid dynamics (CFD) simulation was obtained from MCNP, whereas two nonidentical geometrical models were adopted in the CFD and MCNP simulation, which may make the coupling result unreliable. Further, temperature differences still appear even in a later study in which a geometrical homothetic transformation [ 17 17 ] was employed to adapt the heat source to the real geometry adopted for TH analysis. Laureau et al. further optimized the core geometry of the MSFR, and the hot spot on the curved wall near the fertile blanket was eliminated [ 18 18 ]. However, a new hot spot at 1080 K appeared at the top of the reactor core owing to the relatively low velocity in this region; the hot spot exceeded the upper temperature tolerance (1073 K) of the structural material [ 19 19 ]. For the safety assessment for MSFR, evaluating the THs behavior of the fertile blanket is also of great importance. The small fission fraction of 233 U and the neutron capture of 232 Th in the blanket heats the fertile salt, and the heat cannot be removed by heat convection between the fuel salt and the interface wall of the core and the fertile blanket [ 9 9 , 19 19 , 20 20 ]. Therefore, an external cooling system for the fertile blanket has to be considered in the reactor design. In this work, a new coupling scheme between the open source codes OpenMC and OpenFOAM is developed to physically describe the strong coupling between the neutronics and thermal hydraulics in the MSFR. The distributions of the core power, fuel temperature, and fuel density are exchanged between the two codes to perform the coupling calculation, with the goal of eliminating the hot spot and flow stagnation zones in the core. In addition, a preliminary TH analysis of the fertile blanket is also performed. As mentioned above, the heat deposited in the fertile blanket cannot be evacuated by heat convection, which indicates that the heat transfer between the fuel salt and fertile blanket is relatively small. In addition, the influence of the core heat on the fertile blanket is analyzed in this work by using a conjugate heat transfer calculation for the core and fertile blanket. The results demonstrate that separate modeling of the fertile blanket is reasonable for MSFR TH analysis. This paper is organized as follows. Section 2 describes the methodology and tools. Section 3 presents the core and blanket TH analyses. The conclusions of this work are given in Sect. 4. 2 Methodology and tools 1 2.1 Neutronics model 2 OpenMC is an open source Monte Carlo particle transport simulation code that can perform fixed source, k -eigenvalue, and subcritical multiplication calculations on three-dimensional solid geometry models [ 21 21 ]. It was originally developed by t he Massachusetts Institute of Technology and is widely used for calculations of the criticality, depletion, multigroup cross section generation, etc. The Direct Accelerated Geometry technology has been integrated into OpenMC (DAG-OpenMC hereafter) to perform neutron transport simulation directly on a computer-aided design model [ 22 22 ]. In this work, this technology contributes to the solution of the problem of curved wall modeling for an MSFR owing to its advantages for precise modeling of an arbitrary complex geometry. Based on the ENDF/B-VII.1 nuclear data, a microscopic cross section library between 500 K and 1600 K with 100K internals was produced by NJOY2016 to support the OpenMC continuous-energy transport calculation; a 20 K interval was used between 800 and 1400 K (the operating temperature range of MSFRs). The cross sections are given exactly for only a few points in the table; the temperature values for other materials can be obtained by statistical interpolation [ 23 23 ]. Suppose that cross sections are available at temperatures  T 1 , T 2 ,..., T n , and a material is assigned a temperature  T , where  T i < T < T i +1 . Statistical interpolation is applied as follows. A uniformly distributed random number ( ξ ) of unit intervals (0–1) is sampled. If ξ < ( T − T i )/( T i +1 − T i ), the cross sections at temperature T i +1  are used. Otherwise, the cross sections at  T i  are used. To obtain the power distribution of the MSFR, the calculation model is divided into numerous regions, and each region is labeled with a unique material number. 2.2 Thermal-hydraulic model 2 The open source code OpenFOAM features automatic matrix construction and solution capabilities for scalar and vector equations [ 24 24 ]. Each fluid flow can be represented by a series of partial differential equations; the solvers for these fluid mechanics equations have been developed and implemented in OpenFOAM. The fuel salt flow in the MSFR core without a moderator is considered to be turbulent [ 25 25 ] because the Reynolds number is estimated to be 2.5 × 10 6 in the core region and 5.3 × 10 5 in the out-of-core region [ 11 11 ]. The fuel and fertile salt flows can be assumed to be incompressible because the effect of fuel density variation in the gravity force term is taken into account by using the Boussinesq approximation. Moreover, the ReNormalization Group (RNG) k–ε turbulence model can provide a better description of the flow on the curved wall than the standard k–ε turbulence model [ 26 26 ]. Therefore, considering that the fuel salt flow is incompressible, the Reynolds-averaged Navier–Stokes equations with the Boussinesq closure hypothesis and the RNG k–ε turbulence model are adopted, and buoyantBoussinesqSimpleFoam [ 24 24 ] is selected as the OpenFOAM solver. The steady-state continuity, momentum, and energy conservation equations for the fuel and fertile salts can be written as <math display="block" id="M1"><mrow><mo>∇</mo><mo>⋅</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> <math display="block" id="M2"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mover accent="true"><mi>u</mi><mo>→</mo></mover><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>p</mi><mi>I</mi><mo>+</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mi>k</mi><mi>I</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>ν</mi><mo>+</mo><msub><mi>ν</mi><mtext>t</mtext></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow></mrow><mtext>T</mtext></msup></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>+</mo><mover accent="true"><mi>g</mi><mo>→</mo></mover><mrow><mo>[</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi><mrow><mo>(</mo><mrow><mi>T</mi><mo>−</mo><msub><mi>T</mi><mrow><mtext>ref</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M3"><mrow><mo>∇</mo><mo>⋅</mo><mrow><mo>(</mo><mrow><mover accent="true"><mi>u</mi><mo>→</mo></mover><mi>T</mi></mrow><mo>)</mo></mrow><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mfrac><mi>ν</mi><mrow><mi>P</mi><mi>r</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>ν</mi><mtext>t</mtext></msub></mrow><mrow><mi>P</mi><msub><mi>r</mi><mtext>t</mtext></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>T</mi></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>S</mi><mi>T</mi></msub><mo>,</mo></mrow></math> where <math id="M4"><mover accent="true"><mtext>u</mtext><mo>→</mo></mover></math> represents the Reynolds mean flow velocity vector, m∙s −1 ; p represents the Reynolds average pressure, which is divided by the reference density, m 2 ∙s −2 ; ν and ν t represent the dynamic viscosity and turbulent dynamic viscosity, respectively, m 2 ∙s −1 ; <math id="M5"><mover accent="true"><mtext>g</mtext><mo>→</mo></mover></math> is the gravitational acceleration, m∙s −2 ; T represents the Reynolds-averaged temperature, K ; β is the thermal expansion coefficient in K −1 at a given reference temperature T ref in K; S T is the heat source term, K∙s −1 ; and I , Pr , and Pr t represent the identity matrix, Prandtl number, and turbulent Prandtl number, respectively. The turbulent kinetic energy k and turbulent dissipation rate <math id="M6"><mi>ε</mi></math> are modeled by the RNG k–ε turbulence model, which is given by <math display="block" id="M7"><mrow><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>⋅</mo><mo>∇</mo><mi>k</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>ν</mi><mo>+</mo><mfrac><mrow><msub><mi>ν</mi><mtext>t</mtext></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>k</mi></mrow><mo>]</mo></mrow><mo>+</mo><mi>G</mi><mo>−</mo><mi>ε</mi><mo>,</mo></mrow></math> <math display="block" id="M8"><mrow><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>⋅</mo><mo>∇</mo><mi>ε</mi><mo>=</mo><mo>∇</mo><mo>⋅</mo><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>ν</mi><mo>+</mo><mfrac><mrow><msub><mi>ν</mi><mtext>t</mtext></msub></mrow><mrow><msub><mi>σ</mi><mi>ε</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>∇</mo><mi>ε</mi></mrow><mo>]</mo></mrow><mo>+</mo><mrow><mo>(</mo><mrow><msub><mi>C</mi><mrow><mi>ε</mi><mn>1</mn></mrow></msub><mo>−</mo><mi>R</mi></mrow><mo>)</mo></mrow><mfrac><mi>ε</mi><mi>k</mi></mfrac><mi>G</mi><mo>−</mo><msub><mi>C</mi><mrow><mi>ε</mi><mn>2</mn></mrow></msub><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mi>k</mi></mfrac><mo>,</mo></mrow></math> with <math id="M9"><mrow><mi>R</mi><mo>=</mo><mfrac><mrow><mi>η</mi><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mi>η</mi><mo>/</mo><msub><mi>η</mi><mn>0</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mn>1</mn><mo>+</mo><mi>α</mi><msup><mi>η</mi><mn>3</mn></msup></mrow></mfrac></mrow></math> ; <math id="M10"><mrow><msub><mi>ν</mi><mtext>t</mtext></msub><mo>=</mo><msub><mi>C</mi><mi>μ</mi></msub><mfrac><mrow><msup><mi>k</mi><mn>2</mn></msup></mrow><mi>ε</mi></mfrac></mrow></math> ; <math id="M11"><mrow><mi>G</mi><mo>=</mo><mn>2</mn><msub><mi>ν</mi><mtext>t</mtext></msub><msup><mi>S</mi><mn>2</mn></msup></mrow></math> ; <math id="M12"><mrow><mi>S</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mo>∇</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∇</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow><mo>)</mo></mrow></mrow><mtext>T</mtext></msup></mrow><mo>)</mo></mrow></mrow></math> ; <math id="M13"><mrow><mi>η</mi><mo>=</mo><msqrt><mn>2</mn></msqrt><mi>S</mi><mi>k</mi><mo>/</mo><mi>ε</mi></mrow></math> ; <math id="M14"><mrow><msub><mi>C</mi><mi>μ</mi></msub><mo>=</mo><mn>0.085</mn></mrow></math> ; <math id="M15"><mrow><msub><mi>C</mi><mrow><mi>ε</mi><mn>1</mn></mrow></msub><mo>=</mo><mn>1.42</mn></mrow></math> ; <math id="M16"><mrow><msub><mi>C</mi><mrow><mi>ε</mi><mn>2</mn></mrow></msub><mo>=</mo><mn>1.68</mn></mrow></math> ; <math id="M17"><mrow><msub><mi>σ</mi><mi>k</mi></msub><mo>=</mo><mn>0.71942</mn></mrow></math> ; <math id="M18"><mrow><msub><mi>σ</mi><mi>ε</mi></msub><mo>=</mo><mn>0.71942</mn></mrow></math> ; <math id="M19"><mrow><msub><mi>η</mi><mn>0</mn></msub><mo>=</mo><mn>4.38</mn></mrow></math> ; <math id="M20"><mrow><mi>α</mi><mo>=</mo><mn>0.012</mn></mrow></math> . 2.3 Coupling scheme 2 The coupling calculation flowchart of DAG-OpenMC and OpenFOAM is presented in Fig. 1 Fig. 1 , illustrating the exchange of core power distribution, temperature, and density distribution of the fuel and fertile salt using a Python script (pyMF). As each region of the neutronics calculation model contains many meshes for the TH calculation by OpenFOAM, the average fuel temperature and density of each region are assigned according to the material number (where each region has a unique material number) used in DAG-OpenMC. Similarly, the OpenFOAM meshes located in the same region are assigned the same power density from the neutronics calculation. The key points of the coupling code pyMF are summarized as follows: 1) The temperature and density of the fuel salt and structural material are first assumed to be uniform in the DAG-OpenMC input file to start the neutronics calculation. 2) The power distribution is obtained by transforming the kappa-fission tally of each region from OpenMC to the power density as follows: <math display="block" id="M21"><mrow><msub><mi>P</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><mi>ν</mi><mi>P</mi></mrow><mrow><mi>Q</mi><msub><mi>V</mi><mi>i</mi></msub><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msub><mi>T</mi><mi>i</mi></msub><mo>,</mo></mrow></math> where p i is the power density of region i , J∙m −3 ∙s −1 ; ν represents the neutron yield per fission, neutron/fission; P is the thermal power of the MSFR, W; Q represents the released energy per fission in the MSFR, eV/fission; V i is the volume of region i , m 3 ; k eff represents the effective multiplication factor obtained from the DAG-OpenMC criticality calculation, neutron/source neutron; and T i is the kappa-fission tally result for region i , eV/source neutron. 3) The heat source term in the energy conservation equation [Eq. (3)] is updated from the power distribution, and OpenFOAM is launched to start the TH calculation. The calculation is terminated when the temperature and velocity fields meet the residual convergence criteria set by users; then the temperature and density distributions can be obtained. 4) The above iteration is performed until either k eff and the power density converge or the maximum number of iterations is reached [ 27 27 – 29 29 27–29 ]. In our previous work [ 30 30 ], to verify the reliability of the coupled code, a benchmark model of the MSFR was established, and the effects of the region number in the neutronics calculation and the initial conditions on the k eff , fuel salt velocity, and temperature distributions were studied. The TH results simulated using pyMF were also compared to the reference result for the MSFR [ 31 31 ], and good agreement was obtained for the velocity and temperature fields (Table S1 and Figs. S1 and S2). It can be concluded that the pyMF code can provide reliable results for the steady-state neutronics and TH coupling calculations for an MSFR. 3 Results and Discussion 1 3.1 Description of the molten salt fast reactor 2 A schematic view of the MSFR system is shown in Fig. 2 Fig. 2 . The MSFR is a 3 GWth liquid-fueled reactor with a core fuel salt volume of 9 m 3 [ 9 9 , 19 19 , 20 20 ]. In the primary loop, a molten fluoride salt system functioning as both fuel and coolant works at approximately 923 and 1023 K at the inlet and outlet of the core, respectively. A Hastelloy [ 19 19 , 20 20 ] is selected as the axial reflector, and a radial fertile blanket is adopted to absorb the neutrons escaping from the core and thus improve the breeding capability of the reactor. The primary loop of the MSFR has 16 fuel salt loops, and each loop is equipped with a pump and a heat exchanger. The fuel salt circulates through the core and the 16 external loops, and the nuclear heat generated by fission is transferred to the secondary loop salt by the heat exchangers. The main parameters of the MSFR reactor core are listed in Table 1 Table 1 [ 12 12 , 20 20 , 31 31 ]. 1 1 Parameter 1 1 Value 1 1 Nominal power (MWth) 1 1 3000 1 1 Fuel salt inlet temperature (K) 1 1 923 1 1 Fuel salt outlet temperature (K) 1 1 1023 1 1 Fuel salt composition (mol%) 1 1 77.5LiF-19.886UF4 -2.614ThF4 1 1 Blanket salt composition (mol%) 1 1 77.5LiF-22.5ThF4 1 1 Salt density (kg∙m −3 ) 1 1 4094 − 0.882[ T (K) −1008] 1 1 Salt specific heat (J∙kg −1 ∙K −1 ) 1 1 −1111 + 2.78 T (K) 1 1 Salt thermal conductivity (W∙m −1 ∙K −1 ) 1 1 0.928 + (8.397∙10 −5 ) T (K) 1 1 Salt kinematic viscosity (m 2 ∙s −1 ) 1 1 5.54∙10 −8 exp[3689/ T (K)] 3.2 Core thermal-hydraulic analysis 2 As illustrated in Fig. 3 Fig. 3 , three core geometries are proposed in this work to analyze the TH behavior of the MSFR. Geometry I [ Fig. 3 (a) Fig. 3 (a) ] is described in the literature [ 17 17 , 18 18 ]. Geometry II [ Fig. 3 (b) Fig. 3 (b) ] is obtained by further curving the upper wall of the core, and Geometry III [ Fig. 3 (c) Fig. 3 (c) ] is obtained by curving both the upper and lower walls of the core for comparison. All the arcs of the blanket and absorber in the three geometries are half-ellipses. In the OpenFOAM simulation, only the core, fertile blanket, and primary loop are modeled, and the solid regions, including the reflector and absorber, are not considered, because the coupling calculation is performed only in the above three regions. To obtain a good balance between the calculation efficiency and accuracy, the maximum bulk cell size is set to 3 cm [ 16 16 , 17 17 ]. The RNG k–ε turbulence model and a wall function are adopted to model the flow behavior near the wall. In fluid dynamics, a dimensionless distance to the wall ( y + ) is used as an evaluation criterion of near-wall mesh quality. Therefore, a suitable y + value is needed to obtain accurate results when the k–ε turbulence model and wall function are employed. To ensure that the y + values of all the volume elements adjacent to the wall are distributed within a range of 11.63–300 [ 32 32 ], the cell size is reduced to approximately 3 mm near the core wall. Considering that the MSFR core configuration has good symmetry, only 1/16 of the core geometry ( Fig. 4 Fig. 4 ) is simulated in the TH calculation to conserve computational resources. The pump and heat exchanger are treated as a momentum source in the momentum conservation equation and as a heat sink in the energy conservation equation. The mesh numbers of these three geometries are 124,540 [ Fig. 4 (a) Fig. 4 (a) ], 132,784 [ Fig. 4 (b) Fig. 4 (b) ], and 138,562 [ Fig. 4 (c) Fig. 4 (c) ], respectively. A continuous-energy neutronics simulation is carried out using the DAG-OpenMC code; the calculation models ( Fig. 5 Fig. 5 ) include the core, fertile blanket, and solid regions. All of the components are simulated during the neutron transport simulation; a reflective boundary condition is set for the symmetry planes, and a vacuum boundary condition is set for the walls. No heat exchange is modeled between the core and the fertile blanket or between the fertile blanket and the solid region, and the temperatures of the fertile blanket and solid region are assigned as the average temperature of the fuel salt in the core TH analysis. To obtain the power distribution for the TH calculation, the core is divided into numerous regions, and each region is labeled with a unique material number during the neutronics calculation. Considering the efficiency and accuracy of the neutronics calculation, the size of these regions should be approximately 20 cm, as demonstrated in our previous work [ 30 30 ]. Only a small fraction of the power is released in the external loops, and the strong coupling between the neutronics and thermal hydraulics is significant in the core region. Therefore, no detailed regional division of the external loops is performed. For the neutron transport simulation, 300 generations (the first 50 iterations are skipped) with 50,000 neutrons per generation are adopted, and the obtained standard deviation of k eff is approximately 0.00023. The fuel velocity fields of the three geometries are shown in Fig. 6 Fig. 6 , which shows that a low flow velocity occurs near the outer wall (close to the fertile blanket), whereas the large recirculation region near the fertile blanket in the reference MSFR model (see Fig.S1) is eliminated, and a relatively flat velocity distribution is obtained for all three geometries. Although a small vortex or recirculation pattern still occurs at the core inlet, no hot spot forms owing to the low power density in the area of interest and the low fuel salt temperature at the core inlet. Therefore, further optimization to eliminate this small vortex and recirculation pattern is not needed. The connecting pipe is much smaller than the core; thus, the flow area is widened or narrowed when the fuel salt is flowing into or out of the core. Furthermore, the connecting pipe is curved. Therefore, a large local pressure drop is obtained, which maximizes the flow velocity in the connecting pipe in the three geometries. The maximum fuel salt velocity obtained for Geometry I [ Fig. 6 (a) Fig. 6 (a) ] is approximately 5.05 m∙s −1 , which is slightly higher than those for Geometry II [ Fig. 6 (b) Fig. 6 (b) ] and Geometry III [ Fig. 6 (c) Fig. 6 (c) ]. The reason is that the connecting pipe in Geometry I is smaller than that in the other two geometries. As shown in Fig. 6 (a) Fig. 6 (a) , a small stagnation zone appears near the upper and lower reflectors in the velocity field of Geometry I, whereas a stagnation zone appears in Geometry II only at the core bottom, as shown in Fig. 6 (b) Fig. 6 (b) . For Geometry III, the stagnation zone is eliminated [ Fig. 6 (c) Fig. 6 (c) ], which indicates that Geometry III allows a more uniform velocity distribution. The fuel salt temperature field in each geometry is presented in Fig. 7 Fig. 7 . The fuel salt temperature increases along the axial direction from the bottom to the top for all three geometries, and a small hot spot appears at the top of the central core in Geometry I [ Fig. 7 (a) Fig. 7 (a) ], which might be caused by the low velocity in this region. The temperature profile for Geometry II [ Fig. 7 (b) Fig. 7 (b) ] does not differ significantly from that for Geometry I, although the peak value of the fuel temperature is slightly lower owing to the difference in the geometry. Therefore, the hot spot at the upper core cannot be effectively eliminated by curving only the upper wall of the core. Unlike these two geometries, Geometry III [ Fig. 7 (c) Fig. 7 (c) ] exhibits a more uniform temperature distribution throughout the core and a relatively small radial gradient of the fuel temperature. Moreover, the peak temperature drops sharply, and the difference between the maximum core temperature and average outlet temperature is only approximately 9 K. The turbulent dynamic viscosity field of each geometry is presented in Fig. 8 Fig. 8 . These results are also compared to the result for the reference MSFR model (see Fig.S3). The location of maximum turbulent dynamic viscosity is shifted from the upper core cavity to the lower core cavity. In addition, a zone of high turbulent dynamic viscosity is formed near the entire arc wall of the fertile blanket. Therefore, there is a strong turbulent transport effect in this region, which prevents the formation of the hot spot on the arc wall near the fertile blanket even though the core flow velocity is low in this region (see Figs. 6 Figs. 6 and 7 7 ). The distributions of the turbulent dynamic viscosity for Geometry I [ Fig. 8 (a) Fig. 8 (a) ] and Geometry II [ Fig. 8 (b) Fig. 8 (b) ] are both relatively flat; hence, the heat transfer in the core due to turbulent transport is almost uniform. By contrast, the turbulent dynamic viscosity in the core center in Geometry III [ Fig. 8 (c) Fig. 8 (c) ] is low, and it increases gradually along the radial direction; therefore, the turbulent transport effect also increases gradually. As shown in Fig. 6 Fig. 6 , the mean flow velocity in the core center is higher than that near the fertile blanket; therefore, the cooling capacity resulting from fuel salt flow in the core center is better than that near the fertile blanket. Owing to the competing effects of fuel salt flow and turbulent transport, Geometry III has the smallest radial temperature gradient in the core and the most uniform radial temperature distribution among the three optimized geometries. Table 2 Table 2 lists the average core inlet and outlet temperatures, the temperature difference across the core, the core volumes, the total volume flow rate in the primary loop, and the maximum core temperature of each geometry. The heat exchanges are modeled as suitable heat sinks in the energy conservation equations to obtain the nominal value (923 K) of the core inlet temperature. A different momentum source is added to the momentum conservation equation for each geometry to maintain the temperature difference across the core at the nominal level (100 K), because the core volume and shape of each geometry is different. The temperature difference across the core and the average core inlet and outlet temperatures for Geometry I, II, and III are all approximately 100, 923, and 1023 K, respectively. However, the maximum temperatures differ significantly. The maximum fuel salt temperature in Geometry III is 1032.5 K, which is much lower than those of Geometry I (1070.3 K) and Geometry III (1062.1 K). Geometry III exhibits the best TH performance. 1 1 Geometry 1 1 Inlet temperature (K) 1 1 Outlet temperature (K) 1 1 Temperature difference (K) 1 1 Core volume (m 3 ) 1 1 Flow rate (m 3 ∙s −1 ) 1 1 Maximum temperature (K) 1 1 Geometry I 1 1 923.00 1 1 1023.01 1 1 100.01 1 1 9.03 1 1 4.53 1 1 1070.3 1 1 Geometry II 1 1 923.01 1 1 1023.04 1 1 100.03 1 1 9.27 1 1 4.51 1 1 1062.1 1 1 Geometry III 1 1 923.04 1 1 1023.05 1 1 100.01 1 1 9.51 1 1 4.52 1 1 1032.5 3.3 Preliminary thermal-hydraulic analysis of the fertile blanket 2 The fertile blanket has a much smaller heat source than the core; hence, it is neglected in the TH studies of the core. However, the effect of the core heat on the fertile blanket is unknown, and the feasibility of modeling the fertile blanket separately should be analyzed. In this work, a conjugate heat transfer model of the core and fertile blanket is established; it is shown schematically in Fig. 9 Fig. 9 . The calculation model presented in Fig. 5 (c) Fig. 5 (c) is adopted for the neutronics simulation, and the corresponding TH calculation mesh is shown in Fig. 10 Fig. 10 . A 3 GWth heat source is loaded in the core region, and three different power loading cases for the fertile blanket are analyzed. An adiabatic boundary condition is set for the blanket wall. In the OpenFOAM code, a heat flux balance is maintained on the interface wall between the core and the fertile blanket, which can be described as <math display="block" id="M22"><mrow><mi>q</mi><mo>=</mo><mfrac><mrow><msub><mi>λ</mi><mtext>c</mtext></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>c</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>B</mtext></msub><mo stretchy="false">)</mo></mrow><mrow><msub><mi>δ</mi><mtext>c</mtext></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>λ</mi><mtext>b</mtext></msub><mo stretchy="false">(</mo><msub><mi>T</mi><mtext>B</mtext></msub><mo>−</mo><msub><mi>T</mi><mtext>b</mtext></msub><mo stretchy="false">)</mo></mrow><mrow><msub><mi>δ</mi><mtext>b</mtext></msub></mrow></mfrac><mo>,</mo></mrow></math> where q is the heat flux; T B is the temperature of the interface wall; and λ c and λ b are the thermal conductivity of the fuel salt and fertile salt, respectively. T c and T b are the core and fertile blanket cell center temperature of the first control volume on the interface wall, respectively; δ c and δ b are the distance of the interface wall from the core and the fertile blanket cell center, respectively. The temperature distribution including the conjugate heat transfer of the fuel and fertile salts is illustrated in Fig. 11 Fig. 11 . A total heat source of 0 MWth [ Fig. 11 (a) Fig. 11 (a) ], 0.1 MWth [ Fig. 11 (b) Fig. 11 (b) ], and 1 MWth [ Fig. 11 (c) Fig. 11 (c) ] is loaded in the fertile blanket. Fig. 11 (a) Fig. 11 (a) shows that the temperature of the fertile blanket matches that of the core in the axial direction even when no heat source is loaded in the fertile blanket, which reveals that the nuclear heat of the core affects the temperature of the fertile blanket. Fig. 11 (b) and 11 (c) Fig. 11 (b) and 11 (c) indicate that the nuclear heat generated in the fertile blanket cannot be removed by the fuel salt flow through the interface wall between the core and the fertile blanket. Therefore, an external cooling system for the fertile blanket should be considered. 233 U will be generated in the fertile blanket because the thorium loaded in this region is exposed to the core neutron flux, and the heat source (including fission and neutron capture by thorium) generated in the fertile salt is estimated to be approximately 37 MW [ 19 19 , 20 20 ]. However, the temperature of the fertile blanket greatly exceeds that of the core when a heat source of only 1 MWth is loaded in the fertile blanket [ Fig. 11 (c) Fig. 11 (c) ]. Therefore, from Fig. 11 (a) and 11 (c) Fig. 11 (a) and 11 (c) , one can deduce that the effect of the core heat on the fertile blanket is much smaller than that of the heat source in the fertile blanket itself, which demonstrates that modeling the fertile blanket separately is reasonable for the MSFR TH analysis. Three cases are designed to study the required flow rate for obtaining an acceptable temperature distribution in the fertile blanket; the external circuit is ignored for convenience. A sketch map of the three cases is presented in Fig. 12 Fig. 12 . (case 1 [ Fig. 12 (a) Fig. 12 (a) ], case 2 [ Fig. 12 (b) Fig. 12 (b) ] and case 3 [ Fig. 12 (c) Fig. 12 (c) ]). To reduce the computational demands, a minimal symmetrical unit is adopted for the TH calculation; the CFD meshes for the three cases are shown in Fig. 13 Fig. 13 . According to the grid independence verification, the mesh numbers for case 1 [ Fig. 13 (a) Fig. 13 (a) ], case 2 [ Fig. 13 (b) Fig. 13 (b) ], and case 3 [ Fig. 13 (c) Fig. 13 (c) ] are 108,912, 56,576, and 29,772, respectively. The core region of Geometry III [ Fig. 3 (c) Fig. 3 (c) ] is selected to consider the effect of the core neutronics on the neutronics behavior of the fertile blanket. Consequently, the fertile blanket is also divided into several regions to obtain its power distribution, as illustrated in Fig. 14 Fig. 14 . All of the parts are simulated during the neutron transport simulation; a reflective boundary condition is set for the symmetry planes, and a vacuum boundary condition is set for the walls. The heat source of the fertile blanket is estimated to be approximately 37 MW, as mentioned above. Therefore, the proportion of 233 U loaded in the fertile blanket is set to match the power deposition. Considering the temperature tolerance of the structural material, the temperature difference (temperature gradient) between the core and the fertile blanket should be as small as possible. Therefore, the inlet temperature of the fertile blanket is set to 923 K, which is consistent with the average core inlet temperature, by setting the inlet boundary condition to 923 K. In addition, the volume flow rate is adjusted so that the average outlet temperature of the fertile blanket is consistent with the average core outlet temperature (1023 K). The velocity vector profile of the fertile blanket for the three cases is shown in Fig. 15 Fig. 15 . The inlet jet for case 1 [ Fig. 15 (a) Fig. 15 (a) ] is larger than those for case 2 [ Fig. 15 (b) Fig. 15 (b) ] and case 3 [ Fig. 15 (c) Fig. 15 (c) ]. The main reason that the inlet velocity decreases as the number of external circuits increases is that when there are more external circuits, a lower inlet velocity for each pipe can maintain the same volume flow rate of the fertile blanket. A strong inlet jet causes a large reverse flow phenomenon in the entire fertile blanket, as illustrated in Fig. 15 (a) Fig. 15 (a) . A small vortex appears in Fig. 15 (b) and 15 (c) Fig. 15 (b) and 15 (c) , and no reverse flow occurs near the arc wall owing to the relatively low inlet flow velocity. Fig. 16 Fig. 16 shows the temperature distribution of the fertile blanket for the three cases. The temperature distribution of the fertile blanket becomes more uniform as the number of external circuits increases. To maintain a constant outlet temperature of approximately 1023 K, the total volume flow rates of the fertile blanket for the three cases are almost identical to each other owing to their similar geometries. As illustrated above, the required inlet flow velocity for case 1 is highest among the three cases to match the equivalent volume flow rate. As a result, the inlet jet in case 1 causes a highly nonuniform temperature distribution, and a relatively high temperature appears at the bottom arc, which is higher than the outlet temperature of the fertile blanket. Furthermore, a local high temperature, which is much higher than the outlet temperature, appears in the top cavity of the fertile blanket. By contrast, because no reverse flow occurs near the arc wall, the temperature distribution of the top cavity for cases 2 and 3 is relatively flat, and no hot spot forms. Table 3 Table 3 shows the maximum temperature, average outlet temperature, and required total volume flow rate of the fertile blanket. Case 3 has the lowest maximum temperature and the most uniform temperature distribution. Therefore, it has the best TH performance among the three cases. However, the calculated maximum temperature and temperature distribution for case 2 are also acceptable; this result indicates that eight cooling circuits with a total volume flow rate of approximately 0.0555 m 3 ∙s −1 can also produce an acceptable temperature distribution. Considering the engineering requirements, the geometry model for case 2 is recommended. 1 1   1 1 Maximum temperature (K) 1 1 Outlet temperature (K) 1 1 Flow rate (m 3 ∙s −1 ) 1 1 Case 1 1 1 1097.42 1 1 1023.10 1 1 0.0539 1 1 Case 2 1 1 1026.85 1 1 1023.04 1 1 0.0555 1 1 Case 3 1 1 1026.17 1 1 1023.24 1 1 0.0565 4 Conclusion 1 In the core of the MSFR benchmark model, a large flow recirculation zone appears near the fertile blanket, and the fuel salt is almost stagnant at the core center, which is close to the axial reflectors. The recirculation in this region produces a local hot spot on the wall (near the fertile blanket), and thus results in unnecessary thermal stress and material problems [ 12 12 ]. A preliminary conjugate heat transfer analysis of the core and fertile blanket indicates the feasibility of using a separate model for the MSFR TH analysis and demonstrates that the nuclear heat generated in the fertile blanket cannot be removed by heat convection due to the fuel salt flow. Therefore, an external cooling system for the fertile blanket has to be considered to remove the heat generated by fission and neutron capture by thorium. In this work, a neutronics and TH coupling scheme between the DAG-OpenMC and OpenFOAM codes is adopted to further optimize the core design and conduct a preliminary TH analysis of the fertile blanket. The main conclusions are as follows: 1) By optimizing the core geometry, the large recirculation region near the fertile blanket in the reference MSFR model is eliminated, and a relatively flat velocity distribution is obtained for all three studied geometries. Moreover, the stagnation zones disappear in the velocity field of Geometry III. The hot spot that appears in the upper core cannot be effectively eliminated by curving only the upper wall of the core. However, the peak temperature for Geometry III is much lower than that for the other two geometries. Furthermore, a more uniform distribution and a relatively small radial gradient of the fuel temperature are obtained for Geometry III. In terms of the thermal hydraulics, Geometry III, which has curved upper and lower core walls, has the best performance. 2) For the fertile blanket, the strong inlet jet in case 1 produces a significant reverse flow in the entire fertile blanket, which results in a highly nonuniform temperature distribution and a relatively high temperature at the bottom arc. In addition, a local high temperature appears in the top cavity of the fertile blanket. For cases 2 and 3, no reverse flow occurs near the arc wall, and no hot spot forms owing to the relatively small inlet jet and relatively flat temperature distribution in the top cavity. Case 3 (with 16 circuits) has the most uniform temperature distribution and the lowest maximum temperature; the calculated result for case 2 (with 8 circuits) is also acceptable. Considering the engineering requirements, the geometry model of case 2 is recommended. References Doe U . A technology road map for generation IV nuclear energy system . Proceedings of the Nuclear Energy Research Advisory Committee and the Generation IV International Forum . 239 - 241 ( 2002 ). A technology road map for generation IV nuclear energy system E.S. Bettis , R.W. Schroeder , G.A. Cristy . The aircraft reactor experiment design and construction . Nucl Sci Eng . 2 , 804 - 825 ( 1957 ). The aircraft reactor experiment design and construction R. C. Robertson . MSRE Design and Operations Report Part I: Description of Reactor Design . ORNL, USA ( 1965 ). MSRE Design and Operations Report Part I: Description of Reactor Design X.C. Zhao , D.Y. Cui , X.Z. Cai , et al . Analysis of Th-U breeding capability for an accelerator-driven subcritical molten salt reactor . Nucl. Sci. Tech . 29 , 121 ( 2018 ). https://doi.org/10.1007/s41365-018-0448-3 https://doi.org/10.1007/s41365-018-0448-3 Analysis of Th-U breeding capability for an accelerator-driven subcritical molten salt reactor X.X. Li , D.Y. Cui , Y.W. Ma , et al . Influence of 235U enrichment on moderator temperature coefficient of reactivity (MTC) in a graphite-moderator Molten Salt Reactor (MSR) . Nucl. Sci. Tech . 30 , 166 ( 2019 ). https://doi.org/10.1007/s41365-019-0694-z https://doi.org/10.1007/s41365-019-0694-z Influence of 235U enrichment on moderator temperature coefficient of reactivity (MTC) in a graphite-moderator Molten Salt Reactor (MSR) K.F. Ma , C.G. Yu , X.Z. Cai , et al . Transmutation of 129I in a single-fluid double-zone Thorium Molten Salt Reactor . Nucl. Sci. Tech . 31 , 10 ( 2020 ). https://doi.org/10.1007/s41365-020-0720-1 https://doi.org/10.1007/s41365-020-0720-1 Transmutation of 129I in a single-fluid double-zone Thorium Molten Salt Reactor L.Y. He , G.C. Li , S.P. Xia , et al . Effect of 37Cl enrichment on neutrons in a Molten Chloride Salt Fast Reactor . Nucl. Sci. Tech . 31 , 27 ( 2020 ). https://doi.org/10.1007/s41365-020-0740-x https://doi.org/10.1007/s41365-020-0740-x Effect of 37Cl enrichment on neutrons in a Molten Chloride Salt Fast Reactor V. Ignatiev , O. Feynberg , I. Gnidoi . Progress in development of Li, Be, Na/F molten salt actinide recycler & transmuter concept . Proceedings of ICAPP . 13 - 18 ( 2007 ). Progress in development of Li, Be, Na/F molten salt actinide recycler & transmuter concept Brovchneko , D. Heuer , E. Merle-Lucotte . Preliminary safety calculations to improve the design of Molten Salt Fast Reactor . In: Proc. Int. Conf. PHYSOR 2012 . Knoxville, TN, USA ( 2012 ). Preliminary safety calculations to improve the design of Molten Salt Fast Reactor C. Renault , M. Delpech , C. Le Brun , et al . The MOST project: key-point and challenges for the feasibility of molten salt reactor . Paper , 5208 , 15 - 19 ( 2005 ). The MOST project: key-point and challenges for the feasibility of molten salt reactor E. Linden . Coupled Neutronics and Computational Fluid Dynamics for the Molten Salt Fast Reactor. MSc Thesis . Technical University of Delft , Netherlands ( 2012 ). Coupled Neutronics and Computational Fluid Dynamics for the Molten Salt Fast Reactor. MSc Thesis C. Fiorina , D. Lathouwers , M. Aufiero , et al . Modelling and analysis of the MSFR transient behaviour . Ann. Nucl. Energy . 64 , 485 - 498 ( 2014 ). https://doi.org/10.1016/j.anucene.2013.08.003 https://doi.org/10.1016/j.anucene.2013.08.003 Modelling and analysis of the MSFR transient behaviour Q.Y. Li , B. Xu , C. Zhou , et al . The optimization of flow rate distribution design of 373 MW molten salt reactor-liquid fuel . Nucl. Tech . 42 , 070605 ( 2019 ). https://doi.org/10.11889/j.0253-3219.2019.hjs.42.070605 https://doi.org/10.11889/j.0253-3219.2019.hjs.42.070605 . (in Chinese) The optimization of flow rate distribution design of 373 MW molten salt reactor-liquid fuel L. He , C.G. Yu , W. Guo , et al . The temperature distribution of molten salt reactor under the outer salt layer blocked situation . Nucl. Tech . 41 , 050601 ( 2018 ). https://doi.org/10.11889/j.0253-3219.2018.hjs.41.050601 https://doi.org/10.11889/j.0253-3219.2018.hjs.41.050601 . (in Chinese) The temperature distribution of molten salt reactor under the outer salt layer blocked situation P.R. Rubiolo , A. Laureau , E. Merle-Lucotte , et al . Core thermal-hydraulics analyses of the Molten Salt Fast Reactor (MSFR) in nominal conditions .  Winter Meeting of the American Nuclear Society  ( 2013 ). http://hal.in2p3.fr/in2p3-00907007 http://hal.in2p3.fr/in2p3-00907007 Core thermal-hydraulics analyses of the Molten Salt Fast Reactor (MSFR) in nominal conditions A. Laureau , P.R. Rubiolo , D. Heuer , et al . Coupled neutronics and thermal-hydraulics numerical simulations of a Molten Fast Salt Reactor (MFSR) . Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo ( 2013 ). https://doi.org/10.1051/snamc/201402307 https://doi.org/10.1051/snamc/201402307 Coupled neutronics and thermal-hydraulics numerical simulations of a Molten Fast Salt Reactor (MFSR) H. Rouch , O. Geoffroy , P. Rubiolo , et al . Preliminary thermal-hydraulic core design of the Molten Salt Fast Reactor (MSFR) . Ann. Nucl. Energy . 64 , 449 - 456 ( 2014 ). https://doi.org/10.1016/j.anucene.2013.09.012 https://doi.org/10.1016/j.anucene.2013.09.012 Preliminary thermal-hydraulic core design of the Molten Salt Fast Reactor (MSFR) A. Laureau , D. Heuer , E. Merle-Lucotte , et al . Transient coupled calculations of the Molten Salt Fast Reactor using the Transient Fission Matrix approach . Nucl. Eng. Des . 316 , 112 - 124 ( 2017 ). https://doi.org/10.1016/j.nucengdes.2017.02.022 https://doi.org/10.1016/j.nucengdes.2017.02.022 Transient coupled calculations of the Molten Salt Fast Reactor using the Transient Fission Matrix approach M. Brovchenko , J.L. Kloosterman , L. Luzzi , et al . Neutronic benchmark of the molten salt fast reactor in the frame of the EVOL and MARS collaborative projects . EPJ Nucl. Sci. Technol . 5 , 2 ( 2019 ). https://doi.org/10.1051/epjn/2018052 https://doi.org/10.1051/epjn/2018052 Neutronic benchmark of the molten salt fast reactor in the frame of the EVOL and MARS collaborative projects M. Brovchenko , E. Merle-Lucotte , H. Rouch , et al . Evaluation and Viability of Liquid Fuel Fast Reactor System EVOL ( 2013 ). https://www.li2c.upmc.fr https://www.li2c.upmc.fr Evaluation and Viability of Liquid Fuel Fast Reactor System EVOL K. Paul , E. Nicholas , R. Bryan , et al . OpenMC: a state-ofthe-art Monte Carlo code for research and development . Ann. Nucl. Energy . 82 , 90 - 97 ( 2015 ). https://doi.org/10.1016/j.anucene.2014.07.048 https://doi.org/10.1016/j.anucene.2014.07.048 OpenMC: a state-ofthe-art Monte Carlo code for research and development DAGMC: Direct Accelerated Geometry Monte Carlo . https://svalinn.github.com/DAGMC https://svalinn.github.com/DAGMC DAGMC: Direct Accelerated Geometry Monte Carlo K. Paul , E. Nicholas , R. Bryan , et al . The OpenMC Monte Carlo code - OpenMC documentation . https://docs.openmc. org/en/v0.10.0/methods/cross_sections.html#temperature-treatment https://docs.openmc. org/en/v0.10.0/methods/cross_sections.html#temperature-treatment . The OpenMC Monte Carlo code - OpenMC documentation OpenFOAM Doumentation . https://www.openfoam.org/docs https://www.openfoam.org/docs OpenFOAM Doumentation M. Cheng , Z. Dai . Development of a three dimension multi-physics code for molten salt fast reactor . Nucl. Sci. Tech . 25 , 010601 ( 2014 ). https://doi.org/10.13538/j.1001-8042/nst.25.010601 https://doi.org/10.13538/j.1001-8042/nst.25.010601 Development of a three dimension multi-physics code for molten salt fast reactor F.J. Wang . Computational Fluid Dynamics analysis--Principle and application of CFD software . Tsinghua University Press ( 2004 ). (in Chinese) Computational Fluid Dynamics analysis--Principle and application of CFD software L. He , C.G. Yu , W. Guo , et al . Development of a coupled code for steady-state analysis of the graphite-moderated channel type molten salt reactor . Sci. Tech. Nucl. Instal . 2018 , 1 - 10 ( 2018 ). https://doi.org/10.1155/2018/4053254 https://doi.org/10.1155/2018/4053254 Development of a coupled code for steady-state analysis of the graphite-moderated channel type molten salt reactor L. He , C.G. Yu , R.M. Ji , et al . Development of a dynamics model for graphite-moderated channel-type molten salt reactor . Nucl. Sci. Tech . 30 , 18 ( 2019 ). https://doi.org/10.1007/s41365-018-0541-7 https://doi.org/10.1007/s41365-018-0541-7 Development of a dynamics model for graphite-moderated channel-type molten salt reactor J.H. CHEN , H.Q. ZHANG , Z.Y. ZHU . Coupled neutronic and thermal-hydraulic analysis of TMSR-SF1 at steady state . Nucl. Tech . 40 , 080603 ( 2017 ). https://doi.org/10.11889/j.0253-3219.2017.hjs.40.080603 https://doi.org/10.11889/j.0253-3219.2017.hjs.40.080603 . (in Chinese) Coupled neutronic and thermal-hydraulic analysis of TMSR-SF1 at steady state B. DENG , J.G. CHEN , L. HE , et al . Development and verification of a coupled code of steady-state neutronics and thermal-hydraulics for molten salt fast reactor . Nucl. Tech . 42 , 110602 ( 2019 ). https://doi.org/10.11889/j.0253-3219.2019.hjs.42.110602 https://doi.org/10.11889/j.0253-3219.2019.hjs.42.110602 . (in Chinese) Development and verification of a coupled code of steady-state neutronics and thermal-hydraulics for molten salt fast reactor L.L.W. Frima . Burnup in a molten salt fast reactor. MSc Thesis . Technical University of Delft , Netherlands ( 2013 ). Burnup in a molten salt fast reactor. MSc Thesis H. K. Versteeg . An Introduction to Computational Fluid Dynamics: The Finite Volume Method . New York : Wiley ( 1995 ). An Introduction to Computational Fluid Dynamics: The Finite Volume Method 10.1007/s41365-020-00803-9 This work was supported by the Chinese TMSR Strategic Pioneer Science and Technology Project (No. XDA02010000) and the Frontier Science Key Program of the Chinese Academy of Sciences (No. QYZDY-SSW-JSC016). 2020-05-15 2020-07-21 2020-07-22 2020-08-27 2020-09-01 2021-04-10 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2020 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-020-00803-9 http://dx.doi.org/10.1007/s41365-020-00803-9 . 2020 Nuclear Science and Techniques 9 31